EM: EM Algorithm for State Space Models

In astsa: Applied Statistical Time Series Analysis

EM Algorithm for State Space Models

Description

Estimation of the parameters in general linear state space models via the EM algorithm. Missing data may be entered as NA or as zero (0), however, use NAs if zero (0) can be an observation. Inputs in both the state and observation equations are allowed. This script replaces EM0 and EM1.

Usage

EM(y, A, mu0, Sigma0, Phi, Q, R, Ups = NULL, Gam = NULL, input = NULL, 
    max.iter = 100, tol = 1e-04)

Arguments

y	data matrix (n x q), vector or time series, n = number of observations, q = number of series. Use NA or zero (0) for missing data, however, use NAs if zero (0) can be an observation.
A	measurement matrices; can be constant or an array with dimension dim=c(q,p,n) if time varying. Use NA or zero (0) for missing data.
mu0	initial state mean vector (p x 1)
Sigma0	initial state covariance matrix (p x p)
Phi	state transition matrix (p x p)
Q	state error matrix (p x p)
R	observation error matrix (q x q - diagonal only)
Ups	state input matrix (p x r); leave as NULL (default) if not needed
Gam	observation input matrix (q x r); leave as NULL (default) if not needed
input	NULL (default) if not needed or a matrix (n x r) of inputs having the same row dimension (n) as y
max.iter	maximum number of iterations
tol	relative tolerance for determining convergence

Details

This script replaces EM0 and EM1 by combining all cases and allowing inputs in the state and observation equations. It uses version 1 of the new Ksmooth script (hence correlated errors is not allowed).

The states x_t are p-dimensional, the data y_t are q-dimensional, and the inputs u_t are r-dimensional for t=1, \dots, n. The initial state is x_0 \sim N(\mu_0, \Sigma_0).

The general model is

x_t = \Phi x_{t-1} + \Upsilon u_{t} + w_t \quad w_t \sim iid\ N(0, Q)

y_t = A_t x_{t-1} + \Gamma u_{t} + v_t \quad v_t \sim iid\ N(0, R)

where w_t \perp v_t. The observation noise covariance matrix is assumed to be diagonal and it is forced to diagonal otherwise.

The measurement matrices A_t can be constant or time varying. If time varying, they should be entered as an array of dimension dim = c(q,p,n). Otherwise, just enter the constant value making sure it has the appropriate q \times p dimension.

Value

Phi	Estimate of Phi
Q	Estimate of Q
R	Estimate of R
Ups	Estimate of Upsilon (NULL if not used)
Gam	Estimate of Gamma (NULL if not used)
mu0	Estimate of initial state mean
Sigma0	Estimate of initial state covariance matrix
like	-log likelihood at each iteration
niter	number of iterations to convergence
cvg	relative tolerance at convergence

Note

The script does not allow for constrained estimation directly, however, constrained estimation is possible with some extra manipulations. There is an example of constrained estimation using EM at FUN WITH ASTSA, where the fun never stops.

Author(s)

D.S. Stoffer

References

You can find demonstrations of astsa capabilities at FUN WITH ASTSA.

The most recent version of the package can be found at https://github.com/nickpoison/astsa/.

In addition, the News and ChangeLog files are at https://github.com/nickpoison/astsa/blob/master/NEWS.md.

The webpages for the texts and some help on using R for time series analysis can be found at https://nickpoison.github.io/.

See Also

Kfilter, Ksmooth

Examples

# example used for ssm() 
# x[t] = Ups + Phi x[t-1] + w[t]
# y[t] = x[t] + v[t]
y = gtemp_land  
A = 1; Phi = 1; Ups = 0.01
Q = 0.001; R = 0.01
mu0 = -0.6; Sigma0 = 0.02
input = rep(1, length(y))
( em = EM(y, A, mu0, Sigma0, Phi, Q, R, Ups, Gam=NULL, input) ) 

astsa documentation built on May 29, 2024, 10:29 a.m.