This is written for version 7 of RATS. It requires the new version of MARKOV.SRC (separate thread) as well.
Code: Select all
*
* Procedure file for Markov switching VAR's. This includes several functions, as
* described below.
*
* @MSVARSetup( options )
* # list of dependent variables
*
* Options:
* LAGS=# of VAR lags [1]
* STATES=# of states [2]
* SWITCH=[M]/I/MH/IH
* In each of these, M indicates a Hamilton-type model where the mean of the
* series is state-dependent. I indicates the intercept is state-dependent.
* Those are mutually exclusive. H indicates variances are state-dependent.
*
* This creates the following variables which will generally need to be put into a parmset and
* estimated:
*
* P = nstates-1 x nstates Markov transition matrix
* THETA = nstates-1 x nstates logistic indices for transition matrix
* (P and THETA are alternative representations)
*
* MU = VECTOR[VECTOR] of means/intercepts. MU(j) is the mean/intercept vector at state j
* PHI = VECTOR[RECT] of lag coefficients. PHI(j) is the N x N vector of lag
* coefficients for lag j
* SIGMA = SYMMETRIC (covariance matrix) or
* SIGMAV = VECTOR[SYMMETRIC] for heterogeneous covariance matrices
*
* It also sets up the following:
* Y = VECT[SERIES] of dependent variables
* PT_T,PT_T1 and PSMOOTH are SERIES[VECT] used for computing and saving the state probabilities
* MSVARLagState,MSVARTransLookup,MSVARTransProbs,MSVARIntercept and PSTAR are working matrices
*
* Other procedures and functions:
*
* @MSVARInitial start gend
* computes initial guess values for a model already set up
*
* %MSVARMarginal(pexpand,lags)
* returns a VECTOR of marginal probabilities for states or tuples of states given an expanded
* vector of probabilities (pexpand). lags is the number of lags in the tuple. lags=0 is used
* to get just the current state, 1 gives the pairs of state(t),state(t-1), etc.
*
* %MSVARFVec(time)
* returns the vector of likelihoods (*not* in log form) for the (expanded) states at <<time>>.
*
* %MSVARProb(time)
* compute the likelihood (not logged) of the model at <<time>> for the current set of
* parameters. As a side effect, it computes pt_t1(time) and pt_t(time)
*
* %MSVARInit()
* does the calculations needed at the start of each function evaluation
*
* For EM algorithm:
*
* @MSVAREMGeneralSetup
* needs to be called before using EM to set up the work arrays
*
* @MSVAREMStep gstart gend
* @MSVAREStep gstart gend
* @MSVARMStep gstart gend
* EMStep does the combined E and M steps, EStep does only the E, and MStep does only
* the M. EMStep and MStep both have options to restrict the set of parameters that
* are updated, so some can be fixed.
*
* Revision Schedule:
* 05/2007 Written by Tom Doan, Estima
* 07/2008 Revised to include calculation of the "mask" series for computing
* probabilities of the states.
* 09/2008 Add <<start>> and <<end>> parameters to MSVARInitial to control range
* 09/2008 Add logic to allow nlags==0 in MSVARInitial
* 01/2009 Split up MSVARInit into separate tasks
* 03/2009 Write EM code
*
*******************************************************************************
*
source markov.src
*
procedure MSVARSetup
option integer lags 1
option integer states 2
option choice switch 1 M I MH IH
*
* Parameter vectors
*
dec rect p
dec rect theta
dec vect[vect] mu
dec vect[rect] phi
dec vect[symm] sigmav
dec symm sigma
*
* Numbers of various items
*
dec integer nvar
dec integer nstates
dec integer nlags
dec integer nexpand
dec vect pstar
dec series[vect] pt_t
dec series[vect] pt_t1
dec series[vect] psmooth
dec vect[series] y
dec rect[integer] MSVARLagState
dec rect[integer] MSVARTransLookup
dec vect MSVARTransProbs
dec vect[vect] MSVARIntercept
dec vect MSVARLogDetV
dec vect[symm] MSVARSigmaInv
dec vect[rect] MSVARSigmaInvF
dec integer MSVARMeanSwitch
dec integer MSVARSigmaSwitch
dec integer MSVARPhiSwitch
*
local vect[int] yvec
local integer i j
*
enter(varying) yvec
*
compute nvar =%rows(yvec)
compute nstates=states
compute nlags =lags
*
dim y(nvar)
do i=1,nvar
set y(i) = yvec(i){0}
end do i
*
* Parse out the choices for the switch option
*
if switch==1.or.switch==3.or.switch==5.or.switch==7
compute MSVARMeanSwitch=1
else
compute MSVARMeanSwitch=0
if switch==3.or.switch==4.or.switch==7.or.switch==8
compute MSVARSigmaSwitch=1
else
compute MSVARSigmaSwitch=0
if switch==5.or.switch==6.or.switch==7.or.switch==8
compute MSVARPhiSwitch=1
else
compute MSVARPhiSwitch=0
*
* If the means of the series are switching, the number of states for the mean of
* the data is nstates**(nlags+1), since we have to cover every possible
* combination. If we're only switching coefficients in the VAR directly, there
* will only be nstates distinct possibilities.
*
if MSVARMeanSwitch
compute nexpand=nstates**(nlags+1)
else
compute nexpand=nstates
*
* MSVARLagState(i,j) is a lookup table for the state of expanded state "i" at lag j-1.
*
dim MSVARLagState(nexpand,nlags+1)
if MSVARMeanSwitch==1
ewise MSVARLagState(i,j)=1+%mod((i-1)/nstates**(nlags+1-j),nstates)
else
ewise MSVARLagState(i,j)=i
*
* We also create a mapping from state to state which will pick out the transition
* probability. MSVARTransLookup will be a number from 1 to nstates**2+1. 1 is the most
* common value, which represents a zero transition probability. For instance, in a
* two state setting with current + four lags, if we're in (1,1,2,2,1), the only
* states which can be obtained are {1,1,1,2,2} and {2,1,1,2,2}. The transition
* probability to the first of these would be p(1,1), and for the second it would be
* 1-p(1,1). Given our coding, it's easy to tell whether a state can be obtained from
* another, since the final nlags spots of the new state have to be the same as the
* lead nlags spots of the old state.
*
dim MSVARTransLookup(nexpand,nexpand) MSVARTransProbs(1+nstates**2)
if MSVARMeanSwitch==1
ewise MSVARTransLookup(i,j)=fix(%if(%mod(i-1,nstates**nlags)==(j-1)/nstates,$
nstates*(MSVARLagState(i,1)-1)+MSVARLagState(j,1)+1,1))
else
ewise MSVARTransLookup(i,j)=(i-1)*nstates+j+1
*
* The transition probabilities are parameterized as an nstates-1 x nstates matrix.
* Each column gives the transition probabilities out of a particular state to the
* first nstates-1 states. The theta array is used for an alternative
* parameterization using logistics.
*
dim p(nstates-1,nstates) theta(nstates-1,nstates) pstar(nexpand)
*
dim mu(nstates)
ewise mu(i) =%zeros(nvar,1)
*
dim phi(nlags)
ewise phi(i)=%zeros(nvar,nvar)
*
dim MSVARIntercept(nexpand)
ewise MSVARIntercept(i)=%zeros(nvar,1)
*
dim sigmav(nstates)
dim MSVARLogDetV(nstates)
dim MSVARSigmaInv(nstates)
dim MSVARSigmaInvF(nstates)
end
*******************************************************************************
function %MSVARLookupTable nstates nlags
type rect[integer] %MSVARLookupTable
type integer nstates nlags
*
local integer i j
*
* Returns a lookup table(i,j) for the state of expanded state <<i>> at lag <<j>>+1
*
dim %MSVARLookupTable(nstates^(nlags+1),nlags+1)
ewise %MSVARLookupTable(i,j)=1+%mod((i-1)/nstates**(nlags+1-j),nstates)
end
*******************************************************************************
function %MSVARMarginal pexpand lags
type vector %MSVARMarginal
type vector pexpand
type integer lags
*
local integer targetSize perTarget
*
compute targetSize=nstates^(lags+1)
dim %MSVARMarginal(targetSize)
compute perTarget=%rows(pexpand)/targetSize
compute %MSVARMarginal=%fill(1,perTarget,1.0)*%vectorect(pexpand,perTarget)
end
*******************************************************************************
*
* Computes a standard set of initial guess values for the parameters
*
procedure MSVARInitial start end
type integer start end
*
local integer i j l
local model olsmodel
local rect crect
*
do i=1,nvar
stats(noprint,fractiles) y(i) start end
if nstates==2 {
compute mu(1)(i)=%fract90
compute mu(2)(i)=%fract10
}
else
if nstates==3 {
compute mu(1)(i)=%fract90
compute mu(2)(i)=%median
compute mu(3)(i)=%fract10
}
end do i
*
system(model=olsmodel)
variables y
if nlags>0
lags 1 to nlags
det constant
end(system)
estimate(noprint) start end
compute crect=%modelgetcoeffs(olsmodel)
do i=1,nvar
do l=1,nlags
do j=1,nvar
ewise phi(l)(i,j)=crect((j-1)*nlags+l,i)
end do j
end do l
end do i
*
if MSVARMeanSwitch==0
ewise mu(i)=%modellagsums(olsmodel)*mu(i)
*
ewise sigmav(i)=%sigma
compute sigma=%sigma
*
ewise p(i,j)=%if(i==j,.8,.2/(nstates-1))
compute theta=%MSPLogistic(p)
*
end
*******************************************************************************
*
* Returns the vector of likelihoods (*not* in log form) for the (expanded)
* states at <<time>>.
*
function %MSVARFVec time
type vector %MSVARFVec
type integer time
*
local integer i j state
local vect u
*
dim %MSVARFVec(nexpand)
*
compute u=%xt(y,time)
do j=1,nlags
compute u=u-phi(j)*%xt(y,time-j)
end do i
*
do i=1,nexpand
compute state=MSVARLagState(i,1),$
%MSVARFVec(i)=exp(-.5*MSVARLogDetV(state)-.5*%qform(MSVARSigmaInv(state),u-MSVARIntercept(i)))
end do i
end
*******************************************************************************
function %MSVARProb time
type real %MSVARProb
type integer time
*
local integer i j
local vect p f fp
local real sp fpt
*
dim p(nexpand)
*
compute f=%MSVarFVec(time)
*
do i=1,nexpand
compute sp=0.0
do j=1,nexpand
compute sp=sp+MSVARTransProbs(MSVARTransLookup(i,j))*pstar(j)
end do j
compute p(i)=sp
end do i
*
compute pt_t1(time)=p
compute fp=f.*p
compute fpt=%sum(fp)
compute pstar=fp*(1.0/fpt)
compute pt_t(time)=pstar
compute %MSVARProb=fpt
end
*******************************************************************************
function %MSVARInitIntercepts
local integer i j
*
* Initializes the MSVARIntercept array given the current settings for the lag
* coefficients and mu array. This has no usable return value.
*
if nexpand==nstates
ewise MSVARIntercept(i)=mu(i)
else {
do i=1,nexpand
compute MSVARIntercept(i)=mu(MSVARLagState(i,1))
do j=1,nlags
compute MSVARIntercept(i)=MSVARIntercept(i)-phi(j)*mu(MSVARLagState(i,j+1))
end do j
end do i
}
end
*******************************************************************************
function %MSVARInitVariances
*
* Initializes the sigmav array given the current settings for the sigma's
* This has no usable return value
*
local integer i
local vector eigval
local rect eigvect
do i=1,nstates
if MSVARSigmaSwitch==0 {
compute sigmav(i)=sigma
if i>1 {
compute MSVARLogDetV(i) =MSVARLogDetV(i-1)
compute MSVARSigmaInv(i) =MSVARSigmaInv(i-1)
compute MSVARSigmaInvF(i)=MSVARSigmaInvF(i-1)
next
}
}
compute MSVARLogDetV(i)=0.0
eigen sigmav(i) eigval eigvect
do j=1,nvar
if eigval(j)<1.e-10
compute eigval(j)=1.e-10
compute MSVARLogDetV(i)+=(log(eigval(j))+1.837877066409345)
compute eigval(j)=1.0/eigval(j)
end do j
compute MSVARSigmaInv(i) =eigvect*%diag(eigval)*tr(eigvect)
compute MSVARSigmaInvF(i)=%diag(%sqrt(eigval))*tr(eigvect)
end do i
end
*******************************************************************************
function %MSVARInitTransition
type real %MSVARInitTransition
*
local integer i row col
ewise MSVARTransProbs(i)=row=(i-2)/nstates+1,col=%clock(i-1,nstates),$
%if(i==1,0.0,%if(row==nstates,1-%sum(%xcol(p,col)),p(row,col)))
*
* Return "false" (real zero) if any of the transition probabilities are negative
*
ewise MSVARTransProbs(i)=%if(abs(MSVARTransProbs(i))<.00001,0.0,MSVARTransProbs(i))
compute %MSVARInitTransition=(%minvalue(MSVARTransProbs)>=0.0)
end
*******************************************************************************
function %MSVARInit
type vector %MSVARInit
*
local rect a
local integer i j
*
* If any of the transition probabilities are negative, the parameter setting will
* be rejected.
*
if %MSVARInitTransition()==0.0 {
compute %MSVARInit=%zeros(nexpand,1)
return
}
compute %MSVARInitIntercepts()
compute %MSVARInitVariances()
*
dim a(nexpand+1,nexpand)
ewise a(i,j)=%if(i>nexpand,1.0,(i==j)-MSVARTransProbs(MSVARTransLookup(i,j)))
compute %MSVARInit=%xcol(inv(tr(a)*a)*tr(a),nexpand+1)
end
*******************************************************************************
procedure MSVARStaticForecast start end yfore
type integer start end
type vect[series] *yfore
*
option series[vect] phat pt_t1
*
local integer time l k
local vector yvec
*
dim yfore(nvar)
do k=1,nvar
set yfore(k) start end = 0.0
end do k
*
compute %MSVARInitIntercepts()
do time=start,end
compute yvec=%zeros(nvar,1)
do l=1,nlags
compute yvec=yvec+phi(l)*%xt(y,time-1)
end do l
do k=1,nexpand
compute yvec=yvec+phat(time)(k)*MSVARIntercept(k)
end do k
compute %pt(yfore,time,yvec)
end do time
end
*******************************************************************************
*
* EM requires smoothed probabilities of all combinations of {S(t),S(t-1)}. When
* the means are switching, we need to compute the smoothed probabilities of all
* combinations of {S(t),...,S(t-p)} anyway, so we can just compute the above
* probabilities by marginalizing.
*
declare rect[int] EMLagState
declare integer EMsize
declare integer EMCMomSize
dec series[vect] EMPT_T
dec series[vect] EMPT_T1
dec series[vect] EMPTSM
dec rect[rect] EMMBlock
dec vect[vect] EMYBlock
declare vector EMP0
*******************************************************************************
procedure MSVAREMGeneralSetup
*
* We need to do (at minimum) pairs of S(t),S(t-1). If we had to expand states
* anyway, we're all set. If not, we need to set up an nstates^2 array of lag
* states.
*
if nexpand==nstates {
compute EMsize=nstates^2
compute EMLagState=%MSVARLookupTable(nstates,1)
}
else {
compute EMsize=nexpand
compute EMLagState=MSVARLagState
}
dim EMMBlock(nstates,nstates)
dim EMYBlock(nstates)
*
if MSVARSigmaSwitch
compute EMCMomSize=nvar^2*nlags
else
compute EMCMomSize=nvar*(nlags+1)
*
compute EMP0=%fill(nexpand,1,1.0/nexpand)
end
*******************************************************************************
function %MSVARInitEM
type vector %MSVARInitEM
compute %MSVarInit()
compute %MSVARInitEM=EMP0
end
*******************************************************************************
function %MSVARPmat time
type rect %MSVARPmat
type int time
*
local integer i j
local rect pexpand
*
if nexpand==nstates {
dim pexpand(nstates,nstates-1)
ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
compute %MSVARPmat=pexpand*p
ewise %MSVARPmat(i,j)=%if(i==nstates,1.0+%MSVARPmat(i,j),%MSVARPmat(i,j))
}
else {
dim %MSVARPmat(nexpand,nexpand)
ewise %MSVARPmat(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
}
end
*******************************************************************************
*
* Does the "E" step for EM, computing the smoothed probabilities of pairs of
* states at t and t-1.
*
procedure MSVAREStep gstart gend
type integer gstart gend
*
option switch ergodic 0
*
local integer time i j
local rect pfull pexpandx
local vect pstar f
local integer targetstate sourcestate
local real ptt ft
*
***************************************************
* EM calculations for transitions
***************************************************
*
gset EMPT_T1 gstart gend = %zeros(EMSize,1)
gset EMPT_T gstart gend = %zeros(EMSize,1)
gset EMPTSM gstart gend = %zeros(EMSize,1)
*
compute %logl=0.0
if ergodic==0
compute pstar=%MSVARInitEM()
else
compute pstar=%MSVARInit()
compute EMP0=pstar
do time=gstart,gend
*
* This is either dimension <<nstates>> or <<nexpand>> depending upon
* the type of switching.
*
compute f=%MSVARFvec(time)
*
* This is the full <<nstates>> x <<nstates>> or <<nexpand>> x <<nexpand>>
* array of transition probabilities.
*
compute pfull=%MSVARPmat(time)
*
* Compute the probabilities for expanded state i
*
do i=1,EMSize
do j=1,EMSize
compute targetState=EMLagState(i,1)
compute sourceState=EMLagState(j,1)
if nstates==nexpand {
if EMLagState(j,2)<>1.or.EMLagState(i,2)<>sourceState
next
compute ptt=pfull(targetState,sourceState)*pstar(sourceState)
compute ft =f(targetState)
}
else {
compute ptt=pfull(i,j)*pstar(j)
compute ft=f(i)
}
compute EMPT_T1(time)(i)=EMPT_T1(time)(i)+ptt
compute EMPT_T(time)(i) =EMPT_T(time)(i) +ptt*ft
end do j
end do i
compute %logl=%logl+log(%sum(EMPT_T(time)))
compute EMPT_T1(time)=EMPT_T1(time)/%sum(EMPT_T1(time))
compute EMPT_T(time) =EMPT_T(time)/%sum(EMPT_T(time))
*
* Marginalize this to get the current state
*
if nstates==nexpand
compute pstar=%MSVARMarginal(EMPT_T(time),0)
else
compute pstar=EMPT_T(time)
end do time
*
* Compute smoothed probabilities
*
gset EMPTSM gstart gend = %zeros(EMSize,1)
*
compute EMPTSM(gend)=EMPT_T(gend)
compute pfull=%MSVARPmat(gend)
if %rows(pfull)<>EMSize {
dim pexpandx(EMSize,EMSize)
ewise pexpandx(i,j)=%if(EMLagState(i,2)==EMLagState(j,1),pfull(EMLagState(i,1),EMLagState(j,1)),0.0)
}
else
compute pexpandx=pfull
*
do time=gend-1,gstart,-1
compute EMPTSM(time)=%mssmoothcalc(pexpandx,EMPT_T(time),EMPTSM(time+1),EMPT_T1(time+1))
end do time
if nexpand<>nstates
compute EMP0=%mssmoothcalc(pexpandx,EMP0,EMPTSM(gstart),EMPT_T1(gstart))
else
compute EMP0=%vectorect(EMPTSM(gstart),nstates)*%fill(nstates,1,1.0)
end
*******************************************************************************
*
* Does the "M" step of EM, updating the various sets of coefficients.
*
procedure MSVARMStep gstart gend
*
option switch domu 1
option switch dovar 1
option switch dosigma 1
option switch dop 1
*
local integer time i j k s lag
local integer notthere
local rect phimat phisum
local vect qhat mublock
local rect xxem testxx
local rect sxem
local vect xem yvec xvec xyem
local rect EMPTSMsum
local vect qsum
local vect sumwts
local rect phistack
local rect ifactor
if MSVARPhiSwitch {
}
else {
*
* Compute the means or intercepts given the previous VAR coefficients
*
if domu==0
branch skipmu
if MSVARMeanSwitch==0 {
*
* Switching intercepts is easy because there's no multiplicative relationship
* between the VAR coefficients and the intercepts. It's just a weighted average
* of means of the residuals from the VAR
*
ewise EMYBlock(i)=%zeros(nvar,1)
compute sumwts=%zeros(nstates,1)
do time=gstart,gend
compute qhat=%MSVARMarginal(EMPTSM(time),0)
compute yvec=%xt(y,time)
do lag=1,nlags
compute yvec=yvec-phi(lag)*%xt(y,time-lag)
end do lag
ewise EMYBlock(s)=EMYBlock(s)+qhat(s)*yvec
ewise sumwts(s)=sumwts(s)+qhat(s)
end do time
ewise mu(s)=EMYBlock(s)/sumwts(s)
}
else {
ewise EMYBlock(i)=%zeros(nvar,1)
ewise EMMBlock(i,j)=%zeros(nvar,nvar)
*
do s=1,nstates
do k=1,nexpand
compute phisum=%zeros(nvar,nvar)
compute notthere=1
do lag=0,nlags
if lag==0
compute phimat = -1.0*%identity(nvar)
else
compute phimat = phi(lag)
if MSVARLagstate(k,lag+1)==s
compute phisum=phisum+phimat,notthere=0
end do lag
if notthere
next
if MSVARSigmaSwitch
compute phisum=tr(phisum)*MSVARSigmaInv(MSVARLagState(k,1))
else
compute phisum=tr(phisum)*MSVARSigmaInv(1)
do time=gstart,gend
compute qhat=EMPTSM(time)
do lag=0,nlags
if lag==0
compute phimat = -1.0*%identity(nvar)
else
compute phimat = phi(lag)
compute EMYBlock(s)=EMYBlock(s)+qhat(k)*phisum*phimat*%xt(y,time-lag)
compute EMMBlock(s,MSVARLagState(k,lag+1))=EMMBlock(s,MSVARLagState(k,lag+1))+qhat(k)*phisum*phimat
end do lag
end do time
end do k
end do s
compute mublock=%solve(%BLOCKGLUE(EMMBlock),%BLOCKGLUE(EMYBlock))
do i=1,nstates
compute mu(i)=%xsubvec(mublock,nvar*(i-1)+1,nvar*i)
end do i
}
:skipmu
*
* Compute VAR lag coefficients given the new means or intercepts
*
if dovar==0
branch skipvar
if MSVARSigmaSwitch==0 {
*
* With sigma not being state-dependent, this can be done simply
* as a multivariate regression with state-dependent weights on
* the observations.
*
compute xxem=%zeros(EMCMomSize,EMCMomSize)
compute xem =%zeros(EMCMomSize,1)
*
do k=1,nexpand
do time=gstart,gend
if MSVARMeanSwitch
compute qhat=EMPTSM(time)
else
compute qhat=%MSVARMarginal(EMPTSM(time),0)
compute xem=%const(0.0)
compute yvec=%xt(y,time)-mu(MSVARLagState(k,1))
compute %psubvec(xem,1,yvec)
do lag=1,nlags
compute xvec=%xt(y,time-lag)
if MSVARMeanSwitch
compute xvec=xvec-mu(MSVARLagState(k,lag+1))
compute %psubvec(xem,1+nvar*lag,xvec)
end do lag
compute xxem=xxem+qhat(k)*%outerxx(xem)
end do time
end do k
compute testxx=%sweeplist(xxem,%seq(nvar+1,EMCMOMSize))
do lag=1,nlags
compute phi(lag)=tr(%xsubmat(testxx,1+nvar*lag,nvar*(lag+1),1,nvar))
end do lag
}
else {
*
* With state-dependent sigmas, this is a more time-consuming weighted
* SUR.
*
compute xxem=%zeros(EMCMomSize,EMCMomSize)
compute xyem=%zeros(EMCMomSize,1)
compute xem =%zeros(nvar*nlags,1)
*
do k=1,nexpand
compute ifactor=MSVARSigmaInvF(MSVARLagState(k,1))
do time=gstart,gend
if MSVARMeanSwitch
compute qhat=EMPTSM(time)
else
compute qhat=%MSVARMarginal(EMPTSM(time),0)
compute yvec=ifactor*(%xt(y,time)-mu(MSVARLagState(k,1)))
compute xem=%const(0.0)
do lag=1,nlags
compute xvec=%xt(y,time-lag)
if MSVARMeanSwitch
compute xvec=xvec-mu(MSVARLagState(k,lag+1))
compute %psubvec(xem,1+nvar*(lag-1),xvec)
end do lag
compute sxem=%kroneker(ifactor,tr(xem))
compute xxem=xxem+qhat(k)*%innerxx(sxem)
compute xyem=xyem+qhat(k)*tr(sxem)*yvec
end do time
end do k
compute phistack=tr(%vectorect(%solve(xxem,xyem),nvar*nlags))
do lag=1,nlags
compute phi(lag)=%xsubmat(phistack,1,nvar,1+nvar*(lag-1),nvar*lag)
end do lag
}
:skipvar
if dosigma==0
branch skipsigma
if MSVARSigmaSwitch==0.and.MSVARMeanSwitch.and.dovar {
*
* We can do this as a side-effect of the multivariate regression
*
compute sigma=%xsubmat(testxx,1,nvar,1,nvar)/(gend-gstart+1)
}
else {
if MSVARSigmaSwitch
ewise sigmav(i)=%zeros(nvar,nvar)
else
compute sigma=%zeros(nvar,nvar)
compute sumwts=%zeros(nstates,1)
do time=gstart,gend
if MSVARMeanSwitch
compute qhat=EMPTSM(time)
else
compute qhat=%MSVARMarginal(EMPTSM(time),0)
do k=1,nexpand
compute yvec=%xt(y,time)-mu(MSVARLagState(k,1))
do lag=1,nlags
compute xvec=%xt(y,time-lag)
if MSVARMeanSwitch
compute xvec=xvec-mu(MSVARLagState(k,lag+1))
compute yvec=yvec-phi(lag)*xvec
end do lag
if MSVARSigmaSwitch {
compute sigmav(MSVARLagState(k,1))=sigmav(MSVARLagState(k,1))+qhat(k)*%outerxx(yvec)
compute sumwts(MSVARLagState(k,1))+=qhat(k)
}
else
compute sigma=sigma+qhat(k)*%outerxx(yvec)
end do k
end do time
if MSVARSigmaSwitch
ewise sigmav(i)=sigmav(i)/sumwts(i)
else
compute sigma=sigma/(gend-gstart+1)
}
}
:skipsigma
*
* Update P. For an expanded set of states, we need to marginalize to
* the pairs of (S(t),S(t-1))
*
if dop==0
branch skipp
compute EMPTSMsum=%zeros(nstates,nstates)
do time=gstart,gend
do k=1,EMSize
compute EMPTSMsum(EMLagState(k,1),EMLagState(k,2))+=EMPTSM(time)(k)
end do k
end do time
compute qsum=%fill(1,nstates,1.0)*EMPTSMsum
ewise p(i,j)=EMPTSMsum(i,j)/qsum(j)
:skipp
end MSVAREMStep
*******************************************************************************
*
* MSVAREMStep does a single step (both "E" and "M") for a Markov switching VAR
* The "M" step in most cases will be a generalized "M" step, is the sense that
* it doesn't maximize the likelihood given the smoothed probabilities for the
* states, but just increases it.
*
procedure MSVAREMStep gstart gend
*
option switch ergodic 0
option switch domu 1
option switch dovar 1
option switch dosigma 1
option switch dop 1
*
@MSVAREStep(ergodic=ergodic) gstart gend
@MSVARMStep(domu=domu,dovar=dovar,dosigma=dosigma,dop=dop) gstart gend
end
*******************************************************************************
*
* Does a function evaluation for the EM algorithm. The function value is in %logl
*
procedure MSVAREMEval gstart gend
type integer gstart gend
*
option switch ergodic 0
*
local vector f pstar fp
local rect pfull
local integer time
local real sumfp
*
compute %logl=0.0
if ergodic==0
compute pstar=%MSVARInitEM()
else
compute pstar=%MSVARInit()
do time=gstart,gend
*
* This is either dimension <<nstates>> or <<nexpand>> depending upon
* the type of switching.
*
compute f=%MSVARFvec(time)
*
* This is the full <<nstates>> x <<nstates>> or <<nexpand>> x <<nexpand>>
* array of transition probabilities.
*
compute pfull=%MSVARPmat(time)
*
* Compute the probabilities for expanded state i
*
compute fp=%zeros(nexpand,1)
do i=1,nexpand
do j=1,nexpand
compute fp(i)=fp(i)+f(i)*pfull(i,j)*pstar(j)
end do j
end do i
compute sumfp=%sum(fp)
compute pstar=fp/sumfp
compute %logl=%logl+log(sumfp)
end do time
endCode: Select all
*
* MS(2)-VAR(1) model
*
source msvarsetup.src
open data wbc.xls
calendar(q) 1960
*
data(format=xls,org=columns) 1960:1 1991:4 usa can uk frg jap aus d73q3
set dusa = usa-usa{1}
set djap = jap-jap{1}
set dfrg = frg-frg{1}
set duk = uk-uk{1}
set dcan = can-can{1}
set daus = aus-aus{1}
*
compute gstart=1962:1,gend=1991:4
*
@msvarsetup(lags=1,states=2)
# dusa djap dfrg duk dcan daus
*******************************************************************************
*
* EM iterations from the standard guess values. This hits a stationary point with
* a lower log likelihood. This is basically decomposing the sample into recession
* and expansion states.
*
@msvarinitial gstart gend
@msvarEMgeneralsetup
do emits=1,50
@msvaremstep gstart gend
disp "Iteration" emits "Log Likelihood" %logl
end do emits
disp "mu(1)=" mu(1)
disp "mu(2)=" mu(2)
disp "P=" p
*
set smoothp1 gstart gend = pstar=%MSVARMarginal(EMPTSM(t),0),pstar(1)
graph(footer="Recession/Expansion Mode")
# smoothp1
*******************************************************************************
*
* These are Krolzig's values for the mu's. From here, we pretty much reproduce his
* results until the p matrix produces an absorbing state.
*
@msvarinitial gstart gend
input mu(1)
1.0951 2.2973 1.1473 0.8920 1.4705 1.4651
input mu(2)
0.4749 0.9422 0.4920 0.3717 0.6735 0.6097
@msvarEMgeneralsetup
do emits=1,50
@msvaremstep gstart gend
disp "Iteration" emits "Log Likelihood" %logl
end do emits
*****************************************************************************
*
* These are almost imperceptibly different from the ones above, but seem to avoid
* the absorbing state, and produce a higher log likelihood
*
@msvarinitial gstart gend
input mu(1)
1.03093 2.81999 1.23774 1.01214 1.37031 1.57119
input mu(2)
0.58775 1.00435 0.57165 0.42210 0.82514 0.71478
@msvarEMgeneralsetup
do emits=1,50
@msvaremstep gstart gend
disp "Iteration" emits "Log Likelihood" %logl
end do emits
disp "mu(1)=" mu(1)
disp "mu(2)=" mu(2)
disp "P=" p
set smoothp1 gstart gend = pstar=%MSVARMarginal(EMPTSM(t),0),pstar(1)
graph(footer="End of Golden Age Mode")
# smoothp1
*****************************************************************************
*
* Full sample ML
* With some effort (it nearly gets stuck at a lower log L), this converges to same
* basic model as just estimated
*
nonlin(parmset=varparms) mu phi sigma
nonlin(parmset=msparms) theta
gset pt_t gstart gend = %zeros(nexpand,1)
gset pt_t1 gstart gend = %zeros(nexpand,1)
frml msvarf = log(%MSVARProb(t))
@msvarinitial gstart gend
maximize(trace,parmset=varparms+msparms,start=%(p=%mslogisticp(theta),pstar=%MSVARInit()),$
reject=%minvalue(MSVARTransProbs)<0.0,method=bfgs,iters=400) msvarf gstart gend
do emits=1,30
@msvaremstep gstart gend
disp "Iteration" emits "Log Likelihood" %logl
end do emits