KF: Kalman Filter for State Space Models
In KFKSDS: Kalman Filter, Smoother and Disturbance Smoother
Description
Usage
Arguments
Details
Value
State space representation
References
See Also
Examples
Description
These functions run the iterative equations of the Kalman filter
for a state space model.
Usage
1
2
3
4
5
Arguments
y
a numeric time series or vector.
ss
a list containing the matrices of the state space model.
xreg
optional matrix or list of external regressors. See details.
convergence
a numeric vector of length two to control and determine
the convergence of the filter. See details below.
t0
a numeric indicating the index of the first observation
at which the contributions to the likelihood are addep up.
return.all
logical. If TRUE, extended output containing elements
to be used by KS.deriv is returned.
Details
The implementation is a direct transcription of the iterative equations
of the filter that are summarized below. Details can be found in the references given
below and in many other textbooks.
The source code follows the notation used in Durbin and Koopman (2001).
The elements in the argument ss must be named in accordance with
the notation given below for the state space representation.
For those models defined in the package stsm,
a convenient way to create the argument ss is by means of the
function char2numeric.
The contributions to the likelihood function of the first observations
may be omitted by choosing a value of t0 greater than one.
The functions with ‘.deriv’ in the name compute the derivatives of some
of the elements involved in the filter with respect to the parameters of the model.
The functions KF and KF.deriv are fully implemented
in R while KF.deriv.C calls to compiled C code.
Using KF.deriv.C with return.all = FALSE returns minimal output with
the elements necessary to compute the derivative of the log-likelihood function.
Using return.all = TRUE further elements to be used in
KS.deriv are returned.
Missing observations are allowed. If a missing value is observed after the
filter has converged then all operations of the filter are run instead of using
steady state values until convergence is detected again.
Parameters to control the convergence of the filter.
In some models, the Kalman filter may converge to a steady state. Finding the
explicit expression of the steady state values can be cumbersome in some models.
Alternatively, at each iteration of the filter it can be checked whether a steady state
has been reached. For it, some control parameters can be defined in the argument convergence.
It is considered that convergence was reached when the following is observed:
the change in the variance of the prediction error over the last convergence[2]
consecutive iterations of the filter is below the tolerance value convergence[1].
When the steady state is reached, the values from the last iteration are used
in the remaining iteration of the filter.
Thus, the number of matrix operations can be reduced substantially as pointed
in Harvey (1989) Sec. 3.3.4.
External regressors.
A matrix of external regressors can be passed in argument xreg.
If xreg is a matrix then it is assumed that the time series
passed in argument y has been already adjusted for the effect
of these regressors, that is, y_t^{adj} = y_t - X γ.
If y is the observed series, then xreg should be a list containing
the following elements: xreg, the matrix of regressors; and
coefs, the vector of coefficients, γ, related to the regressors.
The coefficients must be placed in coefs in the same order
as the corresponding vectors are arranged by columns in xreg.
The number of rows of the matrix of regressors must be equal to the length
of the series y.
Column names are necessary for KF.deriv and
are optional for KF.deriv.C.
Value
A list containing the following elements:
v
prediction error.
f
variance of the prediction error.
K
Kalman gain.
L
auxiliar matrices to be used by the smoother.
a.pred
one step ahead prediction of the state vector.
P.pred
covariance matrix of a.pred.
a.upd
update of a.pred after the observation at time t
that is predicted in a.pred enters in the recursive filter.
P.upd
update of P.pred.
convit
the iteration at which the filter converged. If convergence
was not observed it is set to NULL.
mll
value of the negative of the log-likelihood function.
The function KF.C is a simplified and faster version that
records and returns only the value of negative of the log-likelihood function.
It is suited to be passed as argument to optim in the
stats package.
The functions that evaluate the derivatives include in the output the derivatem
terms: da.pred, dP.pred, da.upd, dP.upd, dv, df,
dvof (derivative of quotient v over f), dK and dL.
KF.deriv.C does not return a.upd and P.upd and their
derivative terms da.upd and dP.upd.
If return.all = TRUE, this function returns: dvof,
dL, da.pred, dP.pred, which are the derivative terms necessary
to evaluate the gradient of the log-likelihood function. By default they are not returned,
return.all = FALSE. They are in any case computed, the operations that are omitted
in the latter case is the arrangement of the output from the call to compiled C code into
matrices of the proper dimension containing the data in the right order.
Derivatives of the likelihood function are implemented in package stsm.
Although the Kalman filter provides information to evaluate the likelihood function,
it is not its primary objective. That's why the derivatives of the likelihood are
currently part of the package stsm, which is specific to likelihood methods
in structural time series models.
State space representation
The general univariate linear Gaussian state space model is defined
as follows:
y[t] = Za[t] + e[t], e[t] \sim N(0, H)
a[t+1] = Ta[t] + Rw[t], w[t] \sim N(0, V)
for t=1,…,n and a[1] \sim N(a0, P0).
Z is a matrix of dimension 1xm;
H is 1x1;
T is mxm;
R is mxr;
V is rxr;
a0 is mx1 and
P0 is mxm,
where m is the dimension of the state vector a and
r is the number of variance parameters in the state vector.
The Kalman filtering recursions for the model above are:
Prediction
a[t] = T a[t-1]
P[t] = T P[t-1] T' + R V R'
v[t] = y[t] - Z a[t]
F[t] = Z P[t] Z' + H
Updating
K[t] = P[t] Z' F[t]^{-1}
a[t] = a[t] + K[t] v[t]
P[t] = P[t] - K[t] Z P[t]'
for t=2,…,n, starting with a[1] and P[1] equal
to a0 and P0. v[t] is the prediction error at observation
in time t and F[t] is the variance of v[t].
References
Durbin, J. and Koopman, S. J. (2001).
Time Series Analysis by State Space Methods.
Oxford University Press.
Harvey, A. C. (1989).
Forecasting, Structural Time Series Models and the Kalman Filter.
Cambridge University Press.
See Also
char2numeric in package stsm.
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
# local level plus seasonal model with arbitrary parameter values
# for the 'JohnsonJohnson' time series
m <- stsm::stsm.model(model = "llm+seas", y = JohnsonJohnson,
pars = c("var1" = 2, "var2" = 15, "var3" = 30))
ss <- stsm::char2numeric(m)
# run the Kalman filter
kf <- KF(m@y, ss)
plot(kf$a.upd[,1:2], main = "filtered state vector")
# 'KF.C' is a faster version that returns only the
# value of the negative of the likelihood function
kfc <- KF.C(m@y, ss)
all.equal(kf$mll, kfc)
# compute also derivative terms used below
kfd <- KF.deriv(m@y, ss)
all.equal(kfc, kfd$mll)
kfdc <- KF.deriv.C(m@y, ss, return.all = TRUE)
all.equal(kf$mll, kfdc$mll)
# as expected the versions that use compiled C code
# are faster that the versions written fully in R
# e.g. not including derivatives
## Not run:
system.time(for(i in seq_len(10)) kf <- KF(m@y, ss))
system.time(for(i in seq_len(10)) kfc <- KF.C(m@y, ss))
# e.g. including derivatives
system.time(for(i in seq_len(10)) kfd <- KF.deriv(m@y, ss))
system.time(for(i in seq_len(10)) kfdc <- KF.deriv.C(m@y, ss, return.all = TRUE))
## End(Not run)
# compare analytical and numerical derivatives
# they give same results up to a tolerance error
fcn <- function(x, model, type = c("v", "f"))
{
m <- stsm::set.pars(model, x)
ss <- stsm::char2numeric(m)
kf <- KF(m@y, ss)
switch(type, "v" = sum(kf$v), "f" = sum(kf$f))
}
dv <- numDeriv::grad(func = fcn, x = m@pars, model = m, type = "v")
all.equal(dv, colSums(kfd$dv), check.attributes = FALSE)
all.equal(dv, colSums(kfdc$dv), check.attributes = FALSE)
df <- numDeriv::grad(func = fcn, x = m@pars, model = m, type = "f")
all.equal(df, colSums(kfd$df), check.attributes = FALSE)
all.equal(df, colSums(kfdc$df), check.attributes = FALSE)
# compare timings in version written in R with numDeriv::grad
# no calls to compiled C code in either case
# looking at these timings, using analytical derivatives is
# expected to be useful in optimization algorithms
## Not run:
system.time(for (i in seq_len(10))
numdv <- numDeriv::grad(func = fcn, x = m@pars, model = m, type = "v"))
system.time(for(i in seq_len(10)) kfdv <- colSums(KF.deriv(m@y, ss)$dv))
## End(Not run)
# compare timings when convergence is not checked with the case
# when steady state values are used after convergence is observed
# computation time is reduced substantially
## Not run:
n <- length(m@y)
system.time(for(i in seq_len(20)) a <- KF.deriv(m@y, ss, convergence = c(0.001, n)))
system.time(for(i in seq_len(20)) b <- KF.deriv(m@y, ss, convergence = c(0.001, 10)))
# the results are the same up to a tolerance error
all.equal(colSums(a$dv), colSums(b$dv))
## End(Not run)
Example output
[1] TRUE
[1] TRUE
[1] TRUE
user system elapsed
0.090 0.004 0.093
user system elapsed
0.011 0.000 0.012
user system elapsed
0.304 0.005 0.308
user system elapsed
0.127 0.000 0.127
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
user system elapsed
0.738 0.000 0.738
user system elapsed
0.181 0.000 0.180
user system elapsed
0.327 0.000 0.328
user system elapsed
0.216 0.000 0.216
[1] TRUE
KFKSDS documentation built on May 2, 2019, 8:51 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value State space representation References See Also Examples
Description
These functions run the iterative equations of the Kalman filter for a state space model.
Usage
1 2 3 4 5 |
Arguments
y |
a numeric time series or vector. |
ss |
a list containing the matrices of the state space model. |
xreg |
optional matrix or list of external regressors. See details. |
convergence |
a numeric vector of length two to control and determine the convergence of the filter. See details below. |
t0 |
a numeric indicating the index of the first observation at which the contributions to the likelihood are addep up. |
return.all |
logical. If |
Details
The implementation is a direct transcription of the iterative equations of the filter that are summarized below. Details can be found in the references given below and in many other textbooks. The source code follows the notation used in Durbin and Koopman (2001).
The elements in the argument ss must be named in accordance with
the notation given below for the state space representation.
For those models defined in the package stsm,
a convenient way to create the argument ss is by means of the
function char2numeric.
The contributions to the likelihood function of the first observations
may be omitted by choosing a value of t0 greater than one.
The functions with ‘.deriv’ in the name compute the derivatives of some
of the elements involved in the filter with respect to the parameters of the model.
The functions KF and KF.deriv are fully implemented
in R while KF.deriv.C calls to compiled C code.
Using KF.deriv.C with return.all = FALSE returns minimal output with
the elements necessary to compute the derivative of the log-likelihood function.
Using return.all = TRUE further elements to be used in
KS.deriv are returned.
Missing observations are allowed. If a missing value is observed after the filter has converged then all operations of the filter are run instead of using steady state values until convergence is detected again.
Parameters to control the convergence of the filter.
In some models, the Kalman filter may converge to a steady state. Finding the
explicit expression of the steady state values can be cumbersome in some models.
Alternatively, at each iteration of the filter it can be checked whether a steady state
has been reached. For it, some control parameters can be defined in the argument convergence.
It is considered that convergence was reached when the following is observed:
the change in the variance of the prediction error over the last convergence[2]
consecutive iterations of the filter is below the tolerance value convergence[1].
When the steady state is reached, the values from the last iteration are used
in the remaining iteration of the filter.
Thus, the number of matrix operations can be reduced substantially as pointed
in Harvey (1989) Sec. 3.3.4.
External regressors.
A matrix of external regressors can be passed in argument xreg.
If xreg is a matrix then it is assumed that the time series
passed in argument y has been already adjusted for the effect
of these regressors, that is, y_t^{adj} = y_t - X γ.
If y is the observed series, then xreg should be a list containing
the following elements: xreg, the matrix of regressors; and
coefs, the vector of coefficients, γ, related to the regressors.
The coefficients must be placed in coefs in the same order
as the corresponding vectors are arranged by columns in xreg.
The number of rows of the matrix of regressors must be equal to the length
of the series y.
Column names are necessary for KF.deriv and
are optional for KF.deriv.C.
Value
A list containing the following elements:
v |
prediction error. |
f |
variance of the prediction error. |
K |
Kalman gain. |
L |
auxiliar matrices to be used by the smoother. |
a.pred |
one step ahead prediction of the state vector. |
P.pred |
covariance matrix of |
a.upd |
update of |
P.upd |
update of |
convit |
the iteration at which the filter converged. If convergence
was not observed it is set to |
mll |
value of the negative of the log-likelihood function. |
The function KF.C is a simplified and faster version that
records and returns only the value of negative of the log-likelihood function.
It is suited to be passed as argument to optim in the
stats package.
The functions that evaluate the derivatives include in the output the derivatem
terms: da.pred, dP.pred, da.upd, dP.upd, dv, df,
dvof (derivative of quotient v over f), dK and dL.
KF.deriv.C does not return a.upd and P.upd and their
derivative terms da.upd and dP.upd.
If return.all = TRUE, this function returns: dvof,
dL, da.pred, dP.pred, which are the derivative terms necessary
to evaluate the gradient of the log-likelihood function. By default they are not returned,
return.all = FALSE. They are in any case computed, the operations that are omitted
in the latter case is the arrangement of the output from the call to compiled C code into
matrices of the proper dimension containing the data in the right order.
Derivatives of the likelihood function are implemented in package stsm. Although the Kalman filter provides information to evaluate the likelihood function, it is not its primary objective. That's why the derivatives of the likelihood are currently part of the package stsm, which is specific to likelihood methods in structural time series models.
State space representation
The general univariate linear Gaussian state space model is defined as follows:
y[t] = Za[t] + e[t], e[t] \sim N(0, H)
a[t+1] = Ta[t] + Rw[t], w[t] \sim N(0, V)
for t=1,…,n and a[1] \sim N(a0, P0). Z is a matrix of dimension 1xm; H is 1x1; T is mxm; R is mxr; V is rxr; a0 is mx1 and P0 is mxm, where m is the dimension of the state vector a and r is the number of variance parameters in the state vector.
The Kalman filtering recursions for the model above are:
Prediction
a[t] = T a[t-1]
P[t] = T P[t-1] T' + R V R'
v[t] = y[t] - Z a[t]
F[t] = Z P[t] Z' + H
Updating
K[t] = P[t] Z' F[t]^{-1}
a[t] = a[t] + K[t] v[t]
P[t] = P[t] - K[t] Z P[t]'
for t=2,…,n, starting with a[1] and P[1] equal
to a0 and P0. v[t] is the prediction error at observation
in time t and F[t] is the variance of v[t].
References
Durbin, J. and Koopman, S. J. (2001). Time Series Analysis by State Space Methods. Oxford University Press.
Harvey, A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.
See Also
char2numeric in package stsm.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 | # local level plus seasonal model with arbitrary parameter values
# for the 'JohnsonJohnson' time series
m <- stsm::stsm.model(model = "llm+seas", y = JohnsonJohnson,
pars = c("var1" = 2, "var2" = 15, "var3" = 30))
ss <- stsm::char2numeric(m)
# run the Kalman filter
kf <- KF(m@y, ss)
plot(kf$a.upd[,1:2], main = "filtered state vector")
# 'KF.C' is a faster version that returns only the
# value of the negative of the likelihood function
kfc <- KF.C(m@y, ss)
all.equal(kf$mll, kfc)
# compute also derivative terms used below
kfd <- KF.deriv(m@y, ss)
all.equal(kfc, kfd$mll)
kfdc <- KF.deriv.C(m@y, ss, return.all = TRUE)
all.equal(kf$mll, kfdc$mll)
# as expected the versions that use compiled C code
# are faster that the versions written fully in R
# e.g. not including derivatives
## Not run:
system.time(for(i in seq_len(10)) kf <- KF(m@y, ss))
system.time(for(i in seq_len(10)) kfc <- KF.C(m@y, ss))
# e.g. including derivatives
system.time(for(i in seq_len(10)) kfd <- KF.deriv(m@y, ss))
system.time(for(i in seq_len(10)) kfdc <- KF.deriv.C(m@y, ss, return.all = TRUE))
## End(Not run)
# compare analytical and numerical derivatives
# they give same results up to a tolerance error
fcn <- function(x, model, type = c("v", "f"))
{
m <- stsm::set.pars(model, x)
ss <- stsm::char2numeric(m)
kf <- KF(m@y, ss)
switch(type, "v" = sum(kf$v), "f" = sum(kf$f))
}
dv <- numDeriv::grad(func = fcn, x = m@pars, model = m, type = "v")
all.equal(dv, colSums(kfd$dv), check.attributes = FALSE)
all.equal(dv, colSums(kfdc$dv), check.attributes = FALSE)
df <- numDeriv::grad(func = fcn, x = m@pars, model = m, type = "f")
all.equal(df, colSums(kfd$df), check.attributes = FALSE)
all.equal(df, colSums(kfdc$df), check.attributes = FALSE)
# compare timings in version written in R with numDeriv::grad
# no calls to compiled C code in either case
# looking at these timings, using analytical derivatives is
# expected to be useful in optimization algorithms
## Not run:
system.time(for (i in seq_len(10))
numdv <- numDeriv::grad(func = fcn, x = m@pars, model = m, type = "v"))
system.time(for(i in seq_len(10)) kfdv <- colSums(KF.deriv(m@y, ss)$dv))
## End(Not run)
# compare timings when convergence is not checked with the case
# when steady state values are used after convergence is observed
# computation time is reduced substantially
## Not run:
n <- length(m@y)
system.time(for(i in seq_len(20)) a <- KF.deriv(m@y, ss, convergence = c(0.001, n)))
system.time(for(i in seq_len(20)) b <- KF.deriv(m@y, ss, convergence = c(0.001, 10)))
# the results are the same up to a tolerance error
all.equal(colSums(a$dv), colSums(b$dv))
## End(Not run)
|
Example output
[1] TRUE
[1] TRUE
[1] TRUE
user system elapsed
0.090 0.004 0.093
user system elapsed
0.011 0.000 0.012
user system elapsed
0.304 0.005 0.308
user system elapsed
0.127 0.000 0.127
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
user system elapsed
0.738 0.000 0.738
user system elapsed
0.181 0.000 0.180
user system elapsed
0.327 0.000 0.328
user system elapsed
0.216 0.000 0.216
[1] TRUE
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.