Homework 8, 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.numeric

library(tseries)
library(forecast)

library(vars)

## Loading required package: MASS

## Loading required package: strucchange

## Loading required package: sandwich

## Loading required package: urca

## Loading required package: lmtest

library(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':
## 
##     startsWith

library(KFAS)
library(zoo)
library(dlmodeler)

Quandl.api_key('CEFP3eWxEJwr_uUP9a2D')

Problem 3

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 ...

##  done

plot(Nile)

x <- datasets::Nile

x.SSM <- SSModel(x~ SSMtrend(1,Q=list(NA)) + SSMseasonal(period=12,sea.type='dummy',Q=NA), H=NA)
x.SSM

## Call:
## SSModel(formula = x ~ 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.

b)

Checking the state variable

x.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

x.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           0

x.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    0

x.SSM$Q

## , , 1
## 
##      [,1] [,2]
## [1,]   NA    0
## [2,]    0   NA

x.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           0

x.SSM$H

## , , 1
## 
##      [,1]
## [1,]   NA

Define update function for maximum likelihood estimation

updatefn_KSI <- function(pars,model,...){
    model$H[] <- exp(pars[1])
    diag(model$Q[,,1])<- exp(c(pars[2:3]))
    model
}

estimating model parameters using maximum likelihood

x.ML <- fitSSM(inits=log(c(var(x),0.001,0.0001)),
             model=x.SSM,
             updatefn=updatefn_KSI,
             method='BFGS')
str(x.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 = x ~ SSMtrend(1, Q = list(NA)) + SSMseasonal(period = 12,      sea.type = "dummy", Q = NA), H = NA)
##  $ terms       :Classes 'terms', 'formula'  language x ~ SSMtrend(1, Q = list(NA)) + SSMseasonal(period = 12, sea.type = "dummy",      Q = NA)
##   .. ..- attr(*, "variables")= language list(x, 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] "x" "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 parameter, variance of innovations in measurement equation, state level equation, seasonal component equation

res <- data.frame(sigma2.measurement = x.ML$model$H[,,1],
                  sigma2.level = diag(x.ML$model$Q[,,1])[1],
                  sigma2.seasonality = diag(x.ML$model$Q[,,1])[2])

x.ML$model$Q

## , , 1
## 
##             [,1]         [,2]
## [1,] 0.001001128 0.000000e+00
## [2,] 0.000000000 9.999991e-05

x.ML$model$H

## , , 1
## 
##          [,1]
## [1,] 30156.08

Kalman filtering and smoothing

x.KFS <- KFS(x.ML$model, filtering=c('state'), smoothing=c('state','disturbance','mean'))
x.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.683

str(x.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(x.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

x.lvl.KS <- x.KFS$alphahat[,"level"]
x.sea.KS  <- x.KFS$alphahat[,"sea_dummy1"]
x.eps.KS  <- x.KFS$epshat



par(mfrow=c(3,1), cex=0.7)

annual data and smoothed state component

tmp <- ts( cbind(x, x.lvl.KS), start=1895, frequency=12 )
plot.ts(tmp, plot.type="single", xlab="",ylab="log of KSI", col=c("green","black"), lwd=c(1,2))
abline(v=1995:2013, lty= "dotted",col="green")
legend("topright", c("original data","Kalman smoothed level"), col=c("green","black"), lwd=c(1,2), bty="n")

smoothed seasonal component

tmp <- ts( x.sea.KS, start=1895, frequency=12)
plot(tmp, xlab="", ylab="log of KSI", col="green", lwd=2)
abline(h=0,col="black")
abline(v=1895:2013, lty= "dotted",col="lightgrey")
legend("topleft", "Kalman smoothed seasonal component", col="green", lwd=2, bty="n")

smoothed irregular component

tmp <- ts( x.eps.KS, start=1895, frequency=12)
plot(tmp,xlab="",ylab="log of KSI", col="blue", lwd=2)
abline(h=0,col="black")
abline(v=1895:2013, lty= "dotted",col="black")
legend("topleft", "irregular component", col="blue", lwd=2, bty="n")