compute endobs = 5000, fstart=endobs+1, fend=endobs+20



cal(daily) 1990 1 1

all endobs

open data archdata.rat

data(for=rats) / y x



set v = 0.0

nonlin b0 b1 a0 a1

frml regresid = y-b0-b1*x

frml archvar  = a0+a1*regresid(t-1)**2

frml archlogl = (v(t)=archvar(t)),-.5*(log(v)+regresid(t)**2/v)

linreg y

# constant x  

com b0=%beta(1),b1=%beta(2),a0=%seesq,a1=0.05

maximize(met=bhhh,recursive) archlogl 2 *





**** Compute forecasts:

set x endobs+1 endobs+20 = t

frml(identity) yfrml y = b0 + b1*x 

group ymod  yfrml>>yhat



smpl fstart-20 fend



forecast(model=ymod) * 

print / y yhat

graph(key=upleft) 2

# y

# yhat



**** simulate forecasted residuals using variance information, compute 

****   fitted values using those residuals.

compute ndraws=500

clear usim vsim ysim 

set xx 1 endobs+20 = b0 + b1*x

compute sigma=a0 , alpha=a1 

frml(variance=sigma**2) vdef vsim = 0.0

frml(identity) udef usim = vsim*sqrt(1+alpha*usim{1}**2)

frml(identity) ydef ysim = xx + usim

smpl fstart fend

set yfirst  = 0.0

set ysecond = 0.0

group archmod vdef>>vsim udef ydef>>ysim

compute usim(endobs) = y(endobs) - b0 - b1*x(endobs)

*

do draws=1,ndraws

  simulate(model=archmod) 2 

  set yfirst  = ysim + yfirst

  set ysecond = ysim**2  + yfirst

end do draws



set yfirst = (1./ndraws)*yfirst

set ysecond = (1./ndraws)*ysecond

set stderr = sqrt(ysecond)

set upper = ysim + 2*stderr

set lower = ysim - 2*stderr

smpl fstart-20 fend

print



graph(key=upleft) 5

# y

# yhat

# ysim

# upper

# lower

smpl

dedit(new) archsim.rat

store y yhat ysim upper lower xx

save