## markovSwitchingExample.ssc
## Examples of Markov switching models
##
## author: Eric Zivot
## created: November 3, 2005
## updated: November 8, 2005
##
## Comments
## 1. Requires S-PLUS 7 and S+FinMetrics 2.0



## Section 5. Markov Switching State Space Models

# compute ergodic probabilities from 2-state Markov chain
mchain.p0(matrix(c(0.47, 0.05, 0.53, 0.95), 2, 2))

# compute state space form for MS-AR(2)
GetSsfMSAR = function(parm)
{
  mDelta = mDelta.other = rep(0, 3)
  mDelta[1] = parm[1]
  mDelta.other[1] = parm[2]
#
  mPhi = mPhi.other = matrix(0, 3, 2)
  mPhi[1,] = c(parm[3], parm[4])
  mPhi.other[1,] = c(parm[5], parm[6])
  mPhi[2:3,1] = mPhi.other[2:3,1] = 1
#
  mOmega = mOmega.other = matrix(0, 3, 3)
  mOmega[1,1] = parm[7]
  mOmega.other[1,1] = parm[8]
#
  mSigma = matrix(0, 3, 2)
  mSigma[1:2, 1:2] = diag(1e+6, 2)
#
  mTrans = matrix(0, 2, 2)
  mTrans[1,1] = parm[9]
  mTrans[1,2] = 1 - mTrans[1,1]
  mTrans[2,2] = parm[10]
  mTrans[2,1] = 1 - mTrans[2,2]
#
  list(mDelta=mDelta, mDelta.other=mDelta.other, 
       mPhi=mPhi, mPhi.other=mPhi.other, 
       mOmega=mOmega, mOmega.other=mOmega.other,
       mSigma=mSigma, mTrans=mTrans)
}

# test parameters
# mu1 = -0.5, mu2 = 0.5, phi.11 = 0.7, phi.21 = 0.2, phi.12 = 0.5, phi.22 = 0,
# sigma1.sq = 1, sigma2.sq = 2, p11 = 0.9, p22 = 0.9)

msAR2.parm.names = c("mu.1", "mu.2", "phi.11", "phi.21", "phi.12", "phi.22",
                      "sigma2.1", "sigma2.2", "p11", "p22")

msAR2.parm = c(-0.5, 0.5, 0.7, 0.2, 0.5, 0, 1, 2, 0.9, 0.9)
names(msAR2.parm) = msAR2.parm.names

msAR2.ssf = GetSsfMSAR(msAR2.parm)
msAR2.ssf

# compute state space form for Hamilton's MS-AR(4) model
GetSsfHamilton <- function(parm)
# parm = (mu1,mu2,phi1,phi2,phi3,phi4,sigma,p11,p22)
{
	mDelta = mDelta.other = rep(0, 5)
	mDelta[5] = parm[1]
	mDelta.other[5] = parm[2]
#
	mPhi = matrix(0, 5, 4)
	mPhi[1, ] = parm[3:6]
	mPhi[2:4, 1:3] = diag(3)
	mPhi[5, 1] = 1
#
	mOmega = matrix(0, 5, 5)
	mOmega[1,1] = parm[7]
#
	mSigma = matrix(0, 5, 4)
	mSigma[1:4, 1:4] = diag(1e+6, 4)
#
  mTrans = matrix(0, 2, 2)
  mTrans[1,1] = parm[8]
  mTrans[1,2] = 1 - mTrans[1,1]
  mTrans[2,2] = parm[9]
  mTrans[2,1] = 1 - mTrans[2,2]
#
  list(mDelta=mDelta, mDelta.other=mDelta.other, 
       mPhi=mPhi, mOmega=mOmega,
       mSigma=mSigma, mTrans=mTrans)
}

msHamilton.parm.names = c("mu.1", "mu.2", "phi.1", "phi.2",
                          "phi.3", "phi.4", "sigma", "p11",
                          "p22")
