The RATS Software Forum

There is a newer version of this at http://www.estima.com/forum/viewtopic.php?f=7&t=512

This is a revised set of procedures and functions for estimating Markov switching VAR's. This adds an option for homogeneous vs heterogeneous variance matrices and generates a "mask" matrix which maps the augmented states into the current state, so the probabilities at time t of the states can be obtained. (The %MSSMOOTH procedure can't be applied until the state probabilities have been remapped).

Note that this is written for RATS v7 or later.

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=[MEANS]/INTERCEPTS
*     SWITCH=MEANS is the Hamilton type model where the mean of the series is
*       state-dependent.
*     SWITCH=INTERCEPTS has the intercept term in each VAR equation being state
*       dependent.
*   VARIANCE=[HOMOGENEOUS]/HETEROGENEOUS
*     VARIANCE=HOMOGENEOUS is a single variance (covariance matrix).
*     VARIANCE=HETEROGENEOUS means a separate one for each state.
*
*  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
*   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
*    MASK is a RECTANGULAR used to convert expanded probabilities into current period probabilities
*    LAGSTATE,TRANSIT,TP,VARINTERCEPT and PSTAR are working matrices
*    VARHANDLING is an integer flag which indicates whether the variances switch or not.
*
* @MSVARInitial
*   computes initial guess values for a model already set up
*
* @MSVARProb(time)
*   computes the (not logged) likelihood of the model at "time" for the current set of
*   parameters
*
* @MSVARInit()
*   does the calculations needed at the start of each function evaluation
*
* 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.
*
********************************************************
procedure MSVARSetup
option integer  lags       1
option integer  states   2
option choice   switch   1 means intercepts
option choice   variance 1 homogeneous heterogeneous
*
dec rect[integer] lagstate
dec rect          p
dec rect[integer] transit
dec vect          tp
dec vect          pstar
dec vect[vect]    mu
dec vect[rect]    phi
dec vect[symm]    sigmav
dec symm          sigma
dec vect          pstar
dec series[vect]  pt_t pt_t1 psmooth
dec vect[series]  y
dec vect[vect]    varintercept
dec integer       varhandling
dec rect          mask
*
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
*
* 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 switch==1
   compute nexpand=nstates**(nlags+1)
else
   compute nexpand=nstates
*
* lagstate(i,j) is a lookup table for the state of expanded state "i" at lag j-1.
*
dim lagstate(nexpand,nlags+1)
ewise lagstate(i,j)=1+%mod((i-1)/nstates**(nlags+1-j),nstates)
*
* We also create a mapping from state to state which will pick out the transition
* probability. transit 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 transit(nexpand,nexpand) tp(1+nstates**2)
if switch==1
   ewise transit(i,j)=fix(%if(%mod(i-1,nstates**nlags)==(j-1)/nstates,nstates*(lagstate(i,1)-1)+lagstate(j,1)+1,1))
else
   ewise transit(i,j)=(i-1)*nstates+j+1
*
dim mask(nstates,nexpand)
ewise mask(i,j)=lagstate(j,1)==i
*
* 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.
*
dim p(nstates-1,nstates) pstar(nexpand)
*
dim mu(nstates) phi(nlags) sigmav(nstates) varintercept(nexpand)
ewise mu(i) =%zeros(nvar,1)
ewise phi(i)=%zeros(nvar,nvar)
ewise varintercept(i)=%zeros(nvar,1)
*
* Stuff the value of the "variance" option into the global "varhandling"
*
compute varhandling=variance
*
end
*******************************************************************************
procedure MSVARInitial
local integer i j l
local model olsmodel
local rect crect
*
do i=1,nvar
   stats(noprint,fractiles) y(i)
   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
lags 1 to nlags
det constant
end(system)
estimate(noprint)
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 nexpand==nstates
   ewise mu(i)=%modellagsums(olsmodel)*mu(i)
