Code: Select all
*
* 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