*
* Procedure file for Markov Switching univariate linear regressions,
* with either the full coefficient vector switching or part of the
* coefficient vector switching and part fixed. The residual variance can
* be either switching or fixed. Use the MSVARSetup procedures for
* switching autoregressions and VAR's, and MSSysReg for other types of
* multivariate regressions.
*
* @MSRegression( options ) depvar
* # list of regressors (if you don't use the EQUATION option)
*
* Options:
*   REGIMES=# of regimes [2]
*   SWITCH=[C]/CH/H
*    C means that coefficients are switching, H indicates variances are switching.
*   EQUATION=equation for the regressors
*   NFIX=# of explanatory variables which are fixed across regimes [0 or all]
*    If SWITCH=C or CH, the default is 0, if SWITCH=H, it's all. The fixed
*    coefficients must be placed first in the supplementary card for regressors.
*
* 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)
*
*   BETAS  = VECT[VECT] of switching equation coefficients. BETAS(s) is the
*            vector of coefficients for regime <<s>>.
*   GAMMAS = VECTOR of fixed equation coefficients.
*   SIGSQ  = Variance (for homogeneous variances)
*   SIGSQV = VECTOR of heterogeneous variances
*
*  It also sets up the following:
*    YUNI = copy of the dependent variable
*    PT_T,PT_T1 and PSMOOTH are SERIES[VECT] used for computing and saving
*      the state probabilities.
*    PSTAR are working matrices
*
* @MSRegInitial start end
*   generates standard guess values for the coefficients and the MS
*   parameters.
*
* %MSRegFVec(time)
*   returns the vector of likelihoods (*not* in log form) for the
*   regimes at <<time>>.
*
* %MSRegProb(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)
*
* %MSRegInit()
*   does the calculations needed at the start of each function evaluation
*
* %MSRegInitTransition()
*   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.
*
* @MSRegSmoothed 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.
*
* @MSRegResids(regime=SERIES of regime values/specific=specific regime) resids start end
*   computes the equation residuals for specific regimes. Note that these are *not*
*   useful for diagnostics---they're only used as part of larger calculations.
*
* @MSRegStdResiduals resids start end
*   computes a series of standardized one-step residuals, with regimes
*   weighted by the one-step predictive probabilities. These *are* useful
*   for diagnostics, as they should be serially uncorrelated, mean zero
*   and variance one if the model is correct. (They are not, however,
*   Normally distributed even if the assumption of underlying Normality is
*   correct).
*
* @MSRegSetEquation(regime=which regime,define=equation to define)
*   Defines the equation (coefficients and variance) for a particular regime.
*
* @MSRegParmset(PARMSET=parmset to define)
*   Defines a parmset with the free parameters for the switching
*   regression model (all except the transition probability parameters).
*
* For EM algorithm:
*
* @MSRegEMGeneralSetup
*   needs to be called before using EM to set up the work arrays
*
* @MSRegEMStep gstart gend
* @MSRegEStep gstart gend
* @MSRegMStep 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 algorithm:
*
* @MSRegRelabel swaps
*   reorders all switching components based upon the index array swaps.
*
* Revision Schedule:
*   04/2011 Written by Tom Doan, Estima, based upon MSVARSETUP procedures.
*   05/2012 Revised to allow for time-varying probabilities and added
*           MSRegSetEquation to save a regime-specific version of
*   10/2015 Clean up %MSRegProb to use %MSRegPMat rather than fixed P
*   08/2016 Add DEFINE option to MSRegSetEquation
*   10/2016 Update to use function pointers
*   01/2017 Update to use SVD calculation for ergodic regime probabilities
*   04/2017 Add guessregimes option to @MSRegInitial
*   09/2017 Correct calculation of MSRegStdResiduals (used (t), not (time) for
*           one branch).
*   10/2017 Add @MSRegParmSet procedure
*   12/2017 Add REGIMES as a (preferred) synonym for STATES
*   12/2017 Add @MSRegRelabel
*******************************************************************************
*
* Pull in the required lower level procedures
*
source msemsetupstd.src
*
procedure MSRegression depvar
type series depvar
*
option integer  regimes  2
option integer  states
option choice   switch   1  C CH H
option equation equation
option integer  nfix
*
compute MSDebug=0
*
* Parameter vectors
*
dec vect[vect]    betas
dec vect          gammas
dec vect          sigsqv
dec real          sigsq
*
* Numbers of various items
*
dec integer       nvar
dec integer       nreg
*
dec series        yuni
dec vect          MSRegLogDetV
dec vect[symm]    MSRegSigmaInv
dec vect[rect]    MSRegSigmaInvF
dec integer       MSRegSigmaSwitch
dec integer       MSRegNFix
dec integer       MSRegNSwitch
dec equation      MSRegEqn
dec equation      MSRegEqnSwitch
dec equation      MSRegEqnFix
dec series        MSRegU
*
local integer     i
local rect[int]   ptable etable
*
compute nregimes=%if(%defined(states),states,regimes)
compute nstates=nregimes
@MSEMSetupStd(regimes=nregimes)
if %defined(equation) {
   equation MSRegEqn MSRegU
   # %rlfromtable(%eqntable(equation))
}
else
   equation MSRegEqn MSRegU
