# StateSpaceExamples.ssc		script file for examples used in State Space chapter
#
# author: Eric Zivot
# created: January 12, 2001
# revised: May 17, 2001
#
# Examples: 
#
# 1. Follow examples in statistical algorithms for models in state space using
# ssfpack2.2 by Koopman, Shephard and Doornik, Econometric Journal, 1999.
# 2. Utilize updated examples in SsfPack 3.0 beta: Statistical algorithms
# for models in state space by Koopman, Shephard and Doornik
# 3. Utilize updated ssfpack 3.0 functions.
#

# state space form used by ssfpack 3.0
#
# 1. alpha(t+1) = d(t) + T(t)*alpha(t) + H(t)*n(t)		transition equation
#  (m x 1)
# 2. theta(t) = c(t) + Z(t)*alpha(t)						signal
# 3. y(t) = theta(t) + G(t)*e(t)							measurement equation
# (N x 1)
# 4. n(t) ~ N(0,I), e(t) ~ N(0,I), alpha(1) ~ N(a,P)
# 5. E[e(t)n(t)']=0
#
# system matrices: T(t), H(t), Z(t), G(t)
# fixed components: c(t), d(t)
#
# state space form used by ssfpack 3.0
#
# 6. vec(alpha(t+1),y(t)) = delta(t) + Phi(t)*alpha(t) + u(t)
# ((m+N) x 1)                    ((m+N) x m) (m x 1)
# 7. u(t) ~ N(0,Omega(t))
#           ((m+N)x(m+N))
# 8. delta(t) = vec(d(t),c(t)), Phi(t) = (T(t)' Z(t)')', u(t) = (H(t)' G(t)')'e(t)
# 
# 9. Omega(t) = |H(t)H(t)'   0     |
#               |   0     G(t)G(t)'|
# 
# 10. Sigma = |P |
#             |a'|

# 
# Ex. local level model
#
# y(t) = a(t) + e(t), e(t) ~ N(0,s2e)
# a(t+1) = a(t) = n(t), n(t) ~ N(0,s2n)
# a(1) ~ N(a1,P1)
# 


# construct state space form
sigma.e = 1
sigma.n = 0.5
a1 = 0
P1 = -1
ssf.ll.list = list(mPhi=as.matrix(c(1,1)),
mOmega=diag(c(sigma.n^2,sigma.e^2)),
mSigma=as.matrix(c(P1,a1)))
ssf.ll.list

ssf.ll = CheckSsf(ssf.ll.list)
class(ssf.ll)
names(ssf.ll)

# use GetSsfStsm 

ssf.stsm = GetSsfStsm(irregular=1,level=0.5)
class(ssf.stsm)
ssf.stsm

#
# time varying CAPM
#
X.mat = cbind(1,as.matrix(seriesData(excessReturns.ts[,"SP500"])))
Phi.t = rbind(diag(2),rep(0,2))
Omega = diag(c((.01)^2,(.05)^2,(.1)^2))
J.Phi = matrix(-1,3,2)
J.Phi[3,1] = 1
J.Phi[3,2] = 2
Sigma = -Phi.t
ssf.tvp = list(mPhi=Phi.t,
mOmega=Omega,
mJPhi=J.Phi,
mSigma=Sigma,
mX=X.mat)
ssf.tvp


#
# ARMA Models
#

# GetSsfArma creates state space matrices for ARMA model
args(GetSsfArma)

# AR(1)
ssf.ar1 = GetSsfArma(ar=0.75,sigma=.5)
ssf.ar1

# AR(2)
ar2.mod = list(ar=c(1.25,-0.5),sigma=1)
ssf.ar2 = GetSsfArma(model=ar2.mod)
ssf.ar2

# ARMA(2,1)
# y(t) = 0.6y(t-1)+0.2y(t-2)+e(t)-0.2e(t-1), e(t) ~ N(0,0.9)
arma21.mod = list(ar=c(0.6,0.2),ma=c(-0.2),sigma=sqrt(0.9))
ssf.arma21 = GetSsfArma(model=arma21.mod)
ssf.arma21

