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.
EM(y, A, mu0, Sigma0, Phi, Q, R, Ups = NULL, Gam = NULL, input = NULL,
max.iter = 100, tol = 1e-04)Estimate of Phi
Estimate of Q
Estimate of R
Estimate of Upsilon (NULL if not used)
Estimate of Gamma (NULL if not used)
Estimate of initial state mean
Estimate of initial state covariance matrix
-log likelihood at each iteration
number of iterations to convergence
relative tolerance at convergence
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.
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.
initial state mean vector (p x 1)
initial state covariance matrix (p x p)
state transition matrix (p x p)
state error matrix (p x p)
observation error matrix (q x q - diagonal only)
state input matrix (p x r); leave as NULL (default) if not needed
observation input matrix (q x r); leave as NULL (default) if not needed
NULL (default) if not needed or a
matrix (n x r) of inputs having the same row dimension (n) as y
maximum number of iterations
relative tolerance for determining convergence
D.S. Stoffer
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.
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/.
Kfilter, Ksmooth
# 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) )
Run the code above in your browser using DataLab