*
* KIMNP126.RPF
* Application 3 from pp 126-133 of Kim and Nelson, "State-space Models
* with Regime Switching".
*
* Dynamic factor model with Markov switching component
*
open data s&w.xls
cal(m) 1959:1
data(format=xls,org=columns) 1959:01 1995:01 ip gmyxpq mtq lpnag dcoinc
*
set ipgrow    = 100.0*log(ip/ip{1})
set incgrow   = 100.0*log(gmyxpq/gmyxpq{1})
set salesgrow = 100.0*log(mtq/mtq{1})
set empgrow   = 100.0*log(lpnag/lpnag{1})
*
* Parameters for index
*
dec vect phi(2)
compute phi=%zeros(2,1)
dec real sigmac
compute sigmac=1.0
*
* Parameters for the series. sigma is initialized based upon the sample
* means.
*
dec vect beta(4)
dec vect[vect] gamma(4)
dec vect sigma(4)
dec vect[vect] psi(4)
ewise psi(i)=%zeros(2,1)
*
stats(noprint) ipgrow
compute beta(1)=%mean,sigma(1)=%variance*.5
compute gamma(1)=%fill(1,1,0.6)
stats(noprint) incgrow
compute beta(2)=%mean,sigma(2)=%variance*.5
compute gamma(2)=%fill(1,1,0.5)
stats(noprint) salesgrow
compute beta(3)=%mean,sigma(3)=%variance*.5
compute gamma(3)=%fill(1,1,0.4)
stats(noprint) empgrow
compute beta(4)=%mean,sigma(4)=%variance*.5
compute gamma(4)=%fill(4,1,0.2)
*
* Construct fixed parts of the A matrix. We need four lags for the cycle
* (first block) because of the required lags for the loadings on the
* employment series.
*
dec rect ar4(4,4)
dec rect ar2(2,2)
ewise ar4(i,j)=(i==j+1)
ewise ar2(i,j)=(i==j+1)
compute a=ar4~\ar2~\ar2~\ar2~\ar2
compute ndlm=%rows(a)
*
* These are the positions in the state vector of the noise terms for the
* series.
*
dec vect[integer] spos(4)
compute spos=||5,7,9,11||
*
* Construct the F matrix (always fixed). Again, we need a size 4 block
* for the cycle.
*
dec rect f4(4,1)
dec rect f2(2,1)
compute f4=%unitv(4,1)
compute f2=%unitv(2,1)
compute f=f4~\f2~\f2~\f2~\f2
*
* Construct the fixed parts of the C matrix
*
compute c=%zeros(4,4)~~(f2~\f2~\f2~\f2)
*
dec symm sw(5,5)
*********************************************************************
*
* Setup for Kim filter
*
*********************************************************************
*
* This is common to any MS state space model analyzed using the
* Kim filter.
*
@MSSetup(states=2)
*
* These are the collapsed SSM mean and variance from the previous time
* period.
*
dec vect[vect] xstar(nstates)
dec vect[symm] sstar(nstates)
*
* These are the work SSM mean and variance for the current time period
*
dec rect[vect] xwork(nstates,nstates)
dec rect[symm] swork(nstates,nstates)
*
dec rect       fwork(nstates,nstates)
*
dec series[vect] xstates
*
* Transition probabilities (in logistic form)
*
compute theta=%zeros(nstates-1,nstates)

*********************************************************************
*
* Set up the switching component (the "z" shift in the transition
* equation).
*
dec vect[vect] zswitch(nstates)
dec vect mu(nstates)
compute zswitch(1)=%zeros(ndlm,1)
compute zswitch(2)=%zeros(ndlm,1)
*
compute mu(1)=-1.0,mu(2)=1.0
*********************************************************************
function DLMStart
*
* Insert the lag coefficients for the cycle
*
local symm sx0
*
compute %psubmat(a,1,1,tr(phi))
*
* Insert the lag coefficients for the noise terms into the A matrix and
* the loading coefficients into the C matrix.
*
do i=1,4
   compute %psubmat(a,spos(i),spos(i),tr(psi(i)))
   compute %psubmat(c,1,i,gamma(i))