# AR(2)
SS.AR2 = GetSsfArma(ar=c(0.6,0.2), sigma=sqrt(0.9))
SS.AR2

# MA(1)
SS.MA1 = GetSsfArma(ma=c(-0.2), sigma=sqrt(0.9))
SS.MA1

#
# Structural Time Series Models
#
args(GetSsfStsm)

# Example 2 from KSD
ssf.stsm = GetSsfStsm(irregular=1,
level=0.5,
slope=0.1,
seasonalDummy=c(0.2,3.0))
ssf.stsm

#
# Regression models
#
args(GetSsfReg)

# regression on a constant and time trend
ssf.reg = GetSsfReg(cbind(1,1:10))
class(ssf.reg)
names(ssf.reg)
ssf.reg

tmp = CheckSsf(ssf.reg)
tmp

# time varying time trend model
ssf.reg.tvp = ssf.reg
ssf.reg.tvp$mOmega[2,2] = 0.5
ssf.reg.tvp$mOmega

# time trend regression with AR(2) errors
args(GetSsfRegArma)
ssf.reg.ar2 = GetSsfRegArma(cbind(1,1:10),
ar=c(1.25,-0.5))

#
# non-parametric cubic spline models
#

args(GetSsfSpline)
GetSsfSpline()
t.vals = c(2,3,5,9,12,17,20,23,25)
delta.t = diff(t.vals)
GetSsfSpline(snr=0.2,delta=delta.t)
#
# simulate data using SsfSim
#

args(SsfSim)

# simulate local level model

set.seed(112)
ll.sim = SsfSim(ssf.ll.list,n=250)
class(ll.sim)
colIds(ll.sim)

tsplot(ll.sim)
legend(0,4,legend=c("State","Response"),lty=1:2)

# simulate time varying CAPM

nrow(X.mat)

set.seed(444)
tvp.sim = SsfSim(ssf.tvp,n=nrow(X.mat),a1=c(0,1))
colIds(tvp.sim)

par(mfrow=c(3,1))
tsplot(tvp.sim[,"state.1"],main="Time varying intercept",
ylab="alpha(t)")
tsplot(tvp.sim[,"state.2"],main="Time varying slope",
ylab="beta(t)")
tsplot(tvp.sim[,"response"],main="Simulated returns",
ylab="returns")
par(mfrow=c(1,1))

#
# Algorithms: Kalman filter, moment estimation, conditional density estimation
# Kalman smoother, moment smoother
#

# Kalman Filter for local level model
# create data
y.ll = ll.sim[,"response"]

# Kalman filtering
args(KalmanFil)

KalmanFil.ll = KalmanFil(y.ll,ssf.ll,task="STFIL")
class(KalmanFil.ll)
names(KalmanFil.ll)
KalmanFil.ll$mOut
# show filtered estimates
KalmanFil.ll$mEst

# Kalman smoother
args(KalmanSmo)
KalmanSmo.ll = KalmanSmo(KalmanFil.ll,ssf.ll)
class(KalmanSmo.ll)
names(KalmanSmo.ll)
plot(KalmanSmo.ll,layout=c(1,2))

# SsfMomentEst
# valid tasks are
# STSMO	state smoothing
# STPRED	state prediction
# STFIL	state filtering
# DSSMO	disturbance smoothing

# moment filtering with variances
FilteredEst.ll = SsfMomentEst(y.ll,ssf.ll,task="STFIL")
class(FilteredEst.ll)
names(FilteredEst.ll)
FilteredEst.ll$state.moment[1:5]
FilteredEst.ll$response.moment[1:5]

# note: plot does not show std error bands
plot(FilteredEst.ll,layout=c(1,2))

# plot with standard error bands

upper.state = FilteredEst.ll$state.moment + 
2*sqrt(FilteredEst.ll$state.variance)
lower.state = FilteredEst.ll$state.moment - 
2*sqrt(FilteredEst.ll$state.variance)
tsplot(FilteredEst.ll$state.moment,upper.state,lower.state,
lty=c(1,2,2),ylab="filtered state")

