## structuralChange.ssc		script files for structural change 
##								analysis of exchange rate data
##
## data
##
## author: Eric Zivot
## created: March 31, 2003
## revised: 
## April 4, 2003
##		Added sup-F tests (QLR)
## April 3, 2003
##		Added simulated data from structural change models
## April 2, 2003
##		Added rolling regression analysis
##		Nyblom structural change tests
##		Hansen structural change tests
##	April 1, 2003
##		Added estimation of TVP models

###################################################################
##	Data sets
###################################################################

##
## create simulated data
##
nobs = 200
set.seed(33)
x = rnorm(nobs)
set.seed(59)
e = rnorm(nobs,sd=0.5)
set.seed(78)
u = rnorm(nobs,sd=0.1)
bt = 1 + cumsum(u)
idx = 1:nobs

# no structural change
y0 = x + e
# change in mean
y1 = y0 + 2*(idx > 100)
# change in slope
y2 = y0 + 2*(idx > 100)*x
# change in error variance
y3 = x + e*(idx < 100) + sqrt(0.5)*e*(idx > 100)
# RW slope
y4 = bt*x + e

par(mfrow=c(2,2))
tsplot(y0,main="No structural change",ylab="y",xlab="time")
tsplot(y1,main="Structural change: mean shift",ylab="y",xlab="time")
tsplot(y2,main="Structural change: slope shift",ylab="y",xlab="time")
tsplot(y3,main="Structural change: variance shift",ylab="y",xlab="time")
par(mfrow=c(2,1))
tsplot(bt,main="RW slope",ylab="slope",xlab="time")
tsplot(y4,main="Structural change: RW slope",ylab="y",xlab="time")
par(mfrow=c(1,1))

## full sample estimation
##
fit0 = OLS(y0~x)
summary(fit0)
plot(fit0)
fit1 = OLS(y1~x)
summary(fit1)
plot(fit1)
fit2 = OLS(y2~x)
summary(fit2)
plot(fit2)
fit3 = OLS(y3~x)
summary(fit3)
plot(fit3)
fit4 = OLS(y4~x)
summary(fit4)
plot(fit4)

##
## extract US/DM data from FinMetrics timeSeries lexrates.dat
##
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)
dm.dat = seriesMerge(usdm.ds,tslag(usdm.fp))

## plot data on single graph
##
plot(dm.dat,main="US/DM exchange rate data",reference.grid=F)
legend(0,-6,legend=colIds(dm.dat),lty=c(1,1),col=c(2,1))

## show time series plots as separate plots
##
seriesPlot(usdm.ts[,c(3,4)],strip.text=colIds(usdm.ts[,c(3,4)]),
one.plot=F,main="US/DM exchange rate data")

## show autocorrelations of data
##
par(mfrow=c(2,1))
tmp1 = acf(usdm.ts[,"DUSDMS"],plot=F)
tmp2 = acf(usdm.ts[,"USDMFP"],plot=F)
acf.plot(tmp1,main="Change in US/DM spot rate")
acf.plot(tmp2,main="US/DM forward discount")
par(mfrow=c(1,1))

## show differences regression with regression line
##
guiPlot( PlotType = "Linear Fit", DataSet = "dm.dat", Columns = "DUSDMS,USDMFP.lag1")

## full sample OLS estimation
##
ar1.fit = OLS(USDMFP~tslag(USDMFP),data=usdm.ts)
summary(ar1.fit,correction="white")

ols.fit = OLS(DUSDMS~USDMFP.lag1,data=dm.dat)
summary(ols.fit,correction="white")

###################################################################
##	rolling regression analysis
###################################################################

## no-structural change data
##
data0 = data.frame(y0,x)
fit0.roll = rollOLS(y0~x,data=data0,width=25,incr=1)
fit0.roll
plot(fit0.roll,which=2,sub="No structural change data")
plot(fit0.roll,which=3,sub="No structural change data")

