## StateSpaceSurvey.ssc		script file for examples used in survey paper
##
## author: Eric Zivot
## created: August 5, 2002
## revised: 
## October 11, 2002
##		Adjusted comments
##		Modified code to utilize new method functions for objects of class "SsfFit"
##		Added TVP model example using exchange rate regression
##

#############################################################################################

###################################################################
##				Traditional time series models
###################################################################

##
## AR(1) model
##

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

## simulate 250 observations
##
set.seed(443)
ar1.sim = SsfSim(ssf.ar1,n=250)
tsplot(ar1.sim[,1])

##
## Beveridge-Nelson decomposition of ARIMA(1,1,0)
##

## create ARIMA(1,1,0) process with zero initial value
##
dy.ar1 = ar1.sim[,"response"]
y.arima.1 = cumsum(dy.ar1)
tsplot(y.arima.1)

## compute Kalman filter estimation of dy.ar1
##
filteredEst.ar1 = SsfMomentEst(dy.ar1,ssf.ar1,task="STFIL")
xt.t = filteredEst.ar1$state.moment

## compute BN decomposition using formula from Morley (2002)
##
F.mat = 0.75
BN.t = y.arima.1 + (F.mat/(1-F.mat))*xt.t
c.t = y.arima.1 - BN.t
par(mfrow=c(2,1))
tsplot(cbind(BN.t,y.arima.1),main="Data and Beveridge-Nelson Trend")
tsplot(c.t,main="Beveridge-Nelson Cycle")

##
## Beveridge-Nelson decomposition of ARIMA(1,1,0) with non-zero mean
##

## create ARIMA(1,1,0) process with zero initial value
##
dy.ar1 = ar1.sim[,"response"] + 0.25
y.arima.1 = cumsum(dy.ar1)
tsplot(y.arima.1)

## compute Kalman filter estimation of dy.ar1
##
filteredEst.ar1 = SsfMomentEst((dy.ar1-0.25),ssf.ar1,task="STFIL")
xt.t = filteredEst.ar1$state.moment

## compute BN decomposition using formula from Morley (2002)
##
F.mat = 0.75
## note: BN trend here contains deterministic terms
BN.t = y.arima.1 + (F.mat/(1-F.mat))*xt.t
c.t = y.arima.1 - BN.t
par(mfrow=c(2,1))
tsplot(cbind(BN.t,y.arima.1),main="Data and Beveridge-Nelson Trend")
tsplot(c.t,main="Beveridge-Nelson Cycle")

##
## estimation of ARIMA(2,1,2) for log real GDP and BN decomposition
## Ref: Morley, Nelson and Zivot (2002)
##
lny = as.matrix(lngdpq.df)*100
nobs = nrow(lny)
dlny = diff(lny)

## demean data
##
dlny.dm = dlny - mean(dlny)

## create timeSeries object
##
td = timeSeq(from="1/1/1947",to="4/1/1998",by="quarters")
lny.ts = timeSeries(data=lny,positions=td)
lny.ts@title = "Log Postwar Quarterly Real GDP"
dlny.ts = diff(lny.ts)
dlny.ts.dm = dlny.ts - mean(dlny.ts)

## plot data
##
par(mfrow=c(2,1))
plot(lny.ts,main="Log Postwar Quarterly Real GDP")
plot(diff(lny.ts),main="Quarterly Growth Rate")

## create state space model for ARMA(2,2)
##
arma22.mod = function(parm) {
	phi.1 = parm[1]
	phi.2 = parm[2]
	theta.1 = parm[3]
	theta.2 = parm[4]
	sigma2 = exp(parm[5])	# require positive variance
	ssf.mod = GetSsfArma(ar=c(phi.1,phi.2),ma=c(theta.1,theta.2),
							sigma=sqrt(sigma2))
	CheckSsf(ssf.mod)
}

## Estimation by mle
## use arima.mle to get starting values?
##
arma22.start = c(1.34,-0.70,-1.05,0.51,-0.08)
names(arma22.start) = c("phi.1","phi.2","theta.1","theta.2","ln.sigma2")
arma22.mle = SsfFit(arma22.start,dlny.ts.dm,"arma22.mod")
arma22.mle
summary(arma22.mle)
exp(arma22.mle$parameters["ln.sigma2"])

## residual diagnostics from Kalman filter
##	show ACF plot: option 6
##
kf.arma22 = KalmanFil(dlny.ts.dm,
arma22.mod(arma22.mle$parameters))
class(kf.arma22)
names(kf.arma22)
plot(kf.arma22)

## BN decomposition
##
filteredEst.arma22 = SsfMomentEst(dlny.ts.dm,
arma22.mod(arma22.mle$parameters),task="STFIL")
at.t = filteredEst.arma22$state.moment
T.mat = ssf.arma22$mPhi[1:3,1:3]
tmp = t(T.mat%*%solve((diag(3)-T.mat))%*%t(at.t))
BN.t = lny[2:nobs,] + tmp[,1]
c.t = lny[2:nobs,] - BN.t

## create timeSeries objectsand plot
## 
BN.ts = timeSeries(data=BN.t,positions=filteredEst.arma22$positions)
c.ts = lny.ts - BN.ts
par(mfrow=c(2,1))
plot(lny.ts,BN.ts,main="Log Real GDP and BN Trend")
plot(c.ts,main="BN Cycle")

## remark: use Andrew's shadePlot function here?

##
## VAR(1) model
##

A.mat = matrix(c(0.75, 0, 0, 0.5), 2, 2)
Sigma.mat = matrix(c(1, 0.25, 0.25, 1), 2, 2)
Omega.mat = matrix(0, 4, 4)
Omega.mat[1:2, 1:2] = Sigma.mat
## solve for initial variance of stationary part
##
vecS = as.vector(Sigma.mat)
vecP = solve(diag(4) - kronecker(A.mat, A.mat)) %*% vecS
P.VAR1 = matrix(vecP, 2, 2)
Sigma = matrix(0, 3, 2)
Sigma[1:2, 1:2] = P.VAR1
## create state space
##
ssf.VAR1 = list(mPhi=rbind(A.mat,diag(2)),
mOmega=Omega.mat,
mSigma = Sigma)
ssf.VAR1