#  moment smoothing
SmoothedEst.ll = SsfMomentEst(y.ll,ssf.ll,task="STSMO")
SmoothedEst.ll$state.moment[1:5]
SmoothedEst.ll$response.moment[1:5]

plot(SmoothedEst.ll)

# plot with standard error bands

upper.state = SmoothedEst.ll$state.moment + 
2*sqrt(SmoothedEst.ll$state.variance)
lower.state = SmoothedEst.ll$state.moment - 
2*sqrt(SmoothedEst.ll$state.variance)

tsplot(SmoothedEst.ll$state.moment,upper.state,lower.state,
lty=c(1,2,2),ylab="smoothed state")


# note: SmoothedEst.ll2$state.moment is identical to
# SmoothedEst.ll$state


# SsfCondDens
# returns smoothed estimates of state vector and response
# note: only state smoothing (task="STSMO") and disturbance smoothing
# (task="DSSMO") are supported

smoothedEst2.ll = SsfCondDens(y.ll,ssf.ll,task="STSMO")
class(smoothedEst2.ll)
names(smoothedEst2.ll)
smoothedEst2.ll$state[1:5]
smoothedEst2.ll$response[1:5]

plot(smoothedEst2.ll)

# note: in this example 
# smoothedEst.ll$state = smoothedEst2.ll$response
smoothedEst2.ll$state[1:5]
smoothedEst2.ll$response[1:5]

# disturbance smoothing
disturbEst.ll = SsfMomentEst(y.ll,ssf.ll,task="DSSMO")

# moment prediction

# append 5 missing values to the end of y.ll
y.ll.new = c(y.ll,rep(NA,5))
PredictedEst.ll = SsfMomentEst(y.ll.new,ssf.ll,task="STPRED")
y.ll.fcst = PredictedEst.ll$response.moment
fcst.var = PredictedEst.ll$response.variance
upper = y.ll.fcst + 2*sqrt(fcst.var)
lower = y.ll.fcst - 2*sqrt(fcst.var)
upper[1:250] = lower[1:250] = NA
tsplot(y.ll.new[240:255],y.ll.fcst[240:255],
upper[240:255],lower[240:255],lty=c(1,2,2,2))

#
# exact mle of AR(1)
# 

# specify state space form
# simulate 250 observations
ssf.ar1 = GetSsfArma(ar=0.75,sigma=.5)
set.seed(598)
sim.ssf.ar1 = SsfSim(ssf.ar1,n=250)
y.ar1 = sim.ssf.ar1[,"response"]

# do OLS estimation
OLS(y.ar1~tslag(y.ar1)-1)

# Note: variance is defined as var=exp(theta)
# and the likelihood is maximized wrt theta, which
# is defined for all real values
# no transformation is used to ensure stationary AR1
ar1.mod2 = function(parm) {
	phi = parm[1]
	sigma2 = exp(parm[2])
	ssf.mod = GetSsfArma(ar=phi,sigma=sqrt(sigma2))
	CheckSsf(ssf.mod)
}

ar1.start = c(0.5,0)
names(ar1.start) = c("phi","ln(sigma2)")
ar1.mle = SsfFit(ar1.start,y.ar1,"ar1.mod2")
class(ar1.mle)
names(ar1.mle)
ar1.mle$parameters
exp(ar1.mle$parameters["ln(sigma2)"])
sqrt(exp(ar1.mle$parameters["ln(sigma2)"]))

# note: standard errors are not computed unless
# a formula for computing the hessian is provided

# Note: no transformation is used to guarantee a positive
# variance or to force stationarity
ar1.mod = function(parm) {
	phi = parm[1]
	sigma2 = parm[2]
	ssf.mod = GetSsfArma(ar=phi,sigma=sqrt(sigma2))
	CheckSsf(ssf.mod)
}

ar1.start = c(0.5,1)
names(ar1.start) = c("phi","sigma2")
# evaluate log-likelihood at starting values
KalmanFil(y.ar1,ar1.mod(ar1.start),task="KFLIK")

# set lower bound equal to zero in ensure positive variance
# and stationary AR(1) parameter
ar1.mle = SsfFit(ar1.start,y.ar1,"ar1.mod",
lower=c(-.999,0),upper=c(0.999,Inf))
names(ar1.mle)
ar1.mle$parameters
# evaluate log-likelihood at mles
KalmanFil(y.ar1,ar1.mod(ar1.mle$parameters),task="KFLIK")

