# BNdecomposition.ssc				Beveridge-Nelson decomposition
#
# created: January 8, 2005

##
## 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)
## Real GDP data is available from Eric Zivot's website
## http://faculty.washington.edu/ezivot/ezresearch.htm
##
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
##
ssf.arma22 = arma22.mod(arma22.mle$parameters)
filteredEst.arma22 = SsfMomentEst(dlny.ts.dm,
ssf.arma22,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 objects and 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",reference.grid=F)
plot(c.ts,main="BN Cycle",reference.grid=F)

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