*
* 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/2007 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)
end





