# stsmExamples.ssc			structural time series examples
#
# author: Eric Zivot
# created: December 28, 2005
# updated: January 5, 2006
#
# Comments:
#
# 1. Example files for specifying and estimating structural time series models (stsm's)
#    as described in the STAMP 6 manual by Koopman, Harvey, Doornik and Shephard
# 2. Examples require S-PLUS 7 and S+FinMetrics 2.0
# 3. Example scripts are described in detail in the technical document
#    "Analysis of Structural Time Series Models Using SsfPack in S+FinMetrics
#    2.0" by Eric Zivot, which may be downloaded from
#    http://faculty.washington.edu/ezivot/modelingFinancialTimeSeries.htm
#

#
# local level model
#

# create state space list
ssf.ll = GetSsfStsm(irregular=1, level=0.5)

# simulate and plot 250 obvs from state space
ssf.ll.sim = SsfSim(ssf.ll, n = 250, seed = 10)
tsplot(ssf.ll.sim, col=c(1,5), lty=c(1,1),
       main="Simulations from Local Level Model")
legend(0, 3, legend=c("state", "response"), col=c(1,5), lty=c(1,1))

# ML estimation of full likelihood

ssf.ll.mod = function(parm) {
## parm[1] = log(sig.w)
## parm[2] = log(sig.v)
	ssf.ll = GetSsfStsm(irregular=exp(parm[1]), 
	                    level=exp(parm[2]))
	CheckSsf(ssf.ll)
}

parm.ll.start = c(0,0)
names(parm.ll.start) = c("ln sigma.w", "ln signa.v")

ssf.ll.fit = SsfFit(parm=parm.ll.start,
                    data=ssf.ll.sim[,"response"],
                    FUN="ssf.ll.mod")
ssf.ll.fit
summary(ssf.ll.fit)
sd.hat = exp(coef(ssf.ll.fit))
sqrt.q.mle = as.vector(sd.hat[2]/sd.hat[1])

# extract smoothed state vectors
ssf.ll.hat = ssf.ll.mod(coef(ssf.ll.fit))
smoothedEst.ll = SsfCondDens(ssf.ll.sim[,"response"], ssf.ll.hat)
plot(smoothedEst.ll, layout=c(1,2))

# compare actual states with smoothed states
tsplot(ssf.ll.sim[,"state"], smoothedEst.ll$state, 
       main="Actual and Smoothed States",
       col=c(1,5), lty=c(1,1))
legend(0, 2, legend=c("Acutal", "Smoothed"), col=c(1,5), lty=c(1,1))

# ML estimation of concentrated likelihood

ssf.ll.modc = function(parm) {
## parm[1] = sqrt(q)
## 
	ssf.llc = GetSsfStsm(irregular=1, 
	                    level=parm[1])
	CheckSsf(ssf.llc)
}

parm.llc.start = 1
names(parm.llc.start) = c("sqrt.q")

ssf.llc.fit = SsfFit(parm=parm.llc.start,
                    data=ssf.ll.sim[,"response"],
                    FUN="ssf.ll.modc",
                    conc=T)
ssf.llc.fit
summary(ssf.llc.fit)

ssf.llc.hat = ssf.ll.modc(coef(ssf.llc.fit))
smoothedEst.llc = SsfCondDens(ssf.ll.sim[,"response"], ssf.llc.hat)
plot(smoothedEst.llc, layout=c(1,2))

# compare actual states with smoothed states
# they are the same as with the full loglikelihood
tsplot(ssf.ll.sim[,"state"], smoothedEst.llc$state, 
       main="Actual and Smoothed States",
       col=c(1,5), lty=c(1,1))
legend(0, 2, legend=c("Acutal", "Smoothed"), col=c(1,5), lty=c(1,1))


#
# local linear trend with quarterly stochastic dummy variable seasonal
#

ssf.lltdum = GetSsfStsm(irregular=1, 
                        level=0.5,
                        slope=0.1, 
                        seasonalDummy = c(0.2, 4))
# set initial values for the seasonal
ssf.lltdum$mSigma[6,3:5] = c(2.5, 5.0, -2.5)
ssf.lltdum

# create and plot simulated data
lltdum.sim = SsfSim(ssf.lltdum, n=80, seed=10)
seriesPlot(lltdum.sim,one.plot=F)

# Create deterministic seasonal
# must specify the initial value for the deterministic seasonal component
ssf.lltdum2 = GetSsfStsm(irregular=1, 
                        level=0.5,
                        slope=0.1, 
                        seasonalDummy = c(0, 4))
ssf.lltdum2$mSigma[6,3:5] = c(2.5, 5.0, -2.5)
ssf.lltdum2
lltdum2.sim = SsfSim(ssf.lltdum2, n=80, seed=10)
seriesPlot(lltdum2.sim,one.plot=F)