## simulate from VAR1
##
set.seed(112)
VAR1.sim = SsfSim(ssf.VAR1,n=250)

##
## VAR(1) for 1st difference
##
dY.VAR1 = VAR1.sim[,3:4]
Y.VAR1 = cumsum(dY.VAR1)
tsplot(Y.VAR1)

## compute Kalman filter estimation of dy.ar1
##
filteredEst.VAR1 = SsfMomentEst(dY.VAR1,ssf.VAR1,task="STFIL")
Xt.t = filteredEst.VAR1$state.moment

## compute BN decomposition using formula from Morley (2002)
##
F.mat = A.mat
BN.t = t(Y.VAR1) + F.mat%*%solve((diag(2)-F.mat))%*%t(Xt.t)
c.t = t(Y.VAR1) - BN.t
BN.t = t(BN.t)
c.t = t(c.t)
par(mfrow=c(2,1))
tsplot(cbind(BN.t,Y.VAR1),main="Data and Beveridge-Nelson Trend")
tsplot(c.t,main="Beveridge-Nelson Cycle")

#############################################################################################
#					Regression Models
#############################################################################################

##
## CAPM regression with fixed parameters
##
X.mat = cbind(1,as.matrix(seriesData(excessReturns.ts[,"SP500"])))
msft.ret = excessReturns.ts[,"MSFT"]
ssf.reg = GetSsfReg(X.mat)
ssf.reg

## Kalman Filter estimation of state variables
## initial state is zero
## note: SsfFit is not used
##
filteredEst.reg = SsfMomentEst(msft.ret,ssf.reg,task="STFIL")
class(filteredEst.reg)
names(filteredEst.reg)
filteredEst.reg$state.moment
colIds(filteredEst.reg$state.moment) = c("intercept","slope")
plot(filteredEst.reg)

## the state vector gives the recursive least squares estimates
##
rls.coef = timeSeries(pos=filteredEst.reg$positions,
data=filteredEst.reg$state.moment)
seriesPlot(rls.coef,one.plot=F,strip.text=c("alpha","beta"))

## full sample OLS estimates are last row of rls estimates
##
filteredEst.reg$state.moment[131,]

## OLS estimates should match last row of xt.t
##
ols.fit = OLS(MSFT~SP500,data=excessReturns.ts)
coef(ols.fit)

## compare with FinMetrics RLS function
##
rls.fit = RLS(MSFT~SP500,data=excessReturns.ts)

## note: rls matches exactly with the KF after observation 2
## KF can compute estimates for observations 1 and 2 due to the
## exact initialization of the filter. 
##
coef(rls.fit)

## compute recursive residuals
##
kf.reg = KalmanFil(msft.ret,ssf.reg)
class(kf.reg)
names(kf.reg)
w.t = kf.reg$std.innov[-c(1,2)]

## compute cusum test
##
cusum.t = cumsum(w.t)/stdev(w.t)
nobs = length(cusum.t)
tmp = 0.948*sqrt(nobs)
upper = seq(tmp,3*tmp,length=nobs)
lower = seq(-tmp,-3*tmp,length=nobs) 
tsplot(cusum.t,lower,upper,lty=c(1,2,2))

## create time series object and plot
##
tmp.ts = timeSeries(pos=kf.reg$positions[-c(1,2)],
data=cbind(cusum.t,upper,lower))
plot(tmp.ts,reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,2,2)))

## compute smoothed residuals for structural change diagnostics
##
ks.reg = KalmanSmo(kf.reg,ssf.reg)
class(ks.reg)
names(ks.reg)
colIds(ks.reg$aux.residuals) = c("intercept","slope","excess returns")
plot(ks.reg)

t.stats = timeSeries(data=ks.reg$aux.residuals[1:131,1:2],
pos=ks.reg$positions)
plot(t.stats,reference.grid=F,
plot.args=list(lty=c(1,3)))
legend(0,0.2,legend=c("t-alpha","t-beta"),
lty=c(1,3))

## compute ols residuals
##
ssf.reg$mSigma[3,] = filteredEst.reg$state.moment[131,]
ssf.reg$mSigma[1,1]=ssf.reg$mSigma[2,2]=0
kf.ols = KalmanFil(msft.ret,ssf.reg)
res.ols = kf.ols$innov

## change state space parameters to have slope time varying
##
ssf.reg$mOmega[2,2] = 0.05^2


###################################################################
# 						Time Varying Parameter Models
###################################################################

##
## CAPM with time varying parameters
##

## create state space: all parameters time varying
##
tvp.mod = function(parm,mX=NULL) {
	parm = exp(parm)
	ssf.tvp = GetSsfReg(mX=mX)
	diag(ssf.tvp$mOmega) = parm
	CheckSsf(ssf.tvp)
}

## estimate using mle
##
tvp.start = c(0,0,0)
names(tvp.start) = c("ln(s2.alpha)","ln(s2.beta)","ln(s2.y)")
tvp.mle = SsfFit(tvp.start,msft.ret,"tvp.mod",mX=X.mat)
class(tvp.mle)
names(tvp.mle)
tvp.mle
summary(tvp.mle)

## compute state space with estimated parameters
##
ssf.tvp = tvp.mod(tvp.mle$parameters,mX=X.mat)

## compute mles of transformed parameters
## and delta method variances
##
tvp.mle$parameters = exp(tvp.mle$parameters/2)		# std deviations
names(sigma.mle) = c("s.alpha","s.beta","s.y")
dg = diag(tvp.mle$parameters/2)
tvp.mle$vcov = dg%*%tvp.mle$vcov%*%dg
summary(tvp.mle)

## filtered estimates
##
filteredEst.tvp = SsfMomentEst(msft.ret,ssf.tvp,task="STFIL")
class(filteredEst.tvp)
names(filteredEst.tvp)

