k41m4n/ssm

This repository provides code in R reproducing examples of the states space models presented in book "An Introduction to State Space Time Series Analysis" by J.J.F. Commandeur and S.J. Koopman.

This repository provides code in R reproducing examples of the states space models presented in book "An Introduction to State Space Time Series Analysis" by Jacques J.F. Commandeur and Siem Jan Koopman.

The repository uses extensively the KFAS package of Jouni Helske which includes computationally efficient functions for Kalman filtering, smoothing, forecasting, and simulation of multivariate exponential family state space models. Additionally, some own functions has been created to facilitate the calculation and presentation of diagnostics.

The code is provided in file "SSM.R" and split into sections corresponding to the following parts of the book:

• Introduction
• Chapter 2 The Local Level Model
• Chapter 3 The Local Linear Trend Model
• Chapter 4 The Local Level Model with Seasonal
• Chapter 5 The Local Level Model with Explanatory Variable
• Chapter 6 The Local Level Model with Intervention Variable
• Chapter 7 The UK Seat Belt and Inflation Models
• Chapter 8 General Treatment of Univariate State Space Models
• Chapter 9 Multivariate Time Series Analysis

In R Studio, each section of the code can be executed with keys CTRL+ALT+T, after placing a cursor in that section. Please make sure to execute the first section of the code including the own defined functions that are used by the other sections of the code.

Below, the code of the stochastic level and slope model of chapter 3 is shown as an example.

Loading data on UK drivers killed or seriously injured (KSI):

dataUKdriversKSI <- log(read.table("UKdriversKSI.txt")) %>% 
  ts(start = 1969, frequency = 12)
head(dataUKdriversKSI, 24)
#>           Jan      Feb      Mar      Apr      May      Jun      Jul
#> 1969 7.430707 7.318540 7.317876 7.233455 7.397562 7.320527 7.351800
#> 1970 7.468513 7.475906 7.448334 7.351158 7.362011 7.326466 7.498316
#>           Aug      Sep      Oct      Nov      Dec
#> 1969 7.396335 7.364547 7.410347 7.674153 7.672292
#> 1970 7.495542 7.449498 7.604894 7.715124 7.815207
tail(dataUKdriversKSI, 24)
#>           Jan      Feb      Mar      Apr      May      Jun      Jul
#> 1983 7.309212 6.963190 7.104965 7.063048 7.119636 6.981006 7.068172
#> 1984 7.213032 7.060476 7.156177 7.012115 7.167809 7.077498 7.108244
#>           Aug      Sep      Oct      Nov      Dec
#> 1983 7.037906 7.263330 7.304516 7.301822 7.321850
#> 1984 7.157735 7.275172 7.362011 7.459915 7.474772

Defining the model using function SSModel() of the KFAS package:

model <- SSModel(dataUKdriversKSI ~ SSMtrend(degree = 2, 
         Q = list(matrix(NA), matrix(NA))),  H = matrix(NA))
ownupdatefn <- function(pars, model){
  model$H[,,1] <- exp(pars[1])
  diag(model$Q[,,1]) <- exp(pars[2:3])
  model
}
(model)
#> Call:
#> SSModel(formula = dataUKdriversKSI ~ SSMtrend(degree = 2, Q = list(matrix(NA), 
#>     matrix(NA))), H = matrix(NA))
#> 
#> State space model object of class SSModel
#> 
#> Dimensions:
#> [1] Number of time points: 192
#> [1] Number of time series: 1
#> [1] Number of disturbances: 2
#> [1] Number of states: 2
#> Names of the states:
#> [1]  level  slope
#> Distributions of the time series:
#> [1]  gaussian
#> 
#> Object is a valid object of class SSModel.

Providing the number of diffuse initial values in the state:

d <- q <- 2 

Defining the number of estimated hyperparameters (two state disturbance variances + irregular disturbance variance):

w <- 3

Providing the autocorrelation lag l for r-statistic (ACF function):

l <- 12

Defining the first k autocorrelations to be used in Q-statistic:

k <- 15

Providing the number of observations:

n <- 192

Fitting the model using function fitSSM() and extracting the output using function KFS() of the KFAS package:

fit <- fitSSM(model, inits = log(c(0.001, 0001, 0001)), method = "BFGS")
outKFS <- KFS(fit$model, smoothing = c("state", "mean", "disturbance"))

