6.1 Fitting a state-space model with MARSS

The MARSS package fits multivariate auto-regressive models of this form: \[\begin{equation} \begin{gathered} \mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1}+\mathbf{u}+\mathbf{w}_t \text{ where } \mathbf{w}_t \sim \,\text{N}(0,\mathbf{Q}) \\ \mathbf{y}_t = \mathbf{Z}\mathbf{x}_t+\mathbf{a}+\mathbf{v}_t \text{ where } \mathbf{v}_t \sim \,\text{N}(0,\mathbf{R}) \\ \mathbf{x}_0 = \boldsymbol{\mu} \end{gathered} \tag{6.1} \end{equation}\] To fit your time series model with the MARSS package, you need to put your model into the form above. The \(\mathbf{B}\), \(\mathbf{Z}\), \(\mathbf{u}\), \(\mathbf{a}\), \(\mathbf{Q}\), \(\mathbf{R}\) and \(\boldsymbol{\mu}\) are parameters that are (potentially) estimated. The \(\mathbf{y}\) are your data. The \(\mathbf{x}\) are the hidden state(s). Everything in bold is a matrix; if it is a small bolded letter, it is a matrix with 1 column.

Important: In the state-space model equation, \(\mathbf{y}\) is always the data and \(\mathbf{x}\) is a hidden random walk estimated from the data.

A basic MARSS() call looks like fit=MARSS(y, model=list(...)). The argument model tells the function what form the parameters take. The list has the elements with the names: B, U, Q, etc. The names correspond to the parameters with the same names in Equation (6.1) except that \(\boldsymbol{\mu}\) is called x0. tinitx indicates whether the initial \(\mathbf{x}\) is specified at \(t=0\) so \(\mathbf{x}_0\) or \(t=1\) so \(\mathbf{x}_1\).

Here’s an example. Let’s say we want to fit a univariate AR(1) model observed with error. Here is that model: \[\begin{equation} \begin{gathered} x_t = b x_{t-1} + w_t \text{ where } \mathbf{w}_t \sim \,\text{N}(0,q) \\ y_t = x_t+v_t \text{ where } v_t \sim \,\text{N}(0,r) \\ x_0 = \mu \end{gathered} \tag{6.2} \end{equation}\]

To fit this with MARSS(), we need to write Equation (6.2) as Equation (6.1). Equation (6.1) is in MATRIX form. In the model list, the parameters must be written EXACTLY like they would be written for Equation (6.1). For example, 1 is the number 1 in R. It is not a matrix:

class(1)

[1] "numeric"

If you need a 1 (or 0) in your model, you need to pass in the parameter as a \(1 \times 1\) matrix: matrix(1).

With that mind, our model list for Equation (6.2) is:

mod.list <- list(B = matrix(1), U = matrix(0), Q = matrix("q"), 
    Z = matrix(1), A = matrix(0), R = matrix("r"), x0 = matrix("mu"), 
    tinitx = 0)

We can simulate some AR(1) plus error data like so

q <- 0.1
r <- 0.1
n <- 100
y <- cumsum(rnorm(n, 0, sqrt(q))) + rnorm(n, 0, sqrt(r))

And then fit with MARSS() using mod.list above:

fit <- MARSS(y, model = mod.list)

Success! abstol and log-log tests passed at 16 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 16 iterations. 
Log-likelihood: -65.70444 
AIC: 137.4089   AICc: 137.6589   
 
      Estimate
R.r     0.1066
Q.q     0.0578
x0.mu  -0.2024
Initial states (x0) defined at t=0

Standard errors have not been calculated. 
Use MARSSparamCIs to compute CIs and bias estimates.

If we wanted to fix \(q=0.1\), then \(\mathbf{Q}=[0.1]\) (a \(1 \times 1\) matrix with 0.1). We just change mod.list$Q and re-fit:

mod.list$Q <- matrix(0.1)
fit <- MARSS(y, model = mod.list)