*
* Example from section 10.2, pp 318-329
*
open data table10_1.xls
calendar(q) 1973
data(format=xls,org=columns) 1973:1 1984:4 sales
*
graph(footer="Figure 10.1 Agricultural sales series")
# sales
*
@SeasonalDLM(type=fourier) as cs sws
@LocalDLM(type=trend) al cl swl
*
compute a=al~\as
compute c=cl~~cs
compute discounts=||.85,.85,.97,.97,.97||
dlm(y=sales,a=a,c=c,$
  exact,discount=discounts,sv=1.0,var=conditional,$
  vhat=vhat,yhat=yhat,svhat=svhat,sighist=variances,dfhist=df) / xstates vstates
*
set onestep = yhat(t)(1)
set stderr  = sqrt(svhat(t)(1,1)*variances{1})
set upper   = onestep+%invttest(.10,df{1})*stderr
set lower   = onestep-%invttest(.10,df{1})*stderr
graph(footer="Figure 10.2 Agricultural sales with one-step ahead forecasts",overlay=dots,ovsame) 4
# onestep 7 *
# lower   7 * 2
# upper   7 * 2
# sales
*
set onesteperror = vhat(t)(1)/stderr
set upper =  %invttest(.10,df{1})
set lower = -%invttest(.10,df{1})
graph(footer="Figure 10.3 Standardized one-step ahead forecast errors") 3
# onesteperror 7 *
# lower 7 * 2
# upper 7 * 2
*
compute cx=||1.0,0.0,0.0,0.0,0.0||
set trend  = %dot(cx,xstates(t))
set vtrend = %qform(vstates(t),cx)*variances
set upper  = trend+%invttest(.10,df)*sqrt(vtrend)
set lower  = trend-%invttest(.10,df)*sqrt(vtrend)
graph(footer="Figure 10.4 On-line estimated trend in sales",overlay=dots,ovsame) 4
# trend 7 *
# lower 7 *
# upper 7 *
# sales
*
compute cx=||0.0,0.0,1.0,0.0,1.0||
set seasonal = %dot(cx,xstates(t))
set vseasonal = %qform(vstates(t),cx)*variances
set upper = seasonal+%invttest(.10,df)*sqrt(vseasonal)
set lower = seasonal-%invttest(.10,df)*sqrt(vseasonal)
graph(style=vertical,footer="Figure 10.5 On-line estimated seasonal pattern in sales") 2
# lower 7 *
# upper 7 *
*
*  Compute Bayes factors for the discounted vs fixed models
*
dlm(y=sales,a=a,c=c,sv=1.0,exact,var=chisquare,likely=lxd,yhat=yhatd,discount=discounts)
dlm(y=sales,a=a,c=c,sv=1.0,exact,var=chisquare,likely=lxs,yhat=yhats)
set factors = lxd-lxs
set logweights = factors-factors{1}
*
spgraph(footer="Figure 10.7 Log weights and Bayes' factors plotted over time",vfields=2)
graph(header="Log Weights")
# logweights
graph(header="Log Factors")
# factors
spgraph(done)
*
set predictd = yhatd(t)(1)
set predicts = yhats(t)(1)
*
@uforeerrors(title="Dynamic Model") sales predictd
@uforeerrors(title="Static Model") sales predicts