*
compute nvar=1
set yuni = depvar{0}
*
* Parse out the choices for the switch option
*
if switch==2.or.switch==3
   compute MSRegSigmaSwitch=1
else
   compute MSRegSigmaSwitch=0
*
* The transition probabilities are parameterized as an nregimes-1 x
* nregimes matrix. Each column gives the transition probabilities out of
* a particular state to the first nregimes-1 regimes. The theta array is
* used for an alternative parameterization using logistics.
*
dim p(nregimes-1,nregimes) theta(nregimes-1,nregimes) pstar(nregimes)
*
compute nreg=%eqnsize(MSRegEqn)
if .not.%defined(nfix) {
   if switch==1.or.switch==2
      compute MSRegNswitch=nreg
   else
      compute MSRegNswitch=0
}
else {
   compute MSRegNSwitch=%imax(0,nreg-nfix)
}
*
compute MSRegNFix=nreg-MSRegNSwitch
compute etable=%eqntable(MSRegEqn)
dim ptable(2,MSRegNSwitch)
ewise ptable(i,j)=etable(i,j+MSRegNFix)
equation MSRegEqnSwitch MSRegU
# %rlfromtable(ptable)
dim ptable(2,MSRegNFix)
ewise ptable(i,j)=etable(i,j)
equation MSRegEqnFix MSRegU
# %rlfromtable(ptable)
*
dim   betas(nregimes)
ewise betas(i) =%zeros(MSRegNSwitch,1)
compute gammas =%zeros(MSRegNFix,1)
*
dim   sigsqv(nregimes)
dim   MSRegLogDetV(nregimes)
dim   MSRegSigmaInv(nregimes)
dim   MSRegSigmaInvF(nregimes)
end
*********************************************************************
function  %MSRegPmat  time
type rect %MSRegPmat
type int  time
*
compute %MSRegPMat=%MSPMat(time)
end
*********************************************************************
*
* Creates a parmset with the free parameters from the regression part
* of the model.
*
procedure MSRegParmSet
option parmset *parmset
*
if .not.%defined(parmset) {
   disp "###Usage: @MSRegParmSet(PARMSET=parmset to define)"
   return
}

if MSRegSigmaSwitch
   nonlin(parmset=parmset) gammas betas sigsqv
else
   nonlin(parmset=parmset) gammas betas sigsq
*
end
*********************************************************************
*
* Computes a standard set of initial guess values for the parameters
*
procedure MSRegInitial start end
type integer  start end
*
option series[int] guessregimes
*
local integer i j
local real    spacing
*
linreg(noprint,equation=MSRegEqn) yuni start end
if nregimes>1
   ewise sigsqv(i)=%seesq*(nregimes^(2*i-nregimes-1.0))
else
   compute sigsqv(i)=%seesq
compute sigsq=%seesq
*
if MSRegNFix>0
   compute gammas=%xsubvec(%beta,1,MSRegNFix)
