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