## plot filtered moment estimates
##
colIds(filteredEst.tvp$state.moment) = c("intercept","slope")
## names(filteredEst.tvp$response.moment) = "Expected returns"
plot(filteredEst.tvp,strip.text=c("alpha(t)",
"beta(t)","Expected returns"),main="Filtered Estimates")

## plot filtered estimates with 95% error bands
## note: need to exclude some initial observations
##
filteredEst.ts = timeSeries(data=filteredEst.tvp$state.moment,
positions=filteredEst.tvp$positions)
filteredSd.ts = sqrt(timeSeries(data=filteredEst.tvp$state.variance,
positions=filteredEst.tvp$positions))
upper = filteredEst.ts + 2*filteredSd.ts
lower = filteredEst.ts - 2*filteredSd.ts
par(mfrow=c(2,1))
plot(filteredEst.ts[,1],upper[,1],lower[,1],
main="Filtered estimates of alpha(t)",
plot.args=list(lty=c(1,2,2)))
plot(filteredEst.ts[,2],upper[,2],lower[,2],
main="Filtered estimates of beta(t)",
plot.args=list(lty=c(1,2,2)))

## compute smoothed state estimates
##
smoothedEst.tvp = SsfCondDens(msft.ret,ssf.tvp,task="STSMO")
class(smoothedEst.tvp)
names(smoothedEst.tvp)
colIds(smoothedEst.tvp$state) = c("alpha(t)","beta(t)")
plot(smoothedEst.tvp,strip.text=c("alpha(t)",
"beta(t)","Expected returns"),main="Smoothed Estimates")

## compute smoothed state estimates with estimated variances
## and plot with 95% confidence intervals
## note: need to exclude some initial observations
##
SmoothedEst.tvp = SsfMomentEst(msft.ret,ssf.tvp,task="STSMO")
SmoothedEst.ts = timeSeries(data=SmoothedEst.tvp$state.moment,
pos=SmoothedEst.tvp$positions)
SmoothedSd.ts = sqrt(timeSeries(data=SmoothedEst.tvp$state.variance,
pos=SmoothedEst.tvp$positions))
upper = SmoothedEst.ts + 2*SmoothedSd.ts
lower = SmoothedEst.ts -2*SmoothedSd.ts
par(mfrow=c(2,1))
plot(SmoothedEst.ts[,1],upper[,1],lower[,1],
main="Smoothed estimates of alpha(t)",
plot.args=list(lty=c(1,2,2)))
plot(SmoothedEst.ts[,2],upper[,2],lower[,2],
main="Smoothed estimates of beta(t)",
plot.args=list(lty=c(1,2,2)))

## tests for time varying parameters
## note: cannot use simple LR statistic
##

##
## exchange rate regressions with time varying parameters
##

##
## bivariate VAR model with time-varying parameters
## Use Jeff's example here.
##


###############################################################################
##					Unobserved Component Models
###############################################################################

##
## simple signal plus noise model
##

## construct state space form
##
ssf.ll = GetSsfStsm(irregular=1,level=0.5)
class(ssf.ll)
ssf.ll

## simulate 250 observations and plot
##
set.seed(112)
ll.sim = SsfSim(ssf.ll,n=250)
class(ll.sim)
colIds(ll.sim)
tsplot(ll.sim,main="Simulated observations from local level model")
legend(0,4,legend=c("State","Response"),lty=1:2)

##
## Clark model. Ref Clark (1987), JPE
##

## create state space using parameters based on estimates in MNZ (2002)
##
delta = 0.8119
phi1 = 1.5303
phi2 = -0.6097
sigma.v = 0.6893
sigma.w = 0.6199
bigV = diag(c(sigma.v^2,sigma.w^2))
Omega = matrix(0,4,4)
Omega[1:2,1:2] = bigV
a1 = matrix(0,3,1)
## solve for initial variance of stationary part
##
bigF = matrix(c(phi1,1,phi2,0),2,2)
vecV = c(sigma.w^2,0,0,0)
vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
P.ar2 = matrix(vecP,2,2)
Sigma = matrix(0,4,3)
Sigma[1,1] = -1
Sigma[2:3,2:3] = P.ar2
## create state space list
##
ssf.clark = list(mDelta=c(delta,0,0,0),
mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)),
mOmega=Omega,
mSigma = Sigma)
## print out state space
##
ssf.clark
## create state space object
##
ssf.clark = CheckSsf(ssf.clark)

## simulate from model
##
set.seed(569)
clark.sim = SsfSim(ssf.clark,n=250)
par(mfrow=c(2,1))
tsplot(clark.sim[,c(1,4)],main="Simulated response and trend")
tsplot(clark.sim[,2],main="Simulated cycle")

##
## estimate Clark model for U.S. postwar quarterly real GDP
## ref: Morley, Nelson and Zivot (2002)
##

## create state space model for estimation
##
Clark.mod = function(parm) {
	delta = parm[1]
	phi1 = parm[2]
	phi2 = parm[3]
	sigma2.v = exp(parm[4])
	sigma2.w = exp(parm[5])
	bigV = diag(c(sigma2.v,sigma2.w))
	Omega = matrix(0,4,4)
	Omega[1:2,1:2] = bigV
	a1 = matrix(0,3,1)
# solve for initial variance of stationary part
	bigF = matrix(c(phi1,1,phi2,0),2,2)
	vecV = c(sigma2.w,0,0,0)
	vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
	P.ar2 = matrix(vecP,2,2)
	Sigma = matrix(0,4,3)
	Sigma[1,1] = -1
	Sigma[2:3,2:3] = P.ar2
# create state space list
	ssf.mod = list(mDelta=c(delta,0,0,0),
	mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)),
	mOmega=Omega,
	mSigma = Sigma)
	CheckSsf(ssf.mod)
}

## estimate by mle
## 
Clark.start=c(0.81,1.53,-0.61,-0.74,-0.96)
names(Clark.start) = c("mu","phi.1","phi.2",
"ln.sigma2.v","ln.sigma2.w")
Clark.mle = SsfFit(Clark.start,lny.ts,"Clark.mod")

