State space form
The following notation is closest to the one of Koopman et al.
The state space model is represented by the transition equation and
the measurement equation. Let \(m\) be the dimension of the state
variable, \(d\) be the dimension of the observations, and \(n\)
the number of observations. The transition equation and the
measurement equation are given by
$$\alpha_{t + 1} = d_t + T_t \cdot \alpha_t + H_t \cdot \eta_t$$
$$y_t = c_t + Z_t \cdot \alpha_t + G_t \cdot \epsilon_t,$$
where \(\eta_t\) and \(\epsilon_t\) are iid
\(N(0, I_m)\) and iid \(N(0, I_d)\),
respectively, and \(\alpha_t\) denotes the state
variable. The parameters admit the following dimensions:
\(\alpha_{t} \in R^{m}\) | \(d_{t} \in R^m\) | \(\eta_{t} \in R^m\) |
\(T_{t} \in R^{m \times m}\) | \(H_{t} \in R^{m \times m}\) | |
\(y_{t} \in R^d\) | \(c_t \in R^d\) | \(\epsilon_{t} \in R^d\) |
\(Z_{t} \in R^{d \times m}\) | \(G_{t} \in R^{d \times d}\) | |
Note that fkf
takes as input HHt
and GGt
which
corresponds to \(H_t H_t^\prime\) and \(G_t G_t^\prime\).
Iteration:
The filter iterations are implemented using the expected values
$$a_{t} = E[\alpha_t | y_1,\ldots,y_{t-1}]$$
$$a_{t|t} = E[\alpha_t | y_1,\ldots,y_{t}]$$
and variances
$$P_{t} = Var[\alpha_t | y_1,\ldots,y_{t-1}]$$
$$P_{t|t} = Var[\alpha_t | y_1,\ldots,y_{t}]$$
of the state \(\alpha_{t}\) in the following way
(for the case of no NA's):
Initialisation: Set \(t=1\) with \(a_{t} = a0\) and \(P_{t}=P0\)
Updating equations:
$$v_t = y_t - c_t - Z_t a_t$$
$$F_t = Z_t P_t Z_t^{\prime} + G_t G_t^\prime$$
$$K_t = P_t Z_t^{\prime} F_{t}^{-1}$$
$$a_{t|t} = a_t + K_t v_t$$
$$P_{t|t} = P_t - P_t Z_t^\prime K_t^\prime$$
Prediction equations:
$$a_{t+1} = d_{t} + T_{t} a_{t|t}$$
$$P_{t+1} = T_{t} P_{t|t} T_{t}^{\prime} + H_t H_t^\prime$$
Next iteration: Set \(t=t+1\) and goto “Updating equations”.
NA-values:
NA-values in the observation matrix yt
are supported. If
particular observations yt[,i]
contain NAs, the NA-values are
removed and the measurement equation is adjusted accordingly. When
the full vector yt[,i]
is missing the Kalman filter reduces to
a prediction step.
Parameters:
The parameters can either be constant or deterministic
time-varying. Assume the number of observations is \(n\)
(i.e. \(y = (y_t)_{t = 1, \ldots, n}, y_t = (y_{t1}, \ldots,
y_{td})\)). Then, the parameters admit the following
classes and dimensions:
dt | either a \(m \times n\) (time-varying) or a \(m \times 1\) (constant) matrix. |
Tt | either a \(m \times m \times n\) or a \(m \times m \times 1\) array. |
HHt | either a \(m \times m \times n\) or a \(m \times m \times 1\) array. |
ct | either a \(d \times n\) or a \(d \times 1\) matrix. |
Zt | either a \(d \times m \times n\) or a \(d \times m \times 1\) array. |
GGt | either a \(d \times d \times n\) or a \(d \times d \times 1\) array. |
yt | a \(d \times n\) matrix. |
BLAS and LAPACK routines used:
The R function fkf
basically wraps the C
-function
FKF
, which entirely relies on linear algebra subroutines
provided by BLAS and LAPACK. The following functions are used:
BLAS: | dcopy , dgemm , daxpy . |
LAPACK: | dpotri , dpotrf . |
FKF
is called through the .Call
interface. Internally,
FKF
extracts the dimensions, allocates memory, and initializes
the R-objects to be returned. FKF
subsequently calls
cfkf
which performs the Kalman filtering.
The only critical part is to compute the inverse of \(F_t\)
and the determinant of \(F_t\). If the inverse can not be
computed, the filter stops and returns the corresponding message in
status
(see Value). If the computation of the
determinant fails, the filter will continue, but the log-likelihood
(element logLik
) will be “NA”.
The inverse is computed in two steps:
First, the Cholesky factorization of \(F_t\) is
calculated by dpotrf
. Second, dpotri
calculates the
inverse based on the output of dpotrf
.
The determinant of \(F_t\) is computed using again the
Cholesky decomposition.
The first element of both at
and Pt
is filled with the
function arguments a0
and P0
, and the last, i.e. the (n +
1)-th, element of at
and Pt
contains the predictions for the next time step.