## mean-shift data
##
data1 = data.frame(y1,x)
fit1.roll = rollOLS(y1~x,data=data1,width=25,incr=1)
fit1.roll
plot(fit1.roll,which=2,sub="Mean shift data")
plot(fit1.roll,which=3,sub="Mean shift data")

## slope-shift data
##
data2 = data.frame(y2,x)
fit2.roll = rollOLS(y2~x,data=data2,width=25,incr=1)
fit2.roll
plot(fit2.roll,which=2,sub="slope shift data")
plot(fit2.roll,which=3,sub="slope shift data")

## variance-shift data
##
data3 = data.frame(y3,x)
fit3.roll = rollOLS(y3~x,data=data3,width=25,incr=1)
fit3.roll
plot(fit3.roll,which=2,sub="Error variance shift data")
plot(fit3.roll,which=3,sub="Error variance shift data")


## RW slope data
##
data4 = data.frame(y4,x)
fit4.roll = rollOLS(y4~x,data=data4,width=25,incr=1)
fit4.roll
plot(fit4.roll,which=2,sub="RW slope data")
plot(fit4.roll,which=3,sub="RW slope data")


## AR(1) regression
##
ar1.roll.24 = rollOLS(USDMFP~tslag(USDMFP),data=usdm.ts,
width=24,incr=1)
ar1.roll.24
ar1.roll.48 = rollOLS(USDMFP~tslag(USDMFP),data=usdm.ts,
width=48,incr=12)
ar1.roll.48

plot(ar1.roll.24,which=2,sub="US/DM forward discount")
plot(ar1.roll.24,which=3,sub="US/DM forward discount")
plot(ar1.roll.48,which=2,sub="US/DM forward discount")
plot(ar1.roll.48,which=3,sub="US/DM forward discount")

## differences regression
##
dr.roll.24 = rollOLS(DUSDMS~USDMFP.lag1,data=dm.dat,
width=24,incr=1)
dr.roll.24
dr.roll.48 = rollOLS(DUSDMS~USDMFP.lag1,data=dm.dat,
width=48,incr=12)
dr.roll.48
plot(dr.roll.24,which=2,sub="US/DM differences regression")
plot(dr.roll.24,which=3,sub="US/DM differences regression")
plot(dr.roll.48,which=2,sub="US/DM differences regression")
plot(dr.roll.48,which=3,sub="US/DM differences regression")

###################################################################
##	Recursive regression analysis
###################################################################

## no structural change data
##
rls0.fit = RLS(y0~x,data=data0)
plot(rls0.fit,which=1,sub="No structural change")
plot(rls0.fit,which=2,sub="No structural change")
plot(rls0.fit,which=3,sub="No structural change")

## mean shift data
##
rls1.fit = RLS(y1~x,data=data1)
plot(rls1.fit,which=1,sub="Mean shift data")
plot(rls1.fit,which=2,sub="Mean shift data")
plot(rls1.fit,which=3,sub="Mean shift data")

## slope shift data
##
rls2.fit = RLS(y2~x,data=data2)
plot(rls2.fit,which=1,sub="slope shift data")
plot(rls2.fit,which=2,sub="slope shift data")
plot(rls2.fit,which=3,sub="slope shift data")

## variance shift data
##
rls3.fit = RLS(y3~x,data=data3)
plot(rls3.fit,which=1,sub="variance shift data")
plot(rls3.fit,which=2,sub="variance shift data")
plot(rls3.fit,which=3,sub="variance shift data")

## rw slope data
##
rls4.fit = RLS(y4~x,data=data3)
plot(rls4.fit,which=1,sub="RW slope data")
plot(rls4.fit,which=2,sub="RW slope data")
plot(rls4.fit,which=3,sub="RW slope data")

## AR(1) model for forward discount
##
rls.ar1 = RLS(USDMFP~tslag(USDMFP),data=usdm.ts,na.rm=T)
rls.ar1

