*
*  SSBBoot.PRG
*
*

seed 12534
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.99

comp nboot = 1000



*********************************************
*
* Bias-corrected Bootstrap Estimator
*
*********************************************

Procedure BCBE x bstart bend nboot bcar1
type series x
type int nboot
type real *bcar1
type int bstart
type int bend
local series bcar1simu
local series xb
local int iboot
local series residx
local series pathx
local vect[series] bootsample
local real betahat
local real bootbetabar



linreg(noprint,define=eqnx) x bstart bend residx
# x{1} constant

comp betahat = %beta(1)

set bcar1simu 1 nboot = 0.0
do iboot=1, nboot
boot entries 2 useobs bstart+1 bend
set pathx = residx(entries)

forecast(paths,results=bootsample) 1 useobs-1 2
# eqnx
# pathx


set xb = bootsample(1)

linreg(noprint) xb
# xb{1} constant
comp bcar1simu(iboot) = %beta(1)



end do iboot

stat(noprint) bcar1simu
comp bootbetabar = %mean
comp bcar1 = 2.0*betahat - bootbetabar

end BCBE









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
set betasimu_bc 1 ndraws = 0.0

infobox(action=define,progress,lower=1,upper=ndraws) 'Progress'
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)

@BCBE y init+1 endobs nboot bcbar1
comp betasimu_bc(draws) = bcbar1


infobox(current=draws)
end do draws
infobox(action=remove)

density(type=gaussian) betasimu 1 ndraws sgrid sdensity
density(type=gaussian) betasimu_bc 1 ndraws sgrid_bc sdensity_bc


scatter(style=lines,overlay=line,patterns,ovsame,hgrid=||beta||) 2
# sgrid sdensity
# sgrid_bc sdensity_bc



stat betasimu

stat betasimu_bc


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


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

stat bias_bc


/*
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
*/