*
* 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]
*   REGIMES=# of regimes [2]
*   SWITCH=[M]/I/MH/IH/C/CH
*     In each of these, M indicates a Hamilton-type model where the mean of
*     the series is regime-dependent. I indicates the intercept is regime-
*     dependent. C indicates that the full coefficient matrix switches.
*     Those are mutually exclusive. H indicates variances are
*     regime-dependent.
*
* This creates the following variables which will generally need to be
* put into a parmset and estimated:
*
*   P      = nregimes-1 x nregimes Markov transition matrix
*   THETA  = nregimes-1 x nregimes logistic indices for transition matrix
*     (P and THETA are alternative representations)
*
*   MU     = VECTOR[VECTOR] of means/intercepts. MU(s) is the mean/intercept
*            vector at regime s
*   PHI    = VECTOR[RECT] of lag coefficients. PHI(k) is the N x N matrix of
*            coefficients for lag k
*   PHIV   = RECT[RECT] of lag coefficients. PHIV(s,k) is the N x N matrix of
*            coefficients for lag k at regime s.
*   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 regime probabilities.
*    MSVARLagRegime,MSVARTransLookup,MSVARTransProbs,MSVARIntercept and
*      PSTAR are working matrices
*
* Other procedures and functions:
*
* @MSVARInitial(guessregimes=VECT[INT] of regime guesses) start gend
*   computes initial guess values for a model already set up
*
* @MSVARParmSet(parmset=parmset to define with VAR free parameters)
*
* %MSVARMarginal(pexpand,lags)
*   returns a VECTOR of marginal probabilities for regimes or tuples of
*   regimes 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 regime, 1 gives the pairs of regime(t),regime(t-1), etc.
*
* %MSVARFVec(time)
*   returns the vector of likelihoods (*not* in log form) for the
*   (expanded) regimes 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
*
* %MSVARInitTransition()
*   for a model with fixed transitions, expands the transition matrix if
*   required and scans for negative probabilities. Returns 1 if the
*   probabilities are all non-negative, and 0 if not.
*
* @MSVARSmoothed gstart gend psmoothed
*   calculates the smoothed (marginal) probabilities of the regimes, that
*   is psmoothed(t) is the smoothed probability of being at the various
*   regimes at time t. (For models with switching means, the smoothing
*   algorithm has to be applied to the expanded regime vector).
*
* 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.
*
* For MCMC:
*
* @MSVARRelabel swaps
*   Reorders all switching components based upon the index array <<swaps>>
*
* Revision Schedule:
*   05/2007 Written by Tom Doan, Estima
*   07/2008 Revised to include calculation of the "mask" series for computing
*           probabilities of the regimes.
*   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
*   12/2009 Revise MSVAREStep to handle time-varying probabilities in smoothing
*   02/2010 Add MSVARSmoothed procedure
*   04/2011 Change order from lowest to highest on guess values. Add SWITCH=C and
*           SWITCH=CH.
*   10/2013 Add MSVARStdResiduals
*   07/2016 Update to use function pointers for transition matrices
*   08/2016 Check for missing values in range during MSVARInitial
*   01/2017 Change %MSVARInit to use SVD to permit test parameters with more than
*           one stationary point.
*   09/2017 Add GUESSREGIMES option to @MSVARInitial
*   10/2017 Add @MSVarParmSet
*   11/2017 Add NREGIMES as the preferred synonym for nregimes
*   12/2017 Add @MSVARRelabel procedure
*
*******************************************************************************
*
source mssetup.src
*
procedure MSVARSetup
option integer  lags	    1
option integer  regimes  2
option integer  states
option choice   switch   1  M I MH IH C CH
*
* Parameter vectors
*
dec vect[vect]    mu
dec vect[rect]    phi
dec rect[rect]    phiv
dec vect[symm]    sigmav
dec symm          sigma
*
* Numbers of various items
*
dec integer       nvar

dec vect[series]  y
dec model         MSVARModel
dec vect[series]  MSVARU
dec rect[integer] MSVARLagRegime
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
*
* 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] MSVARFilterLagRegime
declare integer   MSVARFilterSize
declare integer   EMCMomSize
dec series[vect]  MSVARPT_T
dec series[vect]  MSVARPT_T1
dec series[vect]  MSVARPTSM
dec rect[rect]    EMMBlock
dec vect[vect]    EMYBlock
declare vector    MSVARP0
*******************************************************************************
*
local vect[int] yvec
local integer i j
*
enter(varying) yvec
@MSSetup(regimes=%if(%defined(states),states,regimes),lags=lags)
*
compute nvar   =%rows(yvec)
*
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
   compute MSVARMeanSwitch=1
else
   compute MSVARMeanSwitch=0

if switch==3.or.switch==4.or.switch==6
   compute MSVARSigmaSwitch=1
else
   compute MSVARSigmaSwitch=0

if switch==5.or.switch==6
   compute MSVARPhiSwitch=1
else
   compute MSVARPhiSwitch=0