# plot(rls.ar1,which=1)
ncoef = nrow(coef(rls.ar1))
tmp = coef(rls.ar1)
colIds(tmp) = c("(Intercept)","Lag1")
seriesPlot(tmp[10:ncoef,],one.plot=F,strip.text=colIds(tmp),
main="RLS estimation of AR(1) for forward discount")
plot(rls.ar1,which=2)
plot(rls.ar1,which=3,sub="AR(1) for forward discount")

## exchange rate differences regression
##
rls.dr = RLS(DUSDMS~USDMFP.lag1,data=dm.dat)
rls.dr

# plot(rls.dr,which=1)
ncoef = nrow(coef(rls.dr))
tmp = coef(rls.dr)
colIds(tmp) = c("(Intercept)","f(t) - s(t)")
seriesPlot(tmp[10:ncoef,],one.plot=F,strip.text=colIds(tmp),
main="RLS estimation of differences regression")
plot(rls.dr,which=2,sub="differences regression")
plot(rls.dr,which=3,sub="differences regression")

###################################################################
##	parameter stability tests
###################################################################

## no structural change data
##
Nyblom.test0 = tvpTest(fit0)
Nyblom.test0
Hansen.test0 = HansenTest(fit0)
Hansen.test0

## mean shift data
##
Nyblom.test1 = tvpTest(fit1)
Nyblom.test1
Hansen.test1 = HansenTest(fit1)
Hansen.test1

## slope shift data
##
Nyblom.test2 = tvpTest(fit2)
Nyblom.test2
Hansen.test2 = HansenTest(fit2)
Hansen.test2

## variance shift data
##
Nyblom.test3 = tvpTest(fit3)
Nyblom.test3
Hansen.test3 = HansenTest(fit3)
Hansen.test3

## RW slope data
##
Nyblom.test4 = tvpTest(fit4)
Nyblom.test4
Hansen.test4 = HansenTest(fit4)
Hansen.test4


## forward difference AR(1)
##
Nyblom.test.ar1 = tvpTest(ar1.fit)
Nyblom.test.ar1

Hansen.test.ar1 = HansenTest(ar1.fit)
Hansen.test.ar1

## differences regression
##
Nyblom.test.dr = tvpTest(ols.fit)
Nyblom.test.dr

Hansen.test.dr = HansenTest(ols.fit)
Hansen.test.dr


###################################################################
#	sup-F tests (QLR)
###################################################################

## function to compute F-test on dummy variables
##

Ftest = function(ols.fit,R,r) {	
	bhat = coef(ols.fit)
	tmp1 = R%*%bhat - r
	tmp2 = R%*%vcov(ols.fit)%*%t(R)
	ans = t(tmp1)%*%solve(tmp2)%*%tmp1/length(r)
	as.numeric(ans)
}

##
## sup-T tests for simulated data
##

idx = 1:nrow(data0)
R = rbind(c(0,0,1,0),c(0,0,0,1))
r = c(0,0)

trim.fraction = 0.15
trim1 = floor(trim.fraction*nrow(data0))
trim2 = floor((1-trim.fraction)*nrow(data0))

## No structural change data
##

