*
* These four functions are fairly generic Markov Chain/Markov Switching
* model support functions. In all cases, the transition matrix P is
* represented by an (N-1) x N matrix, with the final row omitted, as it
* just forces the columns to sum to 1. This makes it easier to estimate
* the free parameters. Note, however, that the vector of probabilities of
* the states is the full N vector.
*
* %MCSTATE(P,PSTAR) returns the updated state probabilities given initial
* probabilities PSTAR and transition matrix P.
*
* %MCERGODIC(P) returns the ergodic state probabilities for transition
* matrix P.
*
* %MSUPDATE(F,PSTAR,FPT) returns the updated state probabilities given
* the vector of likelihoods (levels, not logs) of the states in F and the
* prior probabilities of the states in the vector PSTAR. FPT is a REAL returned
* by the function which gives the sum over i of F x PSTAR. If F is a vector of
* likelihoods given the states, FPT will be the overall likelihood for this
* entry.
*
* @%MSSMOOTH P PT_T PT_T1 PSMOOTH is a procedure which produces the VECT[SERIES]
* PSMOOTH given the coded transition probability matrix P and the SERIES[VECT]
* PT_T and PT_T1, which give the estimated probabilities at t given information
* through t and t-1 respectively.
*
*******************************************************************************
function %mcstate p pstar
type vect %mcstate pstar
type rect p
local vect pstub
local integer i
*
compute pstub=p*pstar
dim %mcstate(%cols(p))
ewise %mcstate(i)=%if(i<=%rows(p),pstub(i),1-%sum(pstub))
end
*******************************************************************************
function %mcergodic p
type vect %mcergodic
type rect p
*
local rect a
local integer i j n
*
compute n=%cols(p)
dim a(n+1,n)
ewise a(i,j)=%if(i>n,1,(i==j)-%if(i<n,p(i,j),1-%sum(%xcol(p,j))))
compute %mcergodic=%xcol(inv(%innerxx(a))*tr(a),n+1)
end
*******************************************************************************
function %msupdate f pstar fpt
type vect %msupdate f pstar
type real *fpt
*
local vector fp
*
compute fp=f.*pstar
compute fpt=%sum(fp)
*
* If the f's come in zero, just keep the old pstar. The parameters
* will get rejected anyway since the log(fpt) will be NA.
*
if %minvalue(pstar)<0.00
   compute %msupdate=pstar,fpt=%na
else
if fpt>0
   compute %msupdate=fp/fpt
else
   compute %msupdate=pstar
end
*******************************************************************************
procedure %mssmooth p pt_t pt_t1 psmooth
type rect p
type series[vect] pt_t pt_t1 *psmooth
*
local rect pexpand pfull
local integer i j nstates time
local vect pstar
*
* Expand p to a full nxn matrix
*
compute nstates=%cols(p)
dim pexpand(nstates,nstates-1)
ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
compute pfull=pexpand*p
ewise pfull(i,j)=%if(i==nstates,1.0+pfull(i,j),pfull(i,j))
*
compute pstar=pt_t(%regstart())
gset psmooth %regstart() %regend() = pstar
*
do time=%regend()-1,%regstart(),-1
   gset psmooth time time = pt_t.*(tr(pfull)*(psmooth{-1}./pt_t1{-1}))
end do time
end





