# StateSpaceExamples.ssc		script file for examples used in State Space chapter
# 									updated for FM 2.0
#
# author: Eric Zivot
# created: January 12, 2001
# revised: May 17, 2001
# revised: January 5, 2005
#	Fixed AR1 with concentrated loglikelihood
#	Added SV example
# revised: January 7, 2004
#	made slight changes to estimation of TVP CAPM
# revised: July 6, 2007
#
#
# 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

# plot method
plot(KalmanFil.ll)

# 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: 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")

# use nlminb control options to
# set lower bound equal to zero to ensure positive variance
# and force 0 < phi < 1 for stationary AR(1) 
# see help(nlminb) and help(nlminb.control) for details
ar1.mle = SsfFit(ar1.start,y.ar1,"ar1.mod",
lower=c(-.999,0),upper=c(0.999,Inf))
class(ar1.mle)
names(ar1.mle)
ar1.mle
summary(ar1.mle)

# State space for reparameterized model
# variance is defined as var=exp(theta), theta is unrestricted
# exponential transformation is used to ensure 0 < phi < 1
# for stationary AR(1). Note: negative phi < 0 is ruled out
# apriori
ar1.mod2 = function(parm) {
	phi = exp(parm[1])/(1 + exp(parm[1]))
	sigma2 = exp(parm[2])
	ssf.mod = GetSsfArma(ar=phi,sigma=sqrt(sigma2))
	CheckSsf(ssf.mod)
}

ar1.start = c(0,0)
names(ar1.start) = c("ln(phi/(1-phi))","ln(sigma2)")
ar1.mle = SsfFit(ar1.start,y.ar1,"ar1.mod2")
summary(ar1.mle)


# compute standard errors using the delta method
ar1.mle2 = ar1.mle
tmp = coef(ar1.mle)
ar1.mle2$parameters[1] = exp(tmp[1])/(1 + exp(tmp[1]))
ar1.mle2$parameters[2] = exp(tmp[2])
dg1 = exp(-tmp[1])/(1 + exp(-tmp[1]))^2
dg2 = exp(tmp[2])
dg = diag(c(dg1,dg2))
ar1.mle2$vcov = dg%*%ar1.mle2$vcov%*%dg
summary(ar1.mle2)

# evaluate log-likelihood at mles
KalmanFil(y.ar1,ar1.mod2(ar1.mle$parameters),task="KFLIK")

# maximization of the concentrated log-likelihood function
# Note: Hessian computation is not correct here
# the state space needs to be re-specified so that only one
# parameter is being estimated. This way the Hessian will be properly
# computed
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

# Function to compute state space as a function of phi only
# Note: by default, GetSsfArma sets sigma=1 which
# is required to compute concentrated log-likelihood
# in SsfPack
ar1.modc = function(parm) {
	phi = parm[1]
	ssf.mod = GetSsfArma(ar=phi)
	CheckSsf(ssf.mod)
}
ar1.start = c(0.7)
names(ar1.start) = c("phi")
ar1.cmle = SsfFit(ar1.start,y.ar1,"ar1.modc",conc=T,
lower=0,upper=0.999)
summary(ar1.cmle)

# run Kalman filter to recover mle of sig2
ar1.KF = KalmanFil(y.ar1,ar1.modc(ar1.cmle$parameters),task="KFLIK")
ar1.KF$dVar

# note: standard error for sig2 may be recovered by evaluating 
# hessian of full log-likelihood at mles.


#
# 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
#

# function to compute state space for TVP regression
tvp.mod = function(parm,mX=NULL) {
	parm = exp(parm)
	ssf.tvp = GetSsfReg(mX=mX)
	diag(ssf.tvp$mOmega) = parm
	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)
summary(tvp.mle)

