EM {astsa}

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) )