*
* If the means of the series are switching, the number of regimes for the
* mean of the data is nregimes^(nlags+1), since we have to cover every
* possible combination. If we're only switching coefficients in the VAR
* directly, there will only be <<nregimes>> distinct possibilities.
*
if MSVARMeanSwitch
   compute nexpand=nregimes^(nlags+1)
else
   compute nexpand=nregimes
*
* MSVARLagRegime(i,j) is a lookup table for the regime of expanded regime
* "i" at lag j-1.
*
dim MSVARLagRegime(nexpand,nlags+1)
if MSVARMeanSwitch==1
   ewise MSVARLagRegime(i,j)=1+%mod((i-1)/nregimes^(nlags+1-j),nregimes)
else
   ewise MSVARLagRegime(i,j)=i
*
* We also create a mapping from regime to regime which will pick out the
* transition probability. MSVARTransLookup will be a number from 1 to
* nregimes^2+1. 1 is the most common value, which represents a zero
* transition probability. For instance, in a two-regime setting with
* current + four lags, if we're in (1,1,2,2,1), the only regimes 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
* regime can be obtained from another, since the final nlags spots of the
* new regime have to be the same as the lead nlags spots of the old regime.
*
dim MSVARTransLookup(nexpand,nexpand) MSVARTransProbs(1+nregimes^2)
if MSVARMeanSwitch==1
   ewise MSVARTransLookup(i,j)=fix(%if(%mod(i-1,nregimes^nlags)==(j-1)/nregimes,$
       nregimes*(MSVARLagRegime(i,1)-1)+MSVARLagRegime(j,1)+1,1))
else
   ewise MSVARTransLookup(i,j)=(i-1)*nregimes+j+1
*
* Set up the model
*
system(model=MSVARModel)
variables yvec
if nlags>0
   lags 1 to nlags
det constant
end(system)
dim MSVARU(nvar)
*
* The transition probabilities are parameterized as an nregimes-1 x
* nregimes matrix. Each column gives the transition probabilities out of
* a particular regime to the first nregimes-1 regimes. The theta array is
* used for an alternative parameterization using logistics.
*
dim pstar(nexpand)
*
dim   mu(nregimes)
ewise mu(i) =%zeros(nvar,1)
*
dim   phi(nlags)
ewise phi(i)=%zeros(nvar,nvar)
*
dim   phiv(nregimes,nlags)
ewise phiv(i,j)=%zeros(nvar,nvar)
*
dim   MSVARIntercept(nexpand)
ewise MSVARIntercept(i)=%zeros(nvar,1)
*
dim   sigmav(nregimes)
dim   MSVARLogDetV(nregimes)
dim   MSVARSigmaInv(nregimes)
dim   MSVARSigmaInvF(nregimes)
end
*******************************************************************************
*
* Creates a parmset with the free parameters from the regression part
* of the model.
*
procedure MSVARParmSet
option parmset *parmset
*
if .not.%defined(parmset) {
   disp "###Usage: @MSVARParmSet(PARMSET=parmset to define)"
   return
}

