kfLGSSM: Kalman forward filtering, prediction, smoothing and...

In ilyaZar/RcppSMCkalman: Kalman Filtering and Backward Simulation via RcppSMC

Kalman forward filtering, prediction, smoothing and likelihood evaluation

Description

Under a Linear Gaussian State Space Model (LGSSM) specification, the function runs an (exact) Kalman forward (marginal) filter, computes the prediction, marginal and joint smoothing densities (the former sometimes referred to as the RST-smoother), and generates an unbiased estimate of the (log-)likelihood.

Usage

kfLGSSM(
  yObs,
  uReg = NULL,
  wReg = NULL,
  A = NULL,
  B = NULL,
  C = NULL,
  D = NULL,
  Q = NULL,
  R = NULL,
  initX = NULL,
  initP = NULL,
  initU = NULL,
  nSimMF = 1,
  nSimPD = 1,
  nSimMS = 1,
  nSimJS = 1,
  computeMFD = TRUE,
  computePDD = FALSE,
  computeMSD = FALSE,
  computeJSD = FALSE,
  computeLLH = FALSE
)

Arguments

yObs	A matrix or vector of measurements (observations): rows: multivariate dimension columns: time series dimension T If Y is a univariate process, yObs can be passed as a vector of length T. If nrow(yObs) = 1, then yObs becomes a vector of length T.
uReg	Matrix (vector) of regressors for the latent state process of dimension ncol(B) x T. For a single regressors uReg is a vector of length T.
wReg	Matrix (vector) of regressors for the measurement process of dimension ncol(D) x T. For a single regressors wReg is a vector of length T.
A	Parameter (or system) matrix of dimension dimX x dimX.
B	Parameter (or system) matrix of dimension dimX x numU.
C	Parameter (or system) matrix of dimension dimY x dimX.
D	Parameter (or system) matrix of dimension dimY x numW.
Q	Error VCM of state process of dimension dimX x dimX.
R	Error VCM of measurement process of dimension dimY x dimY
initX	Initial value for state process. Think of X_{t=0} as a starting/initial condition for a time series running for t=1,...,T. If not specified, then initX is sampled from the prior under the stationary distribution of the latent state process (see Details).
initP	VCM matrix initialization for state process; if not specified, then the prior covariance-matrix under the stationary distribution of the latent states is used (see Details).
initU	Initial values of regressors of the state process specification (cannot be missing).
nSimMF	Number of forward filtering runs; defaults to nSimMF=1.
nSimPD	Number of prediction runs; defaults to nSimPD=1.
nSimMS	Number of marginal smoothing runs; defaults to nSimMS=1.
nSimJS	Number of joint smoothing (backward simulation) runs; defaults to nSimJS=1.
computeMFD	Logical: if TRUE, then the marginal filtering density (means and variances) and nSimMF forward filtering runs are computed and returned.
computePDD	Logical: if TRUE, then the prediction density (means and variances) and nSimPD predictions are
computeMSD	Logical: if TRUE, then the marginal smoothing density (means and variances) and nSimMS marginal smoothing runs are computed and returned.
computeJSD	Logical: if TRUE, then the joint smoothing density (means and variances) and nSimJS joint smoothing (i.e backward simulation) runs are computed and returned.
computeLLH	Logical: if TRUE, then the log-likelihood (data density) is returned.

Details

The amount of computation can be controlled via the computeXXX-flags (e.g. only a forward filter and log-likelihood computations, but skipping prediction the smoothers). If required the nSimXX-type arguments set the number of simulations from above densities; the defaults are set to 1, but if set to 0, then no simulation is performed and only the means and variances of the Gaussian distributions of the computeXXX-quantities are returned.

Model specification

The specification of the LGSSM is of the following form:

x_t = Ax_{t-1} + Bu_t + v_t \\

y_t = Cx_t + Dw_t + \varepsilon_t\;,

with serially uncorrelated state innovations v_t\sim\mathcal{N}(0, Q) and measurement errors \varepsilon_t\sim\mathcal{N}(0,R), for t=1,\ldots,T. The argument names match the previous equations. If an argument is missing, the corresponding component is dropped (i.e. set to zero), and forward filtering, prediction and backward smoothing are run under this specification as usual.

There can be an initial state (vector) value passed via initX as well as an initial "prediction" or VCM matrix at period t=0 via initP for initialization of the algorithm. If these are not provided, the stationary prior is used as a default initializer having the following form:

x_0\sim\mathcal{N}\left(Bu_0(I-A)^{-1},\left(I-A\right)^{-1}Q \left[\left(I-A\right)^{-1}\right]^{\top}\right)\;,

where I is the identity matrix, and u_0 is initU. In the univariate case (length(initX) = 1), the prior reduces to x_0\sim\mathcal{N}\left(\frac{B\cdot u_0}{1-A}, \frac{Q}{(1-A^2)}\right), and B\cdot u_0 is the "dot" or scalar product if the corresponding number of u-type regressors >1.

Note that above specification requires an initial regressor (vector) value initU if initX is missing, but can be dropped if initX is provided.

(Marginal) Filtering and prediction - implemented recursions

1. For t=0: compute \hat{x}_{0|0} and \hat{x}_{0|0} according to:

   \hat{x}_{0|0} = Bu_0\left(I-A\right)^{-1}

   \hat{P}_{0|0} =\left(I-A\right)^{-1}Q \left[\left(I-A\right)^{-1}\right]^{\top}