if MSRegNSwitch>0 {
   *
   * Assuming the first switching coefficient is the intercept, use
   * guess values which spread across a reasonable range for it.
   *
   compute spacing=%stderrs(MSRegNFix+1)*.5*nregimes/(nregimes-1)
   do i=1,nregimes
      compute betas(i)=%xsubvec(%beta,MSRegNFix+1,nreg)
      compute betas(i)(1)=betas(i)(1)+spacing*(2*i-(nregimes+1))
   end do i
}
*
if %defined(guessregimes) {
   do i=1,nregimes
      linreg(noprint,equation=MSRegEqn,smpl=guessregimes{0}==i) yuni start end
      if MSRegNSwitch>0
         compute betas(i)=%xsubvec(%beta,MSRegNFix+1,nreg)
      if MSRegSigmaSwitch
         compute sigsqv(i)=%seesq
   end do i
}

ewise p(i,j)=%if(i==j,.8,.2/(nregimes-1))
compute theta=%MSPLogistic(p)
end
*********************************************************************
procedure MSRegRelabel swaps
type vect[int] swaps
*
local rect tempp
local vect[vect] tempbetas
local vect tempsigsq
*
if MSRegNSwitch>0 {
   compute tempbetas=betas
   ewise betas(i)=tempbetas(swaps(i))
}
if MSRegSigmaSwitch {
   compute tempsigsq=sigsqv
   ewise sigsqv(i)=tempsigsq(swaps(i))
}
*
* Relabel the transitions
*
compute tempp=p
ewise p(i,j)=tempp(swaps(i),swaps(j))
end
*********************************************************************
*
* Returns the series of regime-dependent residuals
*
procedure MSRegResids resids gstart gend
type series  *resids
type integer gstart gend
*
option series[int]  regime     MSRegime
option integer      specific
*
local integer thisRegime
*
set resids gstart gend = $
   thisRegime=%if(%defined(specific),specific,regime(t)),$
   yuni-%if(MSRegNFix>0,%dot(gammas,%eqnxvector(MSRegEqnFix,t)),0.0)-$
        %if(MSRegNSwitch>0,%dot(betas(thisRegime),%eqnxvector(MSRegEqnSwitch,t)),0.0)
end
*********************************************************************
procedure MSRegStdResiduals resids start end
*
* Returns a series of standardized one-step residuals, with regimes weighted
* by the one-step predictive probabilities.
*
type integer start end
type series  *resids
*
option series[vect] phat pt_t1
*
local integer time k
local real    ey eyk eyy
*
set resids start end = 0.0
*
do time=start,end
   compute ey=0.0
   compute eyy=0.0
   do k=1,nregimes
      compute eyk=0.0
      if MSRegNFix>0
         compute eyk=eyk+%dot(gammas,%eqnxvector(MSRegEqnFix,time))
      if MSRegNSwitch>0
         compute eyk=eyk+%dot(betas(k),%eqnxvector(MSRegEqnSwitch,time))
      *
      compute ey =ey +phat(time)(k)*eyk
      compute eyy=eyy+phat(time)(k)*(sigsqv(k)+eyk^2)
   end do k
   *
   * Convert to variance
   *
   compute eyy=eyy-ey^2
   *
   * Convert to residuals
   *
   compute ey=yuni(time)-ey
   *
   * Standardize and save
   *
   compute resids(time)=ey/sqrt(eyy)
end do time
end
*********************************************************************
*
* Returns the series of residuals from the fixed part
*
procedure MSRegFixResids resids gstart gend
type series  *resids
type integer gstart gend
*
set resids gstart gend = yuni-$
   %if(MSRegNFix>0,%dot(gammas,%eqnxvector(MSRegEqnFix,t)),0.0)
end
*********************************************************************
*
* Returns the series of residuals from the switching part
*
procedure MSRegSwitchResids resids gstart gend
type series  *resids
type integer gstart gend
*
option series[int]  regime
*
set resids gstart gend = yuni-$
   %if(MSRegNSwitch>0,%dot(betas(regime),%eqnxvector(MSRegEqnSwitch,t)),0.0)
