*
* Example 12.2 from pp 434-438
* Mixture model
*
open data table2_2.xls
calendar(m) 1975
data(format=xls,org=columns) 1975:1 1984:7 usukxrate
set usukxrate = usukxrate*.01
*
compute [vect] deltas=||1.0,0.9,0.8,0.7||
dec vect[series] lhats(4) phats(4)
*
* Save the information from periods 34 and 35 for table 12.1
*
dec vect mu34(4) c34(4) f35(4) q35(4)
do i=1,4
   dlm(y=usukxrate,a=1.0,c=1.0,sx0=deltas(i),sv=1.0,discount=deltas(i),$
      variance=chisquare,pdf=1.0,pscale=.01,likelihood=lhats(i),$
      sighist=sighist,yhat=yhat,svhat=svhat) / xstates vstates
   compute mu34(i)=xstates(34)(1)
   compute c34(i)=vstates(34)(1,1)*sighist(34)
   compute f35(i)=yhat(35)(1)
   compute q35(i)=svhat(35)(1,1)*sighist(34)
end dofor delta
*
dec vect pt(4) lt(4)
compute pt=%const(.25)
do time=1975:1,1984:7
   *
   *  The lhats are cumulative log likelihoods, so we have to convert these into
   *  log likelihoods for a specific time
   *
   if time==1975:1
      compute lt=%xt(lhats,time)
   else
      compute lt=%xt(lhats,time)-%xt(lhats,time-1)
   *
   *  Compute the relative weights on the four states
   *
   compute pt=pt.*%exp(lt)
   *
   *  And normalize to get probabilities
   *
   compute pt=pt/%sum(pt)
   do i=1,4
      set phats(i) time time = pt(i)
   end do i
end do time
*
spgraph(vfields=4,footer="Figure 12.2 USA/UK index: Posterior model probabilities")
do i=1,4
   graph(min=0.0,max=1.0)
   # phats(i)
end do i
spgraph(done)
*
report(action=define)
report(atrow=1,atcol=1,align=center) "j" "delta" "p34" "m34" "sqrt(C34)" "f35" "sqrt(Q35)"
do i=1,4
   report(atrow=i+1,atcol=1) i deltas(i) phats(i)(34) mu34(i) sqrt(c34(i)) f35(i) sqrt(q35(i))
end do i
report(action=format,atcol=2,tocol=2,picture="#.#")
report(action=format,atcol=3,tocol=3,picture="#.##")
report(action=format,atcol=4,picture="*.###")
report(action=show)
*
*
dec vect[series] d34(4)
set x   1 800 = -.03+(t-1)*.0001
do i=1,4
   set d34(i) 1 800 = %tdensity((x-mu34(i))/sqrt(c34(i)),35)/sqrt(c34(i))
end do i
set mixture 1 800 = %dot(%xt(d34,t),%xt(phats,34))
scatter(footer="Figure 12.3 Posteriors for mu at t=34",style=line) 5
# x d34(1) 1 800 2
# x d34(2) 1 800 2
# x d34(3) 1 800 2
# x d34(4) 1 800 2
# x mixture 1 800 1