*
*  Illustration from section 2.12.3
*  Diagnostics
*
open data nile.dat
calendar 1871
data(format=free,org=columns,skips=1) 1871:1 1970:1 nile
*
nonlin psi
dlm(y=nile,method=bfgs,a=1.0,c=1.0,exact,scale,sv=1.0,sw=exp(psi),$
  vhat=vhat,svhat=svhat) 1871:1 1970:1 xstates vstates
*
* Standardize the one-step prediction error (vhat). With the variance scale factor
* concentrated out, svhat needs to be multiplied by the %variance to convert it
* to the proper variance.
*
set ehat = vhat(t)(1)/sqrt(%variance*svhat(t)(1,1))
*
* The STATISTICS instruction includes separate tests for skewness=0, excess
* kurtosis=0 and the joint test of the two (Jarque-Bera).
*
stats ehat
sstats 2 34 ehat**2>>denom
sstats 1970:1-32 1970:1 ehat**2>>numer
disp 'H(33)=' numer/denom
*
spgraph(vfields=2,hfields=2,footer="Figure 2.7 Diagnostic Plots for Standardized Prediction Errors")
graph(header="Standardized residual")
# ehat
@regcorrs(number=9,method=yule) ehat
*
* P-plot
*
order(ranks=ranks1) ehat
set ranks1 = ranks1/%nobs
stats(noprint) ehat
set testdist = %invnormal(ranks1)*sqrt(%variance)+%mean
scatter(lines=||0.0,1.0||,header="P-plot")
# ehat testdist
*
@histogram ehat
spgraph(done)