*
* Replication from Kim, Chang-Jin & Nelson, Charles R, 2001. "A Bayesian
* Approach to Testing for Markov-Switching in Univariate and Dynamic
* Factor Models,"  International Economic Review, vol. 42(4), 989-1013.
*
* Stock-Watson common factor model without Markov switching
*
* Example with missing observables
*
cal(m) 1952:1
open data fulldta.prn
data(format=free,org=cols) 1 517 z data1 data2 data3 data4
*
* Convert data to growth rates
*
set lndata1 = log(data1)*100
set lndata2 = log(data2)*100
set lndata3 = log(data3)*100
set lndata4 = log(data4)*100
set y1      = lndata1-lndata1{1}
set y2      = lndata2-lndata2{1}
set y3      = lndata3-lndata3{1}
set y4      = lndata4-lndata4{1}
*
* For example, eliminate part of the y4 series
*
set y4 * 1964:12 = %na
*
* Estimates start in 1959:2
*
compute rstart=1959:2,rend=1995:1
*
* This controls the random number generator to allow full replicability
* (within RATS runs).
*
seed 198011
*
* For flexibility, put the observables into a VECT[SERIES]
*
compute n=4
dec vect[series] y(n)
diff(center) y1 / y(1)
diff(center) y2 / y(2)
diff(center) y3 / y(3)
diff(center) y4 / y(4)
*
* Gibbs sampling controls
* (This probably needs a lot more draws than this to get accurate
* results. If you have patience, 5000 and 25000 would be better).
*
compute nburn =1000
compute ndraws=9000
*
* Number of lags in the AR on the common cycle component. The variance
* is pegged to one to identify the scale of the cycle.
*
compute philags=2
dec vect phi
dec real sigsqv
compute sigsqv=1.0
compute phi=%zeros(philags,1)
*
* Number of lags on the loadings onto the cycle for each series.
* (Current is always included).
*
compute [vect[int]] glags=||0,0,0,3||
dec vect[vect] gamma(n)
ewise gamma(i)=%zeros(glags(i)+1,1)
*
* Number of lags in the noise term for each series.
*
compute [vect[int]] psilags=||2,2,2,2||
dec vect[vect] psi(n)
dec vect       sigsqe(n)
ewise psi(i)=%zeros(psilags(i),1)
***********************************************************************
*
* No user input until next ******
*
* Figure out how many lags of the cycle component we'll need in the
* state space representation. This depends upon the phi lags and upon
* the number needed for the loadings.
*
compute clags=philags
do i=1,n
   compute clags=%imax(clags,glags(i)+1)
end do i
*
* Initialize the fixed elements in the A submatrix for the cycle
* component.
*
dec rect ac(clags,clags)
ewise ac(i,j)=(i==j+1)
*
* Set the F matrix for the cycle component
*
dec rect fc(clags,1)
ewise fc(i,j)=(i==1)
*
* Create and zero out the C submatrix for the cycle component. The
* non-zero elements will come from the "gammas".
*
dec rect cc
compute cc=%zeros(clags,n)
*
* Initialize the fixed elements in the A submatrix, set the C submatrix
* for the individual noise terms. The F submatrix is the same as C.
*
dec rect ae(0,0) ce(0,0) apsi
do i=1,n
   dim apsi(psilags(i),psilags(i))
   ewise apsi(i,j)=(i==j+1)
   compute ae=ae~\apsi
   compute ce=ce~\%unitv(psilags(i),1)
end do i
*
* The F matrix for the system is now known
*
compute f=fc~\ce
*
* This is a procedure to put together when needed the working A, C and SW
* matrices given the current settings.
*
function DLMSetup
local integer i j
compute %psubmat(ac,1,1,tr(phi))
compute j=1
do i=1,n
   compute %psubmat(ae,j,j,tr(psi(i)))
   compute j=j+psilags(i)
   compute %psubmat(cc,1,i,gamma(i))