# d.vals = matrix(0,length(trim1:trim2),4)
F.stats0 = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*data0$x
	tmp.data = data.frame(data0,D1,D2)
	fit = OLS(y0~x+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats0[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,data0$y0[trim1:trim2],type="l",
main="No structural change data",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats0,type="l",
main="F-stats for structural change: No structural change data",
ylim = c(0,6),
ylab="F statistics",xlab="break date")
abline(h=5.86)

## find break date
##
qlr0 = max(F.stats0)
m.hat0 = which(F.stats0==max(F.stats0)) + trim1 - 1
lamda.hat0 = m.hat0/200

## mean shift data
##
F.stats1 = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*data1$x
	tmp.data = data.frame(data1,D1,D2)
	fit = OLS(y1~x+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats1[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,data1$y1[trim1:trim2],type="l",
main="Mean shift data",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats1,type="l",
main="F-stats for structural change: Mean shift data",
ylab="F statistics",xlab="break date")
abline(h=5.86)

## compute QLR and find break date
##
qlr1 = max(F.stats1)
m.hat1 = which(F.stats1==max(F.stats1)) + trim1 - 1
lamda.hat1 = m.hat1/200

## slope shift data
##
F.stats2 = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*data2$x
	tmp.data = data.frame(data2,D1,D2)
	fit = OLS(y2~x+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats2[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,data2$y2[trim1:trim2],type="l",
main="Slope shift data",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats2,type="l",
main="F-stats for structural change: Slope shift data",
ylab="F statistics",xlab="break date")
abline(h=5.86)

## find break date
##
qlr2 = max(F.stats2)
m.hat2 = which(F.stats2==max(F.stats2)) + trim1 - 1
lamda.hat2 = m.hat2/200

## var shift data
##
F.stats3 = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*data3$x
	tmp.data = data.frame(data3,D1,D2)
	fit = OLS(y3~x+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats3[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,data3$y3[trim1:trim2],type="l",
main="Error variance shift data",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats3,type="l",
main="F-stats for structural change: Var shift data",
ylab="F statistics",xlab="break date")
abline(h=5.86)

## find break date
##
qlr3 = max(F.stats3)
m.hat3 = which(F.stats3==max(F.stats3)) + trim1 - 1
lamda.hat3 = m.hat3/200

## RW slope data
##
F.stats4 = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*data4$x
	tmp.data = data.frame(data4,D1,D2)
	fit = OLS(y4~x+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats4[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,data4$y4[trim1:trim2],type="l",
main="RW slope data",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats4,type="l",
main="F-stats for structural change: RW slope data",
ylab="F statistics",xlab="break date")
abline(h=5.86)

## find break date
##
qlr4 = max(F.stats4)
m.hat4 = which(F.stats4==max(F.stats4)) + trim1 - 1
lamda.hat4 = m.hat4/200


## AR(1) forward discount data
##

ar1.data = seriesMerge(usdm.ts[,"USDMFP"],tslag(usdm.ts[,"USDMFP"],trim=T))
ar1.data = as.data.frame(ar1.data)

R = rbind(c(0,0,1,0),c(0,0,0,1))
r = c(0,0)

trim.fraction = 0.15
trim1 = floor(trim.fraction*nrow(ar1.data))
trim2 = floor((1-trim.fraction)*nrow(ar1.data))
idx = 1:nrow(ar1.data)

F.stats.ar1 = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*ar1.data$USDMFP.lag1
	tmp.data = data.frame(ar1.data,D1,D2)
	fit = OLS(USDMFP~USDMFP.lag1+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats.ar1[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,ar1.data$USDMFP[trim1:trim2],type="l",
main="US/DM forward discount",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats.ar1,type="l",
main="F-stats for structural change: US/DM forward discount",
ylab="F statistics",xlab="break date")
abline(h=5.86)

## find break date
##
qlr.ar1 = max(F.stats.ar1)
m.hat.ar1 = which(F.stats.ar1==max(F.stats.ar1)) + trim1 - 1
lamda.hat.ar1 = m.hat.ar1/nrow(ar1.data)


## differences regression data
##

dr.data = as.data.frame(seriesData(dm.dat))

R = rbind(c(0,0,1,0),c(0,0,0,1))
r = c(0,0)

trim.fraction = 0.15
trim1 = floor(trim.fraction*nrow(dr.data))
trim2 = floor((1-trim.fraction)*nrow(dr.data))
idx = 1:nrow(dr.data)

F.stats.dr = matrix(0,length(trim1:trim2),1)
for (i in trim1:trim2) {
	D1 = as.integer((idx >= i))
	D2 = D1*dr.data$USDMFP.lag1
	tmp.data = data.frame(dr.data,D1,D2)
	fit = OLS(DUSDMS~USDMFP.lag1+D1+D2,data=tmp.data)
#	d.vals[i-trim1+1,] = as.vector(coef(fit))
	F.stats.dr[i-trim1+1] = Ftest(fit,R,r)
}

## plot F-stats
##
par(mfrow=c(2,1))
plot(trim1:trim2,dr.data$DUSDMS[trim1:trim2],type="l",
main="US/DM spot exchange rate",
ylab="y",xlab="break date")

plot(trim1:trim2,F.stats.dr,type="l",
ylim = c(0,6),
main="F-stats for structural change: US/DM differences reg",
ylab="F statistics",xlab="break date")
abline(h=5.86)

tmp = timeSeries(data=F.stats.dr,
pos=positions(dm.dat[trim1:trim2,]))
par(mfrow=c(2,1))
plot(dm.dat[trim1:trim2,1],reference.grid=F)
plot(tmp,reference.grid=F)

## find break date
##
qlr.dr = max(F.stats.dr)
m.hat.dr = which(F.stats.dr==max(F.stats.dr)) + trim1 - 1
lamda.hat.dr = m.hat.dr/nrow(dr.data)

###################################################################
## multiple break plots
###################################################################

## pure sc model
##
m1 = -0.084
m2 = -0.412
m3 = -0.249
m4 = 0.023
m5 = 0.375
m6 = -0.141

idx = 1:nrow(dm.dat)
mb = m1*(idx < 18) + m2*(idx >= 18 & idx <103) +
m3*(idx >= 103 & idx < 163) + m4*(idx >=163 & idx < 180) +
m5*(idx >= 180 & idx < 218) + m6*(idx >= 218)

mb.ts = timeSeries(data=-1*mb,pos=positions(dm.dat))
plot(mb.ts,dm.dat[,2],reference.grid=F,ylab="forward discount",
main="Mean shift model: US/DM forward discount")

###################################################################
## time varying parameter regression state space model
###################################################################

## state space model: all coefficients vary
##
tvp.mod = function(parm,mX=NULL) {
	parm = exp(parm)		# log variances
	ssf.tvp = GetSsfReg(mX=mX)
	diag(ssf.tvp$mOmega) = parm
	CheckSsf(ssf.tvp)
}

## state space model: only slope varies
##
tvp.mod2 = function(parm,mX=NULL) {
	parm = exp(parm)		# log variances
	ssf.tvp = GetSsfReg(mX=mX)
	diag(ssf.tvp$mOmega[2:3,2:3]) = parm
	CheckSsf(ssf.tvp)
}

##
## analysis of simulated data
##
tmp = as.matrix(data4$x)
X4.mat = cbind(1,tmp)

## time varying parameter model estimation of differences regression
## y(t) = a(t) + b(t)*x(t-1) + e(t)
##
tvp4.start = c(0,0,0)			# log variances
names(tvp4.start) = c("ln(s2.alpha)","ln(s2.beta)","ln(s2.y)")
tvp4.mle = SsfFit(tvp4.start,data4$y4,"tvp.mod",mX=X4.mat)
summary(tvp4.mle)

## create state space form
##
ssf4.tvp = tvp.mod(tvp4.mle$parameters,mX=X4.mat)

## transform parameters and compute delta method variances
##
tvp4.mle$parameters = exp(tvp4.mle$parameters/2)	# std devs
names(tvp4.mle$parameters) = c("s.alpha","s.beta","s.y")
dg = diag(tvp4.mle$parameters/2)
tvp4.mle$vcov = dg%*%tvp4.mle$vcov%*%dg
summary(tvp4.mle)

filteredEst.tvp4 = SsfMomentEst(data4$y4,
ssf4.tvp,task="STFIL")
colIds(filteredEst.tvp4$state.moment) = c("intercept","slope")
plot(filteredEst.tvp4)

SmoothedEst.tvp4 = SsfMomentEst(data4$y4,
ssf4.tvp,task="STSMO")
colIds(SmoothedEst.tvp4$state.moment) = c("intercept","slope")
plot(SmoothedEst.tvp4)

## plot smoothed estimates with SE bars
##
state4 = SmoothedEst.tvp4$state.moment
sd4 = sqrt(SmoothedEst.tvp4$state.variance)
upper4 = state4 + 1*sd4
lower4 = state4 -1*sd4

plot(filteredEst.tvp4)
plot(SmoothedEst.tvp4)
par(mfrow=c(2,1))
tsplot(state4[-(1:2),1],upper4[-(1:2),1],lower4[-(1:2),1],
main="Smoothed estimates of alpha(t)",lty=c(1,2,2),
col=c(1,3,3))
tsplot(state4[-(1:2),2],upper4[-(1:2),2],lower4[-(1:2),2],
main="Smoothed estimates of beta(t)",lty=c(1,2,2),col=c(1,3,3))

## plot actual vs. smoothed beta(t)
##
par(mfrow=c(1,1))
tsplot(bt,state4[,2],main="Actual and smoothed slope",
lty=c(1,1),col=c(1,3))
legend(0,0.5,legend=c("Actual","Smoothed"),lty=c(1,1),col=c(1,3))

##
## analysis of US/DM exchange rates: forward discount
##

##
## data matrix for TVP estimation
tmp = as.matrix(seriesData(tslag(usdm.ts[,"USDMFP"],trim=T)))
X.ar1.mat = cbind(1,tmp)

## time varying parameter model estimation of differences regression
## y(t) = a(t) + b(t)*y(t-1) + e(t)
##
tvp.ar1.start = c(0,0,0)			# log variances
names(tvp.ar1.start) = c("ln(s2.alpha)","ln(s2.beta)","ln(s2.y)")
tvp.ar1.mle = SsfFit(tvp.ar1.start,usdm.ts[-1,"USDMFP"],
"tvp.mod",mX=X.ar1.mat)
summary(tvp.ar1.mle)

## create state space form
##
ssf.ar1.tvp = tvp.mod(tvp.ar1.mle$parameters,mX=X.ar1.mat)

## transform parameters and compute delta method variances
##
tvp.ar1.mle2 = tvp.ar1.mle
tvp.ar1.mle2$parameters = exp(tvp.ar1.mle$parameters/2)	# std devs
names(tvp.ar1.mle2$parameters) = c("s.alpha","s.beta","s.y")
dg = diag(tvp.ar1.mle2$parameters/2)
tvp.ar1.mle2$vcov = dg%*%tvp.ar1.mle2$vcov%*%dg
summary(tvp.ar1.mle2)

## filtered estimates
##
filteredEst.ar1.tvp = SsfMomentEst(usdm.ts[-1,"USDMFP"],
ssf.ar1.tvp,task="STFIL")
colIds(filteredEst.ar1.tvp$state.moment) = c("intercept","slope")

## plot filtered estimates with 95% error bands
##
filteredEst.ar1.ts = timeSeries(data=filteredEst.ar1.tvp$state.moment,
positions=filteredEst.ar1.tvp$positions)
filteredSd.ar1.ts = sqrt(timeSeries(data=filteredEst.ar1.tvp$state.variance,
positions=filteredEst.ar1.tvp$positions))
upper.ar1 = filteredEst.ar1.ts + 2*filteredSd.ar1.ts
lower.ar1 = filteredEst.ar1.ts - 2*filteredSd.ar1.ts

par(mfrow=c(2,1))
plot(filteredEst.ar1.ts[-(1:3),1],upper.ar1[-(1:3),1],lower.ar1[-(1:3),1],
main="Filtered estimates of alpha(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

plot(filteredEst.ar1.ts[-(1:3),2],upper.ar1[-(1:3),2],lower.ar1[-(1:3),2],
main="Filtered estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

## smoothed estimates
##
SmoothedEst.ar1.tvp = SsfMomentEst(usdm.ts[-1,"USDMFP"],
ssf.ar1.tvp,task="STSMO")

## plot smoothed estimates
##
SmoothedEst.ar1.ts = timeSeries(data=SmoothedEst.ar1.tvp$state.moment,
pos=SmoothedEst.ar1.tvp$positions)
SmoothedSd.ar1.ts = sqrt(timeSeries(data=SmoothedEst.ar1.tvp$state.variance,
pos=SmoothedEst.ar1.tvp$positions))
supper.ar1 = SmoothedEst.ar1.ts + 1*SmoothedSd.ar1.ts
slower.ar1 = SmoothedEst.ar1.ts -1*SmoothedSd.ar1.ts


par(mfrow=c(2,1))
plot(filteredEst.ar1.ts[-(1:3),1],upper.ar1[-(1:3),1],lower.ar1[-(1:3),1],
main="Filtered estimates of alpha(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

plot(filteredEst.ar1.ts[-(1:3),2],upper.ar1[-(1:3),2],lower.ar1[-(1:3),2],
main="Filtered estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

par(mfrow=c(2,1))
plot(SmoothedEst.ar1.ts[-(1:2),1],supper.ar1[-(1:2),1],slower.ar1[-(1:2),1],
main="Smoothed estimates of alpha(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

plot(SmoothedEst.ar1.ts[-(1:2),2],supper.ar1[-(1:2),2],slower.ar1[-(1:2),2],
main="Smoothed estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

## smoothed estimate of mean
##
SmoothedMean = SmoothedEst.ts[,1]/(1 - SmoothedEst.ts[,2])
plot(SmoothedMean,reference.grid=F)

##
## analysis of US/DM exchange rates: differences regression
##

## data matrix for TVP estimation
##
X.dr.mat = cbind(1,as.matrix(seriesData(dm.dat[,2])))	

## full sample regression
##
dm.fit = OLS(DUSDMS~USDMFP.lag1,data=dm.dat)

## time varying parameter model estimation of differences regression
## ds(t) = a(t) + b(t)*(f(t-1) - s(t-1)) + e(t)
##
tvp.dr.start = c(0,0,0)			# log variances
names(tvp.dr.start) = c("ln(s2.alpha)","ln(s2.beta)","ln(s2.y)")
tvp.dr.mle = SsfFit(tvp.dr.start,dm.dat[,1],"tvp.mod",mX=X.dr.mat)

##	note: Hessian fails to invert
##
tvp.dr. mle
## sqrt(diag(tvp.mle$vcov))

sigma.dr.mle = sqrt(exp(tvp.dr.mle$parameters))
names(sigma.dr.mle) = c("s.alpha","s.beta","s.y")
sigma.dr.mle
## note: variance for alpha(t) is almost zero
## try estimation to see if variance goes to zero
## delta method variance
## dg = diag(exp(tvp.mle$parameters)/2)
## se.sigma = sqrt(diag(dg%*%tvp.mle$vcov%*%dg))
## names(se.sigma) = names(sigma.mle)
## se.sigma

## filtered estimates
##
filteredEst.dr.tvp = SsfMomentEst(dm.dat[,1],
tvp.mod(tvp.dr.mle$parameters,mX=X.dr.mat),task="STFIL")
colIds(filteredEst.dr.tvp$state.moment) = c("intercept","slope")

## plot filtered moment estimates
## note: first few elements are enormous
##
plot(filteredEst.dr.tvp,strip.text=c("alpha(t)",
"beta(t)","Spot returns"),main="Filtered Estimates")

## plot filtered estimates excluding first 7 obvs
##
seriesPlot(filteredEst.dr.tvp$state.moment[-(1:7),],one.plot=F,
strip.text=c("alpha(t)","beta(t)"),main="Filtered Estimates",
layout=c(1,2))

## plot filtered estimates with 95% error bands
##
filteredEst.dr.ts = timeSeries(data=filteredEst.dr.tvp$state.moment,
positions=filteredEst.dr.tvp$positions)
filteredSd.dr.ts = sqrt(timeSeries(data=filteredEst.dr.tvp$state.variance,
positions=filteredEst.dr.tvp$positions))
fupper.dr = filteredEst.dr.ts + 2*filteredSd.dr.ts
flower.dr = filteredEst.dr.ts - 2*filteredSd.dr.ts

par(mfrow=c(2,1))
plot(filteredEst.dr.ts[-(1:7),1],fupper.dr[-(1:7),1],flower.dr[-(1:7),1],
main="Filtered estimates of alpha(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

plot(filteredEst.dr.ts[-(1:7),2],fupper.dr[-(1:7),2],flower.dr[-(1:7),2],
main="Filtered estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

## alternative plot using seriesPlot
##
tmp = seriesMerge(filteredEst.ts[-(1:7),2],upper[-(1:7),2],lower[-(1:7),2])
colIds(tmp) = c("slope","upper","lower")
seriesPlot(tmp,strip.text="beta(t)")

## smoothed estimates
##
SmoothedEst.dr.tvp = SsfMomentEst(dm.dat[,1],
tvp.mod(tvp.dr.mle$parameters,mX=X.dr.mat),
task="STSMO")

## plot smoothed estimates
##
SmoothedEst.dr.ts = timeSeries(data=SmoothedEst.dr.tvp$state.moment,
pos=SmoothedEst.dr.tvp$positions)
SmoothedSd.dr.ts = sqrt(timeSeries(data=SmoothedEst.dr.tvp$state.variance,
pos=SmoothedEst.dr.tvp$positions))
supper.dr = SmoothedEst.dr.ts + 1*SmoothedSd.dr.ts
slower.dr = SmoothedEst.dr.ts -1*SmoothedSd.dr.ts

par(mfrow=c(2,1))
plot(filteredEst.dr.ts[-(1:7),1],fupper.dr[-(1:7),1],flower.dr[-(1:7),1],
main="Filtered estimates of alpha(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

plot(filteredEst.dr.ts[-(1:7),2],fupper.dr[-(1:7),2],flower.dr[-(1:7),2],
main="Filtered estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))


par(mfrow=c(2,1))
plot(SmoothedEst.dr.ts[-(1:2),1],supper.dr[-(1:2),1],slower.dr[-(1:2),1],
main="Smoothed estimates of alpha(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

plot(SmoothedEst.dr.ts[-(1:2),2],supper.dr[-(1:2),2],slower.dr[-(1:2),2],
main="Smoothed estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))

##
## analysis of differences regression with fixed intercept
##

X.mat = cbind(1,as.matrix(seriesData(dm.dat[,2])))
tvp.start = c(0,0)			# log variances
names(tvp.start) = c("ln(s2.beta)","ln(s2.y)")
tvp.mle2 = SsfFit(tvp.start,dm.dat[,1],"tvp.mod2",mX=X.mat)
summary(tvp.mle2)

## compute state space from estimated parameters
##
ssf.tvp = tvp.mod2(tvp.mle2$parameters,mX=X.mat)

## transform parameters and compute delta method variances
##
tvp.mle2$parameters = exp(tvp.mle2$parameters/2)	# std devs
names(tvp.mle2$parameters) = c("s.beta","s.y")
dg = diag(tvp.mle2$parameters/2)
tvp.mle2$vcov = dg%*%tvp.mle2$vcov%*%dg
summary(tvp.mle2)

## smoothed estimates
##
SmoothedEst.tvp = SsfMomentEst(dm.dat[,1],ssf.tvp,task="STSMO")
colIds(SmoothedEst.tvp$state.moment) = c("Intercept","Slope")

## smoothed estimate of intercept with standard error
##
n = nrow(SmoothedEst.tvp$state.moment)
SmoothedEst.tvp$state.moment[n,1]
sqrt(SmoothedEst.tvp$state.variance[n,1])

## plot smoothed estimates
##
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 + 1*SmoothedSd.ts
lower = SmoothedEst.ts -1*SmoothedSd.ts
par(mfrow=c(1,1))
plot(SmoothedEst.ts[-(1:2),2],upper[-(1:2),2],lower[-(1:2),2],
main="US/DM: Smoothed estimates of beta(t)",reference.grid=F,
plot.args=list(lty=c(1,2,2),col=c(1,3,3)))