*
* Chapter 8. UK data. Diagnostics.
*
open data ukdriversksi.txt
calendar(m) 1969
data(format=free,org=columns,skips=1) 1969:01 1984:12 ksi
set logksi = log(ksi)
open data logukpetrolprice.txt
data(format=free,org=columns,skips=1) 1969:01 1984:12 logukpetrol
set seatbelt = t>=1983:2
*
* Fixed coefficients regression done with linreg
*
seasonal seasons
linreg logksi
# constant logukpetrol seatbelt seasons{0 to -10}
*
compute sigsqeps=%sigmasq,beta=%beta(2),lambda=%beta(3)
*
* Create the standard structural model for local level + seasonal
*
@LocalDLM(type=level,a=al,c=cl,f=fl)
@SeasonalDLM(type=additive,a=as,c=cs,f=fs)
compute a=al~\as,f=fl~\fs,c=cl~~cs
*
frml explan = beta*logukpetrol+lambda*seatbelt
*
* Model with stochastic level, deterministic seasonal
*
compute sigsqxi=sigsqeps*.01
nonlin sigsqeps sigsqxi sigsqomega=0.00 lambda beta
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  method=bfgs,vhat=vhat,svhat=svhat)
*
set resids = %scalar(vhat)/sqrt(%scalar(svhat))
*
graph(footer="Figure 8.8 Standardized one-step prediction errors")
# resids
*
* RATS uses slightly different (bias-corrected) formulas for the
* skewness and kurtosis, hence the normality statistic is a bit
* different.
*
@STAMPDiags(ncorrs=10) resids
*
@histogram(maxgrid=20,distrib=normal,$
 footer="Figure 8.10 Histogram of standardized one-step prediction errors") resids
*
* Section 4.3 does the stochastic level, deterministic seasonal, but
* without the explanatory variables. The smoothed disturbances are
* computed using Kalman smoothing.
*
nonlin sigsqeps sigsqxi sigsqomega=0.00
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi,$
  method=bfgs,type=smooth,vhat=vhat,what=what)
set outlier = %scalar(vhat)
diff(standardize) outlier
set breaks  = %scalar(what)
diff(standardize) breaks
spgraph(vfields=2,samesize,$
 footer="Figure 8.11 Standardized smoothed level and observation disturbances\\"+$
 "for model without explanatory variables")
graph(hlabel="Structural level break t-tests",vgrid=||-2.0,2.0||)
# breaks
graph(hlabel="Outlier t-tests",vgrid=||-2.0,2.0||)
# outlier
spgraph(done)
*
* Same thing, but with the model with the explanatory variables
*
nonlin sigsqeps sigsqxi sigsqomega=0.00 lambda beta
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  method=bfgs,type=smooth,vhat=vhat,what=what)
set outlier = %scalar(vhat)
diff(standardize) outlier
set breaks  = %scalar(what)
diff(standardize) breaks
spgraph(vfields=2,samesize,$
 footer="Figure 8.12 Standardized smoothed level and observation disturbances\\"+$
 "for model with explanatory variables")
graph(hlabel="Structural level break t-tests",vgrid=||-2.0,2.0||)
# breaks
graph(hlabel="Outlier t-tests",vgrid=||-2.0,2.0||)
# outlier
spgraph(done)