# maximization of the concentrated log-likelihood function

ar1.start = c(0.5,1)
names(ar1.start) = c("phi","sigma2")
ar1.cmle = SsfFit(ar1.start,y.ar1,"ar1.mod",conc=T)
names(ar1.cmle)
ar1.cmle$parameters
ar1.KF = KalmanFil(y.ar1,ar1.mod(ar1.cmle$parameters),task="KFLIK")
ar1.KF$dVar
#
# mle of local level model
#

ll.mod = function(parm) {
	parm = exp(parm)
	ssf.mod = GetSsfStsm(irregular=sqrt(parm[1]),
							level=sqrt(parm[2]))
	CheckSsf(ssf.mod)
}

ll.start = c(0,0)
names(ll.start) = c("ln(sigma2.n)","ln(sigma2.e)")
ll.mle = SsfFit(ll.start,y.ll,"ll.mod")
ll.mle$parameters
exp(ll.mle$parameters)
sqrt(exp(ll.mle$parameters))

#
# mle of time varying CAPM
#
tvp.mod = function(parm,mX=NULL) {
	parm = exp(parm)
	Phi.t = rbind(diag(2),rep(0,2))
   Omega = diag(parm)
   J.Phi = matrix(-1,3,2)
   J.Phi[3,1] = 1
   J.Phi[3,2] = 2
   Sigma = -Phi.t
   ssf.tvp = list(mPhi=Phi.t,
   mOmega=Omega,
   mJPhi=J.Phi,
   mSigma=Sigma,
   mX=mX)
   CheckSsf(ssf.tvp)
}

tvp.start = c(0,0,0)
names(tvp.start) = c("ln(s2.alpha)","ln(s2.beta)","ln(s2.y)")

y.capm = as.matrix(seriesData(excessReturns.ts[,"MSFT"]))
X.mat = cbind(1,as.matrix(seriesData(excessReturns.ts[,"SP500"])))
tvp.mle = SsfFit(tvp.start,y.capm,"tvp.mod",mX=X.mat)
tvp.mle$parameters
sigma.mle = sqrt(exp(tvp.mle$parameters))
names(sigma.mle) = c("s.alpha","s.beta","s.y")
sigma.mle

# compute smoothed state estimates
smoothedEst.tvp = SsfCondDens(y.capm,
tvp.mod(tvp.mle$parameters,mX=X.mat),
task="STSMO")
plot(smoothedEst.tvp,strip.text=c("alpha(t)",
"beta(t)","Expected returns"),main="Smoothed Estimates")

# compute filtered and smoothed estimates and plot results

FilteredEst.tvp = SsfMomentEst(y.capm,
tvp.mod(tvp.mle$parameters,mX=X.mat),
task="STFIL")

SmoothedEst.tvp = SsfMomentEst(y.capm,
tvp.mod(tvp.mle$parameters,mX=X.mat),
task="STSMO")

upper.state = SmoothedEst.tvp$state.moment + 
2*sqrt(SmoothedEst.tvp$state.variance)
lower.state = SmoothedEst.tvp$state.moment - 
2*sqrt(SmoothedEst.tvp$state.variance)

par(mfrow=c(2,1))
tsplot(SmoothedEst.tvp$state.moment[,1],upper.state[,1],
lower.state[,1],
lty=c(1,2,2),ylab="smoothed state")
tsplot(SmoothedEst.tvp$state.moment[,2],upper.state[,2],
lower.state[,2],
lty=c(1,2,2),ylab="smoothed state")
par(mfrow=c(1,1))

#
# simulation smoothing
#

args(SimSmoDraw)
KalmanFil.ll = KalmanFil(y.ll,ssf.ll,task="STSIM")
ll.state.sim = SimSmoDraw(KalmanFil.ll,ssf.ll.list,
task="STSIM")
class(ll.state.sim)
names(ll.state.sim)
plot(ll.state.sim,layout=c(1,2))