HW 7
A S M Shakil Haider
22 May 2016
Problem 1
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(zoo)
library(tseries)
library(forecast)## Loading required package: timeDate## This is forecast 6.2library(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.sizelibrary(KFAS)## Warning: package 'KFAS' was built under R version 3.2.5library(dlmodeler)## Warning: package 'dlmodeler' was built under R version 3.2.51(A)
cpi <-Quandl("FRED/CPIAUCNS",api_key= "5NyVfHfxdScCYsskDP59", type = "zoo")
lcpi<- log(cpi)
plot(lcpi, xlab="Years", ylab="",main="CPI")
y<-lcpi
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")
1(B)
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$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-051(C)
Kalman Filtering:
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"KFdlyt<- dlyt.out$a %*% as.vector(dlyt.out$model$Z)
KSdlyt<- dlyt.out$alphahat %*% as.vector(dlyt.out$model$Z)ilv <- which(colnames(dlyt.out$a)=="level")
ise <- which(colnames(dlyt.out$a)=="se_dummy1")KSdlyt.CI.lvl <- as.vector(dlyt.out$alphahat[,ilv]) + qnorm(0.90)* sqrt(cbind(dlyt.out$V[ilv,ilv,])) %*% cbind(-1,1)
KSdlyt.CI.sea <- as.vector(dlyt.out$alphahat[,ise]) + qnorm(0.90)* sqrt(cbind(dlyt.out$V[ise,ise,])) %*% cbind(-1,1)KFdlyt.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)
KFdlyt.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], KFdlyt.CI.lvl), col=c("green","yellow","red","blue"), lty=c(1,1,3,3), plot.type="single", xlab="", ylab="", main="Filtered Values (real/actual)" )
1(D)
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)cppi <- Quandl("FRED/CPIAUCNS", api_key="5NyVfHfxdScCYsskDP59", type = "zoo",start_date="1990-01-01", end_date="2016-12-31")
lcppi <- log(cppi)
dy <- diff(lcppi)
lcolor <- c(1,2,2,2)
lwidth <- c(2,2,1,1)
ltypes <- c(1,1,2,2)
plot.ts( cbind(dy, dlyt.fcst), plot.type="single", col=lcolor, lwd=lwidth, lty=ltypes, xlab="", ylab="", main=" Predicting rate of Inflation")## Warning in cbind(dy, dlyt.fcst): number of rows of result is not a multiple
## of vector length (arg 1)
Problem 2
2(A)
library("Quandl")
library("tseries")
library("urca")## Warning: package 'urca' was built under R version 3.2.4library("KFAS")y <- datasets::JohnsonJohnson
y.model <- SSModel(y ~ SSMtrend(degree=2, Q=list(NA, NA)) + SSMseasonal(period=4, sea.type="dummy", Q=NA), H=NA)
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 0Updating for MLE:
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 0y.model$Z## , , 1
##
## level slope sea_dummy1 sea_dummy2 sea_dummy3
## [1,] 1 0 1 0 0y.model$Q## , , 1
##
## [,1] [,2] [,3]
## [1,] NA 0 0
## [2,] 0 NA 0
## [3,] 0 0 NAy.model$H## , , 1
##
## [,1]
## [1,] NAy.updatefn <- function(pars, model) {
model$Q[,,1] <- diag(exp(pars[1:3]))
model$H[,,1] <- exp(pars[3])
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] -12.081664 -6.964867 -3.297725y.fit$model$Q## , , 1
##
## [,1] [,2] [,3]
## [1,] 5.662394e-06 0.0000000000 0.00000000
## [2,] 0.000000e+00 0.0009444887 0.00000000
## [3,] 0.000000e+00 0.0000000000 0.03696716y.fit$model$H## , , 1
##
## [,1]
## [1,] 0.03696716Kalman Filter:
y.out <- KFS(y.fit$model, smoothing=c("state","disturbance"))
colnames(y.out$alphahat)## [1] "level" "slope" "sea_dummy1" "sea_dummy2" "sea_dummy3"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 )
2(B)
t.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,] NAt.out <- KFS(y.fit$model, smoothing=c("state","disturbance"))
colnames(t.out$alphahat)## [1] "level" "slope" "sea_dummy1" "sea_dummy2" "sea_dummy3"t.KF <- t.out$a %*% as.vector(t.out$model$Z)
t.KS <- t.out$alphahat %*% as.vector(t.out$model$Z)
ilvl <- which(colnames(t.out$a)=="level")
isea <- which(colnames(t.out$a)=="sea_dummy1")
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)
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)2(C)
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)
2(D)
h <- 16
t.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 )
Problem 3
3(A)
library(dummy)## dummy 0.1.3## dummyNews()data(Nile)
plot(Nile)
T<-length(Nile)year<-seq(1871,1970)
dummy<-ifelse(year>1889,1,0)
y<-datasets::Nile
y<-cbind(y,dummy)
y.LLM <- SSModel(y ~ SSMtrend(1, Q=list(NA))+dummy, data=y, H=NA )
y.LLM## Call:
## SSModel(formula = y ~ SSMtrend(1, Q = list(NA)) + dummy, data = y,
## 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: 1
## [1] Number of states: 2
## Names of the states:
## [1] dummy level
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.y.LLM$a## [,1]
## dummy 0
## level 0y.LLM$Z## , , 1
##
## dummy level
## [1,] 0 1
##
## , , 2
##
## dummy level
## [1,] 0 1
##
## , , 3
##
## dummy level
## [1,] 0 1
##
## , , 4
##
## dummy level
## [1,] 0 1
##
## , , 5
##
## dummy level
## [1,] 0 1
##
## , , 6
##
## dummy level
## [1,] 0 1
##
## , , 7
##
## dummy level
## [1,] 0 1
##
## , , 8
##
## dummy level
## [1,] 0 1
##
## , , 9
##
## dummy level
## [1,] 0 1
##
## , , 10
##
## dummy level
## [1,] 0 1
##
## , , 11
##
## dummy level
## [1,] 0 1
##
## , , 12
##
## dummy level
## [1,] 0 1
##
## , , 13
##
## dummy level
## [1,] 0 1
##
## , , 14
##
## dummy level
## [1,] 0 1
##
## , , 15
##
## dummy level
## [1,] 0 1
##
## , , 16
##
## dummy level
## [1,] 0 1
##
## , , 17
##
## dummy level
## [1,] 0 1
##
## , , 18
##
## dummy level
## [1,] 0 1
##
## , , 19
##
## dummy level
## [1,] 0 1
##
## , , 20
##
## dummy level
## [1,] 1 1
##
## , , 21
##
## dummy level
## [1,] 1 1
##
## , , 22
##
## dummy level
## [1,] 1 1
##
## , , 23
##
## dummy level
## [1,] 1 1
##
## , , 24
##
## dummy level
## [1,] 1 1
##
## , , 25
##
## dummy level
## [1,] 1 1
##
## , , 26
##
## dummy level
## [1,] 1 1
##
## , , 27
##
## dummy level
## [1,] 1 1
##
## , , 28
##
## dummy level
## [1,] 1 1
##
## , , 29
##
## dummy level
## [1,] 1 1
##
## , , 30
##
## dummy level
## [1,] 1 1
##
## , , 31
##
## dummy level
## [1,] 1 1
##
## , , 32
##
## dummy level
## [1,] 1 1
##
## , , 33
##
## dummy level
## [1,] 1 1
##
## , , 34
##
## dummy level
## [1,] 1 1
##
## , , 35
##
## dummy level
## [1,] 1 1
##
## , , 36
##
## dummy level
## [1,] 1 1
##
## , , 37
##
## dummy level
## [1,] 1 1
##
## , , 38
##
## dummy level
## [1,] 1 1
##
## , , 39
##
## dummy level
## [1,] 1 1
##
## , , 40
##
## dummy level
## [1,] 1 1
##
## , , 41
##
## dummy level
## [1,] 1 1
##
## , , 42
##
## dummy level
## [1,] 1 1
##
## , , 43
##
## dummy level
## [1,] 1 1
##
## , , 44
##
## dummy level
## [1,] 1 1
##
## , , 45
##
## dummy level
## [1,] 1 1
##
## , , 46
##
## dummy level
## [1,] 1 1
##
## , , 47
##
## dummy level
## [1,] 1 1
##
## , , 48
##
## dummy level
## [1,] 1 1
##
## , , 49
##
## dummy level
## [1,] 1 1
##
## , , 50
##
## dummy level
## [1,] 1 1
##
## , , 51
##
## dummy level
## [1,] 1 1
##
## , , 52
##
## dummy level
## [1,] 1 1
##
## , , 53
##
## dummy level
## [1,] 1 1
##
## , , 54
##
## dummy level
## [1,] 1 1
##
## , , 55
##
## dummy level
## [1,] 1 1
##
## , , 56
##
## dummy level
## [1,] 1 1
##
## , , 57
##
## dummy level
## [1,] 1 1
##
## , , 58
##
## dummy level
## [1,] 1 1
##
## , , 59
##
## dummy level
## [1,] 1 1
##
## , , 60
##
## dummy level
## [1,] 1 1
##
## , , 61
##
## dummy level
## [1,] 1 1
##
## , , 62
##
## dummy level
## [1,] 1 1
##
## , , 63
##
## dummy level
## [1,] 1 1
##
## , , 64
##
## dummy level
## [1,] 1 1
##
## , , 65
##
## dummy level
## [1,] 1 1
##
## , , 66
##
## dummy level
## [1,] 1 1
##
## , , 67
##
## dummy level
## [1,] 1 1
##
## , , 68
##
## dummy level
## [1,] 1 1
##
## , , 69
##
## dummy level
## [1,] 1 1
##
## , , 70
##
## dummy level
## [1,] 1 1
##
## , , 71
##
## dummy level
## [1,] 1 1
##
## , , 72
##
## dummy level
## [1,] 1 1
##
## , , 73
##
## dummy level
## [1,] 1 1
##
## , , 74
##
## dummy level
## [1,] 1 1
##
## , , 75
##
## dummy level
## [1,] 1 1
##
## , , 76
##
## dummy level
## [1,] 1 1
##
## , , 77
##
## dummy level
## [1,] 1 1
##
## , , 78
##
## dummy level
## [1,] 1 1
##
## , , 79
##
## dummy level
## [1,] 1 1
##
## , , 80
##
## dummy level
## [1,] 1 1
##
## , , 81
##
## dummy level
## [1,] 1 1
##
## , , 82
##
## dummy level
## [1,] 1 1
##
## , , 83
##
## dummy level
## [1,] 1 1
##
## , , 84
##
## dummy level
## [1,] 1 1
##
## , , 85
##
## dummy level
## [1,] 1 1
##
## , , 86
##
## dummy level
## [1,] 1 1
##
## , , 87
##
## dummy level
## [1,] 1 1
##
## , , 88
##
## dummy level
## [1,] 1 1
##
## , , 89
##
## dummy level
## [1,] 1 1
##
## , , 90
##
## dummy level
## [1,] 1 1
##
## , , 91
##
## dummy level
## [1,] 1 1
##
## , , 92
##
## dummy level
## [1,] 1 1
##
## , , 93
##
## dummy level
## [1,] 1 1
##
## , , 94
##
## dummy level
## [1,] 1 1
##
## , , 95
##
## dummy level
## [1,] 1 1
##
## , , 96
##
## dummy level
## [1,] 1 1
##
## , , 97
##
## dummy level
## [1,] 1 1
##
## , , 98
##
## dummy level
## [1,] 1 1
##
## , , 99
##
## dummy level
## [1,] 1 1
##
## , , 100
##
## dummy level
## [1,] 1 1y.LLM$T## , , 1
##
## dummy level
## dummy 1 0
## level 0 1y.LLM$R## , , 1
##
## [,1]
## dummy 0
## level 1y.LLM$H## , , 1
##
## [,1]
## [1,] NAy.LLM$Q## , , 1
##
## [,1]
## [1,] NApars <- log( rep(var(Nile, na.rm=TRUE), 2) )
y.fit <- fitSSM( model=y.LLM, inits=pars, method='BFGS' )
y.fit## $optim.out
## $optim.out$par
## [1] 7.536911 9.556134
##
## $optim.out$value
## [1] 625.7254
##
## $optim.out$counts
## function gradient
## 29 9
##
## $optim.out$convergence
## [1] 0
##
## $optim.out$message
## NULL
##
##
## $model
## Call:
## SSModel(formula = y ~ SSMtrend(1, Q = list(NA)) + dummy, data = y,
## 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: 1
## [1] Number of states: 2
## Names of the states:
## [1] dummy level
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.y.out <- KFS(y.fit$model)
y.KF.CI <- predict(y.fit$model, interval="confidence", level=0.9, filtered=TRUE)
y.KS.CI <- predict(y.fit$model, interval="confidence", level=0.9)
y.KF.CI[1,] <- NA
par(mfrow=c(2,2), cex=0.6)
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend = c("data","filtered state","90% confidence interval"),
col = c("black","blue","blue"), lty = c(1,1,2), lwd = c(1,1,1),bty="n", cex=0.9 )
plot( ts( c(y.out$v[-1]), start=1871), col="blue", lwd=2, xlab="", ylab="", main="forecast error")
abline(h=0, lty=3)
plot( ts( c(y.out$P)[-1], start=1871), col="blue", lwd=2, xlab="", ylab="", main="variance of state")
plot( ts( c(y.out$F)[-1], start=1871), col="blue", lwd=2, xlab="", ylab="", main="variance of forecast error")
par(mfrow=c(1,2), mar=c(3,3,2,1), cex=0.8)
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","filtered","90% confidence interval"), col=c(1,4,4), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
plot.ts( cbind(y, y.KS.CI, Nile), plot.type="single",
col=c(1,2,2,2,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","smoothed","90% confidence interval"), col=c(1,2,2), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
3(B)
par(mfrow=c(1,2), mar=c(3,3,2,1), cex=0.8)
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","filtered","90% confidence interval"), col=c(1,4,4), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
plot.ts( cbind(y, y.KS.CI, Nile), plot.type="single",
col=c(1,2,2,2,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","smoothed","90% confidence interval"), col=c(1,2,2), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
3(C)
data(Nile)
plot(Nile)
T<-length(Nile)
year<-seq(1871,1970)
dummy<-ifelse(year>1889,1,0)
y<-datasets::Nile
y<-cbind(y,dummy)
y[21:50] <- NA
y[71:80] <- NAy.LLM <- SSModel(y ~ SSMtrend(1, Q=list(NA))+dummy, data=y, H=NA )
y.LLM## Call:
## SSModel(formula = y ~ SSMtrend(1, Q = list(NA)) + dummy, data = y,
## 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: 1
## [1] Number of states: 2
## Names of the states:
## [1] dummy level
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.y.LLM$a## [,1]
## dummy 0
## level 0y.LLM$Z## , , 1
##
## dummy level
## [1,] 0 1
##
## , , 2
##
## dummy level
## [1,] 0 1
##
## , , 3
##
## dummy level
## [1,] 0 1
##
## , , 4
##
## dummy level
## [1,] 0 1
##
## , , 5
##
## dummy level
## [1,] 0 1
##
## , , 6
##
## dummy level
## [1,] 0 1
##
## , , 7
##
## dummy level
## [1,] 0 1
##
## , , 8
##
## dummy level
## [1,] 0 1
##
## , , 9
##
## dummy level
## [1,] 0 1
##
## , , 10
##
## dummy level
## [1,] 0 1
##
## , , 11
##
## dummy level
## [1,] 0 1
##
## , , 12
##
## dummy level
## [1,] 0 1
##
## , , 13
##
## dummy level
## [1,] 0 1
##
## , , 14
##
## dummy level
## [1,] 0 1
##
## , , 15
##
## dummy level
## [1,] 0 1
##
## , , 16
##
## dummy level
## [1,] 0 1
##
## , , 17
##
## dummy level
## [1,] 0 1
##
## , , 18
##
## dummy level
## [1,] 0 1
##
## , , 19
##
## dummy level
## [1,] 0 1
##
## , , 20
##
## dummy level
## [1,] 1 1
##
## , , 21
##
## dummy level
## [1,] 1 1
##
## , , 22
##
## dummy level
## [1,] 1 1
##
## , , 23
##
## dummy level
## [1,] 1 1
##
## , , 24
##
## dummy level
## [1,] 1 1
##
## , , 25
##
## dummy level
## [1,] 1 1
##
## , , 26
##
## dummy level
## [1,] 1 1
##
## , , 27
##
## dummy level
## [1,] 1 1
##
## , , 28
##
## dummy level
## [1,] 1 1
##
## , , 29
##
## dummy level
## [1,] 1 1
##
## , , 30
##
## dummy level
## [1,] 1 1
##
## , , 31
##
## dummy level
## [1,] 1 1
##
## , , 32
##
## dummy level
## [1,] 1 1
##
## , , 33
##
## dummy level
## [1,] 1 1
##
## , , 34
##
## dummy level
## [1,] 1 1
##
## , , 35
##
## dummy level
## [1,] 1 1
##
## , , 36
##
## dummy level
## [1,] 1 1
##
## , , 37
##
## dummy level
## [1,] 1 1
##
## , , 38
##
## dummy level
## [1,] 1 1
##
## , , 39
##
## dummy level
## [1,] 1 1
##
## , , 40
##
## dummy level
## [1,] 1 1
##
## , , 41
##
## dummy level
## [1,] 1 1
##
## , , 42
##
## dummy level
## [1,] 1 1
##
## , , 43
##
## dummy level
## [1,] 1 1
##
## , , 44
##
## dummy level
## [1,] 1 1
##
## , , 45
##
## dummy level
## [1,] 1 1
##
## , , 46
##
## dummy level
## [1,] 1 1
##
## , , 47
##
## dummy level
## [1,] 1 1
##
## , , 48
##
## dummy level
## [1,] 1 1
##
## , , 49
##
## dummy level
## [1,] 1 1
##
## , , 50
##
## dummy level
## [1,] 1 1
##
## , , 51
##
## dummy level
## [1,] 1 1
##
## , , 52
##
## dummy level
## [1,] 1 1
##
## , , 53
##
## dummy level
## [1,] 1 1
##
## , , 54
##
## dummy level
## [1,] 1 1
##
## , , 55
##
## dummy level
## [1,] 1 1
##
## , , 56
##
## dummy level
## [1,] 1 1
##
## , , 57
##
## dummy level
## [1,] 1 1
##
## , , 58
##
## dummy level
## [1,] 1 1
##
## , , 59
##
## dummy level
## [1,] 1 1
##
## , , 60
##
## dummy level
## [1,] 1 1
##
## , , 61
##
## dummy level
## [1,] 1 1
##
## , , 62
##
## dummy level
## [1,] 1 1
##
## , , 63
##
## dummy level
## [1,] 1 1
##
## , , 64
##
## dummy level
## [1,] 1 1
##
## , , 65
##
## dummy level
## [1,] 1 1
##
## , , 66
##
## dummy level
## [1,] 1 1
##
## , , 67
##
## dummy level
## [1,] 1 1
##
## , , 68
##
## dummy level
## [1,] 1 1
##
## , , 69
##
## dummy level
## [1,] 1 1
##
## , , 70
##
## dummy level
## [1,] 1 1
##
## , , 71
##
## dummy level
## [1,] 1 1
##
## , , 72
##
## dummy level
## [1,] 1 1
##
## , , 73
##
## dummy level
## [1,] 1 1
##
## , , 74
##
## dummy level
## [1,] 1 1
##
## , , 75
##
## dummy level
## [1,] 1 1
##
## , , 76
##
## dummy level
## [1,] 1 1
##
## , , 77
##
## dummy level
## [1,] 1 1
##
## , , 78
##
## dummy level
## [1,] 1 1
##
## , , 79
##
## dummy level
## [1,] 1 1
##
## , , 80
##
## dummy level
## [1,] 1 1
##
## , , 81
##
## dummy level
## [1,] 1 1
##
## , , 82
##
## dummy level
## [1,] 1 1
##
## , , 83
##
## dummy level
## [1,] 1 1
##
## , , 84
##
## dummy level
## [1,] 1 1
##
## , , 85
##
## dummy level
## [1,] 1 1
##
## , , 86
##
## dummy level
## [1,] 1 1
##
## , , 87
##
## dummy level
## [1,] 1 1
##
## , , 88
##
## dummy level
## [1,] 1 1
##
## , , 89
##
## dummy level
## [1,] 1 1
##
## , , 90
##
## dummy level
## [1,] 1 1
##
## , , 91
##
## dummy level
## [1,] 1 1
##
## , , 92
##
## dummy level
## [1,] 1 1
##
## , , 93
##
## dummy level
## [1,] 1 1
##
## , , 94
##
## dummy level
## [1,] 1 1
##
## , , 95
##
## dummy level
## [1,] 1 1
##
## , , 96
##
## dummy level
## [1,] 1 1
##
## , , 97
##
## dummy level
## [1,] 1 1
##
## , , 98
##
## dummy level
## [1,] 1 1
##
## , , 99
##
## dummy level
## [1,] 1 1
##
## , , 100
##
## dummy level
## [1,] 1 1y.LLM$T## , , 1
##
## dummy level
## dummy 1 0
## level 0 1y.LLM$R## , , 1
##
## [,1]
## dummy 0
## level 1y.LLM$H## , , 1
##
## [,1]
## [1,] NAy.LLM$Q## , , 1
##
## [,1]
## [1,] NApars <- log( rep(var(Nile, na.rm=TRUE), 2) )
y.fit <- fitSSM( model=y.LLM, inits=pars, method='BFGS' )
y.fit## $optim.out
## $optim.out$par
## [1] 6.453453 9.540353
##
## $optim.out$value
## [1] 366.9446
##
## $optim.out$counts
## function gradient
## 21 9
##
## $optim.out$convergence
## [1] 0
##
## $optim.out$message
## NULL
##
##
## $model
## Call:
## SSModel(formula = y ~ SSMtrend(1, Q = list(NA)) + dummy, data = y,
## 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: 1
## [1] Number of states: 2
## Names of the states:
## [1] dummy level
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.y.out <- KFS(y.fit$model)
# 90% confidence intervals for filtered state
y.KF.CI <- predict(y.fit$model, interval="confidence", level=0.9, filtered=TRUE)
y.KS.CI <- predict(y.fit$model, interval="confidence", level=0.9)
y.KF.CI[1,] <- NA
par(mfrow=c(2,2), cex=0.6)
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend = c("data","filtered state","90% confidence interval"),
col = c("black","blue","blue"), lty = c(1,1,2), lwd = c(1,1,1),bty="n", cex=0.9 )
plot( ts( c(y.out$v[-1]), start=1871), col="blue", lwd=2, xlab="", ylab="", main="forecast error")
abline(h=0, lty=3)
plot( ts( c(y.out$P)[-1], start=1871), col="blue", lwd=2, xlab="", ylab="", main="variance of state")
plot( ts( c(y.out$F)[-1], start=1871), col="blue", lwd=2, xlab="", ylab="", main="variance of forecast error")
MLE & Kalman Filtering
par(mfrow=c(1,2), mar=c(3,3,2,1), cex=0.8)
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","filtered","90% confidence interval"), col=c(1,4,4), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
plot.ts( cbind(y, y.KS.CI, Nile), plot.type="single",
col=c(1,2,2,2,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","smoothed","90% confidence interval"), col=c(1,2,2), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )