*
* Example 12.6
* SWARCH.PRG
* Revision Schedule:
* 11/2004 Altered to save filtered probabilities and to compute smoothed
* probabilities.
*
all 6237
open data g10xrate.xls
data(format=xls,org=columns) / usxjpn
*
* Convert to percent daily returns
*
set x = 100.0*log(usxjpn/usxjpn{1})
*
* These set the number of states, and the number of ARCH lags
*
source markov.src
compute nstates=3
compute q=2
*
* P is the matrix of transition probabilities
*
dec rect p(nstates-1,nstates)
*
* HV will be the relative variances in the states. This will be normalized to a
* relative variance of 1 in the first state, so it's dimensioned N-1.
*
dec vect hv(nstates-1)
*
* This will be the vector of ARCH parameters
*
dec vect a(q)
*
* We have three parts of the parameter set: the mean equation parameters (here,
* just an intercept), the ARCH model parameters, and the Markov switching parameters.
*
nonlin(parmset=meanparms) mu
nonlin(parmset=archparms) a0 a
nonlin(parmset=msparms) p hv
*
dec vector pstar
*
dec series[vect] pt_t pt_t1 psmooth
gset pt_t    = %zeros(nstates,1)
gset pt_t1   = %zeros(nstates,1)
gset psmooth = %zeros(nstates,1)
*
* uu and u are used for the series of squared residuals and the series
* of residuals.
*
clear uu u
*
* ARCHStateF returns a vector of likelihoods for the various states at
* time given residual e. The likelihoods differ in the states based upon
* the values of hv, where the intercept in the ARCH equation is scaled up
* by hv. Again, the elements of hv are offset by 1 since they're normalized
* with the first state at 1.
*
function ARCHStateF time e
type vector ARCHStateF
type real e
type integer time
*
local integer i j
local real vi
dim ARCHStateF(nstates)
do i=1,nstates
   compute vi=a0*%if(i>1,hv(i-1),1)
   do j=1,q
      compute vi=vi+a(j)*uu(time-j)
   end do i
   compute ARCHStateF(i)=%if(vi>0,%density(e/sqrt(vi))/sqrt(vi),0.0)
end do i
end
*
* As is typically the case with Markov switching models, there is no global
* identification of the states. By defining a fairly wide spread for hv, we'll
* hope that we'll stay in the zone of the likelihood where state 1 is low variance,
* state 2 is medium and state 3 is high.
*
compute hv=||10,100||
compute p=||.8,.2,.0|.2,.6,.4||
stats x
compute mu=%mean,a0=%variance,a=%const(0.05)
set uu = %variance
*
* The log likelihood function recursively calculates the vector pstar of estimated
* state probabilities. The first step at each t is to compute the mean and update
* the u and uu series. The likelihoods in the states are computed by the ARCHStateF
* function, and then %mcstate and %msupdate update the pstar vector and compute the
* likelihood. We use a start formula on maximize to start pstar at the ergodic
* probabilities.
*
* Because we need 2 lags of uu, the estimation starts at 3.
*
frml logl = u(t)=(x-mu),uu(t)=u(t)**2,f=ARCHStateF(t,u(t)),$
   pt_t1=%mcstate(p,pstar),pt_t=pstar=%msupdate(f,pt_t1,fpt),log(fpt)
frml init = (pstar=%mcergodic(p)),0
maximize(start=init,method=bhhh,iters=100,pmethod=simplex,piters=5,parmset=meanparms+archparms+msparms) logl 3 *
*
* Expand p to a full nxn matrix
*
dec rect pexpand(nstates,nstates-1)
ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
compute pfull=pexpand*p
ewise pfull(i,j)=%if(i==nstates,1.0+pfull(i,j),pfull(i,j))
*
gset psmooth 1 %regend() = pstar
*
do time=%regend()-1,3,-1
   gset psmooth time time = pt_t.*(tr(pfull)*(psmooth{-1}./pt_t1{-1}))
end do time