*
*  SmallSampleBias.PRG
*
*

seed 125388
comp ndraws=1000
comp useobs=100
comp init=50
comp endobs=init+useobs
all endobs

comp sigma = 1.0
comp alpha = 0.1


comp beta = 0.50

frml(variance=sigma**2) vdef v = 0.0
frml(identity) ydef y = alpha + beta*y{1} + v
group ar1model vdef>>v ydef>>y

set v 1 1 = 0.0
set y 1 1 = 0.0
set betasimu 1 ndraws = 0.0

do draws=1,ndraws
simulate(model=ar1model) * endobs-1 2
linreg(noprint) y init+1 endobs resids
# y{1} constant 


comp betasimu(draws) = %beta(1)

end do draws

density(type=gaussian) betasimu 1 ndraws sgrid sdensity
scatter(style=lines,hgrid=||beta||) 1
# sgrid sdensity

stat betasimu

set bias 1 ndraws = 0.0
do jj=1,ndraws
comp bias(jj) = %if(betasimu(jj)<beta,1.0,0.0) 
end do draws

stat bias

comp c01=fix(0.01*ndraws)
comp c05=fix(0.05*ndraws)
comp c10=fix(0.10*ndraws)
comp c50=fix(0.50*ndraws)
comp c90=fix(0.90*ndraws)
comp c95=fix(0.95*ndraws)
comp c975=fix(0.975*ndraws)
comp c99=fix(0.99*ndraws)

comp sortbetasimu=%sort(betasimu)

comp beta01=sortbetasimu(c01,1) 
comp beta05=sortbetasimu(c05,1) 
comp beta10=sortbetasimu(c10,1) 
comp beta50=sortbetasimu(c50,1) 
comp beta90=sortbetasimu(c90,1) 
comp beta95=sortbetasimu(c95,1) 
comp beta99=sortbetasimu(c99,1) 


disp "beta"
disp $
"      1%            5%            10%           50%           90%          95%           99%"
disp beta01 beta05 beta10 beta50 beta90 beta95 beta99