(A newer version of this was posted on 16 July 2008 under a different thread).
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.
*
* 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)
*
* 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
* LAGSTATE,TRANSIT,TP,VARINTERCEPT and PSTAR are working matrices
*
*
* @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/07 Written by Tom Doan, Estima
*
********************************************************
procedure MSVARSetup
option integer lags 1
option integer states 2
option choice switch 1 means intercepts
*
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
*
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
*
* 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)
*
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 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
}
*
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)
endExample
Code: Select all
source msvarsetup.src
frml msvarf = log(MSVARProb(t))
*
open data wbc.xls
calendar(q) 1960
*
data(format=xls,org=columns) 1960:01 1991:04 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}
*
@msvarsetup(lags=1)
# dusa djap dfrg duk dcan daus
*
* With the mu's and phi's, you have to list the elements explicitly when defining
* the parmset. (NONLIN won't accept a VECT[VECT] or VECT[RECT]).
*
nonlin(parmset=msparms) p
nonlin(parmset=varparms) mu(1) mu(2) phi(1) sigma
@msvarinitial
compute gstart=1962:1,gend=1991:4
gset pt_t gstart gend = %zeros(nexpand,1)
gset pt_t1 gstart gend = %zeros(nexpand,1)
compute t1=%cputime()
maximize(parmset=varparms+msparms,start=(pstar=MSVARInit()),reject=%minvalue(tp)<0.0,method=bfgs,iters=400) msvarf gstart gend
disp %cputime()-t1