*
* BOOTARMA.RPF
* Bootstraps an ARMA model.
* Reference: David S. Stoffer and Kent D. Wall(1991), "Bootstrapping
* State-Space Models: Gaussian Maximum Likelihood Estimation and the
* Kalman Filter",  Journal of the American Statistical Association, Vol.
* 86, No. 416, pp. 1024-1033
*
open data income.asc
cal(q) 1971
all 1985:2
data(format=free,org=columns) / income consum
*
set logy = log(income)
@bjident(diffs=1) logy
*
* Estimate the ARMA model in differenced form. It's simplest to handle
* the mean (if any) by extracting it first, rather than estimating it as
* a coefficient.
*
set dy = logy-logy{1}
boxjenk(ma=||2,4||,demean,maxl,define=dyeq) dy
compute bstart=%regstart(),bend=%regend()
*
* Save the extracted mean
*
compute dymean = %mean
*
* Convert the ARMA equation to a state-space representation
*
@armadlm(a=adlm,f=fdlm) dyeq
***********************************************************************************
*
* Random simulation can be done in a straightforward fashion using DLM
* with TYPE=SIMULATE.
*
dlm(type=simulate,presample=ergodic,a=adlm,f=fdlm,sw=%seesq) bstart bend xstates
*
* Add the extracted mean back to give the same mean as the original series
*
set dysim = xstates(t)(1)+dymean
***********************************************************************************
* Bootstrap
* First, need to get the residuals and the predictive variances
*
dlm(method=solve,presample=ergodic,a=adlm,f=fdlm,y=dy-dymean,c=%unitv(%rows(adlm),1),sw=1.0,$
   type=filter,vhat=vhat,svhat=svhat) bstart bend xstates vstates
*
* Get standardized residuals
*
set svscale   bstart bend = sqrt(%scalar(svhat))
set stdresids bstart bend = %scalar(vhat)/svscale
*
* The following (until the "end do time") is repeated for each bootstrap
* simulation.
*
boot entries bstart bend
*
compute [vector] x0=%zeros(%rows(adlm),1)
compute [symm]  sx0=%psdinit(adlm,%outerxx(fdlm))
clear dyboot
do time=bstart,bend
   *
   * Compute the Kalman gain for this period's update
   *
   compute sx0 =adlm*sx0*tr(adlm)+%outerxx(fdlm)
   compute gain=%xcol(sx0,1)*(1.0/sx0(1,1))
   *
   * Adjust the state vector by using the Kalman gain applied to the
   * reshuffled standardized residual reflated as appropriate to this
   * period's predictive variance.
   *
   compute x0  =adlm*x0+gain*svscale(time)*stdresids(entries(time))
   compute sx0 =vstates(time)
   compute dyboot(time)=x0(1)+dymean
end do time