## create state space from fitted parameters
##
ssf.Clark = Clark.mod(Clark.mle$parameters)

## estimated parameters and asymptotic standard errors
##
summary(Clark.mle)

## estimated standard deviations and delta method standard errors
##
sqrt(exp(Clark.mle$parameters[4:5]))

## compute filtered trend and cycle
## use plot method to show filtered estimates
##

filteredEst.Clark = SsfMomentEst(lny.ts,ssf.Clark,task="STFIL")
names(filteredEst.Clark)
colIds(filteredEst.Clark$state.moment) = 
c("trend(t)","cycle(t)","cycle(t-1)")
plot(filteredEst.Clark,
strip.text=c("trend(t)","cycle(t)","cycle(t-1)","y(t)"),
main="Filtered estimates")

## create timeSeries and plot
## 
trend.filter = timeSeries(data=filteredEst.Clark$state.moment[,1],
positions=filteredEst.Clark$positions)
cycle.filter = timeSeries(data=filteredEst.Clark$state.moment[,2],
positions=filteredEst.Clark$positions)
par(mfrow=c(2,1))
plot(trend.filter,lny.ts,
main="Log Real GDP and Filtered Trend from Clark Model")
plot(cycle.filter,main="Filtered Cycle from Clark Model")

## compute smoothed trend and cycle and plot
##
smoothedEst.Clark = SsfMomentEst(lny.ts,ssf.Clark,task="STSMO")
trend.smooth = timeSeries(data=smoothedEst.Clark$state.moment[,1],
positions=smoothedEst.Clark$positions)
cycle.smooth = timeSeries(data=smoothedEst.Clark$state.moment[,2],
positions=smoothedEst.Clark$positions)
par(mfrow=c(2,1))
plot(trend.smooth,lny.ts,
main="Log Real GDP and Smoothed Trend from Clark Model")
plot(cycle.smooth,main="Smoothed Cycle from Clark Model")

## plot smoothed cycle with 95% error bands
##
trend.sd = sqrt(timeSeries(data=smoothedEst.Clark$state.variance[,1],
positions=smoothedEst.Clark$positions))
cycle.sd= sqrt(timeSeries(data=smoothedEst.Clark$state.variance[,2],
positions=smoothedEst.Clark$positions))
upper = cycle.smooth + 2*cycle.sd
lower = cycle.smooth - 2*cycle.sd
par(mfrow=c(1,1))
plot(cycle.smooth,upper,lower,
plot.args=list(lty=c(1,2,2)))

## Do diagnostics on standardized innovations from Kalman filter
##
kf.Clark = KalmanFil(lny.ts,ssf.Clark)
plot(kf.Clark)
summaryStats(kf.Clark$std.innov)
autocorrTest(kf.Clark$std.innov)

## KPSS specification test
##
stationaryTest(lny.ts,trend="ct")

## Nyblom specification test
##
n = nrow(msft.ret)
ssf.Clark$mOmega[1,1] = 0
ssf.Clark$mSigma[1,1] = 0
ssf.Clark$mSigma[4,1] = smoothedEst.Clark$state.moment[n,1]
kf.Clark = KalmanFil(msft.ret[n:1],ssf.Clark)
test.stat = sum(cumsum(kf.Clark$innov)^2)

##
## MNZ model: Ref. Morley, Nelson and Zivot (2002), "Why are
## Beveridge-Nelson and Unobserved Components Decompositions of GDP
## so Different?", forthcoming in Review of Economics and Statistics
##

## Note: parameters are taken from Table 4 of 
## Morley, Nelson and Zivot (2002)
##
delta = 0.8156
phi1 = 1.3419
phi2 = -0.7060
sigma.v = 1.2368
sigma.w = 0.7485
sigma.vw = -0.8389
bigV = matrix(c(sigma.v^2,sigma.vw,sigma.vw,sigma.w^2),2,2)
Omega = matrix(0,4,4)
Omega[1:2,1:2] = bigV
a1 = matrix(0,3,1)
# solve for initial variance of stationary part
bigF = matrix(c(phi1,1,phi2,0),2,2)
vecV = c(sigma.w^2,0,0,0)
vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
P.ar2 = matrix(vecP,2,2)
Sigma = matrix(0,4,3)
Sigma[1,1] = -1
Sigma[2:3,2:3] = P.ar2
ssf.mnz= list(mDelta=c(delta,0,0,0),
mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)),
mOmega=Omega,
mSigma = Sigma)
ssf.mnz

## simulate from model
##
set.seed(569)
mnz.sim = SsfSim(ssf.mnz,n=250)
par(mfrow=c(2,1))
tsplot(mnz.sim[,c(1,4)],main="Simulated response and trend")
tsplot(mnz.sim[,2],main="Simulated cycle")

## estimate model using US postwar quarterly real GDP
##
MNZ.mod = function(parm) {
	delta = parm[1]
	phi1 = parm[2]
	phi2 = parm[3]
	sigma.v = exp(parm[4])
	sigma.w = exp(parm[5])
	rho.vw = parm[6]
	sigma.vw = sigma.v*sigma.w*rho.vw
	bigV = matrix(c(sigma.v^2,sigma.vw,sigma.vw,sigma.w^2),2,2)
	Omega = matrix(0,4,4)
	Omega[1:2,1:2] = bigV
	a1 = matrix(0,3,1)
# solve for initial variance of stationary part
	bigF = matrix(c(phi1,1,phi2,0),2,2)
	vecV = c(sigma.w^2,0,0,0)
	vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
	P.ar2 = matrix(vecP,2,2)
	Sigma = matrix(0,4,3)
	Sigma[1,1] = -1
	Sigma[2:3,2:3] = P.ar2
	ssf.mod= list(mDelta=c(delta,0,0,0),
	mPhi=rbind(c(1,0,0),c(0,phi1,phi2),c(0,1,0),c(1,1,0)),
	mOmega=Omega,
	mSigma = Sigma)
	CheckSsf(ssf.mod)
}

