Homework 7 Problem 1
Niraj Prasad Koirala
May 19, 2016
#PROBLEM ONE#Objective
#a) Plotting the time series of logarithmically transormed data of Consumer Price Index.
#b) Setting up and estimating local level model and seasonal component.#
#c) Plotting the smoothed and smoothed seasonal components at 90 percentage interval and estimating irregular disturbance terms.formula.#
#d) Create a 36 period ahead forecast.##SOURCES OF DATA
#Consumer Price Index of All Consumers: [FRED/CPIAUCNS](https://www.quandl.com/data/FRED/CPIAUCNS-Consumer-Price-Index-for-All-Urban-Consumers-All-Items)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(tseries)
library(forecast)## Loading required package: timeDate## This is forecast 7.1library(vars)## Loading required package: MASS## Loading required package: strucchange## Loading required package: sandwich## Loading required package: urca## Loading required package: lmtestlibrary(gdata)## gdata: Unable to locate valid perl interpreter
## gdata:
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata:
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.## ## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.## ## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.##
## Attaching package: 'gdata'## The following objects are masked from 'package:xts':
##
## first, last## The following object is masked from 'package:stats':
##
## nobs## The following object is masked from 'package:utils':
##
## object.size## The following object is masked from 'package:base':
##
## startsWithlibrary(KFAS)
library(zoo)
library(dlmodeler)cpi <- Quandl("FRED/CPIAUCNS",api_key="T8gDpM8iFjXR9pfbXBmx",type = "zoo")y <- log(cpi)
plot(y, xlab="Years", ylab="",main="Consumer Price Index for All Consumers")
# (A)
# Restricting the Data from January/1955 to December/2016 and plotting time series.#yt <-window(cpi,start="Jan 1955", end = "Dec 2016")
lyt <- log(yt)
dlyt <- diff(lyt)
plot(dlyt, xlab="Years", ylab="",main=" Consumer Index from 1955 to 2013")
# (B)
#Setting up and estimating local level model with seasonal component#Estimating local level model with seasonal component.
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(735, 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] -259.75012 -104.64014 -11.36904dlyt.fit$model$Q## , , 1
##
## [,1] [,2]
## [1,] 1.55581e-113 0.000000e+00
## [2,] 0.00000e+00 3.592241e-46dlyt.fit$model$H## , , 1
##
## [,1]
## [1,] 1.154748e-05- Kalman Filter
dlyt.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)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" )
- Forecasting
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/CPIAUCNS",api_key="T8gDpM8iFjXR9pfbXBmx",type = "zoo",start_date="1990-01-01", end_date="2016-12-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)
Problem 2
#Objective
#a) Examining the 16 period ahead forecast.
#b)Modifiying the model from local level model to local linear model.
#c) Plotting the smoothed level component and smoothed seasonal component.
#d) Creating a plot of 16 period ahead forecast.#Sources of Data
#Johnson and Johnson's earnings per share from the book of Tsay.#Methodology#Question Ajnj <- ts(scan(file="http://faculty.chicagobooth.edu/ruey.tsay/teaching/fts3/q-jnj.txt"), start=c(1960,1), frequency=4)y <- log(jnj)T <- length(y)y.model <- SSModel(y ~ SSMtrend(degree=1, Q=NA) + SSMseasonal(period=4, sea.type="dummy", Q=NA), H=NA)#Checking System Matrices
y.model$T## , , 1
##
## level sea_dummy1 sea_dummy2 sea_dummy3
## level 1 0 0 0
## sea_dummy1 0 -1 -1 -1
## sea_dummy2 0 1 0 0
## sea_dummy3 0 0 1 0y.model$R## , , 1
##
## [,1] [,2]
## level 1 0
## sea_dummy1 0 1
## sea_dummy2 0 0
## sea_dummy3 0 0y.model$Z## , , 1
##
## level sea_dummy1 sea_dummy2 sea_dummy3
## [1,] 1 1 0 0y.model$Q## , , 1
##
## [,1] [,2]
## [1,] NA 0
## [2,] 0 NAy.model$H## , , 1
##
## [,1]
## [1,] NA#Defining 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
}#Estimating the Parameters of the model.
pars <- log(c(0.1,0.01,0.001))
y.fit <- fitSSM(y.model, inits=pars, updatefn=y.updatefn, method="BFGS")
y.fit$optim.out$par## [1] -5.243001 -7.059233 -14.832445y.fit$model$Q## , , 1
##
## [,1] [,2]
## [1,] 0.005284377 0.0000000000
## [2,] 0.000000000 0.0008594368y.fit$model$H## , , 1
##
## [,1]
## [1,] 3.617019e-07#Kalman Filtering
y.out <- KFS(y.fit$model, smoothing=c("state","disturbance"))
colnames(y.out$alphahat)## [1] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3"#Forecasting
h<-16
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)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.5 )
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.5 )
#Question Bt.model <- SSModel(y ~ SSMtrend(degree=2) + SSMseasonal(period=4, sea.type="dummy"), H=NA)t.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 0t.model$R## , , 1
##
## [,1]
## level 0
## slope 0
## sea_dummy1 0
## sea_dummy2 0
## sea_dummy3 0t.model$Z## , , 1
##
## level slope sea_dummy1 sea_dummy2 sea_dummy3
## [1,] 1 0 1 0 0t.model$Q## , , 1
##
## [,1]
## [1,] 0t.model$H## , , 1
##
## [,1]
## [1,] NA#Kalman Filtering
t.out <- KFS(y.fit$model, smoothing=c("state","disturbance"))
colnames(t.out$alphahat)## [1] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3"#Constructing Smoothed Time Series
t.KF <- t.out$a %*% as.vector(t.out$model$Z)
t.KS <- t.out$alphahat %*% as.vector(t.out$model$Z)#Finding the position of level and first seasonal dummy.
ilvl <- which(colnames(t.out$a)=="level")
isea <- which(colnames(t.out$a)=="sea_dummy1")#Constructing 90% intervals for filtered state
t.KF.CI.lvl <- as.vector(t.out$a[-(T+1),ilvl]) + qnorm(0.95)* sqrt(cbind(t.out$P[ilvl,ilvl,-(T+1)])) %*% cbind(-1,1)
t.KF.CI.sea <- as.vector(t.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
t.KS.CI.lvl <- as.vector(t.out$alphahat[,ilvl]) + qnorm(0.95)* sqrt(cbind(t.out$V[ilvl,ilvl,])) %*% cbind(-1,1)
t.KS.CI.sea <- as.vector(t.out$alphahat[,isea]) + qnorm(0.95)* sqrt(cbind(y.out$V[isea,isea,])) %*% cbind(-1,1)#Question C# log transformed earnigns per share - actual, filtered and seasonal component
par(mfrow=c(2,1))
plot.ts( cbind(y, t.out$a[-(T+1),ilvl], t.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(t.out$a[-(T+1),isea], t.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, t.out$alphahat[,ilvl], t.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], t.KS.CI.sea), col="red", lty=c(1,3,3), plot.type="single", xlab="", ylab="", main="seasonal component" )
abline(h=0, lty=3)
#D Forecastingh <- 16t.fcst <- predict(t.out$model, interval="confidence", level=0.9, n.ahead=h)
t.fcst <- ts(t.fcst, start=(index(y)[length(y)]+0.25), frequency=4)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, t.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.5 )
plot.ts( exp( cbind(y, t.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.5 )
# The last forecast is improvement over the first one. Problem 3
#Objectives
#a) Modifying and Estimating the Model using Dummy Variables.
#b)Plotting with Actual Vs Kalman Filetered Data and Actual Vs Kalman Smoothed Data.
#c)Introducing Missing Data.library(dummy)## dummy 0.1.3## dummyNews()#Annual Flow of the Nile from 1871 to 1970.
data(Nile)
help(Nile)## starting httpd help server ...## doneplot(Nile)
f <- datasets::Nilef.SSM <- SSModel(f~ SSMtrend(1,Q=list(NA)) + SSMseasonal(period=12,sea.type='dummy',Q=NA), H=NA)
f.SSM## Call:
## SSModel(formula = f ~ SSMtrend(1, Q = list(NA)) + SSMseasonal(period = 12,
## sea.type = "dummy", Q = NA), H = NA)
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 100
## [1] Number of time series: 1
## [1] Number of disturbances: 2
## [1] Number of states: 12
## Names of the states:
## [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
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.#Problem b#Checking the State Variable
f.SSM$a## [,1]
## level 0
## sea_dummy1 0
## 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 0# examine system matrices
f.SSM$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 0f.SSM$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 0f.SSM$Q## , , 1
##
## [,1] [,2]
## [1,] NA 0
## [2,] 0 NAf.SSM$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 0f.SSM$H## , , 1
##
## [,1]
## [1,] NA# define update function for maximum likelihhod estimation
updatefn_KSI <- function(pars,model,...){
model$H[] <- exp(pars[1])
diag(model$Q[,,1])<- exp(c(pars[2:3]))
model
}# estimate model parameters using maximum likelihhod
f.ML <- fitSSM(inits=log(c(var(f),0.001,0.0001)),
model=f.SSM,
updatefn=updatefn_KSI,
method='BFGS')
str(f.ML$model)## List of 14
## $ y : Time-Series [1:100, 1] from 1871 to 1970: 1120 1160 963 1210 1160 1160 813 1230 1370 1140 ...
## $ Z : num [1, 1:12, 1] 1 1 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : NULL
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : NULL
## $ H : num [1, 1, 1] 30156
## $ T : num [1:12, 1:12, 1] 1 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : NULL
## $ R : num [1:12, 1:2, 1] 1 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 3
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : NULL
## .. ..$ : NULL
## $ Q : num [1:2, 1:2, 1] 1e-03 0e+00 0e+00 1e-04
## $ a1 : num [1:12, 1] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : NULL
## $ P1 : num [1:12, 1:12] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## $ P1inf : num [1:12, 1:12] 1 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## .. ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...
## $ u : chr "Omitted"
## $ distribution: chr "gaussian"
## $ tol : num 1.49e-08
## $ call : language SSModel(formula = f ~ SSMtrend(1, Q = list(NA)) + SSMseasonal(period = 12, sea.type = "dummy", Q = NA), H = NA)
## $ terms :Classes 'terms', 'formula' language f ~ SSMtrend(1, Q = list(NA)) + SSMseasonal(period = 12, sea.type = "dummy", Q = NA)
## .. ..- attr(*, "variables")= language list(f, SSMtrend(1, Q = list(NA)), SSMseasonal(period = 12, sea.type = "dummy", Q = NA))
## .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "f" "SSMtrend(1, Q = list(NA))" "SSMseasonal(period = 12, sea.type = \"dummy\", Q = NA)"
## .. .. .. ..$ : chr [1:2] "SSMtrend(1, Q = list(NA))" "SSMseasonal(period = 12, sea.type = \"dummy\", Q = NA)"
## .. ..- attr(*, "term.labels")= chr [1:2] "SSMtrend(1, Q = list(NA))" "SSMseasonal(period = 12, sea.type = \"dummy\", Q = NA)"
## .. ..- attr(*, "specials")=Dotted pair list of 6
## .. .. ..$ SSMregression: NULL
## .. .. ..$ SSMtrend : int 2
## .. .. ..$ SSMseasonal : int 3
## .. .. ..$ SSMcycle : NULL
## .. .. ..$ SSMarima : NULL
## .. .. ..$ SSMcustom : NULL
## .. ..- attr(*, "order")= int [1:2] 1 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "class")= chr "SSModel"
## - attr(*, "p")= int 1
## - attr(*, "m")= int 12
## - attr(*, "k")= int 2
## - attr(*, "n")= int 100
## - attr(*, "tv")= int [1:5] 0 0 0 0 0
## - attr(*, "state_types")= chr [1:12] "level" "seasonal" "seasonal" "seasonal" ...
## - attr(*, "eta_types")= chr [1:2] "level" "seasonal"# estimated parameters - variance of innovations in measurement equation, state level equation, seasonal component equation
res <- data.frame(sigma2.measurement = f.ML$model$H[,,1],
sigma2.level = diag(f.ML$model$Q[,,1])[1],
sigma2.seasonality = diag(f.ML$model$Q[,,1])[2])
f.ML$model$Q## , , 1
##
## [,1] [,2]
## [1,] 0.001001128 0.000000e+00
## [2,] 0.000000000 9.999991e-05f.ML$model$H## , , 1
##
## [,1]
## [1,] 30156.08# Kalman filtering and smoothing
f.KFS <- KFS(f.ML$model, filtering=c('state'), smoothing=c('state','disturbance','mean'))
f.KFS## Smoothed values of states and standard errors at time n = 100:
## Estimate Std. Error
## level 918.740 17.393
## sea_dummy1 51.368 55.630
## sea_dummy2 -31.188 55.630
## sea_dummy3 28.590 55.630
## sea_dummy4 11.923 55.630
## sea_dummy5 -1.993 58.683
## sea_dummy6 4.382 58.683
## sea_dummy7 43.007 58.683
## sea_dummy8 -9.743 58.683
## sea_dummy9 23.632 58.683
## sea_dummy10 -107.618 58.683
## sea_dummy11 -47.118 58.683str(f.KFS$alphahat)## Time-Series [1:100, 1:12] from 1871 to 1970: 919 919 919 919 919 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:12] "level" "sea_dummy1" "sea_dummy2" "sea_dummy3" ...dimnames(f.KFS$alphahat)## [[1]]
## NULL
##
## [[2]]
## [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"# extract smoothed level, seasonal, and irregular component
f.lvl.KS <- f.KFS$alphahat[,"level"]
f.sea.KS <- f.KFS$alphahat[,"sea_dummy1"]
f.eps.KS <- f.KFS$epshat
par(mfrow=c(3,1), cex=0.7)# actual data and smoothed state component
tmp <- ts( cbind(f, f.lvl.KS), start=1895, frequency=12 )
plot.ts(tmp, plot.type="single", xlab="",ylab="log KSI", col=c("darkgrey","red"), lwd=c(1,2))
abline(v=1995:2013, lty= "dotted",col="lightgrey")
legend("topright", c("actual data","Kalman smoothed level"), col=c("darkgrey","red"), lwd=c(1,2), bty="n")
# smoothed seasonal component
tmp <- ts( f.sea.KS, start=1895, frequency=12)
plot(tmp, xlab="", ylab="log KSI", col="red", lwd=2)
abline(h=0,col="grey")
abline(v=1895:2013, lty= "dotted",col="lightgrey")
legend("topright", "Kalman smoothed seasonal component", col="red", lwd=2, bty="n")
# smoothed irregular component
tmp <- ts( f.eps.KS, start=1895, frequency=12)
plot(tmp,xlab="",ylab="log KSI", col="red", lwd=2)
abline(h=0,col="grey")
abline(v=1895:2013, lty= "dotted",col="lightgrey")
legend("topright", "irregular component", col="red", lwd=2, bty="n")