Q2

Examine the 16 period ahead forecast from this model. After running the initial settings and forecasts, we look at the graphs of the output.

Looking at the output, the upper bound of the confidence interval does look realistic, but not the forecast itself. Earnings per share has been growing over time in a trend, but the forecast does not reflect this trend and instead looks flat.

B) We modify from the local level to the local linear trend

# 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)
                   
# check system matrices
y.model$T
## , , 1
## 
##            level slope sea_dummy1 sea_dummy2 sea_dummy3
## level          1     1          0          0          0
## slope          0     1          0          0          0
## sea_dummy1     0     0         -1         -1         -1
## sea_dummy2     0     0          1          0          0
## sea_dummy3     0     0          0          1          0
y.model$R
## , , 1
## 
##            [,1] [,2] [,3]
## level         1    0    0
## slope         0    1    0
## sea_dummy1    0    0    1
## sea_dummy2    0    0    0
## sea_dummy3    0    0    0
y.model$Z
## , , 1
## 
##      level slope sea_dummy1 sea_dummy2 sea_dummy3
## [1,]     1     0          1          0          0
y.model$Q
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]   NA    0    0
## [2,]    0   NA    0
## [3,]    0    0   NA
y.model$H
## , , 1
## 
##      [,1]
## [1,]   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
}

The updates to the linear model is above.

C Plot Smoothed Level Component and Seasonal Component.

# log transformed earnigns per share - actual, filtered and seasonal component
par(mfrow=c(2,1))
plot.ts( cbind(y, y.out$a[-(T+1),ilvl], y.KF.CI.lvl), col=c("black","blue","blue","blue"), lty=c(1,1,3,3), plot.type="single", xlab="", ylab="", main="actual and filtered values" )
plot.ts( cbind(y.out$a[-(T+1),isea], y.KS.CI.sea), col="blue", lty=c(1,3,3), plot.type="single", xlab="", ylab="", main="seasonal component" )
abline(h=0, lty=3)

# log transformed earnigns per share - actual, smoothed and seasonal component
par(mfrow=c(2,1))
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)

# earnigns per share - actual, smoothed and seasonal component
par(mfrow=c(3,1))
plot.ts( exp( 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( exp( 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=1, lty=3)
plot.ts( exp( y.out$epshat ), col="red", plot.type="single", xlab="", ylab="", main="irregular component" )
abline(h=1, lty=3)

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

D) Create forecast

# plot the forecast
par(mfcol=c(2,1), cex=0.9)
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 )

Yes, this does look more plausible than the one in part A)