if MSVARSigmaSwitch {
   if MSVARPhiSwitch
      nonlin(parmset=parmset) mu phiv sigmav
   else
      nonlin(parmset=parmset) mu phi sigmav
}
else {
   if MSVARPhiSwitch
      nonlin(parmset=parmset) mu phiv sigma
   else
      nonlin(parmset=parmset) mu phi sigma
}
*
end
*******************************************************************************
procedure MSVARRelabel swaps
type vect[integer] swaps
*
local integer i
local vect[symm] tempsigma
local rect[rect] tempphi
local vect[vect] tempmu
local rect       tempp
*
if MSVARSigmaSwitch {
   compute tempsigma=sigmav
   ewise sigmav(i)=tempsigma(swaps(i))
}
if MSVARPhiSwitch {
   compute tempphi=phiv
   ewise phiv(j,i)=tempphi(j,swaps(i))
}
if MSVARMeanSwitch {
   compute tempmu=mu
   ewise mu(i)=tempmu(swaps(i))
}
*
compute tempp=p
ewise p(i,j)=tempp(swaps(i),swaps(j))
end
*******************************************************************************
function %MSVARLookupTable nregimes nlags
type rect[integer] %MSVARLookupTable
type integer nregimes nlags
*
local integer i j
*
* Returns a lookup table(i,j) for the regime of expanded regime <<i>> at
* lag <<j>>+1.
*
dim   %MSVARLookupTable(nregimes^(nlags+1),nlags+1)
ewise %MSVARLookupTable(i,j)=1+%mod((i-1)/nregimes^(nlags+1-j),nregimes)
end
*******************************************************************************
function %MSVARMarginal pexpand lags
type     vector  %MSVARMarginal
type     vector  pexpand
type     integer lags
*
local integer targetSize perTarget
*
compute targetSize=nregimes^(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
option series[int]   guessregimes
*
local integer i j k l startl endl
local rect    crect
local real    spacing
*
* By default, this makes all means increase together.
*
do i=1,nvar
   stats(noprint) y(i) start end
   compute spacing=sqrt(%variance)*.5*nregimes/(nregimes-1)
   do k=1,nregimes
      compute mu(k)(i)=%mean+spacing*(2*k-(nregimes+1))
   end do k
end do i
*
* Estimate the coefficients, and split them out by lags.
*
estimate(model=msvarmodel,noprint) start end
compute startl=%regstart(),endl=%regend()
compute %nmiss=(endl-startl+1)-%nobs
if %nmiss>0 {
   disp "## @MSVARInitial range of " %datelabel(start) "to" %datelabel(end) "includes" %nmiss "missing observations."
   disp "## Correct the range before estimating model."
}
compute crect=%modelgetcoeffs(MSVARModel)
do i=1,nvar
   do l=1,nlags
      do j=1,nvar
         compute phi(l)(i,j)=crect((j-1)*nlags+l,i)
      end do j
   end do l
end do i
*
* Compute the coefficients to the phiv in case we need that.
*
ewise phiv(i,j)=phi(j)
*
* If we're switching intercepts, not means, back out the intercepts from
* the means.
*
if MSVARMeanSwitch==0
   ewise mu(i)=%modellagsums(MSVARModel)*mu(i)
*
* Make all the variance matrices equal
*
ewise sigmav(i)=%sigma
compute sigma=%sigma
*
* Use moderately persistent regimes
*
ewise p(i,j)=%if(i==j,.8,.2/(nregimes-1))
compute theta=%MSPLogistic(p)
*
* Override if we have guessregimes used
*
gset pt_t      = %zeros(nexpand,1)
gset pt_t1     = %zeros(nexpand,1)
gset psmooth   = %zeros(nexpand,1)
do i=1,nvar
   set MSVARU(i) = 0.0
end do i
*
if .not.%defined(guessregimes)
   return
*
* Use user-specified guess values for regimes
*
do k=1,nregimes
   if MSVARMeanSwitch {
      do i=1,nvar
         stats(noprint,smpl=guessregimes{0}==k) y(i) startl endl
         compute mu(k)(i)=%mean
      end do i
   }
   if MSVARPhiSwitch.or.MSVARSigmaSwitch {
      estimate(noprint,model=MSVARModel,smpl=guessregimes{0}==k) startl endl
      if MSVARPhiSwitch {
         compute crect=%modelgetcoeffs(MSVARModel)
         do i=1,nvar
            do l=1,nlags
               do j=1,nvar
                  compute phiv(k,l)(i,j)=crect((j-1)*nlags+l,i)
               end do j
            end do l
         end do i
      }
      if MSVARSigmaSwitch {
         compute sigmav(k)=%sigma
      }
   }
end do k
end
*********************************************************************
function  %MSVARPmat  time
type rect %MSVARPmat
type int  time
*
compute %MSVARPmat=%MSPMat(time)
end
*******************************************************************************
*
* Returns the vector of likelihoods (*not* in log form) for the
* (expanded) regimes at <<time>>.
*
function %MSVARFVec time
type vector  %MSVARFVec
type integer time
*
local integer i j regime
local vect    u
local real    iconst
*
dim %MSVARFVec(nexpand)
*
compute iconst=(2*%pi)^(-.5*nvar)
if MSVARPhiSwitch {
   *
   * Entire coefficient vector switches, so the calculation needs to be
   * done inside the loop over regimes.
   *
   do i=1,nexpand
      *
      * regime should always be = i.
      *
      compute regime=MSVARLagRegime(i,1)
      compute u=%xt(y,time)
      do j=1,nlags
         compute u=u-phiv(regime,j)*%xt(y,time-j)
      end do j
      compute %MSVARFVec(i)=iconst*exp(-.5*MSVARLogDetV(regime)-$
         .5*%qform(MSVARSigmaInv(regime),u-MSVARIntercept(regime)))
   end do i
}
else {
   *
   * Lag coefficients are fixed. The raw residuals (excepting the
   * mean/intercepts) are common to all regimes.
   *
   compute u=%xt(y,time)
   do j=1,nlags
      compute u=u-phi(j)*%xt(y,time-j)
   end do j
   *
   do i=1,nexpand
      compute regime=MSVARLagRegime(i,1),$
         %MSVARFVec(i)=iconst*exp(-.5*MSVARLogDetV(regime)-$
           .5*%qform(MSVARSigmaInv(regime),u-MSVARIntercept(i)))
   end do i
}
end
*********************************************************************
*********************************************************************
*
* Functions/procedures needed for internal calculations
*
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==nregimes
   ewise MSVARIntercept(i)=mu(i)
else {
   do i=1,nexpand
      compute MSVARIntercept(i)=mu(MSVARLagRegime(i,1))
      if MSVARPhiSwitch {
         do j=1,nlags
            compute MSVARIntercept(i)=MSVARIntercept(i)-phiv(i,j)*mu(MSVARLagRegime(i,j+1))
         end do j
      }
      else {
         do j=1,nlags
            compute MSVARIntercept(i)=MSVARIntercept(i)-phi(j)*mu(MSVARLagRegime(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,nregimes
   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))
      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
*******************************************************************************
*
* Returns the vector of residuals at <<time>> using only the
* coefficients which are fixed across regimes.
*
function %MSVARFixResids time
type vector %MSVARFixResids
type integer time
*
local vector  u
local integer j
*
compute u=%xt(y,time)
*
* If the full coefficient vector switches, just return the dependent
* variables.
*
if MSVARPhiSwitch {
}
else {
   do j=1,nlags
      compute u=u-phi(j)*%xt(y,time-j)
   end do j
}
compute %MSVARFixResids=u
end
*******************************************************************************
*
* Returns the vector of residuals at <<time>> using only the lag
* coefficients. This shouldn't be used if the lag coefficients switch.
*
function %MSVARRawResids time
type vector %MSVARRawResids
type integer time
*
local vector  u
local integer j
*
compute u=%xt(y,time)
*
* If the full coefficient vector switches, just return the dependent
* variables.
*
if MSVARPhiSwitch {
}
else {
   do j=1,nlags
      compute u=u-phi(j)*%xt(y,time-j)
   end do j
}
compute %MSVARRawResids=u
end
*******************************************************************************
*******************************************************************************
*
* Returns the residuals at <<time>> for MS VAR given actual values of
* the current and lagged regimes.
*
procedure MSVARResids v start end
type vect[series]    *v
type integer         start end
*
option series[int]   regime MSRegime
option integer       specific
*
local vector  u x
local integer j time
local rect    phix
local integer thisRegime
*
dim v(nvar)
*
do time=start,end
   if %defined(specific)
      compute thisRegime=specific
   else
      compute thisRegime=regime(time)
   *
   compute u=%xt(y,time)-mu(thisRegime)
   do j=1,nlags
      if MSVARPhiSwitch
         compute phix=phiv(thisRegime,j)
      else
         compute phix=phi(j)
      compute x=%xt(y,time-j)
      if MSVARMeanSwitch
         if %defined(specific)
            compute x=x-mu(specific)
         else
            compute x=x-mu(regime(time-j))
      compute u=u-phix*x
   end do j
   compute %pt(v,time,u)
end do time
end
*******************************************************************************
*
* Returns a VECT[SERIES] of standardized one-step residuals, with
* regimes weighted by the one-step predictive probabilities.
*
procedure MSVARStdResiduals start end resids
*
type integer start end
type vect[series]  *resids
*
option series[vect] phat pt_t1
*
local integer time k
local integer regime
local vector  ey eyk eycommon
local symm    eyy
*
dim resids(nvar)
do k=1,nvar
   set resids(k) start end = 0.0
end do k
*
compute %MSVARInitIntercepts()
do time=start,end
   if MSVARPhiSwitch {
      *
      * Entire coefficient vector switches, so the calculation needs to be
      * done inside the loop over regimes.
      *
      compute ey =%zeros(nvar,1)
      compute eyy=%zeros(nvar,nvar)
      do k=1,nexpand
         compute eyk=%zeros(nvar,1)
         compute regime=MSVARLagRegime(k,1)
         do j=1,nlags
            compute eyk=eyk+phiv(regime,j)*%xt(y,time-j)
         end do j
         compute eyk=eyk+MSVARIntercept(k)
         compute ey =ey+phat(time)(k)*eyk
         compute eyy=eyy+phat(time)(k)*(sigmav(regime)+%outerxx(eyk))
      end do k
   }
   else {
      *
      * Lag coefficients are fixed. The raw residuals (excepting the
      * mean/intercepts) are common to all regimes.
      *
      compute eycommon=%zeros(nvar,1)
      do j=1,nlags
         compute eycommon=eycommon+phi(j)*%xt(y,time-j)
      end do j
      *
      compute ey =%zeros(nvar,1)
      compute eyy=%zeros(nvar,nvar)
      do k=1,nexpand
         compute regime=MSVARLagRegime(k,1)
         compute eyk=eycommon+MSVARIntercept(k)
         compute ey =ey+phat(time)(k)*eyk
         compute eyy=eyy+phat(time)(k)*(sigmav(regime)+%outerxx(eyk))
      end do k
   }
   *
   * Convert to variance
   *
   compute eyy=eyy-%outerxx(ey)
   *
   * Convert to residuals
   *
   compute ey=%xt(y,time)-ey
   *
   * Standardize and save
   *
   compute %pt(resids,time,inv(%decomp(eyy))*ey)
end do time
end
********************************************************************************
*
* Defines the model (coefficients and covariance matrix) for a
* particular regime.
*
procedure MSVARSetModel
option integer regime 1
*
local rect crect
local int  i l j cpos
*
compute crect=%ModelGetCoeffs(MSVARModel)
do i=1,nvar
   do l=1,nlags
      do j=1,nvar
         compute crect((j-1)*nlags+l,i)=%if(MSVARPhiSwitch,phiv(regime,l)(i,j),phi(l)(i,j))
      end do j
   end do l
end do i
compute cpos=nlags*nvar+1,%MSVARInitIntercepts()
do i=1,nvar
   compute crect(cpos,i)=MSVARIntercept(regime)(i)
end do i
*
compute %ModelSetCoeffs(MSVARModel,crect)
compute %ModelSetVCV(MSVARModel,sigmav(regime))
end
*******************************************************************************
function %MSVARProb time
type real    %MSVARProb
type integer time
*
local vect    pw f fp
local real    fpt
local rect    pfull
*
* Get the full nexpand x nexpand transition matrix
*
compute pfull=%MSVARPMat(time)
*
* Get the full nexpand vector of likelihoods
*
compute f=%MSVarFVec(time)
*
* Compute the predicted probabilities
*
compute pw=pfull*pstar
*
* Save them
*
compute pt_t1(time)=pw
*
* Do the update given the likelihoods.
*
compute fp=f.*pw
compute fpt=%sum(fp)
*
* Save the filtered probabilities, both in the global pstar (used for
* the next step) and in the pt_t series (for general use).
*
compute pstar=fp*(1.0/fpt)
compute pt_t(time)=pstar
*
* Return the likelihood of the observation
*
compute %MSVARProb=fpt
end
*******************************************************************************
function  %MSVARInitTransition
type real %MSVARInitTransition
*
local integer i row col
if %rows(p)==%cols(p) {
   ewise MSVARTransProbs(i)=row=(i-2)/nregimes+1,col=%clock(i-1,nregimes),$
     %if(i==1,0.0,p(row,col))
}
else {
   ewise MSVARTransProbs(i)=row=(i-2)/nregimes+1,col=%clock(i-1,nregimes),$
     %if(i==1,0.0,%if(row==nregimes,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 %MSVARPInit
type  rect     %MSVARPInit
local integer  i j
*
dim %MSVARPInit(nexpand,nexpand)
ewise %MSVARPInit(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
end
*******************************************************************************
function    %MSVARInit
type vector %MSVARInit
*
* 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()
*
compute %MSVARInit=%MCErgodic(%MSVARPInit())
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
*******************************************************************************
procedure MSVARFilterSetup
option switch EM	0
*
* For EM, we need to do (at minimum) pairs of S(t),S(t-1). If we had to
* expand regimes anyway, we're all set. If not, we need to set up an
* nregimes^2 array of lag regimes.
*
if nexpand==nregimes.and.EM {
   compute MSVARFilterSize=nregimes^2
   compute MSVARFilterLagRegime=%MSVARLookupTable(nregimes,1)
}
else {
   compute MSVARFilterSize=nexpand
   compute MSVARFilterLagRegime=MSVARLagRegime
}
compute MSVARP0=%fill(nexpand,1,1.0/nexpand)

if MSVARPhiSwitch
   compute EMCMomSize=nvar*nlags+1
else
if MSVARSigmaSwitch
   compute EMCMomSize=nvar^2*nlags
else
   compute EMCMomSize=nvar*(nlags+1)

end
*******************************************************************************
procedure MSVAREMGeneralSetup
*
* We need to do (at minimum) pairs of S(t),S(t-1). If we had to expand
* regime anyway, we're all set. If not, we need to set up an nregimes^2
* array of lag regimes.
*
@MSVARFilterSetup(EM)
dim EMMBlock(nregimes,nregimes)
dim EMYBlock(nregimes)
*
end
*******************************************************************************
function %MSVARInitEM
type vector %MSVARInitEM
compute %MSVarInit()
compute %MSVARInitEM=MSVARP0
end
*******************************************************************************
procedure MSVARFilterInit gstart gend
type integer gstart gend
*
option vector p0
*
gset MSVARPT_T1 gstart gend = %zeros(MSVARFilterSize,1)
gset MSVARPT_T  gstart gend = %zeros(MSVARFilterSize,1)
gset MSVARPTSM  gstart gend = %zeros(MSVARFilterSize,1)
*
compute %logl=0.0
compute pstar=%MSVarInit()
if %defined(p0)
   compute pstar=p0
end
*********************************************************************
procedure MSVARFilterStep time f
type  integer time
type  vector  f
*
local integer i j
local rect    pfull
local integer targetregime sourceregime
local real    ptt ft
*
* This is the full <<nregimes>> x <<nregimes>> or <<nexpand>> x <<nexpand>>
* array of transition probabilities.
*
compute pfull=%MSVARPmat(time)
*
* Compute the probabilities for expanded state i
*
do i=1,MSVARFilterSize
   do j=1,MSVARFilterSize
      compute targetregime=MSVARFilterLagRegime(i,1)
      compute sourceregime=MSVARFilterLagRegime(j,1)
      if MSVARFilterSize<>nexpand {
         if MSVARFilterLagRegime(j,2)<>1.or.MSVARFilterLagRegime(i,2)<>sourceregime
            next
         compute ptt=pfull(targetregime,sourceregime)*pstar(sourceregime)
         compute ft =f(targetregime)
      }
      else {
         compute ptt=pfull(i,j)*pstar(j)
         compute ft=f(i)
      }
      compute MSVARPT_T1(time)(i)=MSVARPT_T1(time)(i)+ptt
      compute MSVARPT_T(time)(i) =MSVARPT_T(time)(i) +ptt*ft
   end do j
end do i
compute %logl=%logl+log(%sum(MSVARPT_T(time)))
compute MSVARPT_T1(time)=MSVARPT_T1(time)/%sum(MSVARPT_T1(time))
compute MSVARPT_T(time) =MSVARPT_T(time)/%sum(MSVARPT_T(time))
*
* Marginalize this to get the current regime
*
if MSVARFilterSize==nexpand
   compute pstar=MSVARPT_T(time)
else
   compute pstar=%sumc(%reshape(MSVARPT_T(time),MSVARFilterSize/nregimes,nregimes))
end
*********************************************************************
procedure MSVARFilter gstart gend
type integer gstart gend
*
option switch ergodic 1
*
local vect 		f
local integer 	time
*
if ergodic==0
   @MSVARFilterInit(p0=MSVARP0) gstart gend
else
   @MSVARFilterInit gstart gend

do time=gstart,gend
   compute f=%MSVARFvec(time)
   @MSVARFilterStep time f
end do time
end
*********************************************************************
*
* Compute smoothed marginal probabilities
*
procedure MSVARSmoothed gstart gend psmooth
type series[vect] *psmooth
*
local integer time
local vector  fullp
local rect    pfull
*
gset psmooth gstart gend = %zeros(nregimes,1)
*
* Start with the full sample filter (therefore smoothed) probabilities
* of the expanded regimes.
*
compute fullp=pt_t(gend)
compute psmooth(gend)=%MSVARMarginal(fullp,0);
*
do time=gend-1,gstart,-1
   *
   * Pull up the full transition matrix from time+1
   *
   compute pfull=%MSVARPmat(time+1)
   *
   * Do the smoothing update
   *
   compute fullp=%mssmoothcalc(pfull,pt_t(time),fullp,pt_t1(time+1))
   *
   * Save the marginalized probabilities
   *
   compute psmooth(time)=%MSVARMarginal(fullp,0)
end do time
end
*********************************************************************
procedure MSVARSample gstart gend regimes
type integer gstart gend
type series[int] *regimes
*
option vector  *counts
*
local rect 		  pfull pexpandx
local integer 	i j time exregime
local vect 		  pstar
*
* Draw the final (expanded) regime and determine the time T regime.
*
compute exregime=%ranbranch(MSVARPT_T(gend))
compute regimes(gend)=MSVARLagRegime(exregime,1)
*
compute pstar=%unitv(MSVARFilterSize,exregime)
do time=gend-1,gstart,-1
   compute pfull=%MSVARPmat(time+1)
   if %rows(pfull)<>MSVARFilterSize {
      dim pexpandx(MSVARFilterSize,MSVARFilterSize)
      ewise pexpandx(i,j)=%if(MSVARFilterLagRegime(i,2)==MSVARFilterLagRegime(j,1),$
         pfull(MSVARFilterLagRegime(i,1),MSVARFilterLagRegime(j,1)),0.0)
   }
   else
      compute pexpandx=pfull
   compute pstar=%mssmoothcalc(pexpandx,MSVARPT_T(time),pstar,MSVARPT_T1(time+1))
   compute exregime=%ranbranch(pstar)
   compute regimes(time)=MSVARFilterLagRegime(exregime,1)
   compute pstar=%unitv(MSVARFilterSize,exregime)
end do time
*
* If needed, pull out the pre-sample regimes
*
if MSVARFilterSize<>nregimes {
   do time=gstart-1,gstart-nlags,-1
      compute regimes(time)=MSVARFilterLagRegime(exregime,gstart-time+1)
   end do time
}
*
* If requested, compute the counts of the regimes
*
if %defined(counts) {
   dim counts(nregimes)
   do i=1,nregimes
      sstats gstart gend regimes==i>>counts(i)
   end do i
}
end
*********************************************************************
procedure MSVARPhiCMOM gstart gend xxem
type 	 integer gstart gend
type   symm    *xxem
*
option series[int] regime MSRegime
*
local vect    xem yvec xvec
local integer time lag
local rect    ifactor sxem
*
* This allows for regime-dependent sigmas
*
compute xxem=%zeros(nvar^2*nlags+1,nvar^2*nlags+1)
compute xem =%zeros(nvar*nlags,1)
*
do time=gstart,gend
   compute ifactor=MSVARSigmaInvF(regime(time))
   compute yvec=ifactor*(%xt(y,time)-mu(regime(time)))
   compute xem=%const(0.0)
   do lag=1,nlags
      compute xvec=%xt(y,time-lag)
      if MSVARMeanSwitch
         compute xvec=xvec-mu(regime(time-lag))
      compute %psubvec(xem,1+nvar*(lag-1),xvec)
   end do lag
   compute sxem=%kroneker(ifactor,tr(xem))~tr(yvec)
   compute xxem=xxem+%innerxx(sxem)
end do time
end
*********************************************************************
procedure MSVARPhiUnstack phiIn phiOut
type rect phiIn
type vect[rect] *phiOut
*
local integer lag
local rect    crect
*
compute crect=%reshape(phiIn,nvar*nvar,nlags)
do lag=1,nlags
   compute phiOut(lag)=%reshape(%xsubmat(crect,1,nvar*nvar,lag,lag),nvar,nvar)
end do lag
end
*********************************************************************
procedure MSVARMuCMOM gstart gend xxem
type integer gstart gend
type symm *xxem
*
option series[int] regime MSRegime
*
local vect[rect] phiblock(nregimes)
local vect    yvec xvec
local integer time lag s
local rect    phimat phisum
*
compute xxem=%zeros(nvar*nregimes+1,nvar*nregimes+1)
*
do time=gstart,gend
   *
   * Compute the (regime-dependent) raw residuals
   *
   compute yvec=%xt(y,time)
   do lag=1,nlags
      if MSVARPhiSwitch
         compute phimat=phiv(regime(time),lag)
      else
         compute phimat=phi(lag)
      compute yvec=yvec-phimat*%xt(y,time-lag)
   end do lag
   *
   * Figure out the phi multiplier for each regime
   *
   do s=1,nregimes
      compute phisum=%zeros(nvar,nvar)
      do lag=0,nlags
         if lag==0
            compute phimat = -1.0*%identity(nvar)
         else
         if MSVARPhiSwitch
            compute phimat = phiv(regime(time),lag)
         else
            compute phimat = phi(lag)
         if MSVARMeanSwitch==1.and.regime(time-lag)==s.or.$
            MSVARMeanSwitch==0.and.regime(time)==s
            compute phisum=phisum-phimat
      end do lag
      compute phisum=tr(phisum)*MSVARSigmaInv(regime(time))
      compute phiblock(s)=phisum
   end do s
   compute xvec=%blockglue(phiblock)~~tr(yvec)*MSVARSigmaInv(regime(time))
   compute xxem=xxem+%outerxx(xvec)
end do time
end
*********************************************************************
procedure MSVARMuUnstack muIn muOut
type rect muIn
type vect[vect] *muOut
*
local integer i
*
dim muOut(nregimes)
ewise muOut(i)=%xsubvec(muIn,(i-1)*nvar+1,i*nvar)
end
*********************************************************************
*
* Does the "E" step for EM, computing the smoothed probabilities of
* pairs of regimes at t and t-1. As a side effect, it computes the log
* likelihood.
*
procedure MSVAREStep gstart gend
type integer gstart gend
*
option switch ergodic 0
*
local integer time i j
local rect    pfull pexpandx
*
@MSVARFilter(ergodic=ergodic) gstart gend
*
* Compute smoothed probabilities
*
gset MSVARPTSM gstart gend = %zeros(MSVARFilterSize,1)
compute MSVARPTSM(gend)=MSVARPT_T(gend)
*
do time=gend-1,gstart-1,-1
   compute pfull=%MSVARPmat(time+1)
   if %rows(pfull)<>MSVARFilterSize {
      dim pexpandx(MSVARFilterSize,MSVARFilterSize)
      ewise pexpandx(i,j)=%if(MSVARFilterLagRegime(i,2)==MSVARFilterLagRegime(j,1),$
         pfull(MSVARFilterLagRegime(i,1),MSVARFilterLagRegime(j,1)),0.0)
   }
   else
      compute pexpandx=pfull
   if time>=gstart
      compute MSVARPTSM(time)=%mssmoothcalc(pexpandx,MSVARPT_T(time),MSVARPTSM(time+1),MSVARPT_T1(time+1))
   *
   *  Pre-sample smoothed probabilities
   *
   else
   if nexpand<>nregimes
      compute MSVARP0=%mssmoothcalc(pexpandx,MSVARP0,MSVARPTSM(gstart),MSVARPT_T1(gstart))
   else
      compute MSVARP0=%vectorect(MSVARPTSM(gstart),nregimes)*%fill(nregimes,1,1.0)
end do time
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    MSVARPTSMsum
local vect    qsum
local vect    sumwts
local rect    phistack
local rect    ifactor
local series  wt
local real    sumwt
local symm    rawsigma

if MSVARPhiSwitch {
   *
   * Lag coefficients switch as well.
   *
   compute sigma=%zeros(nvar,nvar)
   do s=1,nregimes
      set wt gstart gend = qhat=%MSVARMarginal(MSVARPTSM(t),0),qhat(s)
      cmom(weight=wt,model=msvarmodel) gstart gend
      sstats gstart gend wt>>sumwt
      compute phistack=%sweeptop(%cmom,EMCMomSize)
      do lag=1,nlags
         ewise phiv(s,lag)(i,j)=phistack((j-1)*nlags+lag,EMCMomSize+i)
      end do lag
      ewise mu(s)(i)=phistack(EMCMOMSize,EMCMOMSize+i)
      compute rawsigma=%xsubmat(phistack,EMCMomSize+1,EMCMOMSize+nvar,EMCMomSize+1,EMCMOMSize+nvar)
      if MSVARSigmaSwitch
         compute sigmav(s)=rawsigma/sumwt
      else
         compute sigma=sigma+rawsigma
   end do s
   if MSVARSigmaSwitch==0
      compute sigma=sigma/%nobs
   goto skipsigma
}
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(nregimes,1)
      do time=gstart,gend
         compute qhat=%MSVARMarginal(MSVARPTSM(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,nregimes
         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 MSVARLagRegime(k,lag+1)==s
                  compute phisum=phisum+phimat,notthere=0
            end do lag
            if notthere
               next
            if MSVARSigmaSwitch
               compute phisum=tr(phisum)*MSVARSigmaInv(MSVARLagRegime(k,1))
            else
               compute phisum=tr(phisum)*MSVARSigmaInv(1)
            do time=gstart,gend
               if MSVARMeanSwitch.and.nlags>0
                  compute qhat=MSVARPTSM(time)
               else
                  compute qhat=%MSVARMarginal(MSVARPTSM(time),0)
               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,MSVARLagRegime(k,lag+1))=EMMBlock(s,MSVARLagRegime(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,nregimes
         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.or.nlags==0
      branch skipvar

   if MSVARSigmaSwitch==0 {
      *
      * With sigma not being regime-dependent, this can be done simply as
      * a multivariate regression with regime-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.and.nlags>0
               compute qhat=MSVARPTSM(time)
            else
               compute qhat=%MSVARMarginal(MSVARPTSM(time),0)
            compute xem=%const(0.0)
            compute yvec=%xt(y,time)-mu(MSVARLagRegime(k,1))
            compute %psubvec(xem,1,yvec)
            do lag=1,nlags
               compute xvec=%xt(y,time-lag)
               if MSVARMeanSwitch
                  compute xvec=xvec-mu(MSVARLagRegime(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 regime-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(MSVARLagRegime(k,1))
         do time=gstart,gend
            if MSVARMeanSwitch.and.nlags>0
               compute qhat=MSVARPTSM(time)
            else
               compute qhat=%MSVARMarginal(MSVARPTSM(time),0)
            compute yvec=ifactor*(%xt(y,time)-mu(MSVARLagRegime(k,1)))
            compute xem=%const(0.0)
            do lag=1,nlags
               compute xvec=%xt(y,time-lag)
               if MSVARMeanSwitch
                  compute xvec=xvec-mu(MSVARLagRegime(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
      goto skipsigma

   if MSVARSigmaSwitch==0.and.MSVARMeanSwitch.and.dovar.and.nlags>0 {
      *
      * 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(nregimes,1)
      do time=gstart,gend
         if MSVARMeanSwitch.and.nlags>0
            compute qhat=MSVARPTSM(time)
         else
            compute qhat=%MSVARMarginal(MSVARPTSM(time),0)
         do k=1,nexpand
            compute yvec=%xt(y,time)-mu(MSVARLagRegime(k,1))
            do lag=1,nlags
               compute xvec=%xt(y,time-lag)
               if MSVARMeanSwitch
                  compute xvec=xvec-mu(MSVARLagRegime(k,lag+1))
               compute yvec=yvec-phi(lag)*xvec
            end do lag
            if MSVARSigmaSwitch {
               compute sigmav(MSVARLagRegime(k,1))=sigmav(MSVARLagRegime(k,1))+qhat(k)*%outerxx(yvec)
               compute sumwts(MSVARLagRegime(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 regimes, we need to marginalize to the
* pairs of (S(t),S(t-1))
*
if dop==0
   goto skipp

compute MSVARPTSMsum=%zeros(nregimes,nregimes)
do time=gstart,gend
   do k=1,MSVARFilterSize
      compute MSVARPTSMsum(MSVARFilterLagRegime(k,1),MSVARFilterLagRegime(k,2))+=MSVARPTSM(time)(k)
   end do k
end do time
compute qsum=%fill(1,nregimes,1.0)*MSVARPTSMsum
ewise p(i,j)=MSVARPTSMsum(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 regimes, 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 <<nregimes>> or <<nexpand>> depending upon
   * the type of switching.
   *
   compute f=%MSVARFvec(time)
   *
   * This is the full <<nregimes>> x <<nregimes>> or <<nexpand>> x
   * <<nexpand>> array of transition probabilities.
   *
   compute pfull=%MSVARPmat(time)
   *
   * Compute the probabilities for expanded regime 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
end
*********************************************************************
*
* Glues together a vector of parameters for testing convergence
*
procedure MSVARTestVector v
type   vector *v
option switch  dop       0
option switch  dotheta   0
option switch  domu      0
option switch  dovar     0
local  integer i
*
dim v(0)
if dop
   compute v=v~~%vec(p)
if dotheta
   compute v=v~~%vec(theta)
if domu {
   do i=1,nregimes
      compute v=v~~mu(i)
   end do i
}
if dovar {
   do i=1,nlags
      compute v=v~~%vec(phi(i))
   end do i
}
end