end do i
*
compute a=ac~\ae
compute c=cc~~ce
compute sw=%diag(sigsqv~~sigsqe)
end
*
* This is used for the patched over y data.
*
dec vect[series] yhat(n)
*
**********************************************************************
**********************************************************************
*
* Prior parameters
*
**********************************************************************
*
* For the psi's
*
dec vect 		nupsi(n)
dec vect 		s2psi(n)
dec vect[vect] bpsi(n)
dec vect[symm]	hpsi(n)
*
* Prior for precision
*
ewise nupsi(i)=4.0
ewise s2psi(i)=4.0
*
* Zero mean with fairly low precision on the coefficients
*
ewise bpsi(i)=%zeros(psilags(i),1)
ewise hpsi(i)=%identity(psilags(i))*1.0
**********************************************************************
*
* For the phi
*
dec vect      bphi
dec symm      hphi
*
* Zero mean with fairly low precision on the coefficients
*
compute bphi=%zeros(philags,1)
compute hphi=%identity(philags)*1.0
**********************************************************************
*
* For the gammas
*
dec vect[vect] bgamma(n)
dec vect[symm] hgamma(n)
ewise bgamma(i)=%zeros(glags(i)+1,1)
ewise hgamma(i)=%identity(glags(i)+1)*1.0
*
* Initial values for Gibbs sampler. Positive loadings on the gamma's can
* help to fix the sign on the cycle component.
*
compute sigsqe=%const(.2)
do i=1,n
   compute gamma(i)(1)=.5
end do i
*
**********************************************************************
nonlin(parmset=fixparms) phi psi gamma sigsqe
dec series[vect] bgibbs
gset bgibbs 1 ndraws = %zeros(%rows(%parmspeek(fixparms)),1)
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampling"
do draw=-nburn,ndraws
   *
   * Draw the state space components (cycle and e's)
   *
   compute DLMSetup()
   dlm(type=csimulate,presample=ergodic,a=a,c=c,f=f,sw=sw,y=%xt(y,t)) rstart rend xstates
   *
   * Put the conditionally simulated values of y into yhat. These will
   * match wherever y exists.
   *
   do i=1,n
      set yhat(i) rstart rend = [vector] ypatch=tr(c)*xstates,ypatch(i)
   end do i
   *
   * Draw the psi's and sigsqe's given the e's
   *
   compute pos=clags+1
   do i=1,n
      set ei rstart rend = xstates(t)(pos)
      compute pos=pos+psilags(i)
      cmom
      # ei{1 to psilags(i)} ei
      compute psi(i)=%ranmvpostcmom(%cmom,1.0/sigsqe(i),hpsi(i),bpsi(i))
   	compute rssplus=nupsi(i)*s2psi(i)+%rsscmom(%cmom,bpsi(i))
   	compute sigsqe(i)=rssplus/%ranchisqr(nupsi(i)+%nobs)
   end do i
   *
   * Draw the phi given the c's
   *
   set cycle rstart rend = xstates(t)(1)
   cmom
   # cycle{1 to philags} cycle
   compute phi=%ranmvpostcmom(%cmom,1.0/sigsqv,hphi,bphi)
   *
   * Draw the gammas. This uses GLS on the (patched) data and cycles.
   *
	do i=1,n
      filter(weights=1.0~~-1.0*psi(i),type=lagging) yhat(i)  / fy
      filter(weights=1.0~~-1.0*psi(i),type=lagging) cycle / fcycle
      cmom
      # fcycle{0 to glags(i)} fy
   	compute gamma(i)=%ranmvpostcmom(%cmom,1.0/sigsqe(i),hgamma(i),bgamma(i))
   end do i
   infobox(current=draw)
   if draw<=0
      next
   *
   *   Do the bookkeeping here.
   *
   compute bgibbs(draw)=%parmspeek(fixparms)
end do draw
infobox(action=remove)
@mcmcpostproc(ndraws=ndraws,mean=bmean,stderrs=bstderrs) bgibbs
report(action=define,title="Gibbs Sampling Estimates")
report(atrow=1,atcol=1,align=center) "Mean" "StdErr"
report(atrow=2,atcol=1,fillby=cols) bmean
report(atrow=2,atcol=2,fillby=cols) bstderrs
report(atcol=1,tocol=2,action=format,picture="*.###")
report(action=show)
*
* Maximum likelihood
* This works exactly as before --- DLM will adjust automatically for the
* missing observables.
*
dlm(method=bfgs,parmset=fixparms,start=DLMSetup(),presample=ergodic,$
  a=a,c=c,f=f,sw=sw,y=%xt(y,t)) rstart rend xstates