Extracting the maximum likelihood using function logLik() of the KFAS package:

(maxLik <- logLik(fit$model)/n)
#> [1] 0.6247902

Calculating the Akaike information criterion (AIC):

(AIC <- (-2*logLik(fit$model)+2*(w+q))/n)
#> [1] -1.197497

Extracting the maximum likelihood estimate of the irregular variance:

(H <- fit$model$H)
#> , , 1
#> 
#>            [,1]
#> [1,] 0.00211807

Extracting the maximum likelihood estimate of the state disturbance variances for level and slope:

(Q <- fit$model$Q)
#> , , 1
#> 
#>           [,1]         [,2]
#> [1,] 0.0121285 0.000000e+00
#> [2,] 0.0000000 2.774437e-09

Extracting the initial values of the smoothed estimates of states using function coef() of the KFAS package:

smoothEstStat <- coef(outKFS)
(initSmoothEstStat <- smoothEstStat[1,])
#>        level        slope 
#> 7.4157359290 0.0002893677

Extracting the values for trend (stochastic level + slope) using function signal() of the KFAS package:

trend <-signal(outKFS, states = "trend")$signal
head(trend, 24)
#>           Jan      Feb      Mar      Apr      May      Jun      Jul
#> 1969 7.415736 7.330297 7.312187 7.261499 7.371391 7.331428 7.353890
#> 1970 7.491953 7.473422 7.440664 7.363989 7.360783 7.350545 7.478187
#>           Aug      Sep      Oct      Nov      Dec
#> 1969 7.388317 7.376828 7.435665 7.639474 7.644703
#> 1970 7.490569 7.474472 7.601383 7.708188 7.775278
tail(trend, 24)
#>           Jan      Feb      Mar      Apr      May      Jun      Jul
#> 1983 7.308694 7.024342 7.090157 7.071173 7.098717 7.006473 7.060062
#> 1984 7.208837 7.089071 7.133044 7.044559 7.141850 7.090491 7.113531
#>           Aug      Sep      Oct      Nov      Dec
#> 1983 7.067213 7.242181 7.296049 7.301432 7.304580
#> 1984 7.166847 7.272340 7.361614 7.448617 7.470927

Showing Figure 3.1. of the book for trend of stochastic linear trend model:

plot(dataUKdriversKSI , xlab = "", ylab = "", lty = 1)
lines(trend, lty = 3)
title(main = "Figure 3.1. Trend of stochastic linear trend model", cex.main = 0.8)
legend("topright",leg = c("log UK drivers KSI", "stochastic level and slope"), 
       cex = 0.5, lty = c(1, 3), horiz = T)

Showing Figure 3.2. of the book for slope of stochastic linear trend model:

plot(smoothEstStat[, "slope"], xlab = "", ylab = "", lty = 1)
title(main = "Figure 3.2. Slope of stochastic linear trend model", 
      cex.main = 0.8)
legend("topleft",leg = "stochastic slope", 
       cex = 0.5, lty = 1, horiz = T)

Extracting auxiliary irregular residuals (non-standardised) using function residuals() of the KFAS package:

irregResid <- residuals(outKFS, "pearson") 
head(irregResid, 24)
#>               Jan          Feb          Mar          Apr          May
#> 1969  0.014971154 -0.011757877  0.005689260 -0.028043196  0.026170167
#> 1970 -0.023439520  0.002484397  0.007669682 -0.012830394  0.001228027
#>               Jun          Jul          Aug          Sep          Oct
#> 1969 -0.010901472 -0.002089700  0.008018528 -0.012281231 -0.025317464
#> 1970 -0.024078894  0.020128715  0.004973270 -0.024974219  0.003511262
#>               Nov          Dec
#> 1969  0.034679100  0.027589003
#> 1970  0.006935629  0.039929223
tail(irregResid, 24)
#>                Jan           Feb           Mar           Apr           May
#> 1983  0.0005184131 -0.0611517047  0.0148088188 -0.0081251822  0.0209190398
#> 1984  0.0041951480 -0.0285946618  0.0231321567 -0.0324432815  0.0259595873
#>                Jun           Jul           Aug           Sep           Oct
#> 1983 -0.0254675168  0.0081097971 -0.0293069304  0.0211484965  0.0084671400
#> 1984 -0.0129927620 -0.0052871789 -0.0091118808  0.0028323560  0.0003965910
#>                Nov           Dec
#> 1983  0.0003903769  0.0172699294
#> 1984  0.0112977453  0.0038452836

