HW8
Mai AlSwayan
5/18/2017
Question8: 1.Obtain monthly data for Consumer Price Index for All Urban Consumers: All Items Less Food and Energy,Not Seasonally Adjusted FRED/CPILFENS. 2.Plot the time series for inflation. 3.et up and estimate local level model with seasonal component, using sample 1955M1-2014M12. 4.Plot the smoothed level and smoothed seasonal components. 5.Create a sequence of 1 step ahead forecasts for period 2015M1-2017M5,and Plot the forecast together with 90% confidence intervals and actual data for the period 2005M1-2016M12. ##Answers #Part 1 obtaining the data , and data summury
library(Quandl)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2. http://CRAN.R-project.org/package=stargazerlibrary(vars)## Loading required package: MASS## Loading required package: strucchange## Loading required package: sandwich## Loading required package: urca## Loading required package: lmtestlibrary(zoo)
library(xts)
library(dygraphs)
library(knitr)
library(forecast)
library(urca)
library(KFAS)
library(g.data)cpi <- Quandl("FRED/CPILFENS", type="zoo")
str(cpi)## 'zooreg' series from Jan 1957 to Apr 2017
## Data: num [1:724] 28.5 28.5 28.7 28.8 28.8 28.9 28.9 29 29.1 29.2 ...
## Index: Class 'yearmon' num [1:724] 1957 1957 1957 1957 1957 ...
## Frequency: 12summary(cpi)## Index cpi
## Min. :1957 Min. : 28.50
## 1st Qu.:1972 1st Qu.: 43.48
## Median :1987 Median :116.45
## Mean :1987 Mean :121.30
## 3rd Qu.:2002 3rd Qu.:189.88
## Max. :2017 Max. :251.64part2
y <- log(cpi)
plot(y, xlab="Years", ylab="",main="Consumer Price Index for All Consumers")
yt <-window(cpi,start="Jan 1955", end = "Dec 2014")
lyt <- log(yt)
dlyt <- diff(lyt)
plot(dlyt, xlab="Years", ylab="",main=" Consumer Index from 1955 to 2014")
#part3
T <- length(dlyt)dlyt.model <- SSModel(dlyt~SSMtrend(degree=1,Q=NA)+SSMseasonal(period=12,sea.type="dummy",Q=NA), H=NA)## Warning in `dim<-.zoo`(`*tmp*`, value = length(x)): setting this dimension
## may lead to an invalid zoo object## Warning in `dim<-.zoo`(`*tmp*`, value = c(695, 1)): setting this dimension
## may lead to an invalid zoo objectdlyt.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 0dlyt.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 0dlyt.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 0dlyt.model$Q## , , 1
##
## [,1] [,2]
## [1,] NA 0
## [2,] 0 NAdlyt.model$H## , , 1
##
## [,1]
## [1,] NAdlyt.updatefn <- function(pars, model) {
model$Q[,,1] <- diag(exp(pars[1:2]))
model$H[,,1] <- exp(pars[3])
model
}dlyt.fit <- fitSSM(dlyt.model, inits=log(c(0.1,0.01,0.001)), updatefn=dlyt.updatefn, method="BFGS")
dlyt.fit$optim.out$par## [1] -246.27582 -99.28775 -11.82852dlyt.fit$model$Q## , , 1
##
## [,1] [,2]
## [1,] 1.106034e-107 0.00000e+00
## [2,] 0.000000e+00 7.58368e-44dlyt.fit$model$H## , , 1
##
## [,1]
## [1,] 7.293567e-06dlyt.out <- KFS(dlyt.fit$model)
colnames(dlyt.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"dlyt.KF <- dlyt.out$a %*% as.vector(dlyt.out$model$Z)
dlyt.KS <- dlyt.out$alphahat %*% as.vector(dlyt.out$model$Z)ilv <- which(colnames(dlyt.out$a)=="level")
ise <- which(colnames(dlyt.out$a)=="se_dummy1")dlyt.KS.CI.lvl <- as.vector(dlyt.out$alphahat[,ilv]) + qnorm(0.90)* sqrt(cbind(dlyt.out$V[ilv,ilv,])) %*% cbind(-1,1)
dlyt.KS.CI.sea <- as.vector(dlyt.out$alphahat[,ise]) + qnorm(0.90)* sqrt(cbind(dlyt.out$V[ise,ise,])) %*% cbind(-1,1)dlyt.KF.CI.lvl <- as.vector(dlyt.out$a[-(T+1),ilv]) + qnorm(0.90)* sqrt(cbind(dlyt.out$P[ilv,ilv,-(T+1)])) %*% cbind(-1,1)
dlyt.KF.CI.se <- as.vector(dlyt.out$a[-(T+1),ise]) + qnorm(0.90)* sqrt(cbind(dlyt.out$P[ise,ise,-(T+1)])) %*% cbind(-1,1)part4
plot.ts( cbind(dlyt, dlyt.out$a[-(T+1),ilv], dlyt.KF.CI.lvl), col=c("green","yellow","red","blue"), lty=c(1,1,3,3), plot.type="single", xlab="", ylab="", main="actual filtered values" )
#part5
h <- 36
dlyt.fcst <- predict(dlyt.out$model, interval="confidence", level=0.9, n.ahead=h)
dlyt.fcst <- ts(dlyt.fcst, start=(index(dlyt)[length(dlyt)]+0.25), frequency=12)cpii <- Quandl("FRED/CPILFENS", type = "zoo",start_date="2015-01-01", end_date="2017-05-31")
lcpii <- log(cpii)
yo <- diff(lcpii)
lcolor <- c(1,2,2,2)
lwidth <- c(2,2,1,1)
ltypes <- c(1,1,2,2)
plot.ts( cbind(yo, dlyt.fcst), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main=" inflation rate prediction")## Warning in cbind(yo, dlyt.fcst): number of rows of result is not a multiple
## of vector length (arg 1)