fks.SP: Fast Kalman Smoother through Sequential Processing

In FKF.SP: Fast Kalman Filtering Through Sequential Processing

Fast Kalman Smoother through Sequential Processing

Description

\loadmathjax

The Kalman smoother is a backwards algorithm that is run after the Kalman filter that allows the user to refine estimates of previous states to produce "smoothed" estimates of state variables. This function performs the "Kalman smoother" algorithm using sequential processing, an approach that can substantially improve processing time over the traditional Kalman filtering/smoothing algorithms. The primary application of Kalman smoothing is in conjunction with expectation-maximization to estimate the parameters of a state space model. This function is called after running fkf.SP. fks.SP wraps the C-function fks_SP which relies upon the linear algebra subroutines of BLAS (Basic Linear Algebra Subprograms).

Usage

fks.SP(FKF.SP_obj)

Arguments

FKF.SP_obj

An S3-object of class "fkf.SP", returned by fkf.SP when using the argument verbose = TRUE.

Details

fks.SP is typically called after the fkf.SP function to calculate "smoothed" estimates of state variables and their corresponding variances. Smoothed estimates are used when utilizing expectation-maximization (EM) to efficiently estimate the parameters of a state space model.

Sequential Processing Kalman smoother solution:

The fks.SP function uses the solution to the Kalman smoother through sequential processing provided in the textbook of Durbin and Koopman (2001).

Given a state space model has been filtered through the sequential processing Kalman filter algorithm described in fkf.SP, the smoother can be reformulated for the univariate series:

\mjdeqn

y_t'=(y_(1,1),y_(1,2),\cdots,y_(1,p_1),y_(2,1),\cdots,y_(t,p_t))y_t'=(y_(1,1),y_(1,2),...y_(1,p_1),y_(2,1),...,y_(t,p_t))

The sequential processing Kalman smoother approach iterates backwards through both observations and time, i.e.: \mjeqni=p_t, \cdots, 1i=p[t], ... ,1 and \mjeqnt=n,\cdots,1t=n, ... ,1, where \mjeqnp_tp_t is the number of observations at time \mjeqntt and \mjeqnnn is the total number of observations.

The initialisations are:

\mjdeqn

r_(n,p_n) = 0r(n,p_n)=0 \mjdeqnN_(n,p_n)=0N(n,p_n)=0

Then, \mjeqnrr and \mjeqnNN are recursively calculated through:

\mjdeqn

L_t,i = I_m - K_t,i Z_t,iL(t,i) = I_m - K(t,i) Z(t,i)

\mjdeqn

r_(t,i-1) = Z_t,i' F_t,i^-1 v_t,i + L_t,i' r_t,ir(t,i-1) = Z(t,i)' F(t,i)^-1 v(t,i) + L(t,i)' r(t,i)

\mjdeqn

N_t,i-1 = Z_t,i' F_t,i^-1 Z_t,i + L_t,i' N_t,i L_t,iN(t,i-1) = Z(t,i)' F(t,i)^-1 Z(t,i) + L(t,i)' N(t,i) L(t,i)

\mjdeqn

r_t-1,p_t = T_t-1' r_t,0r(t-1,p_t) = T(t-1)' r(t,0)

\mjdeqn

N_t-1,p_t = T_t-1' N_t,0 T_t-1N(t-1,p_t) = T(t-1)' N(t,0) T(t-1)

for \mjeqni=p_t,\cdots,1i=p_t, ..., 1 and \mjeqnt=n,\cdots,1t=n, ..., 1

The equations for \mjeqnr_t-1,p_tr(t-1,p_t) and \mjeqnN_t-1,p_tN(t-1,p_t) do not apply for \mjeqnt=1t=1

Under this formulation, the values for \mjeqnr_t,0r(t,0) and \mjeqnN_t,0N(t,0) are the same as the values for the smoothing quantities of \mjeqnr_t-1r(t-1) and \mjeqnN_t-1N(t-1) of the standard smoothing equations, respectively.

The standard smoothing equations for \mjeqn\hata_tahat(t) and \mjeqnV_tV(t) are used:

\mjdeqn\hat

a_t = a_t + P_t r_t-1ahat(t) = a(t) + P(t) r(t-1)

\mjdeqn

V_t = P_t - P_t N_t-1 P_tV(t) = P(t) - P(t) N(t-1) P(t)

Where:

\mjdeqn

a_t=a_t,1a(t) = a(t,1)

\mjdeqn

P_t = P_t,1P(t) = P(t,1)

In the equations above, \mjeqnr_t,ir(t,i) is an \mjeqnm \times 1m X 1 vector, \mjeqnI_mI_m is an \mjeqnm \times mm X m identity matrix, \mjeqnK_t,iK(t,i) is an \mjeqnm \times 1m X 1 column vector, \mjeqnZ_t,iZ(t,i) is a \mjeqn1 \times m1 X m row vector, and both \mjeqnF_t,i^-1F(t,i)^-1 and \mjeqnv_t,iv(t,i) are scalars. The reduced dimensionality of many of the variables in this formulation compared to traditional Kalman smoothing can result in increased computational efficiency.

Finally, in the formulation described above, \mjeqna_ta(t) and \mjeqnP_tP(t) correspond to the values of att and ptt returned from the fkf.SP function, respectively.

Value

An S3-object of class fks.SP, which is a list with the following elements:

ahatt	A m * n-matrix containing the smoothed state variables, i.e. ahatt[,t] = \mjeqna_t|na(t|n)
Vt	A m * m * n-array containing the variances of ahatt, i.e. Vt[,,t] = \mjeqnP_t|nP(t|n)

References

Aspinall, T. W., Harris, G., Gepp, A., Kelly, S., Southam, C., and Vanstone, B. (2022). The Estimation of Commodity Pricing Models with Applications in Capital Investments. Available Online.

Durbin, James, and Siem Jan Koopman (2001). Time series analysis by state space methods. Oxford university press.

Examples

### Perform Kalman Filtering and Smoothing through sequential processing:
#Nile's annual flow:
yt <- Nile

# Incomplete Nile Data - two NA's are present:
yt[c(3, 10)] <- NA

dt <- ct <- matrix(0)
Zt <- Tt <- matrix(1)
a0 <- yt[1]   # Estimation of the first year flow
P0 <- matrix(100)       # Variance of 'a0'

# Parameter estimation - maximum likelihood estimation:
# Unknown parameters initial estimates:
GGt <- HHt <- var(yt, na.rm = TRUE) * .5
HHt = matrix(HHt)
GGt = matrix(GGt)
yt = rbind(yt)
# Filter through the Kalman filter - sequential processing:
Nile_filtered <- fkf.SP(HHt = matrix(HHt), GGt = matrix(GGt), a0 = a0, P0 = P0, dt = dt, ct = ct,
                  Zt = Zt, Tt = Tt, yt = rbind(yt), verbose = TRUE)
# Smooth filtered values through the Kalman smoother - sequential processing:
Smoothed_Estimates <- fks.SP(Nile_filtered)

FKF.SP documentation built on Oct. 10, 2022, 9:06 a.m.