Vectorized Multivariate AR-1 State-space Model
MARSS_vectorized.RdThe EM algorithm (MARSSkem) in the MARSS package works by converting the more familiar MARSS model in matrix form into the vectorized form which allows general linear constraints (Holmes 2012).
The vectorized form is:
$$\mathbf{x}(t) = (\mathbf{x}(t-1)^\top \otimes \mathbf{I}_m)(\mathbf{f}_b(t)+\mathbf{D}_b(t)\beta) + (\mathbf{f}_u(t)+\mathbf{D}_u(t)\upsilon) + \mathbf{w}(t), \textrm{ where } \mathbf{W}(t) \sim \textrm{MVN}(0,\textbf{Q}(t))$$
$$\mathbf{y}(t) = (\mathbf{x}(t)^\top \otimes \mathbf{I}_n)(\mathbf{f}_z(t)+\mathbf{D}_z(t)\zeta) + (\mathbf{f}_a(t)+\mathbf{D}_a(t)\alpha) + \mathbf{v}(t), \textrm{ where } \mathbf{V}(t) \sim \textrm{MVN}(0,\textbf{R}(t))$$
$$\mathbf{x}(1) \sim \textrm{MVN}(x0, V0) \textrm{ or } \mathbf{x}(0) \sim \textrm{MVN}(x0, V0)$$
where \(\beta\), \(\upsilon\), \(\zeta\), and \(\alpha\) are column vectors of estimated values, the \(\mathbf{f}\) are column vectors of inputs (fixed values), and the \(\mathbf{D}\) are perturbation matrices that align the estimated values into the right rows. The \(\mathbf{f}\) and \(\mathbf{D}\) are potentially time-varying. \(\otimes\) means kronecker product and \(\mathbf{I}_p\) is a p x p identity matrix.
Normally the user will specify their model in "marxss" form, perhaps with text short-cuts. The "marxss" form is then converted to "marss" form using the conversion function marxss_to_marss(). In "marss" form, the D, d, C, and c information is put in A and U respectively. If there are inputs (d and c), then this will make A and U time-varying. This is unfortunate, because this slows down the EM algorithm considerably due to the unfortunate decision (early on) to store time-varying parameters as 3-dimensional. The functions for the "marss" form (in the file MARSS_marss.R) convert the "marss" form model into vectorized form and prepares the f (fixed) and D (free) matrices that are at the heart of the model specification.
Note, "marss" is a model form. A model form is defined by a collection of form functions discussed in marssMODEL. These functions are not exported to the user, but are called by MARSS() using the argument form. These internal functions convert the users model list into the vectorized form of a MARSS model and do extensive error-checking. "marxss" is also a model form and these models are also stored in vectorized form (See examples below).
References
Holmes, E. E. (2012). Derivation of the EM algorithm for constrained and unconstrained multivariate autoregressive state-space (MARSS) models. Technical Report. arXiv:1302.3919 [stat.ME]
Examples
dat <- t(harborSealWA)
dat <- dat[2:4, ]
MLEobj <- MARSS(dat)
#> Success! abstol and log-log tests passed at 44 iterations.
#> Alert: conv.test.slope.tol is 0.5.
#> Test with smaller values (<0.1) to ensure convergence.
#>
#> MARSS fit is
#> Estimation method: kem
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> Estimation converged in 44 iterations.
#> Log-likelihood: 17.84491
#> AIC: -15.68982 AICc: -10.45173
#>
#> Estimate
#> R.diag 0.00582
#> U.X.SJF 0.06833
#> U.X.SJI 0.07084
#> U.X.EBays 0.04221
#> Q.(X.SJF,X.SJF) 0.04150
#> Q.(X.SJI,X.SJI) 0.01271
#> Q.(X.EBays,X.EBays) 0.00807
#> x0.X.SJF 5.97602
#> x0.X.SJI 6.70656
#> x0.X.EBays 6.63306
#> Initial states (x0) defined at t=0
#>
#> Standard errors have not been calculated.
#> Use MARSSparamCIs to compute CIs and bias estimates.
#>
# free (D) and fixed (f) matrices
names(MLEobj$model$free)
#> [1] "Z" "A" "R" "B" "U" "Q" "x0" "V0" "D" "C" "d" "c" "G" "H" "L"
names(MLEobj$model$fixed)
#> [1] "Z" "A" "R" "B" "U" "Q" "x0" "V0" "D" "C" "d" "c" "G" "H" "L"
# In marss form, the D, C, d, and c matrices are found in A and U
# If there are inputs, this makes U time-varying
names(MLEobj$marss$free)
#> [1] "Z" "A" "R" "B" "U" "Q" "x0" "V0" "G" "H" "L"
names(MLEobj$marss$fixed)
#> [1] "Z" "A" "R" "B" "U" "Q" "x0" "V0" "G" "H" "L"
# par is in marss form so does not have values for D, C, d, or c
names(MLEobj$par)
#> [1] "Z" "A" "R" "B" "U" "Q" "x0" "V0" "G" "H" "L"
# if you need the par in marxss form, you can use print
tmp <- print(MLEobj, what="par", form="marxss", silent=TRUE)
names(tmp)
#> [1] "Z" "A" "R" "B" "U" "Q" "x0" "V0" "G" "H" "L" "C" "D" "c" "d"