Problem 2

For this task we need to consider local level model for Johnson and Johnson’s earnings per share from Tsay (q-jnj.txt). This exercise is based on codes of lecture 22 (lec22KF_3JNJ.R).

(a)

library(KFAS)
library(Quandl)
epsJNJ <- ts(scan(file="http://faculty.chicagobooth.edu/ruey.tsay/teaching/fts3/q-jnj.txt"), start=c(1960,1), frequency=4)
y <- log(epsJNJ)

Let’s examine a 16 period ahead forecast from this model. Here \(y_t=\log(epsJNJ)\).

par(mfcol=c(2,1), cex=0.9)
plot(epsJNJ, xlab="", ylab="epsJNJ", main="Johnson and Johnson: earnings per share")
plot(y, xlab="", ylab=expression(""*y[t]*""), main="Johnson and Johnson: log transformed earnings per share")

In the clas we derived following forecast for the next 16 quarters:

# forecast horizon
h <- 16
# y.fcst <- predict(y.out$model, interval="confidence", level=0.9, n.ahead=h, filtered=TRUE)
y.fcst <- predict(y.out$model, interval="confidence", level=0.9, n.ahead=h)
y.fcst <- ts(y.fcst, start=(index(y)[length(y)]+0.25), frequency=4)
# plot the forecast
lcolor <- c(1,2,2,2)
lwidth <- c(2,2,1,1)
ltypes <- c(1,1,2,2)
plot.ts( cbind(y, y.fcst), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="Johnson and Johnson: log transformed earnings per share")
legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
       col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.8 )

plot.ts( exp( cbind(y, y.fcst) ), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="Johnson and Johnson: earnings per share")
legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
       col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.8 )

Looking at both graphs and taking into account we can say that it might be little misleading. Since it is clear that earnings per share are stricltly growing over time, but in the forecast it is poorly reflected. That’s why in this paricular case local linear trend model is more appropriate.

(b)

Since we concluded that local linear trend model is more appropriate than local level model in the above mentioned case, we need to modify the model from local level to local linear trend.

T <- length(y)
# define a local level model with seasonal component
y.model <- SSModel(y ~ SSMtrend(degree=2, Q=list(NA, NA)) + SSMseasonal(period=4, sea.type="dummy", Q=NA), H=NA)
# define update function for maximum likelihood estimation
y.updatefn <- function(pars, model) {
    model$Q[,,1] <- diag(exp(pars[1:3]))
    model$H[,,1] <- exp(pars[3])
    model
}
# estimate the model parameters using maximum likelihood
pars <- log(c(0.1,0.01,0.001))
y.fit <- fitSSM(y.model, inits=pars, updatefn=y.updatefn, method="BFGS")
# Kalman filtering and smoothing
y.out <- KFS(y.fit$model, smoothing=c("state","disturbance"))
colnames(y.out$alphahat)
## [1] "level"      "slope"      "sea_dummy1" "sea_dummy2" "sea_dummy3"
# construct smoothed time series - add smmothed level and smoothed seasonal components
y.KF <- y.out$a %*% as.vector(y.out$model$Z)
y.KS <- y.out$alphahat %*% as.vector(y.out$model$Z)
# find position of level and first seasonal dummy
ilvl <- which(colnames(y.out$a)=="level")
isea <- which(colnames(y.out$a)=="sea_dummy1")
# construct 90% confidence intervals for filtered state
y.KF.CI.lvl <- as.vector(y.out$a[-(T+1),ilvl]) + qnorm(0.95)* sqrt(cbind(y.out$P[ilvl,ilvl,-(T+1)])) %*% cbind(-1,1)
y.KF.CI.sea <- as.vector(y.out$a[-(T+1),isea]) + qnorm(0.95)* sqrt(cbind(y.out$P[isea,isea,-(T+1)])) %*% cbind(-1,1)
# construct 90% confidence intervals for smoothed state
y.KS.CI.lvl <- as.vector(y.out$alphahat[,ilvl]) + qnorm(0.95)* sqrt(cbind(y.out$V[ilvl,ilvl,])) %*% cbind(-1,1)
y.KS.CI.sea <- as.vector(y.out$alphahat[,isea]) + qnorm(0.95)* sqrt(cbind(y.out$V[isea,isea,])) %*% cbind(-1,1)

(c)

plot.ts( cbind(y, y.out$alphahat[,ilvl], y.KS.CI.lvl), col=c("black","red","red","red"), lty=c(1,1,3,3), plot.type="single", xlab="", ylab="", main="actual and smoothed values" )

plot.ts( cbind(y.out$alphahat[,isea], y.KS.CI.sea), col="red", lty=c(1,3,3), plot.type="single", xlab="", ylab="", main="seasonal component" )
abline(h=0, lty=3)

(d)

# forecast horizon
h <- 16
# y.fcst <- predict(y.out$model, interval="confidence", level=0.9, n.ahead=h, filtered=TRUE)
y.fcst <- predict(y.out$model, interval="confidence", level=0.9, n.ahead=h)
y.fcst <- ts(y.fcst, start=(index(y)[length(y)]+0.25), frequency=4)
# plot the forecast
lcolor <- c(1,2,2,2)
lwidth <- c(2,2,1,1)
ltypes <- c(1,1,2,2)
plot.ts( cbind(y, y.fcst), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="Johnson and Johnson: log transformed earnings per share")
legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
       col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.8 )

plot.ts( exp( cbind(y, y.fcst) ), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="Johnson and Johnson: earnings per share")
legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
       col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.8 )

Now after applying the local linear trend we can see more or less realistic scenario. Unlike the forecast we observed in part (a) here forecast takes into account the fact that share prices were growing overtime.