JJ=ts(scan("P:/jj.dat"),start=1960,freq=4)

A = t(c(1, 1, 0, 0))
R = .0086
theta = 1.035
Phi = c(theta, 0 , 0, 0, 0, -1, -1, -1, 0, 1, 0, 0, 0, 0, 1, 0)
Phi = matrix(data=Phi, nrow=4, ncol=4, byrow=TRUE)
Q = c(.0169, 0 , 0, 0, 0, .0497, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
Q = matrix(data=Q, nrow=4, ncol=4, byrow=TRUE)
mu = c(.55, .21, .15, .06)
eye = c(1, 0 , 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1)
eye = matrix(data=eye, nrow=4, ncol=4, byrow=TRUE)
Sigma = eye*.01
N = length(JJ)
xf = matrix(0,nrow=4,ncol=N)
xu = matrix(0,nrow=4,ncol=N)
low.bd.fore = vector("numeric",length=N) # To draw upper and lower bounds; not part of the Kalman Żlter.
upp.bd.fore = vector("numeric",length=N)
low.bd.data = vector("numeric",length=N)
upp.bd.data = vector("numeric",length=N)
low.bd.upd = vector("numeric",length=N)
upp.bd.upd = vector("numeric",length=N)

# The Kalman Filter recursions:

# The forecast and update cycle for t=1

xf[,1] = Phi%*%mu # The forecast
K = Sigma%*%t(A)%*%solve(A%*%Sigma%*%t(A) + R) # The Kalman Gain
xu[,1] = xf[,1] + K%*%(JJ[1] - A%*%xf[,1]) # The update
Pu = (eye - K%*%A)%*%Sigma # The posterior error covariance
for(j in 2: N)
{
# The forecast
xf[,j] = Phi%*% xu[,j-1]
Pf = Phi%*%Pu%*%t(Phi) + Q
low.bd.fore[j] = A%*%xf[,j] - 1.96*sqrt(A%*%Pf%*%t(A) ) # State forecast bounds
upp.bd.fore[j] = A%*%xf[,j] + 1.96*sqrt(A%*%Pf%*%t(A) )
low.bd.data[j] = A%*%xf[,j] - 1.96*sqrt(A%*%Pf%*%t(A) + R) # Data forecast bounds
upp.bd.data[j] = A%*%xf[,j] + 1.96*sqrt(A%*%Pf%*%t(A) + R)
# The update
K = Pf%*%t(A)%*%solve(A%*%Pf%*%t(A) + R)
xu[,j] = xf[,j] + K%*%(JJ[j] - A%*%xf[,j])
Pu = (eye - K%*%A)%*%Pf
low.bd.upd[j] = A%*%xu[,j] - 1.96*sqrt(A%*%Pu%*%t(A)) #State Żlter bounds
upp.bd.upd[j] = A%*%xu[,j] + 1.96*sqrt(A%*%Pu%*%t(A))
}

# Plotting
par(mfrow=c(2,1))
quarter = 69:N
plot(quarter,JJ[quarter], type="p", cex = .75, ylim=c(7.5, 17.5), main="Data and Data Forecast 95% PI", xlab="Quarter",
ylab="JJ")
lines(quarter, t(forecast[quarter]))
lines(quarter, low.bd.data[quarter], lty=("dashed"))
lines(quarter, upp.bd.data[quarter], lty=("dashed"))
plot(quarter,JJ[quarter], type="p", cex = .75, , ylim=c(7.5, 17.5), main="Data, State Forecast 95% PI, and Update
95% PI", xlab="Quarter", ylab="JJ")
lines(quarter, t(update[quarter]))
lines(quarter, low.bd.upd[quarter], lty=("dotted"))
lines(quarter, upp.bd.upd[quarter], lty=("dotted"))
lines(quarter, low.bd.fore[quarter], lty=("dashed"))
lines(quarter, upp.bd.fore[quarter], lty=("dashed"))