## set starting values
##
MNZ.start=c(0.81,1.34,-0.70,0.21,-0.30,-0.9)
names(MNZ.start) = c("delta","phi.1","phi.2",
"ln.sigma.v","ln.sigma.w","rho")

## estimate model by MLE
##
low.vals = c(0,0,-2,-Inf,-Inf,-0.999)
up.vals = c(2,2,2,Inf,Inf,0.999)
MNZ.mle = SsfFit(MNZ.start,lny.ts,"MNZ.mod",
lower=low.vals,upper=up.vals)

## estimated parameters and asymptotic standard errors
##
MNZ.mle$parameters
sqrt(diag(MNZ.mle$vcov))
# estimated standard deviations
exp(MNZ.mle$parameters[4:5])
# note: use delta method to get standard errors for standard deviations

# compute filtered trend and cycle
filteredEst.MNZ = SsfMomentEst(lny.ts,
MNZ.mod(MNZ.mle$parameters),task="STFIL")
trend.filter = timeSeries(data=filteredEst.MNZ$state.moment[,1],
positions=filteredEst.MNZ$positions)
cycle.filter = timeSeries(data=filteredEst.MNZ$state.moment[,2],
positions=filteredEst.MNZ$positions)

par(mfrow=c(2,1))
plot(trend.filter,lny.ts,
main="Log Real GDP and Filtered Trend from MNZ Model")
plot(cycle.filter,main="Filtered Cycle from MNZ Model")

# compute smoothed trend and cycle
smoothedEst.MNZ = SsfMomentEst(lny.ts,
MNZ.mod(MNZ.mle$parameters),task="STSMO")
trend.smooth = timeSeries(data=smoothedEst.MNZ$state.moment[,1],
positions=smoothedEst.MNZ$positions)
cycle.smooth = timeSeries(data=smoothedEst.MNZ$state.moment[,2],
positions=smoothedEst.MNZ$positions)

par(mfrow=c(2,1))
plot(trend.smooth,lny.ts,
main="Log Real GDP and Smoothed Trend from MNZ Model")
plot(cycle.smooth,main="Smoothed Cycle from MNZ Model")

# plot smoothed cycle with 95% error bands
cycle.sd= sqrt(timeSeries(data=smoothedEst.MNZ$state.variance[,2],
positions=smoothedEst.MNZ$positions))
upper = cycle.smooth + 2*cycle.sd
lower = cycle.smooth - 2*cycle.sd
par(mfrow=c(1,1))
plot(cycle.smooth,upper,lower,
plot.args=list(lty=c(1,3,3)))

# plot filtered and smoothed cycle together
plot(cycle.filter,cycle.smooth,
plot.args=list(lty=c(1,2)))
legend(0.75,-3,legend=c("filtered","smoothed"),lty=c(1,2))

# plot filtered and smoothed trend together
plot(diff(trend.filter),diff(trend.smooth))
#
# Common trend model. Ref. Hai, Mark and Wu (1997), "Understanding Spot
# and Forward Exchange Rate Regressions", Journal of Applied Econometrics.
#
# remark: The intercept parameters in HMW's model are most likely
# not identified

# Note: parameters calibrated to Pound rate in Table VI of Hai, Mark
# and Wu.
#
# drifts
mu.1 = 0.0827			# these are probably not identified
mu.2 = 0.0388
# variances and covariances
sigma.tau = 3.13		# variance of permanent trend
sigma.1 = 0.97		# variance of transitory spot innovations
sigma.2 = 0.88		# variance of transitory forward innovations
rho.12 = 0.98			# correlation between spot and forward innovations
sigma.12 = rho.12*sigma.1*sigma.2
# VAR(1) dynamics
phi.11 = 0.9
phi.12 = -0.05
phi.21 = 0.14
phi.22 = 0.71
# state space matrices
bigV = matrix(c(sigma.1^2,sigma.12,sigma.12,sigma.2^2),2,2)
Omega = matrix(0,5,5)
Omega[1,1] = sigma.tau^2
Omega[2:3,2:3] = bigV
a1 = matrix(0,3,1)
# solve for initial variance of stationary part
bigF = matrix(c(phi.11,phi.21,phi.12,phi.22),2,2)
vecV = as.vector(bigV)
vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
P.VAR1 = matrix(vecP,2,2)
Sigma = matrix(0,4,3)
Sigma[1,1] = -1
Sigma[2:3,2:3] = P.VAR1
ssf.ct2= list(mDelta=c(0,mu.1,mu.2,0,0),
mPhi=rbind(c(1,0,0),c(0,phi.11,phi.12),c(0,phi.21,phi.22),
c(1,1,0),c(1,0,1)),
mOmega=Omega,
mSigma = Sigma)
ssf.ct2

# simulate from model
set.seed(569)
ct2.sim = SsfSim(ssf.ct2,n=250)
par(mfrow=c(2,2))
tsplot(ct2.sim[,c(4,5)],main="Simulated responses")
tsplot(ct2.sim[,1],main="Simulated trend")
tsplot(ct2.sim[,2],main="Simulated cycle 1")
tsplot(ct2.sim[,3],main="Simulated cycle 2")


#
# revised parameterization for survey paper
# intercepts removed from state equation
# data taken from lexrates.dat
#

# extract data
usdm.s = lexrates.dat[,"USDMS"]*100
usdm.f = lexrates.dat[,"USDMF"]*100
usdm.fp = usdm.f - usdm.s
colIds(usdm.fp) = "USDMFP"
usdm.ds = diff(usdm.s)
colIds(usdm.ds) = "DUSDMS"
usdm.df = diff(usdm.f)
colIds(usdm.df) = "DUSDMF"
usdm.ts = seriesMerge(usdm.s,usdm.f,usdm.fp,usdm.ds,usdm.df)