2. For t=1,\ldots,T: compute:

   • K_t and L_t

     K_t = \hat{P}_{t|t-1} C^{\top}

     L_t = \left[C\hat{P}_{t|t-1} C^{\top} + R\right]^{-1}

     with \hat{x}_{t|t-1} and \hat{P}_{t|t-1} as specified below under prediction.

   • for marginal filtering:

     \hat{x}_{t|t} = \hat{x}_{t|t-1} + K_tL_t \left(y_t - C\hat{x}_{t|t-1}-Dw_t\right)

     \hat{P}_{t|t} = \hat{P}_{t|t-1} - K_t L_t K_t^{\top}

     with \hat{x}_{t|t-1} and \hat{P}_{t|t-1} as specified below under prediction.

   • for prediction:

     \hat{P}_{t|t-1} = A\hat{P}_{t-1|t-1} A^{\top} + Q

     \hat{x}_{t|t-1} = A\hat{x}_{t-1|t-1} + Bu_{t-1}

These recursions are implemented in kfMFPD.

Marginal smoothing - implemented recursions

To obtain the marginal smoothing densities p\left(x_t|y_{1:T}\right)\forall t=1,\ldots,T a complete Kalman forward filtering and prediction pass is needed. The marginal smoothing density for period T equals the marginal filtering density for the same period, i.e. after the Kalman filtering and prediction pass \hat{x}_{T|T} and \hat{P}_{T|T} for p\left(x_T|y_{1:T}\right) are already given. To compute the other marginals, the following recursion is run backwards in time:

For t=T-1, \ldots,1 compute:

1. J_t=\hat{P}_{t|t}A^{\top}\hat{P}_{t+1|t}^{-1}

2. \hat{x}_{t|T}=\hat{x}_{t|t} + J_t \left(\hat{x}_{t+1|T}- \hat{x}_{t+1|t}\right)

3. \hat{P}_{t|T}=\hat{P}_{t|t} + J_t\left(\hat{P}_{t+1|T} - \hat{P}_{t+1|t}\right)J_t^{\top}

Each marginal smoothing density is then given via p\left(x_t|y_{t:T}\right)=\mathcal{N}\left(x_t|\hat{x}_{t|T}, \hat{P}_{t|T}\right).

These recursions are implemented in kfMSD.

Joint smoothing density- implemented backward simulation recursions

Backward simulation is a strategy for generating realizations of a whole trajectory \tilde{x}_{0:T} from the joint smoothing density p\left(x_{0:T}|y_{1:}\right):

• Draw nSimJSD samples from

  \tilde{x}_T\sim p\left(x_T|y_{1:T}\right)=\mathcal{N}\left(x_T|\hat{x}_{T|T}|P_{T:T}\right)

  i.e. from the last iteration of the forward filter.

• Draw nSimJSD samples from backwards in time for t=T-1,\ldots,0:

  \tilde{x}_t\sim p\left(x_t|\tilde{x}_{t+1},y_{1:t}\right) =\mathcal{N}\left(x_t|\hat{\mu}_{t},L_t\right)

with means and variances obtained via the following recursions:

• J_t=\hat{P}_{t|t}A^{\top}\left(A\hat{P}_{t|t}A^{\top} + Q\right)^{-1}

• \mu_t=\hat{x}_{t|t} + J_t \left(x_{t+1} - Bu_t -A\hat{x}_{t|t} \right)

• L_t=\hat{P}_{t|t} - J_t A\hat{P}_{t|t}

After a complete backward sweep, the nSimJSD backward trajectories \tilde{x}_{0:T}^{nSimJS} are (by construction) a realization from the above joint smoothing density. These recursions are implemented in kfJSD.

Data distribution - observed likelihood computation

The data likelihood (if logarithmized also sometimes simply referred to as the log-likelihood), with \theta=A, B, C, D, P, Q, is defined as

p\left(y_{1:T}|\theta\right)\prod_{t=1}^{T} \int p\left(y_t| \theta,x_t\right)p\left(x_t|y_{t:(t-1)}\right)dx_t=\mathcal{N} \left(y_t|C\hat{x}_{t|t-1}+Dw_t, C\hat{P}_{t|t-1}C^{\top} + R\right)

with the logarithmic version given as

\log p\left(y_{1:T}\right) = \sum_{t=1}^T\log\left(\mathcal{N} \left(y_t|C\hat{x}_{t|t-1}+Dw_t, C\hat{P}_{t|t-1}C^{\top} + R\right)\right) = -\frac{Tn_y}{2}\log\left(2\pi\right) + \sum_{t=1}^T \left(-\frac{1}{2} \log\det\lbrace\Lambda_t\rbrace - \frac{1}{2} \left(y_t-C\hat{x}_{t|t-1}-Dw_t\right) \Lambda_t^{-1} \left(y_t-C\hat{x}_{t|t-1}-Dw_t\right)^{\top} \right)

with \Lambda_t = C\hat{P}_{t|t-1}C^{\top} + R. These computations are implemented in kfLLH.

Value

A named list of 4:

• kfMFDout: if computeMFD=TRUE the output produced by kfMFPD with PDSTORE=FALSE (or NULL if computeMFD=FALSE)

• kfMPDout: if computePDD=TRUE the output produced by kfMFPD with PDSTORE=TRUE (or NULL if computePDD=FALSE)

• kfMSDout: if computeMSD=TRUE the output produced by kfMSD (or NULL if computeMSD=FALSE)

• kfJSDout: if computeJSD=TRUE the output produced by kfJSD (or NULL if computeJSD=FALSE)

• kfLLHout: if computeMFD=TRUE the output produced by kfLLH (or NULL if computeMFD=FALSE)

See Also

If simulated data is required, the function dataGenLGSSM can be used to simulate from a LGSSM.

To use real data that comes with this package: data(IBMdata), see IBMdata.

ilyaZar/RcppSMCkalman documentation built on Oct. 19, 2023, 11 a.m.