# compute mles of transformed parameters
# and delta method standard errors
tvp2.mle = tvp.mle
tvp2.mle$parameters = exp(tvp.mle$parameters/2)		
names(tvp2.mle$parameters) = c("s.alpha","s.beta","s.y")
dg = diag(tvp2.mle$parameters/2)
tvp2.mle$vcov = dg%*%tvp.mle$vcov%*%dg
summary(tvp2.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", scales="free")


# 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))

#
# Stochastic Volatility Model
# ref: Harvey, Ruiz and Shephard (1994). "Multivariate Stochastic
# Variance Models," Review of Economic Studies, 61, 247-264.
#

# define linear state space for SV model
# note: could use logit transformation to impose 0 < phi < 1
sv.mod = function(parm) {
	g = parm[1]
	sigma2.n = exp(parm[2])
	phi = parm[3]
	ssf.mod = list(mDelta=c(g,-1.27),
	mPhi=as.matrix(c(phi,1)),
	mOmega=matrix(c(sigma2.n,0,0,0.5*pi^2),2,2),
	mSigma=as.matrix(c((sigma2.n/(1-phi^2)),g/(1-phi))))
	CheckSsf(ssf.mod)
}

# simulate from SV model using parameters calibrated from HRS
# estimates for $/DM weekly exchange rate
parm.hrs = c(-0.3556,log(0.0312),0.9646)
nobs = 1000
set.seed(179)
e = rnorm(nobs)
xi = log(e^2)+1.27
eta = rnorm(nobs,sd=sqrt(0.0312))
sv.sim = SsfSim(sv.mod(parm.hrs),
mRan=cbind(eta,xi),a1=(-0.3556/(1-0.9646)))

# plot simulated values

tsplot(exp(sv.sim[1:250,]),main="Simulated values from SV model",
lty=c(1,1), lwd=c(2,3), col=c(2,5))
legend(0,0.0005,legend=c("volatility","squared returns"),
lty=c(1,1), col=c(2,5), lwd=c(2,3))

# quasi maximum likelihood estimation
# note: sandwhich covariance matrix is not computed with SsfFit
#
sv.start = c(-0.3,log(0.03),0.9)
names(sv.start) = c("g","ln.sigma2","phi")
low.vals = c(-Inf,-Inf,-0.999)
up.vals = c(Inf,Inf,0.999)
sv.mle = SsfFit(sv.start,sv.sim[,2],"sv.mod",
lower=low.vals,upper=up.vals)
sv.mle
summary(sv.mle)
exp(sv.mle$parameters[2])

# compute filtered and smoothed estimates of volatility
# and plot with actual volatility

ssf.sv = sv.mod(sv.mle$parameters)
filteredEst.sv = SsfMomentEst(sv.sim[,2],ssf.sv,task="STFIL")
smoothedEst.sv = SsfCondDens(sv.sim[,2],ssf.sv,task="STSMO")

tsplot(sv.sim[1:250,1],filteredEst.sv$state.moment[1:250],
smoothedEst.sv$state[1:250],main="Log-Volatility",
lty=c(1,1,1),col=c(1,2,5), lwd=c(2,2,2))
legend(0.25,-11,legend=c("actual","filtered","smoothed"),
lty=c(1,1,1),col=c(1,2,5), lwd=c(2,2,2))

# quasi maximum likelihood estimation
# note: sandwhich covariance matrix can be computed with SsfFitFast
# note: SsfFitFast minimizes -1*(average log likelihood) so results
# may be slightly different from SsfFit

sv.start = c(-0.3,log(0.03),0.9)
names(sv.start) = c("g","ln.sigma2","phi")
low.vals = c(-Inf,-Inf,-0.999)
up.vals = c(Inf,Inf,0.999)
sv.qmle = SsfFitFast(sv.start,sv.sim[,2],"sv.mod",
lower=low.vals,upper=up.vals)
sv.qmle
names(sv.qmle)
summary(sv.qmle,method="qmle")
# qmle for sig2
sig2.qmle = exp(coef(sv.qmle)[2])
# delta method se for sig2
dg = exp(coef(sv.qmle)[2])
se.sig2 = sqrt(dg%*%vcov(sv.qmle,method="qmle")[2,2]%*%dg)
rbind(exp(coef(sv.qmle)[2]),se.sig2)

# compute filtered and smoothed estimates of volatility
# and plot with actual volatility

ssf.sv = sv.mod(sv.qmle$parameters)
filteredEst.sv = SsfMomentEst(sv.sim[,2],ssf.sv,task="STFIL")
smoothedEst.sv = SsfCondDens(sv.sim[,2],ssf.sv,task="STSMO")

tsplot(sv.sim[1:250,1],filteredEst.sv$state.moment[1:250],
smoothedEst.sv$state[1:250],main="Log-Volatility",
lty=c(1,2,3),col=c(1,2,5), lwd=c(2,2,2))
legend(0.25,-11,legend=c("actual","filtered","smoothed"),
lty=c(1,2,3),col=c(1,2,5), lwd=2)



#
# 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))