These functions run the iterative equations of the Kalman filter
for a state space model.
Usage
## S3 method for class 'stsmSS'
predict(object, y, n.ahead = 12L, ...)
Arguments
object
a list containing the matrices of the state space model.
y
a numeric time series.
n.ahead
a numeric. The number of steps ahead to predict.
...
further arguments. Currently omitted.
Details
This function computes the same values as the
function predict.StructTS from the stats package but
the predictions of the components are also returned.
Value
A list containing the following elements:
itempreda time series containing n.ahead predictions.
itemsea time series containing the standard errors of pred.
itemaa univariate or multivariate time series object containing n.ahead
predictions for the state vector.
itemPa univariate or multivariate time series object containing the square of
the standard errors of a.
References
Harvey, A. C. (1989).
Forecasting, Structural Time Series Models and the Kalman Filter.
Cambridge University Press.
Examples
## local level model
## Nile time series
require("stsm")
y <- Nile
m <- stsm::stsm.model(model = "local-level", y = y, transPars = "StructTS")
fit <- StructTS(y, "level")
m <- stsm::set.pars(m, as.vector(fit$coef[c(2,1)]) * 100 / var(y))
ss <- stsm::char2numeric(m, P0cov = TRUE)
res <- predict(ss, y, 5)
# display forecasts and confidence intervals
plot(cbind(y, res$pred), type = "n", plot.type = "single")
lines(y)
lines(res$pred, col = "blue")
lines(res$pred + 2 * res$se, col = "red", lty = 2)
lines(res$pred - 2 * res$se, col = "red", lty = 2)
# for the whole series, the above is the same as "predict.StructTS"
all.equal(res$pred, predict(fit, 5)$pred)
all.equal(res$se, predict(fit, 5)$se)
## basic Structural model
## AirPassengers time series (in logarithms)
y <- log(AirPassengers)
m <- stsm::stsm.model(model = "BSM", y = y, transPars = "StructTS")
fit <- StructTS(y, "BSM")
m <- stsm::set.pars(m, as.vector(fit$coef[c(4,1:3)]) * 100 / var(y))
ss <- stsm::char2numeric(m, P0cov = TRUE)
res <- predict(ss, y, 12)
all.equal(res$pred, predict(fit, 12)$pred)
all.equal(res$se, predict(fit, 12)$se)
# forecasts and confidence intervals for the series
# scaled back to original scale
expy <- exp(y)
plot(cbind(expy, exp(res$pred + 2 * res$se)), type = "n", plot.type = "single")
lines(expy)
lines(exp(res$pred), col = "blue")
lines(exp(res$pred + 2 * res$se), col = "red", lty = 2)
lines(exp(res$pred - 2 * res$se), col = "red", lty = 2)
# forecasts for the trend component
# the aproach in StructTS() seems to seasonal fluctuations in the trend
# see the "stsm" package for a more flexible interface for maximum likelihood
# procedures to fit a structural time series model
trend <- exp(fitted(fit)[,1])
plot(cbind(trend, AirPassengers), type = "n", plot.type = "single")
lines(AirPassengers, col = "gray")
lines(trend)
lines(exp(res$a[,1]), col = "blue")
lines(exp(res$a[,1] + 2 * sqrt(res$P[,1])), col = "red", lty = 2)
lines(exp(res$a[,1] - 2 * sqrt(res$P[,1])), col = "red", lty = 2)
# forecasts for the seasonal component
seas <- exp(fitted(fit)[,3])
plot(cbind(seas, exp(res$a[,3]) + 2 * sqrt(res$P[,3])),
type = "n", plot.type = "single")
lines(seas)
lines(exp(res$a[,3]), col = "blue")
lines(exp(res$a[,3] + 2 * sqrt(res$P[,3])), col = "red", lty = 2)
lines(exp(res$a[,3] - 2 * sqrt(res$P[,3])), col = "red", lty = 2)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(KFKSDS)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/KFKSDS/predict-stsmSS.Rd_%03d_medium.png", width=480, height=480)
> ### Name: predict.stsmSS
> ### Title: Kalman Filter for State Space Models
> ### Aliases: predict.stsmSS
> ### Keywords: ts, model
>
> ### ** Examples
>
> ## local level model
> ## Nile time series
> require("stsm")
Loading required package: stsm
> y <- Nile
> m <- stsm::stsm.model(model = "local-level", y = y, transPars = "StructTS")
> fit <- StructTS(y, "level")
> m <- stsm::set.pars(m, as.vector(fit$coef[c(2,1)]) * 100 / var(y))
> ss <- stsm::char2numeric(m, P0cov = TRUE)
> res <- predict(ss, y, 5)
>
> # display forecasts and confidence intervals
> plot(cbind(y, res$pred), type = "n", plot.type = "single")
> lines(y)
> lines(res$pred, col = "blue")
> lines(res$pred + 2 * res$se, col = "red", lty = 2)
> lines(res$pred - 2 * res$se, col = "red", lty = 2)
>
> # for the whole series, the above is the same as "predict.StructTS"
> all.equal(res$pred, predict(fit, 5)$pred)
[1] TRUE
> all.equal(res$se, predict(fit, 5)$se)
[1] TRUE
>
> ## basic Structural model
> ## AirPassengers time series (in logarithms)
> y <- log(AirPassengers)
> m <- stsm::stsm.model(model = "BSM", y = y, transPars = "StructTS")
> fit <- StructTS(y, "BSM")
> m <- stsm::set.pars(m, as.vector(fit$coef[c(4,1:3)]) * 100 / var(y))
> ss <- stsm::char2numeric(m, P0cov = TRUE)
> res <- predict(ss, y, 12)
>
> all.equal(res$pred, predict(fit, 12)$pred)
[1] TRUE
> all.equal(res$se, predict(fit, 12)$se)
[1] TRUE
>
> # forecasts and confidence intervals for the series
> # scaled back to original scale
> expy <- exp(y)
> plot(cbind(expy, exp(res$pred + 2 * res$se)), type = "n", plot.type = "single")
> lines(expy)
> lines(exp(res$pred), col = "blue")
> lines(exp(res$pred + 2 * res$se), col = "red", lty = 2)
> lines(exp(res$pred - 2 * res$se), col = "red", lty = 2)
>
> # forecasts for the trend component
> # the aproach in StructTS() seems to seasonal fluctuations in the trend
> # see the "stsm" package for a more flexible interface for maximum likelihood
> # procedures to fit a structural time series model
> trend <- exp(fitted(fit)[,1])
> plot(cbind(trend, AirPassengers), type = "n", plot.type = "single")
> lines(AirPassengers, col = "gray")
> lines(trend)
> lines(exp(res$a[,1]), col = "blue")
> lines(exp(res$a[,1] + 2 * sqrt(res$P[,1])), col = "red", lty = 2)
> lines(exp(res$a[,1] - 2 * sqrt(res$P[,1])), col = "red", lty = 2)
>
> # forecasts for the seasonal component
> seas <- exp(fitted(fit)[,3])
> plot(cbind(seas, exp(res$a[,3]) + 2 * sqrt(res$P[,3])),
+ type = "n", plot.type = "single")
> lines(seas)
> lines(exp(res$a[,3]), col = "blue")
> lines(exp(res$a[,3] + 2 * sqrt(res$P[,3])), col = "red", lty = 2)
> lines(exp(res$a[,3] - 2 * sqrt(res$P[,3])), col = "red", lty = 2)
>
>
>
>
>
> dev.off()
null device
1
>