#
# drifts
mu = 0.16			# mean of interest rate differential
# variances and covariances
sigma.tau = 3.40	# sd of common trend
sigma.s = 0.30	# sd of spot innov
sigma.f = 0.30	# sd of forward innov
rho.sf = 0.98		# cor bw spot and forward innov
sigma.sf = rho.sf*sigma.s*sigma.f
# VAR(1) dynamics
phi.ss = 0.75
phi.sf = 0
phi.fs = 0
phi.ff = 0.75
# state space matrices
bigV = matrix(c(sigma.s^2,sigma.sf,sigma.sf,sigma.f^2),2,2)
Omega = matrix(0,5,5)
Omega[1,1] = sigma.tau^2
Omega[2:3,2:3] = bigV
a1 = matrix(0,3,1)
# solve for initial variance of stationary part
bigF = matrix(c(phi.ss,phi.fs,phi.sf,phi.ff),2,2)
vecV = as.vector(bigV)
vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
P.VAR1 = matrix(vecP,2,2)
Sigma = matrix(0,4,3)
Sigma[1,1] = -1
Sigma[2:3,2:3] = P.VAR1
# create state space form
ssf.exrate= list(mDelta=c(0,0,0,0,mu),
mPhi=rbind(c(1,0,0),c(0,phi.ss,phi.sf),c(0,phi.fs,phi.ff),
c(1,1,0),c(1,0,1)),
mOmega=Omega,
mSigma = Sigma)
ssf.exrate

# simulate from model
set.seed(569)
exrate.sim = SsfSim(ssf.exrate,n=250)
par(mfrow=c(2,2))
tsplot(ct2.sim[,c(4,5)],main="Simulated responses")
tsplot(ct2.sim[,1],main="Simulated trend")
tsplot(ct2.sim[,2],main="Simulated spot cycle")
tsplot(ct2.sim[,3],main="Simulated forward cycle")

# run Kalman Filter

kf.exrate = KalmanFil(exrate.sim[,4:5],ssf.exrate)

# create function for estimate of model
exrate.mod = function(parm) {
	mu = parm[1]					# mean of interest rate differential
	sigma.tau = exp(parm[2])	# sd of common trend
	sigma.s = exp(parm[3])	 	# sd of spot innov
	sigma.f = exp(parm[4])		# sd of forward innov
	rho.sf = parm[5]				# cor bw spot and forward innov
	sigma.sf = rho.sf*sigma.s*sigma.f
	phi.ss = parm[6]
	phi.sf = parm[7]
	phi.fs = parm[8]
	phi.ff = parm[9]
# state space matrices
	bigV = matrix(c(sigma.s^2,sigma.sf,sigma.sf,sigma.f^2),2,2)
	Omega = matrix(0,5,5)
	Omega[1,1] = sigma.tau^2
	Omega[2:3,2:3] = bigV
	a1 = matrix(0,3,1)
# solve for initial variance of stationary part
	bigF = matrix(c(phi.ss,phi.fs,phi.sf,phi.ff),2,2)
	vecV = as.vector(bigV)
	vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
	P.VAR1 = matrix(vecP,2,2)
	Sigma = matrix(0,4,3)
	Sigma[1,1] = -1
	Sigma[2:3,2:3] = P.VAR1
# create state space form
	ssf.mod = list(mDelta=c(0,0,0,0,mu),
	mPhi=rbind(c(1,0,0),c(0,phi.ss,phi.sf),c(0,phi.fs,phi.ff),
	c(1,1,0),c(1,0,1)),
	mOmega=Omega,
	mSigma = Sigma)
	CheckSsf(ssf.mod)
}
exrate.start=c(0.16,log(3.4),log(0.4),log(0.3),0.98,0.75,0,0,0.75)
names(exrate.start) = c("mu","ln.sigma.z","ln.sigma.s","ln.sigma.f",
"rho","phi.ss","phi.sf","phi.fs","phi.ff")

# estimate model by MLE
low.vals = c(-Inf,-Inf,-Inf,-Inf,-0.999,-2,-2,-2,-2)
up.vals = c(Inf,Inf,Inf,Inf,0.999,2,2,2,2)
exrate.mle = SsfFit(exrate.start,usdm.ts[,1:2],"exrate.mod",
lower=low.vals,upper=up.vals)

exrate.mle$parameters
sqrt(diag(exrate.mle$vcov))
# estimated standard deviations
exp(exrate.mle$parameters[2:4])
# note: use delta method to get standard errors for standard deviations

# compute filtered trend and cycle
filteredEst.exrate = SsfMomentEst(usdm.ts[,1:2],
exrate.mod(exrate.mle$parameters),task="STFIL")
trend.filter = timeSeries(data=filteredEst.exrate$state.moment[,1],
positions=filteredEst.exrate$positions)
cycle.filter = timeSeries(data=filteredEst.exrate$state.moment[,2:3],
positions=filteredEst.exrate$positions)
exrate.filter = timeSeries(data=filteredEst.exrate$response.moment,
positions=filteredEst.exrate$positions)
fp.filter = exrate.filter[,2]-exrate.filter[,1]

# extract risk premium
predictedEst.exrate = SsfMomentEst(usdm.ts[,1:2],
exrate.mod(exrate.mle$parameters),task="STPRED")
predicted.spot = timeSeries(data=predictedEst.exrate$response.moment[,1],
positions=predictedEst.exrate$positions)
rpEst = usdm.ts[,2]-tslag(predicted.spot,-1)
fe = usdm.ts[,1] - tslag(usdm.ts[,2])

par(mfrow=c(2,1))
plot(trend.filter,usdm.ts[,1],
main="Spot Rate and Filtered Trend")
plot(cycle.filter,main="Filtered Cycle")

par(mfrow=c(2,1))

# compute smoothed trend and cycle
smoothedEst.exrate = SsfMomentEst(usdm.ts[,1:2],
exrate.mod(exrate.mle$parameters),task="STSMO")
trend.smooth = timeSeries(data=smoothedEst.exrate$state.moment[,1],
positions=smoothedEst.exrate$positions)
cycle.smooth = timeSeries(data=smoothedEst.exrate$state.moment[,2:3],
positions=smoothedEst.exrate$positions)

par(mfrow=c(2,1))
plot(trend.smooth,usdm.ts[,1],
main="Spot Rate and smoothed Trend")
plot(cycle.smooth,main="smoothed Cycle")



