astsa (version 2.4)

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

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

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

Author

D.S. Stoffer

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.

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

y = ar1miss   # an 'astsa' data set with missing values
A = 1; Phi = .1; Ups = 1
Q = 0.1; R = 0.1; Gam=1
mu0 = 0; Sigma0 = 1
input = rep(1, length(y))
( em = EM(y, A, mu0, Sigma0, Phi, Q, R, Ups, Gam, input) ) 

# run Kalman smoother at the final estimates
ks = Ksmooth(y, A=1, em$mu0, em$Sigma0, em$Phi, sQ=sqrt(em$Q), sR=sqrt(em$R),
              Ups=em$Ups, Gam=em$Gam, input=input) 

# admire your work
tsplot(cbind(y, drop(ks$Xs)), type='o', pch=c(19,NA), col=2*2:3, lwd=1:2, 
        gg=TRUE, spag=TRUE, addLegend=TRUE, legend=c('data', 'smooth'))
miss = cbind(which(is.na(y)), min(y, na.rm=TRUE)-.2)   # show yourself
points(miss, pch=18)