HW 7
Behzod Ahundjanov
May 20, 2016
Problem 1
In this problem we have to construct the time series for month-over-month inflation rate (\(y_t=\Delta \log CPI_t\)) based on Consumer Price Index for All Urban Consumers (FRED/CPIAUCNS). This exercise is based on codes of lecture 22 (lec22KF_3JNJ.R).
library(Quandl)
library(KFAS)(a)
Let’s restrict the the sample to 1955M1-2013M12 and plot the time series for inflation \(y_t\).
cpiauc <- Quandl("FRED/CPIAUCNS", type="zoo", start_date="1955-01-01", end_date="2013-12-31")
str(cpiauc)## 'zooreg' series from Jan 1955 to Dec 2013
## Data: num [1:708] 26.7 26.7 26.7 26.7 26.7 26.7 26.8 26.8 26.9 26.9 ...
## Index: Class 'yearmon' num [1:708] 1955 1955 1955 1955 1955 ...
## Frequency: 12Now let’s express the date in \(y_t\) form:
lcpiauc <- log(cpiauc)
y <- diff(lcpiauc)
plot(y, xlab="", ylab=expression(""*y[t]*""), main=expression("Month-over-month inflation rate "*y[t]*""))
(b)
Now we set up and estimate local level model with seasonal component:
T <- length(y)Now we define a local level model with seasonal component:
y.model <- SSModel(y~SSMtrend(degree=1,Q=NA)+SSMseasonal(period=12,sea.type="dummy",Q=NA), H=NA)Now we can check system matrices:
y.model$T## , , 1
##
## level sea_dummy1 sea_dummy2 sea_dummy3 sea_dummy4 sea_dummy5
## level 1 0 0 0 0 0
## sea_dummy1 0 -1 -1 -1 -1 -1
## sea_dummy2 0 1 0 0 0 0
## sea_dummy3 0 0 1 0 0 0
## sea_dummy4 0 0 0 1 0 0
## sea_dummy5 0 0 0 0 1 0
## sea_dummy6 0 0 0 0 0 1
## sea_dummy7 0 0 0 0 0 0
## sea_dummy8 0 0 0 0 0 0
## sea_dummy9 0 0 0 0 0 0
## sea_dummy10 0 0 0 0 0 0
## sea_dummy11 0 0 0 0 0 0
## sea_dummy6 sea_dummy7 sea_dummy8 sea_dummy9 sea_dummy10
## level 0 0 0 0 0
## sea_dummy1 -1 -1 -1 -1 -1
## sea_dummy2 0 0 0 0 0
## sea_dummy3 0 0 0 0 0
## sea_dummy4 0 0 0 0 0
## sea_dummy5 0 0 0 0 0
## sea_dummy6 0 0 0 0 0
## sea_dummy7 1 0 0 0 0
## sea_dummy8 0 1 0 0 0
## sea_dummy9 0 0 1 0 0
## sea_dummy10 0 0 0 1 0
## sea_dummy11 0 0 0 0 1
## sea_dummy11
## level 0
## sea_dummy1 -1
## sea_dummy2 0
## sea_dummy3 0
## sea_dummy4 0
## sea_dummy5 0
## sea_dummy6 0
## sea_dummy7 0
## sea_dummy8 0
## sea_dummy9 0
## sea_dummy10 0
## sea_dummy11 0y.model$R## , , 1
##
## [,1] [,2]
## level 1 0
## sea_dummy1 0 1
## sea_dummy2 0 0
## sea_dummy3 0 0
## sea_dummy4 0 0
## sea_dummy5 0 0
## sea_dummy6 0 0
## sea_dummy7 0 0
## sea_dummy8 0 0
## sea_dummy9 0 0
## sea_dummy10 0 0
## sea_dummy11 0 0y.model$Z## , , 1
##
## level sea_dummy1 sea_dummy2 sea_dummy3 sea_dummy4 sea_dummy5
## [1,] 1 1 0 0 0 0
## sea_dummy6 sea_dummy7 sea_dummy8 sea_dummy9 sea_dummy10 sea_dummy11
## [1,] 0 0 0 0 0 0y.model$Q## , , 1
##
## [,1] [,2]
## [1,] NA 0
## [2,] 0 NAy.model$H## , , 1
##
## [,1]
## [1,] NAAnd define update function for maximum likelihood estimation:
y.updatefn <- function(pars, model) {
model$Q[,,1] <- diag(exp(pars[1:2]))
model$H[,,1] <- exp(pars[3])
model
}Now we estimate model parameters using maximum likelihood:
y.fit <- fitSSM(y.model, inits=log(c(0.1,0.01,0.001)), updatefn=y.updatefn, method="BFGS")
y.fit$optim.out$par## [1] -249.95806 -100.75446 -11.37012y.fit$model$Q## , , 1
##
## [,1] [,2]
## [1,] 2.783526e-109 0.000000e+00
## [2,] 0.000000e+00 1.749413e-44y.fit$model$H## , , 1
##
## [,1]
## [1,] 1.153504e-05(c)
Applying Kalman filtering and smoothing:
y.out <- KFS(y.fit$model)
colnames(y.out$alphahat)## [1] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" "sea_dummy4"
## [6] "sea_dummy5" "sea_dummy6" "sea_dummy7" "sea_dummy8" "sea_dummy9"
## [11] "sea_dummy10" "sea_dummy11"construct smoothed time series - add smoothed 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 smoothed state
y.KS.CI.lvl <- as.vector(y.out$alphahat[,ilvl]) + qnorm(0.90)* sqrt(cbind(y.out$V[ilvl,ilvl,])) %*% cbind(-1,1)
y.KS.CI.sea <- as.vector(y.out$alphahat[,isea]) + qnorm(0.90)* sqrt(cbind(y.out$V[isea,isea,])) %*% cbind(-1,1)construct 90% confidence intervals for filtered state
y.KF.CI.lvl <- as.vector(y.out$a[-(T+1),ilvl]) + qnorm(0.90)* 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.90)* sqrt(cbind(y.out$P[isea,isea,-(T+1)])) %*% cbind(-1,1)log transformed earnings per share - actual, filtered and seasonal component
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 earnings per share - actual, smoothed and seasonal component
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)
Now we have to create a 36 period ahead forecast i.e. a forecast for period 2014M1-2016M12 and test how well does the model perform when it comes to its forecasting ability.
We set the forecast horizon as 36 and use the forecast code:
h <- 36
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=12)And the plots of the forecast are as follows:
cpiauc2 <- Quandl("FRED/CPIAUCNS", type="zoo", start_date="2000-01-01", end_date="2016-12-31")
lcpiauc2 <- log(cpiauc2)
y2 <- diff(lcpiauc2)
lcolor <- c(1,2,2,2)
lwidth <- c(2,2,1,1)
ltypes <- c(1,1,2,2)
plot.ts( cbind(y2, y.fcst), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="Month-over-month inflation rate prediction")## Warning in cbind(y2, y.fcst): number of rows of result is not a multiple of
## vector length (arg 1)legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.8 )