*
* Example 11.1 fom pp 492-508
*
open data aa-3rv.txt
data(format=free,org=columns) 1 340 aa5m aa10m aa20m
*
set aa10m = log(aa10m)
*
graph(footer="Figure 11.1 Time plot of logarithms of intraday realized volatility")
# aa10m
*
* Estimate by maximum likelihood as an IMA(1,1)
*
boxjenk(maxl,diff=1,ma=1) aa10m
corr(qstats,number=12,dfc=1) %resids
set usq = %resids**2
corr(qstats,number=12,dfc=0) usq
*
* Estimate the equivalent model in state space form
*
nonlin sv sw
compute sv=sw=1.0
dlm(y=aa10m,c=1.0,a=1.0,sv=sv,sw=sw,method=bfgs,exact)
*
* Analysis of fitted model
* (continuation page 496)
*
* The DLM instructions below are basically the same, with minor changes to get
* different pieces of information. The one which does produce generally different
* results is with TYPE=SMOOTHED. That returns smoothed (rather than filtered)
* estimates of the states, variances and errors.
*
*
dlm(y=aa10m,c=1.0,a=1.0,sv=sv,sw=sw,method=solve,exact,vhat=vhat) / xstates
set filtered = xstates(t)(1)
set error    = vhat(t)(1)
*
spgraph(footer="Figure 11.2 Time plots of output of the Kalman filter",vfields=2)
graph(header="(a) Filtered state variable")
# filtered
graph(header="(b) Prediction error")
# error
spgraph(done)
*
* Further analysis of fitted model
* (continuation page 501)
*
dlm(y=aa10m,c=1.0,a=1.0,sv=sv,sw=sw,method=solve,exact) / xstates vstates
set filtered  = xstates(t)(1)
set fvariance = vstates(t)(1,1)
set upper     = filtered+%invnormal(.975)*sqrt(fvariance)
set lower     = filtered+%invnormal(.025)*sqrt(fvariance)
graph(footer="Figure 11.3 Filtered state variable and 95% confidence interval") 3
# filtered
# upper / 2
# lower / 2
dlm(y=aa10m,c=1.0,a=1.0,sv=sv,sw=sw,method=solve,exact,type=smoothed) / xstates vstates
set smoothed  = xstates(t)(1)
set svariance = vstates(t)(1,1)
set upper     = smoothed+%invnormal(.975)*sqrt(svariance)
set lower     = smoothed+%invnormal(.025)*sqrt(svariance)
graph(footer="Figure 11.4 Smoothed state variable and 95% confidence interval") 3
# smoothed
# upper / 2
# lower / 2
*
* Diagnostics (described on page 508)
* Get the series of standardized residuals and run tests on them
*
dlm(y=aa10m,c=1.0,a=1.0,sv=sv,sw=sw,method=solve,exact,vhat=vhat,svhat=svhat) / xstates
set zerrors  = vhat(t)(1)/sqrt(svhat(t)(1,1))
*
@archtest(lags=25) zerrors
@regcorrs(number=25) zerrors