#
# Bivariate Clark model: Ref. Kim and Nelson (1999), pags 37-44
#
sigma.tau = 0.0049
sigma.gamma = 0.0003
sigma.k = 0.0015
sigma.x = 0.0067
sigma.c = 0.0003
Omega = diag(c(sigma.tau^2,sigma.gamma^2,sigma.k^2,sigma.x^2,
0,0,0,sigma.c^2))

phi.1 = 1.4386
phi.2 = -0.5174
alpha.0 = -0.3368
alpha.1 = -0.1635
alpha.2 = -0.0720

a1 = matrix(0,6,1)

# solve for initial variance of stationary part
bigF = matrix(c(phi.1,1,phi.2,0),2,2)
vecV = as.vector(diag(c(sigma.c^2,0)))
vecP = solve(diag(4)-kronecker(bigF,bigF))%*%vecV
P0.stat = matrix(vecP,2,2)
P1.stat = bigF%*%P.stat
P.stat= matrix(0,3,3)
P.stat[1:2,1:2]=P0.stat
P.stat[1:2,3] = P1.stat[,2]
P.stat[3,1] = P.stat[1,3]
P.stat[3,2] = P.stat[2,3]
P.stat[3,3]=P.stat[1,1]
Sigma = rbind(-diag(6),matrix(0,1,6))
Sigma[4:6,4:6] = P.stat
ssf.bclark= list(mPhi=rbind(c(1,1,0,0,0,0),c(0,1,0,0,0,0),
c(0,0,1,0,0,0),c(0,0,0,phi.1,phi.2,0),c(0,0,0,1,0,0),
c(0,0,0,0,1,0),c(1,0,0,1,0,0),c(0,0,1,alpha.0,alpha.1,alpha.2)),
mOmega=Omega,
mSigma = Sigma)
ssf.bclark

set.seed(569)
bclark.sim = SsfSim(ssf.bclark,n=250)

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

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)))
tsplot(exp(sv.sim),main="Simulated values from SV model")
legend(0,0.001,legend=c("volatility","squared returns"),
lty=c(1,2))

sv.start = c(-0.3,log(0.03),0.9)
names(sv.start) = c("g","ln.sigma2","phi")
# estimate model by MLE
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)

##
## Term Structure Examples
##

# show fama.bliss data

plot(fama.bliss,plot.args=list(lty=c(1,2,3,4)))
legend(0,0.16,legend=colIds(fama.bliss),lty=c(1,2,3,4))

# estimate Vasicek model using Fama-Bliss data

vasicek.ssf = function(param, tau=NULL, freq=1/52)
{
## 1. Check for valid input.
  if (length(param) < 5)
    stop("param must have length greater than 4.")
  N = length(param) - 4
  if (length(tau) != N)
    stop("Length of tau is inconsistent with param.")
  names(param) = c("kappa", "theta", "sigma", "lambda", 
                    paste("s",1:N,sep=""))
## 2. Extract parameters and impose constraints.
  Kappa = exp(param[1])  ## Kappa > 0
  Theta = exp(param[2])  ## Theta > 0
  Sigma = exp(param[3])  ## Sigma > 0
  Lamda = param[4]
  Var   = exp(param[1:N+4])
## 3. Compute Gamma, A, and B.
  Gamma = Theta + Sigma * Lamda / Kappa - Sigma^2 / (2 * Kappa^2)
  B = (1 - exp(-Kappa * tau)) / Kappa
  A = Gamma * (B - tau) - Sigma^2 * B^2 / (4 * Kappa)
## 4. Compute a, b, and Phi.
  a = Theta * (1 - exp(-Kappa * freq))
  b = exp(-Kappa * freq)
  Phi = (Sigma^2 / (2 * Kappa)) * (1 - exp(-2 * Kappa * freq))
## 5. Compute the state space form.
  mDelta = matrix(c(a, -A/tau), ncol=1)
  mPhi = matrix(c(b, B/tau), ncol=1)
  mOmega = diag(c(Phi, Var^2))
## 6. Duan and Simonato used this initial setting.
  A0 = Theta
  P0 = Sigma * Sigma / (2*Kappa)
  mSigma = matrix(c(P0, A0), ncol=1)
## 7. Return state space form.
  ssf.mod = list(mDelta=mDelta, mPhi=mPhi, mOmega=mOmega, mSigma=mSigma)
  CheckSsf(ssf.mod)
}

##
## Estimate model using Fama-Bliss Data
##
start.vasicek = c(log(0.1), log(0.06), log(0.02), 0.3, log(0.003), 
log(0.001), log(0.003), log(0.01))
names(start.vasicek) = c("ln.kappa","ln.theta","ln.sigma","lamda",
							"ln.sig.3M","ln.sig.6M","ln.sig.12M","ln.sig.60M")					
start.tau = c(0.25, 0.5, 1, 5)	
ans.vasicek = SsfFit(start.vasicek, fama.bliss, vasicek.ssf, 
               tau=start.tau, freq=1/12, trace=T,
               control=nlminb.control(abs.tol=1e-6, rel.tol=1e-6, 
               x.tol=1e-6, eval.max=1000, iter.max=500))
ssf.fit = vasicek.ssf(ans.vasicek$parameters,tau=start.tau,freq=1/12)
# tmp = ans.vasicek$parameters
# ans.vasicek$parameters = tmp
ans.vasicek$parameters[-4] = exp(ans.vasicek$parameters[-4])
names(ans.vasicek$parameters) = c("kappa","theta","sigma","lamda",
							"sig.3M","sig.6M","sig.12M","sig.60M")
dg = ans.vasicek$parameters; dg[4] = 1
ans.vasicek$vcov = diag(dg) %*% ans.vasicek$vcov %*% diag(dg)
summary(ans.vasicek)

ssf.fit = vasicek.ssf(ans.vasicek$parameters,tau=start.tau,freq=1/12)

