version 10.0
capture log close
log using appx_c, replace

* examine hip data
use hip, clear
histogram y, percent
more

summarize y, detail

* estimate mean
mean y

* generate several normal variables
clear
matrix C = (1, .5 \ .5, 1)
drawnorm x y, n(1000) corr(C) seed(12345)
summarize x y
tabstat x y, statistics (mean median variance semean)
corr x y
twoway scatter y x
more

* central limit theorem
clear
set obs 1000
set seed 12345
gen y1 = sqrt(uniform())
histogram y1
more

gen y2 = sqrt(uniform())
gen y3 = sqrt(uniform())
gen y4 = sqrt(uniform())
gen y5 = sqrt(uniform())
gen y6 = sqrt(uniform())
gen y7 = sqrt(uniform())
gen y8 = sqrt(uniform())
gen y9 = sqrt(uniform())
gen y10 = sqrt(uniform())
gen y11 = sqrt(uniform())
gen y12 = sqrt(uniform())
gen ybar3 = (y1+y3+y3)/3
gen ybar7 = (y1+y3+y3+y4+y5+y6+y7)/7
gen ybar12 = (y1+y3+y3+y4+y5+y6+y7+y8+y9+y10+y11+y12)/12
gen z3 = (ybar3 - 2/3)/sqrt((1/18)/3)
gen z7 = (ybar7 - 2/3)/sqrt((1/18)/7)
gen z12 = (ybar12 - 2/3)/sqrt((1/18)/12)

histogram z3, normal
more

summarize z3, detail

histogram z7, normal
more

summarize z7, detail

histogram z12, normal
more

summarize z12, detail

* interval estimates
* simulated data
clear
set obs 30
set seed 12345
drawnorm x1-x10
gen y1 = 10 + sqrt(10)*x1
gen y2 = 10 + sqrt(10)*x2
gen y3 = 10 + sqrt(10)*x3
gen y4 = 10 + sqrt(10)*x4
gen y5 = 10 + sqrt(10)*x5
gen y6 = 10 + sqrt(10)*x6
gen y7 = 10 + sqrt(10)*x7
gen y8 = 10 + sqrt(10)*x8
gen y9 = 10 + sqrt(10)*x9
gen y10 = 10 + sqrt(10)*x10
ci y1-y10

* hip data
use hip, clear
ci y

quietly summarize y, detail
return list
scalar ybar = r(mean)
scalar nobs = r(N)
scalar df = nobs - 1
scalar tc975 = invttail(df,.025)
scalar sighat = r(sd)
scalar se = sighat/sqrt(nobs)
scalar lb = ybar - tc975*se
scalar ub = ybar + tc975*se

di "lb of 95% confidence interval " lb
di "ub of 95% confidence interval " ub

* hypothesis testing

* right tail test mu = 16.5
* automatic version
ttest y==16.5

* details
quietly summarize y, detail
scalar t1 = (ybar - 16.5)/se
scalar tc95 = invttail(df,.05)
scalar p1 = ttail(df,t1)
di "right tail test"
di "tstat = " t1
di "tc95  = " tc95
di "pval  = " p1

* two tail test mu = 17
* automatic version
ttest y==17

* details
quietly summarize y, detail
scalar t2 = (ybar - 17)/se
scalar p2 = 2*ttail(df,abs(t2))
di "two tail test"
di "tstat = " t2
di "tc975  = " tc975
di "pval  = " p2

* Testing variance = 4

* automatic test
sdtest y == 2

* details
quietly summarize y, detail
scalar s0 = 4
scalar sighat2 = r(Var)
scalar df = r(N)-1
scalar v = df*sighat2/s0
scalar chi2_95 = invchi2(df,.95)
scalar chi2_05 = invchi2(df,.05)
scalar p = 2*chi2(df,v)
di "Chi square test stat " v
di "5th percentile chisquare(49) " chi2_05
di "95th percentile chisquare(49) " chi2_95 
di "2 times p value " p

* testing population means
clear
drawnorm x1 x2, n(50) means(1 2) seed(12345)
summarize

* assuming variances are equal
ttest x1 == x2, unpaired

* assuming variances unequal
drawnorm x3 x4, n(50) means(1 2) sds(1 2) seed(12345)
ttest x3 == x4, unpaired unequal

* testing population variances
sdtest x3 == x4

* test normality
use hip, clear

* automatic test
sktest y

* Jarque_Bera test
quietly summarize y, detail
scalar nobs = r(N)
scalar s = r(skewness)
scalar k = r(kurtosis)
scalar jb = (nobs/6)*(s^2 + ((k-3)^2)/4)
scalar chi2_95 = invchi2(2,.95)
scalar pval = 1 - chi2(2,jb)
di "jb test statistic " jb
di "95th percentile chi2(2) " chi2_95
di "pvalue " pval

log close
translate appx_c.smcl appx_c.txt, replace