# take parameters from Kim and Nelson (1999) Table 4.1
msHamilton.parm = c(-0.2132, 1.1283, 0.0898, -0.0186, 
                    -0.1743, -0.0839, 0.7962, 0.7606, 0.9008)
names(msHamilton.parm) = msHamilton.parm.names
msHamilton.ssf = GetSsfHamilton(msHamilton.parm)
msHamilton.ssf

#
# time varying transition probabilities
#

#
# The time varying transition probability models are estimated as:
#
#        e^5                 e^(a0 + a1*x)
# ---------------------  ---------------------
#  e^5 + e^(a0 + a1*x)    e^5 + e^(a0 + a1*x)
#
#      e^(b0 + b1*x)             e^5
# ---------------------  ---------------------
#  e^5 + e^(b0 + b1*x)    e^5 + e^(b0 + b1*x)
#
# so the coefficients have to specified as:
# 
#   [ a0  a1 ... ]
#   [ b0  b1 ...]
#

Ssf.AR1 = function(parm, trans.data=NULL) { 
# parm = c(mu1, mu2, phi, sigma, p11, p22, 
    mDelta = c(parm[1], 0) 
    mDelta.other = c(parm[2], 0) 
    mPhi = c(parm[3], 1) 
    mOmega = diag(c(parm[4]^2, 0)) 
    mSigma = c(1e+6, 0) 
#
#    mTrans = matrix(0, 2, 2) 
#    mTrans[1,2] = parm[5] 
#    mTrans[1,1] = 1 - mTrans[1,2] 
#    mTrans[2,1] = parm[6] 
#    mTrans[2,2] = 1 - mTrans[2,1] 

		mTrExo = trans.data
		cTrExo = 1
		cftrExo = c(
    list(mDelta=mDelta, 
         mDelta.other=mDelta.other, 
         mPhi=mPhi, 
         mOmega=mOmega, 
         mSigma=mSigma, 
         mTrans=mTrans, mTrExo=matrix(rnorm(151), ncol=1), cTrExo=1, cfTrExo=c(100,100,100,100) 
                        ) 
  } 
temp <- SsfFitMS(c(-0.5, 1, 0.2, 0.3, 0.5, 0.5), kim.ret, 
Ssf.AR1, lower=c(rep(-Inf,2), -1, -Inf, rep(0,2)), 
upper= c(rep(Inf,2), 1, Inf, rep(1,2)), l.start=7) 


#
# Simulations from MS state space models
#

# simulate from MSAR(2) model
args(SsfSimMS)
# function(ssf, n = 100, mRan = NULL, seed = NULL)
msAR2.sim = SsfSimMS(msAR2.ssf, n = 250, seed = 10)
class(msAR2.sim)
names(msAR2.sim)
names(msAR2.sim$states)

# plot simulated data
par(mfrow=c(2,1))
tsplot(msAR2.sim$states[,"response"])
tsplot(msAR2.sim$regimes)
par(mfrow=c(1,1))

# simulation from Hamilton's model
msHamilton.sim = SsfSimMS(msHamilton.ssf, n = 150, seed = 10)
class(msHamilton.sim)
names(msHamilton.sim)

# plot simulated data
par(mfrow=c(2,1))
tsplot(msHamilton.sim$states[,"response"], 
       main="Simulated real GDP growth")
abline(h=0)
tsplot(msHamilton.sim$regimes, main = "Simulated states")
par(mfrow=c(1,1))




#
# MLE of MS state space models
#

GetSsfMSAR2 = function(parm)
{
  parm = as.vector(parm)
#
  mDelta = mDelta.other = rep(0, 3)
  mDelta[1] = parm[1]
  mDelta.other[1] = parm[1] + exp(parm[2])
#
  AR11 = 2/(1+exp(-parm[3])) - 1
  AR12 = 2/(1+exp(-(parm[3]+exp(parm[4])))) - 1
  AR21 = 2/(1+exp(-parm[5])) - 1
  AR22 = 2/(1+exp(-(parm[5]+exp(parm[6])))) - 1
#
  mPhi = mPhi.other = matrix(0, 3, 2)
  mPhi[1,] = c(AR11+AR12, -AR11*AR12)
  mPhi.other[1,] = c(AR21+AR22, -AR21*AR22)
  mPhi[2:3,1] = mPhi.other[2:3,1] = 1
#
  mOmega = matrix(0, 3, 3)
  mOmega[1,1] = exp(parm[7])
#
  mSigma = matrix(0, 3, 2)
  mSigma[1:2, 1:2] = diag(1e+6, 2)
#
  mTrans = matrix(0, 2, 2)
  mTrans[1,2] = 1/(1+exp(-parm[8]))
  mTrans[1,1] = 1 - mTrans[1,2]
  mTrans[2,1] = 1/(1+exp(-parm[9]))
  mTrans[2,2] = 1 - mTrans[2,1]
#
  ssf = list(mDelta=mDelta, mDelta.other=mDelta.other, 
        mPhi=mPhi, mPhi.other=mPhi.other, mOmega=mOmega, 
        mTrans=mTrans, mSigma=mSigma)
  CheckSsf(ssf)
}

ndx.start = c(-2, -0.7, -0.7, 0.7, -0.7, 0.7, -2, -2, -3)
names(ndx.start) = c("mu1", "mu2", "phi11", "phi12", "phi21", "phi22", "sigma", "p", "q") 
ndx.msar = SsfFitMS(ndx.start, log(ndx.rvol), GetSsfMSAR2, l.start=11)
class(ndx.msar)
summary(ndx.msar)

ndx.ssf = GetSsfMSAR2(ndx.msar$parameters)
cbind(ndx.ssf$mDelta, ndx.ssf$mDelta.other)
ndx.ssf$mPhi
ndx.ssf$mPhi.other
ndx.ssf$mOmega
ndx.ssf$mTrans

ndx.f = SsfLoglikeMS(log(ndx.rvol), ndx.ssf, save.rgm=T)
par(mfrow=c(2,1))
plot(timeSeries(ndx.f$regimes[,1], pos=positions(ndx.rvol)), 
     reference.grid=F, main="Filtered Low Vol Regime Prob")
plot(timeSeries(ndx.f$regimes[,3], pos=positions(ndx.rvol)), 
     reference.grid=F, main="Smoothed Low Vol Regime Prob")

## Section 6. An Extended Example

GetSsfCoinIndex = function(parm) {
  parm = as.vector(parm)
  mDelta = mDelta.other = rep(0, 11)
  mDelta[1] = parm[1]
  mDelta.other[1] = parm[1] + exp(parm[2])
#
  AR.C1 = 2/(1+exp(-parm[3])) - 1
  AR.C2 = 2/(1+exp(-(parm[3]+exp(parm[4])))) - 1
#
  AR.e1 = 2/(1+exp(-parm[5])) - 1
  AR.e2 = 2/(1+exp(-parm[6])) - 1
  AR.e3 = 2/(1+exp(-parm[7])) - 1
  AR.e4 = 2/(1+exp(-parm[8])) - 1
#
  mPhi = matrix(0, 11, 7)
  mPhi[1,1:2] = c(AR.C1+AR.C2, -AR.C1*AR.C2)
  mPhi[2,1] = 1
  mPhi[3,3] = AR.e1
  mPhi[4,4] = AR.e2
  mPhi[5,5] = AR.e3
  mPhi[6,6] = AR.e4
  mPhi[7,1] = mPhi[7,7] = 1
#
  mPhi[8:10,1] = parm[9:11]
  mPhi[11,1:2] = parm[12:13]
  mPhi[8,3] = mPhi[9,4] = mPhi[10,5] = mPhi[11,6] = 1
#
  mOmega = matrix(0, 11, 11)
  mOmega[1,1] = 1
  mOmega[3,3] = exp(parm[14])
  mOmega[4,4] = exp(parm[15])
  mOmega[5,5] = exp(parm[16])
  mOmega[6,6] = exp(parm[17])
#
  mTrans = matrix(0, 2, 2)
  mTrans[1,2] = 1/(1+exp(-parm[18]))
  mTrans[1,1] = 1-mTrans[1,2]
  mTrans[2,1] = 1/(1+exp(-parm[19]))
  mTrans[2,2] = 1-mTrans[2,1]
#
  mSigma = matrix(0, 8, 7)
  mSigma[1:7, 1:7] = diag(1e+6, 7)
  ans = list(mDelta=mDelta, mDelta.other=mDelta.other, 
             mSigma=mSigma, mOmega=mOmega, 
             mPhi=mPhi, mTrans=mTrans)
  CheckSsf(ans)
}

DOC.dat = getReturns(DOC.ts[,1:4], percentage=T)
DOC.dat@data = t(t(DOC.dat@data) - colMeans(DOC.dat@data))
DOC.dat@data = t(t(DOC.dat@data) / colStdevs(DOC.dat@data))

DOC.start = c(-1.5, 0.6, 0.3, 0.1, .1, .1, .1, .1, 0.3, 
               0.3, 0.3, 0.3, 0.1, -.5, -.5, -.5, -.5, -1.5, -3)
names(DOC.start) = c("mu1", "mu2", "phi1", "phi2", "psi1", 
       "psi2", "psi3", "psi4", "L1", "L2", "L3", "L41", 
       "L42", "s1", "s2", "s3", "s4", "p", "q")
DOC.fit = SsfFitMS(DOC.start, DOC.dat, GetSsfCoinIndex, 
          l.start=13, trace=T)
summary(DOC.fit)

DOC.ssf = GetSsfCoinIndex(DOC.fit$parameters)
c(DOC.ssf$mDelta[1], DOC.ssf$mDelta.other[1])
print(DOC.ssf$mPhi, digits=3)

DOC.f = SsfLoglikeMS(DOC.dat, DOC.ssf, save.rgm=T, l.start=13)
names(DOC.f)

DOC.dates = positions(DOC.dat)[-(1:12)]
filt.p = timeSeries(DOC.f$regimes[,1], pos=DOC.dates)
smoo.p = timeSeries(DOC.f$regimes[,3], pos=DOC.dates)
par(mfrow=c(2,1))
plot(filt.p, reference.grid=F, main="Filtered Recession Probability")
plot(smoo.p, reference.grid=F, main="Smoothed Recession Probability")

DOC.index = lm(DOC.ts@data[,5]~I(1:433))$residuals[-(1:13)]
filt.ci = rowSums(DOC.f$state[,c(7,14)]*DOC.f$regime[,1:2])
filt.ci = timeSeries(filt.ci, pos=DOC.dates)
plot(filt.ci, reference.grid=F, main="Filtered MS Coincident Index")
doc.ci = timeSeries(DOC.index, pos=DOC.dates)
plot(doc.ci, reference.grid=F, main="DOC Coincident Index")

# Use Andrew's shadePlot to add business cycle dates to plot
NBERShade = list(x=data.frame(timeDate(c("04/01/60","10/01/69","11/01/73","01/01/80","07/01/81","07/01/90")), 
                              timeDate(c("02/01/61","11/01/70","03/01/75","07/01/80","11/01/82","03/01/91"))),
									col=12)

par(mfrow=c(2,1))
shadePlot(filt.p, reference.grid=F, main="Filtered Recession Probability",
		   shade.regions=list(NBERShade))
shadePlot(smoo.p, reference.grid=F, main="Smoothed Recession Probability", 
          shade.regions=list(NBERShade))