' bootstrap sample mean and median
' series proc version (slower than matrix version)
' version 4 
' 11/2/99 h
' revised and last checked 9/27/2000 h


' set size of sample
!n = 50

' create workfile
workfile bootmed1 u 1 !n

' set seed of random number generator
rndseed 456

' generate pseudo-draws from a standard normal
series x = nrnd

' set number of bootstrap replications
!reps = 1000

' store means and medians in matrix
matrix(!reps,2) out

' bootstrap loop
tic
for !i = 1 to !reps
	' generate bootstrap sample
	x.resample x_b
	' store bootstrap median
	out(!i,1) = @median(x_b)
	' store bootstrap mean
	out(!i,2) = @mean(x_b)
next
scalar elp = @toc

' display descriptive statistics
show out.stats

' theoretical standard error of median from normal
' se(med) = 1 / ( 2*m*(obs)^0.5 )
' where m is the median value
' (Kendall-Stuart, 1977, 4th ed, vol 1, p.252)
!pi = @acos(-1)
scalar se_med = @sqrt(0.5*!pi/!n)
scalar se_mea = @sqrt(1/!n)

' convert bootstrap sample into series
expand 1 !reps
smpl 1 !reps

series out_medi
series out_mean
group gout out_medi out_mean
mtos(out, gout)

' display results in table
table tab
setcolwidth(tab,1,20)

tab(1,2) = "median"
tab(1,3) = "mean"

setline(tab,2)

tab(3,1) = "bootstrap mean:"
tab(3,2) = @mean(out_medi)
tab(3,3) = @mean(out_mean)

tab(4,1) = "bootstrap s.d.:"
tab(4,2) = @stdev(out_medi)
tab(4,3) = @stdev(out_mean)

tab(5,1) = "theoretical s.d.:"
tab(5,2) = se_med
tab(5,3) = se_mea

setline(tab,6)

tab(7,1) = "sample size = "
tab(7,2) = @str(!n)

tab(8,1) = "bootstrap replications = " 
tab(8,2) = @str(!reps)

tab(9,1) = "elapsed time = " 
tab(9,2) = @str(elp) + " seconds"
show tab

'get kernel density estimates
!n=100	'points to evaluate
do out_mean.kdensity(!n,o=dist_mea,b=0.04)
do out_medi.kdensity(!n,o=dist_med,b=0.2)

smpl 1 !n

'convert to series
series x_mean
series f_mean
group g1 x_mean f_mean
mtos(dist_mea, g1)

series x_medi
series f_medi
group g2 x_medi f_medi
mtos(dist_med, g2)

'compare bootstrap distribution
group g3 g1 g2 ' f_mean x_mean f_medi x_medi
freeze(gra1) g3.xyline(b)
'gra1.setelem(1) legend(mean)
'gra1.setelem(2) legend(median)
gra1.addtext(t) Kernel density estimates of bootstrap distribution
gra1.legend -display
gra1.addtext(1.3,0.2) Mean
gra1.addtext(0.2,2.2) Median
show gra1

