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