end
*********************************************************************
*
* Returns the vector of likelihoods (*not* in log form) for the
* (expanded) regimes at <<time>>.
*
function %MSRegFVec time
type vector  %MSRegFVec
type integer time
*
local integer i
local real    u ufix
*
dim %MSRegFVec(nregimes)
*
* Integrating constant
*
compute ufix=yuni(time)
if MSRegNFix>0
   compute ufix=ufix-%dot(gammas,%eqnxvector(MSRegEqnFix,time))
do i=1,nregimes
   compute u=%if(MSRegNSwitch==0,ufix,$
        ufix-%dot(betas(i),%eqnxvector(MSRegEqnSwitch,time)))
   compute %MSRegFVec(i)=exp(%logdensity(sigsqv(i),u))
end do i
end
*********************************************************************
function %MSRegProb time
type real      %MSRegProb
type integer   time
*
local vector   f
local rect     pmat
local vect     f fp
local real     fpt
*
compute f=%MSRegFVec(time)
compute pmat=%MSPMat(time)
compute pt_t1(time)=pmat*pstar
compute fp=f.*pt_t1(time)
compute fpt=%sum(fp)
compute pstar=fp*(1.0/fpt)
compute pt_t(time)=pstar
compute %MSRegProb=fpt
end
*********************************************************************
function %MSRegInitVariances
*
* Initializes the sigsqv array given the current settings for the
* sigsq's, returning the smallest variance.
*
local integer i

do i=1,nregimes
   if MSRegSigmaSwitch==0 {
      compute sigsqv(i)=sigsq
      if i>1 {
         compute MSRegLogDetV(i)  =MSRegLogDetV(i-1)
         compute MSRegSigmaInv(i) =MSRegSigmaInv(i-1)
         compute MSRegSigmaInvF(i)=MSRegSigmaInvF(i-1)
         next
      }
   }
   compute MSRegLogDetV(i)  =log(sigsqv(i))
   compute MSRegSigmaInv(i) =1.0/sigsqv(i)
   compute MSRegSigmaInvF(i)=1.0/sqrt(sigsqv(i))
end do i
compute %MSRegInitVariances=%minvalue(sigsqv)
end
*********************************************************************
function  %MSRegInitTransition
type real %MSRegInitTransition
*
compute MSPFull=%MSPmat(%RegStart())
compute %MSRegInitTransition=1
end
*********************************************************************
function    %MSRegInit
type vector %MSRegInit
*
* If any of the transition probabilities are negative, the parameter
* setting will be rejected.
*
if %MSRegInitTransition()==0.0 {
   compute %MSRegInit=%zeros(nregimes,1)
   return
}
compute %MSRegInitVariances()
compute %MSRegInit=%MCErgodic(MSPFull)
end
***************************************************
procedure MSRegFilter gstart gend
type integer gstart gend
*
option vector p0
*
local vect 		f
local integer 	time
*
@MSFilterInit(p0=p0)
do time=gstart,gend
   compute f=%MSRegFvec(time)
   @MSFilterStep time f
end do time
end
***************************************************
* EM calculations
***************************************************
procedure MSRegEMGeneralSetup
option integer regimes nregimes
*
@MSEMSetupStd(regimes=regimes)
dim p(regimes,regimes)
end
**********************************************************************
function %MSRegInitEM
type vector %MSRegInitEM
compute %MSRegInit()
*compute %MSRegInitEM=EMP0
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 MSRegEStep gstart gend
type integer gstart gend
*
option switch ergodic 0
*
local integer time
local vect    f
local rect    pt
*
***************************************************
* EM calculations for transitions
***************************************************
*
@MSEMFilterInit gstart gend
compute %MSRegInitVariances()
do time=gstart,gend
   compute f=%MSRegFvec(time)
   compute pt=%MSPMat(time)
   @MSEMFilterStep(pmat=pt) time f
