' bootstrap sample mean and median
' matrix function version (faster than series version)
' version 4 
' 11/2/99 h
' last checked 9/19/2000 h

' create workfile
workfile bootmedian2 u 1 1

' set seed of random number generator
rndseed 456

' set size of sample
!n = 50

' generate pseudo-draws from a standard normal
matrix(!n,1) x
nrnd(x)

' set number of bootstrap replications
!reps = 1000

' store means and medians in vectors
vector(!reps,1) out_medi
vector(!reps,1) out_mean

' declare matrix to hold bootstrap sample
matrix x_b

' bootstrap loop
' reset timer
tic
for !i = 1 to !reps
	' generate bootstrap sample
	x_b = @resample(x)
	' store bootstrap median
	out_medi(!i) = @median(x_b)
	' store bootstrap mean
	out_mean(!i) = @mean(x_b)
next
' store elapsed time in seconds
scalar elp = @toc

' 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)

' 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