end do i
compute sw=f*%diag(sigmac~~sigma)*tr(f)
*
compute zswitch(1)(1)=mu(1)
compute zswitch(2)(1)=mu(2)
*
* Kim filter initialization
*
* Transform the theta to transition probabilities.
*
compute p=%MSLogisticP(theta)
compute pstar=%mcergodic(p)
compute p=%mspexpand(p)
*
compute sx0=%psdinit(a,sw)
ewise xstar(i)=%solve(%identity(ndlm)-a,zswitch(i))
ewise sstar(i)=sx0
end DLMStart
*
* The observables are the growth rates minus their estimated means. The
* paper uses the sample means rather than estimating these, so we leave
* beta out of the parameter set.
*
diff(standardize) ipgrow    / cip
diff(standardize) incgrow   / cinc
diff(standardize) salesgrow / csales
diff(standardize) empgrow   / cemp
*
* The observables are the growth rates minus their estimated means. The
* paper uses the sample means rather than estimating these, so we leave
* beta out of the parameter set.
*
dec frml[vect] yf
frml yf = ||ipgrow-beta(1),incgrow-beta(2),$
            salesgrow-beta(3),empgrow-beta(4)||
*********************************************************************
*
* This does a single step of the Kim (approximate) filter
*
function KimFilter time
type integer time
*
local integer i j
local vect    yerr
local real    likely
local symm    vhat
local rect    gain
local rect    phat(nstates,nstates)
local rect    fwork(nstates,nstates)
*
* Any "switch" code will be in this loop
*
do i=1,nstates
   do j=1,nstates
      *
      * Do the SSM predictive step. In this model, the shift is in the
      * "Z" component.
      *
      compute xwork(i,j)=a*xstar(j)+zswitch(i)
      compute swork(i,j)=a*sstar(j)*tr(a)+sw
      *
      * Do the prediction error for y under state i, and compute the
      * density function for the prediction error.
      *
      compute yerr=yf(time)-tr(c)*xwork(i,j)
      compute vhat=tr(c)*swork(i,j)*c
      compute gain=swork(i,j)*c*inv(vhat)
      compute fwork(i,j)=exp(%logdensity(vhat,yerr))
      *
      * Do the SSM update step
      *
      compute xwork(i,j)=xwork(i,j)+gain*yerr
      compute swork(i,j)=swork(i,j)-gain*tr(c)*swork(i,j)
   end do j
end do i
*
* From here until the end of the function is the same for all models
*
* Compute the Hamilton filter likelihood
*
compute likely=0.0
do i=1,nstates
   do j=1,nstates
      compute phat(i,j)=p(i,j)*pstar(j)*fwork(i,j)
      compute likely=likely+phat(i,j)
   end do j
end do i
*
* Compute the updated probabilities of the state combinations
*
compute phat=phat/likely
compute pstar=%sumr(phat)
compute pt_t(time)=pstar
*
* Collapse the SSM matrices down to one per state
*
compute xstates(time)=%zeros(ndlm,1)
do i=1,nstates
   compute xstar(i)=%zeros(ndlm,1)
   do j=1,nstates
      compute xstar(i)=xstar(i)+xwork(i,j)*phat(i,j)/pstar(i)
   end do j
   compute xstates(time)=xstates(time)+xstar(i)*pstar(i)
   compute sstar(i)=%zeros(ndlm,ndlm)
   do j=1,nstates
      compute sstar(i)=sstar(i)+phat(i,j)/pstar(i)*$
         (swork(i,j)+%outerxx(xstar(i)-xwork(i,j)))
   end do j
end do i
compute KimFilter=likely
end
*
* Estimate the parameters by (approximate) maximum likelihood. This has
* multiple modes, one of which has a substantially higher likelihood
* value (once corrected for integrating constants) than the one cited,
* but with very odd parameter values. This might indicate that the
* approximation has failed to work adequately.
*
compute mu(1)=-1.0,mu(2)=.25
nonlin(parmset=dlmparms) gamma psi sigma phi mu
nonlin(parmset=msparms)  theta
compute theta(1,1)=log(.8/.2),theta(1,2)=log(.2/.8)
*
frml kimf = kf=KimFilter(t),log(kf)
*********************************************************************
gset xstates = %zeros(ndlm,1)
@MSFilterInit
maximize(parmset=dlmparms+msparms,start=DLMStart(),method=bfgs) kimf 1959:2 *