##
## 9. smoothed estimates of term-structure
##
m.s = SsfCondDens(fama.bliss, ssf.fit)
r.s = timeSeries(data=m.s$state,pos=m.s$positions)
y.s = timeSeries(data=m.s$response,pos=m.s$positions)

plot(r.s,main="Smoothed Estimate of Short Rate")
par(mfrow=c(2,2))
plot(fama.bliss[,1],y.s[,1],main="Actual vs. Smoothed: 3M")
plot(fama.bliss[,2],y.s[,2],main="Actual vs. Smoothed: 6M")
plot(fama.bliss[,3],y.s[,3],main="Actual vs. Smoothed: 12M")
plot(fama.bliss[,4],y.s[,4],main="Actual vs. Smoothed: 60M")

##
## 10. Diagnostics on standardized innovations
##

summaryStats(KalmanFil(fama.bliss,ssf.fit)$std.innov)
autocorTest(KalmanFil(fama.bliss,ssf.fit)$std.innov)

# estimate Chen-Scott two factor model

chen.ssf =  
#
# Gaussian Affine Term Structure Model
# See: Duan and Simonato (1999), "Estimating Exponential-Affine Term 
#      Structure Models by Kalman Filter", Review of Quantitative 
#      Finance and Accounting, 13 (2), 111-135.
#
function(param, tau=NULL, freq=1/52, n.factors=2)
{
## 1. Check for valid input.
  if (length(param) <= n.factors*4)
    stop("param must have length greater than ", 4*n.factors)
  N = length(param) - 4*n.factors
  if (length(tau) != N)
    stop("Length of tau is inconsistent with param.")
## 2. Extract parameters and impose constraints.
  Kappa = exp(param[1:n.factors])               ## Kappa > 0
  Theta = exp(param[n.factors+1:n.factors])     ## Theta > 0
  Sigma = exp(param[2*n.factors+1:n.factors])   ## Sigma > 0
  Lambda = param[3*n.factors+1:n.factors]       ## Lambda < 0
  Var = exp(param[1:N+n.factors*4])
## 3. Compute Gamma, A, and B.
  Gamma = sqrt((Kappa + Lambda)^2 + 2 * Sigma^2)
  A = B = matrix(0, N, n.factors) 
  for (i in 1:n.factors) {
    denom = (Kappa[i] + Lambda[i] + Gamma[i])*(exp(Gamma[i]*tau) - 1) + 2 * Gamma[i]
    A[,i] = (2*Gamma[i]*exp((Kappa[i]+Lambda[i]+Gamma[i])*tau/2)/denom)^(2*Kappa[i]*Theta[i]/(Sigma[i]^2))
    B[,i] = 2*(exp(Gamma[i]*tau)-1)/denom
  }
  AA = rep(1, length(tau)) 
  for (i in 1:n.factors) {
    AA = AA * A[,i]
  } 
## 4. Compute a, b, Phi, and vS.
  a = Theta * (1 - exp(-Kappa * freq))
  b = diag(exp(-Kappa * freq))
  Phi = Theta * (Sigma^2) / (2 * Kappa) * (1 - exp(-Kappa*freq))^2
  mAffine = Sigma^2 / Kappa * (exp(-Kappa*freq) - exp(-2*Kappa*freq))
## 5. Computer the state space form.
  mDelta = matrix(c(a, -log(AA)/tau), ncol=1)
  mPhi = rbind(b, B/tau)
  mOmega = diag(c(Phi, Var^2))
## 6. Duan and Simonato used this initial setting.
  A0 = Theta
  P0 = diag(Sigma^2 * Theta / (2*Kappa))
  mSigma = rbind(P0, A0)
## 7. Return state space form.
  ssf.mod = list(mDelta=mDelta, mPhi=mPhi, mOmega=mOmega, 
             mSigma=mSigma, mAffine=mAffine)
  CheckSsf(ssf.mod)
}

## 8. compute quasi-mles
start.chen = c(log(0.3), log(0.2), log(0.06), log(0.03), log(0.04),
log(0.02), -0.5, -0.5, log(0.001), log(0.001), log(0.001), log(0.001))
names(start.chen) = c(paste("kappa", 1:2, sep="."),
                    paste("theta", 1:2, sep="."),
                    paste("sigma", 1:2, sep="."),
                    paste("lamda", 1:2, sep="."),
                    paste("s",1:4,sep=""))
start.tau = c(0.25, 0.5, 1, 5)
ans.chen = SsfFitFast(start.chen, fama.bliss, chen.ssf, 
            n.factors=2, tau=start.tau, freq=1/12, trace=T,
            control=nlminb.control(abs.tol=1e-6, rel.tol=1e-6, x.tol=1e-6))
ssf.fit = chen.ssf(ans.chen$parameters, tau=start.tau, freq=1/12, n.factors=2)
ans.chen$parameters[-c(7,8)] = exp(ans.chen$parameters[-c(7,8)])

# delta method covariances - not sandwhich covariance
tmp = ans.chen$parameters
tmp[7:8] = 1
ans.chen$vcov = diag(tmp) %*% ans.chen$vcov %*% diag(tmp)
summary(ans.chen,digits=4)

## 9. smoothed estimates of term-structure

m.s = SsfCondDens(fama.bliss, ssf.fit)
r.s = timeSeries(data=m.s$state,pos=m.s$positions)
y.s = timeSeries(data=m.s$response,pos=m.s$positions)

par(mfrow=c(2,1))
plot(r.s[,1],main="Smoothed Estimate of 1st Factor")
plot(r.s[,2],main="Smoothed Estimate of 2nd Factor")

par(mfrow=c(2,2))
plot(fama.bliss[,1],y.s[,1],main="Actual vs. Smoothed: 3M")
plot(fama.bliss[,2],y.s[,2],main="Actual vs. Smoothed: 6M")
plot(fama.bliss[,3],y.s[,3],main="Actual vs. Smoothed: 12M")
plot(fama.bliss[,4],y.s[,4],main="Actual vs. Smoothed: 60M")

##
## 10. Kalman Filter diagnostics
##

autocorTest(KalmanFil(fama.bliss,ssf.fit)$std.innov)