ECO 5316 H7 Q1
Raymond Kim
Q1
We consider the monthly time series for the Consumer Price Index for all Urban Consumers: All Items Not Seasonally Adjusted the FRED Database.
First we load the “Quandl” packages as well as others.
library("Quandl")
library("tseries")
library("urca")
library("KFAS")A) Restrict sample and plot inflation
The symbol is FRED/CPIAUCNS. First lets find the month over month inflation rate, \(y_t= Delta log CPI_t\).
cpia<-Quandl("FRED/CPIAUCNS", type="timeSeries")
cpia<-window(cpia,start="1955-01-01", end="2013-12-01")
cpia.l<-log(cpia)
cpia.ld<-diff(cpia.l,1)
par(mfrow=c(1,2))
plot(cpia.ld,xlab="", ylab="",main=expression("Yearly Inflation Yt"))
plot(cpia.l,xlab="", ylab="",main=expression("Log CPI"))
B) Set up and estimate local level model with seasonal component
T<-length(cpia.ld)
cpi.model<-SSModel(cpia.ld~ SSMtrend(degree=1,Q=NA)+SSMseasonal(period=12,sea.type="dummy",Q=NA),H=NA)
cpi.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 0cpi.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 0cpi.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 0cpi.model$Q## , , 1
##
## [,1] [,2]
## [1,] NA 0
## [2,] 0 NAcpi.model$H## , , 1
##
## [,1]
## [1,] NANow that we defined the local level model with a seasonal component and checked the system matrices, now we define update function for maximum likelihood estimation, and estimate the model parameters using maximum likelihood.
cpi.updatefn<-function(pars, model) {
model$Q[,,1] <- diag(exp(pars[1:2]))
model$H[,,1] <- exp(pars[3])
model}
pars<-log(c(.1,.01,.001))
cpi.fit<-fitSSM(cpi.model,inits=pars,updatefn=cpi.updatefn,method="BFGS")
cpi.fit$optim.out$par## [1] -249.95806 -100.75446 -11.37012cpi.fit$model$Q## , , 1
##
## [,1] [,2]
## [1,] 2.783526e-109 0.000000e+00
## [2,] 0.000000e+00 1.749413e-44cpi.fit$model$H## , , 1
##
## [,1]
## [1,] 1.153504e-05Plot the Smoothed level and Smoothed Seasonal Component \(\mu_{t|T}\) and \(\gamma_{t|T}\) with 90% confidence intervals. Plot the Smoothed Irregular component \(e_{t|T}\)
# Kalman filtering and smoothing
cpi.out <- KFS(cpi.fit$model, smoothing=c("state","disturbance"))
colnames(cpi.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 smmothed level and smoothed seasonal components
cpi.KF <- cpi.out$a %*% as.vector(cpi.out$model$Z)
cpi.KS <- cpi.out$alphahat %*% as.vector(cpi.out$model$Z)
# find position of level and first seasonal dummy
ilvl <- which(colnames(cpi.out$a)=="level")
isea <- which(colnames(cpi.out$a)=="sea_dummy1")
# construct 90% confidence intervals for filtered state
cpi.KF.CI.lvl <- as.vector(cpi.out$a[-(T+1),ilvl]) + qnorm(0.95)* sqrt(cbind(cpi.out$P[ilvl,ilvl,-(T+1)])) %*% cbind(-1,1)
cpi.KF.CI.sea <- as.vector(cpi.out$a[-(T+1),isea]) + qnorm(0.95)* sqrt(cbind(cpi.out$P[isea,isea,-(T+1)])) %*% cbind(-1,1)
# construct 90% confidence intervals for smoothed state
cpi.KS.CI.lvl <- as.vector(cpi.out$alphahat[,ilvl]) + qnorm(0.95)* sqrt(cbind(cpi.out$V[ilvl,ilvl,])) %*% cbind(-1,1)
cpi.KS.CI.sea <- as.vector(cpi.out$alphahat[,isea]) + qnorm(0.95)* sqrt(cbind(cpi.out$V[isea,isea,])) %*% cbind(-1,1)Now plot the series
# actual, filtered and seasonal component
par(mfrow=c(2,1))
plot.ts( cbind(cpia.ld, cpi.out$a[-(T+1),ilvl], cpi.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(cpi.out$a[-(T+1),isea], cpi.KS.CI.sea), col="blue", lty=c(1,3,3), plot.type="single", xlab="", ylab="", main="seasonal component" )
abline(h=0, lty=3)
# actual, smoothed and seasonal component
par(mfrow=c(2,1))
plot.ts( cbind(cpia.ld, cpi.out$alphahat[,ilvl], cpi.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(cpi.out$alphahat[,isea], cpi.KS.CI.sea), col="red", lty=c(1,3,3), plot.type="single", xlab="", ylab="", main="seasonal component" )
abline(h=0, lty=3)
# Irregular component actual, smoothed and seasonal component
par(mfrow=c(3,1))
plot.ts( exp( cbind(cpia.ld, cpi.out$alphahat[,ilvl], cpi.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(cpi.out$alphahat[,isea], cpi.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( cpi.out$epshat ), col="red", plot.type="single", xlab="", ylab="", main="irregular component" )
abline(h=1, lty=3)
D) Create a 36 period ahead forecast together with the 90% confidence intervals and actual data. How well does the model perform when it comes to forecasting ability?
# forecast horizon
h <- 36
cpi.fcst <- predict(cpi.out$model, interval="confidence", level=0.9, n.ahead=h, filtered=TRUE)
cpi.fcst <- predict(cpi.out$model, interval="confidence", level=0.9, n.ahead=h)
cpi.fcst <- ts(cpi.fcst, start=(index(cpia.ld)[length(cpia.ld)]+0.25), frequency=12)
# plot the forecast
lcolor <- c(1,2,2,2)
lwidth <- c(2,2,1,1)
ltypes <- c(1,1,2,2)
plot.ts( cbind(cpia.ld, cpi.fcst), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="CPI Not Seasonally Adjusted Log Transformed")
legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.5 )
plot.ts( exp( cbind(cpia.ld, cpi.fcst) ), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main="CPI Not Seasonally Adjusted")
legend("topleft", legend = c("actual data","forecast","90% confidence interval"),
col=lcolor, lwd=lwidth, lty=ltypes, bty="n", cex=0.5 )
Yes, it does seem as though the predicted values are within the range of the actual values.