Showing Figure 3.3. of the book for irregular component of stochastic trend model:

plot(irregResid  , xlab = "", ylab = "", lty = 2)
abline(h = 0, lty = 1)
title(main = "Figure 3.3. Irregular component of stochastic trend model", cex.main = 0.8)
legend("topright",leg = "irregular",cex = 0.5, lty = 2, horiz = T)

Extracting one-step-ahead prediction residuals (standardised) using function rstandard() of the KFAS package and calculating diagnostic for these residuals using own defined functions qStatistic(), rStatistic(), hStatistic() and nStatistic():

predResid <- rstandard(outKFS) 
qStat <- qStatistic(predResid, k, w)
rStat <- rStatistic(predResid, d, l)
hStat <- hStatistic(predResid, d)
nStat <- nStatistic(predResid, d)

Showing Table 3.2 of the book for diagnostic tests for the local linear trend model applied to the log of the UK drivers KSI using own defined function dTable():

title = "Table 3.2 Diagnostic tests for the local linear trend model applied to \n
the log of the UK drivers KSI"
dTable(qStat, rStat, hStat, nStat, title)
#> Table 3.2 Diagnostic tests for the local linear trend model applied to 
#> 
#> the log of the UK drivers KSI
#> -----------------------------------------------------------------------------
#>                     statistic    value   critical value   asumption satisfied
#> -----------------------------------------------------------------------------
#> independence           Q(15)   100.609            22.36        -
#>                         r(1)     0.005           +-0.15        +
#>                        r(12)     0.532           +-0.15        -
#> homoscedasticity       H(63)     1.058             1.65        +
#> normality                  N    14.946             5.99        -
#> -----------------------------------------------------------------------------

ssm's Issues

Regression Multivariate Local Linear Trend model

Dear,

My dataset consists 4 time series. I already fitted a Local Linear Trend on this data but now I want to add a regression to this model. The only problem is that I want to regress index_new[, 1] on covid_new[, 1], index_new[, 2] on covid_new[, 2] and so on. Is that possible using SSMregression?

Fitting a Multivariate Local Linear Trend model with KFAS

Dear,

I want to fit a Multivariate Local Linear Trend model to my data. My dataset consist of 4 time series but I really don't know how to fit a Multivariate Local Linear Trend model to my data. I already defined something but it is completely wrong. Can you please help me?

covid_new <- ts(data = covid_new, frequency = 1, start = c(1))
ts.plot(window(covid_new[, 1:4]), col = 1:4,
ylab = "Log confirmed cases",
xlab = "Time")
legend("topright",col = 1:4, lty = 1, legend = colnames(covid_new)[1:4])

covid_new <- window(covid_new, start = 1, end = 145)
model <- SSModel(covid_new[, 1:4] ~ SSMtrend(2, Q = list(matrix(NA, 4, 4), matrix(NA, 4, 4))), H = matrix(NA, 4, 4))

updatefn <- function(pars, model, ...) {
diag(H) <- exp(pars[1:4])
Q <- diag(exp(pars[1:4]))
Q[upper.tri(Q)] <- pars[5:10]
model["Q", etas = "level"] <- crossprod(Q)
Q <- diag(exp(pars[11:14]))
Q[upper.tri(Q)] <- pars[15:20]
model["Q", etas = "slope"] <- crossprod(Q)
model
}

init <- chol(cov(covid_new[, 1:4]))
fitinit <- fitSSM(model, updatefn = updatefn,
inits = rep(c(diag(init), init[upper.tri(init)], init[upper.tri(init)]), 2),
method = "BFGS")
-fitinit$optim.out$val

fit <- fitSSM(model, updatefn = updatefn, inits = fitinit$optim.out$par, method = "BFGS", nsim = 250)
-fit$optim.out$val

varcor <- fit$model["Q", etas = "level"]
varcor[upper.tri(varcor)] <- cov2cor(varcor)[upper.tri(varcor)]
print(varcor, digits = 2)

varcor <- fit$model["Q", etas = "slope"]
varcor[upper.tri(varcor)] <- cov2cor(varcor)[upper.tri(varcor)]
print(varcor, digits = 2)

(out <- KFS(fit$model, nsim = 1000))