*
ewise sigmav(i)=%sigma
compute sigma=%sigma
*
ewise p(i,j)=%if(i==j,.8,.1/(nstates-1))
*
end
*******************************************************************************
function MSVARProb time
type real MSVARProb
type integer time
*
local integer i j
local vect u p f fp
local real sp fpt
*
dim f(nexpand) p(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 f(i)=exp(%logdensity(sigmav(lagstate(i,1)),u-varintercept(i)))
end do i
*
do i=1,nexpand
   compute sp=0.0
   do j=1,nexpand
      compute sp=sp+tp(transit(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 MSVARProbEM time
type real MSVARProbEM
type integer time
*
local integer i j
local vect u f
local real fpt
*
dim f(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 f(i)=exp(%logdensity(sigmav(lagstate(i,1)),u-varintercept(i)))
end do i
*
compute fp=f.*psmooth(time)
compute MSVARProbEM=%sum(fp)
end
*******************************************************************************
function MSVARInit
type vector MSVARInit
*
local rect a
local integer i j row col
*
ewise tp(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)))
*
* If any of the transition probabilities are negative, the parameter setting will
* be rejected.
*
if %minvalue(tp)<0.0 {
   compute MSVARInit=%zeros(nexpand,1)
   return
}
*
if nexpand==nstates
   ewise varintercept(i)=mu(i)
else {
   do i=1,nexpand
      compute varintercept(i)=mu(lagstate(i,1))
      do j=1,nlags
         compute varintercept(i)=varintercept(i)-phi(j)*mu(lagstate(i,j+1))
      end do j
   end do i
}
*
if varhandling==1
   ewise sigmav(i)=sigma
*
dim a(nexpand+1,nexpand)
ewise a(i,j)=%if(i>nexpand,1.0,(i==j)-tp(transit(i,j)))
compute MSVARInit=%xcol(inv(tr(a)*a)*tr(a),nexpand+1)
end

Code: Select all

*
* This does the univariate Hamilton model on GDP, using msvarsetup instead of the
* custom code in the example from the RATS manual.
*
cal(q) 1951:1
open data gnpdata.prn
data(format=prn,org=columns) 1951:1 1984:4
*
set g = 100*log(gnp/gnp{1})
source msvarsetup.src
source markov.src
frml msvarf = log(MSVARProb(t))
@msvarsetup(lags=4)
# g
compute gstart=1952:2,gend=1984:4
*
gset pt_t  gstart gend = %zeros(nexpand,1)
gset pt_t1 gstart gend = %zeros(nexpand,1)
nonlin(parmset=msparms) p
nonlin(parmset=varparms) mu phi sigma
@msvarinitial
maximize(trace,parmset=varparms+msparms,start=(pstar=MSVARInit()),$
  reject=%minvalue(tp)<0.0,pmethod=simplex,piters=1,method=bfgs,iters=300) msvarf gstart gend
gset pt_t  gstart gend = mask*pt_t
gset pt_t1 gstart gend = mask*pt_t1
@%mssmooth p pt_t pt_t1 psmooth
*
* To create the shading marking the recessions, create a dummy series which is 1
* when the recessq series is 1, and 0 otherwise. (recessq is 1 for NBER recessions
* and -1 for expansions).
*
set contract = recessq==1
*
set prob2 = 1-psmooth(t)(1)
graph(style=polygon,header="Probability of Economy Being in Contraction",shade=contract)
# prob2 gstart gend

I have a number of queries with respect to the revised MS VAR procedure:

1) Can RATS 6.3 be used to execute the procedure?

I tried to run the code which does the univariate Hamilton model on GDP, using MSVARSETUP but I get the following error messages:

## SX22. Expected Type REAL, Got VECTOR(VECTOR) Instead
>>>>rmset=varparms) mu <<<<


2) As far as I have understood, the SWITCH option allows for either the MEAN or the INTERCEPT to switch. Would it be possible to estimate a specification where also the VAR autoregressive coefficients are regime-dependent?

3) The variance/covariance matrix: when "HETEROGENEOUS" is chosen for the VARIANCE option, the variances are allowed to be regime-dependent. What about the covariances? Is the var/cov matrix diagonal?

4) I am not clear about the meaning of the "mask" matrix

5) Does the MSVARSETUP allow to calculate impulse response functions? If yes, could you indicate how? If no, do you have any suggestion about how to implement that?

Thank you very much in advance for your kind help!

Dear Doan,
Thanks for MSVAR code allowing for switching mean and intercept. But as far as I know determining number of regimes is very tricky in this kind of models and often subject to harsh criticisim by referees.
Krolzig and the other authors having MSVAR code and programs are not clear about this. How can we determine number of states in this code? is there any formal test for or do we have to calculate some information criterias like AIC or BIC? is it possible to conduct some diagnostics on the residuals of MSVAR model by which allows us to reach a decision on the true model.
Thanks again.
Nazif

At least in the case of a VAR, it's highly unlikely that an information criterion will be helpful. The extra state pulls in so many extra parameters that it's almost impossible for the improved fit to be worth the "cost". If two states are adequate, a three state model has a huge number of redundant parameters. The sheer impossibility of getting the expanded model to converge is likely to be a good sign that it isn't necessary.

Dear Doan
Thanks for your answer.
I know It's to much but I have another. In case of unit root do we have to first difference of the variables when we're running this code? I've read lots of paper but couldn't get the clear answer.
Thanks again

Regime switching behaviour itself might introduce unit root. Performance of the regular unit root tests are not promising in the case of regime switching. You can find a nice working paper on the subject here: http://research.stlouisfed.org/wp/2001/2001-013.pdf. I think, you can first build your model as is, and then look for tracks of unit root afterwards. If you find any evidence on unit root, then you can try the usual remedies.

nazif wrote:Dear Doan
Thanks for your answer.
I know It's to much but I have another. In case of unit root do we have to first difference of the variables when we're running this code? I've read lots of paper but couldn't get the clear answer.
Thanks again

Because the VAR part is all conditional on past y's, there's nothing about this that is sensitive to unit roots. The biggest question would be the treatment of any exogenous variables, since the "constant" switches its meaning when you have unit roots.