*   if MSDebug
*      disp "MSRegEStep" time f MSEMMarginal(MSEMPT_T(time))
end do time
@MSEMSmooth
gset psmooth gstart gend = MSEMMarginal(MSEMpt_sm(t))
end
*******************************************************************************
*
* Does the "M" step of EM, updating the various sets of coefficients.
*
procedure MSRegMStep gstart gend
*
option switch dobeta 	1
option switch dosigma 	1
option switch dop     	1
*
local integer i s
local vect    qhat
local rect    phistack
local series  wt
local symm    cmomEM
local series[int] iregime
local real    sigsum wtsum thissumsq thissumwt
*
* Do time-varying given fixed
*
if MSRegNSwitch>0 {
   @MSRegFixResids MSRegU gstart gend
   do s=1,nregimes
      set wt gstart gend = qhat=MSEMMarginal(MSEMpt_sm(t)),qhat(s)
      cmom(weight=wt,equation=MSRegEqnSwitch,depvar) gstart gend
      compute phistack=%sweeptop(%cmom,%ncmom-1)
      ewise betas(s)(i)=phistack(i,%ncmom)
   end do s
}
*
* Do fixed given switching
*
if MSRegNFix>0 {
   compute cmomEM=%zeros(MSRegNFix+1,MSRegNFix+1)
   do s=1,nregimes
      gset iregime gstart gend = s
      @MSRegSwitchResids(regime=iregime) MSRegU gstart gend
      set wt gstart gend = qhat=MSEMMarginal(MSEMpt_sm(t)),qhat(s)
      cmom(weight=wt,spread=sigsqv(s),equation=MSRegEqnFix,depvar) gstart gend
      compute cmomEM=cmomEM+%cmom
   end do s
   compute phistack=%sweeptop(cmomEM,%ncmom-1)
   ewise gammas(i)=phistack(i,%ncmom)
}
if dosigma==0
   branch skipsigma

compute sigsum=0.0
compute wtsum =0.0
do s=1,nregimes
   set wt gstart gend = qhat=MSEMMarginal(MSEMpt_sm(t)),qhat(s)
   gset iregime gstart gend = s
   @MSRegResids(regime=iregime) MSRegU gstart gend
   sstats gstart gend MSRegU^2*wt>>thisSumSq wt>>thisSumWt
   if MSRegSigmaSwitch
      compute sigsqv(s)=thisSumSq/thisSumWt
   else
      compute sigsum=sigsum+thisSumSq,wtsum=wtsum+thisSumWt
end do s
if MSRegSigmaSwitch==0
   compute sigsq=sigsum/wtsum

:skipsigma

*
* Update P.
*
if dop
   @MSEMDoPMatrix
end MSRegMStep
********************************************************************************
*
* MSRegEMStep does a single step (both "E" and "M") for a Markov
* switching regression.
*
procedure MSRegEMStep gstart gend
*
option switch ergodic 0
option switch dobeta  1
option switch dosigma 1
option switch dop     1
*
@MSRegEStep(ergodic=ergodic) gstart gend
@MSRegMStep(dobeta=dobeta,dosigma=dosigma,dop=dop) gstart gend
end
********************************************************************************
*
* This takes an input set of AR coefficients (starting with lag 1) and
* returns 1 if they create an unstable AR and 0 if they are stable.
*
function %MSRegARIsUnstable arc
type integer %MSRegARIsUnstable
type vector  arc
*
local cvector roots
*
compute roots=%polycxroots(-1.0~~arc)
compute %MSRegARIsUnstable=fix(%cabs(roots(%size(arc)))<=1.0)
end
********************************************************************************
*
* Defines the equation (coefficients and variance) for a particular
* regime.
*
procedure MSRegSetEquation
option integer    regime   1
option equation   *define  MSRegEqn
*
compute eqn=MSRegEqn
if MSRegNFix==0
   compute %EqnSetCoeffs(eqn,betas(regime))
else
   compute %EqnSetCoeffs(eqn,gammas~~betas(regime))
if MSRegSigmaSwitch
   compute %EqnSetVariance(eqn,sigsqv(regime))
else
   compute %EqnSetVariance(eqn,sigsq)
end
