*
* These 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.
*
* %MSSMOOTHCALC(pfull,pt_t,pt1_sm,pt1_t) does the smoothing calculation
* for one time period. pt_t is the estimate of the probabilities for
* time period t given data through t, pt1_t is for time period t+1 given
* data through t, pt1_sm is for time period t+1 given all data, and
* pfull is the (full) transition matrix for states at t to t+1. This is
* safeguarded against divides by zero when pt1_t has underflowed down to
* a zero value.
*
* %MSLOGISTICP(theta) returns the result from transforming an N-1 x N
* matrix of logistic indexes into an N-1 x N matrix of transition
* probabilities. (The parameteriziation in terms of logistic indexes can
* be useful in doing ML estimation).
*
* %MSPLOGISTIC(p) is the inverse of %MSLOGISTICP - it takes the N-1 x N
* matrix of transition probabilities and returns the implied logistic
* indexes.
*
* Revision Schedule:
*  09/2008 Correct initialization of psmooth(end) in %MSSMOOTH
*  01/2009 Add P0 option to %MSSMOOTH to compute pre-sample expected value
*  03/2009 Added %MSSMOOTHCALC which does the smoothing update safeguarded
*          against divide by zero.
*  03/2009 Added %MSLOGISTICP which maps a matrix of logistic index values into
*          an N-1 x N matrix of transition probabilities
*
*******************************************************************************
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
*********************************************************************
*
* Transforms an N-1 x N matrix of logistic indexes into an N-1 x N
* matrix of transition probabilities
*
function  %mslogisticp theta
type rect %mslogisticp
type rect theta
*
local vect    totals
local integer nstates i j
*
compute nstates=%cols(theta)
dim %mslogisticp(nstates-1,nstates) totals(nstates)
do j=1,nstates
   compute totals(j)=1.0
   do i=1,nstates-1
      compute %mslogisticp(i,j)=exp(theta(i,j)),totals(j)+=%mslogisticp(i,j)
   end do i
end do j
ewise %mslogisticp(i,j)=%mslogisticp(i,j)/totals(j)
end
*********************************************************************
*
* Transforms an N-1 x N matrix of transition probabilities into an
* N-1 x N matrix of logistic indexes.
*
function  %msplogistic p
type rect %msplogistic
type rect p
*
local vect    totals
local integer nstates i j
*
compute nstates=%cols(p)
dim %msplogistic(nstates-1,nstates) totals(nstates)
do j=1,nstates
   compute totals(j)=1.0
   do i=1,nstates-1
      compute %msplogistic(i,j)=%if(p(i,j)<=0.0,-20.0,log(p(i,j)))
      compute totals(j)-=p(i,j)
   end do i
   compute totals(j)=%if(totals(j)<=0.0,-20.0,log(totals(j)))
end do j
ewise %msplogistic(i,j)=%msplogistic(i,j)-totals(j)
end
*********************************************************************
*
* This does the smoothing update for a single period, safeguarded
* against divides by zero.
*
function    %mssmoothcalc pfull pt_t pt1_sm pt1_t
type vector %mssmoothcalc pt_t pt1_sm pt1_t
type rect pfull
*
local vect pratio(%rows(pfull))
*
ewise pratio(i) = %if(pt1_t(i)==0.0,0.0,pt1_sm(i)/pt1_t(i))
compute %mssmoothcalc = pt_t.*(tr(pfull)*pratio)
end
*********************************************************************
procedure %mssmooth p pt_t pt_t1 psmooth
type rect p
type series[vect] pt_t pt_t1 *psmooth
*
option vect *p0
*
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(%regend())
gset psmooth %regstart() %regend() = pstar
*
do time=%regend()-1,%regstart(),-1
   gset psmooth time time = %mssmoothcalc(pfull,pt_t,psmooth{-1},pt_t1{-1})
end do time
*
* If requested, get smoothed pre-sample value as well
*
if %defined(p0)
   compute p0=%mssmoothcalc(pfull,p0,psmooth(%regstart()),pt_t1(%regstart()))
end


