Model specification

Model specification is via a list with the structure of each of the model parameters: B, U, C, Q, Z, A, D, R, ξ, andΛ. There is a one-to-one correspondence between the model in and the math version of the model. The model written in matrix form is translated into equivalent matrices (or arrays if time-varying) in code. Matrices that combine fixed and estimated values are specified using a list matrix with numerical values for fixed values and character names for the estimated values.

For the MARSS model above, this would be:

B1 <- matrix(list("b",0,0,"b"),2,2)
U1 <- matrix(0,2,1)
Q1 <- matrix(c("q11","q12","q12","q22"),2,2)
Z1 <- matrix(c(1,0,1,1,1,0),3,2)
A1 <- matrix(list("a1",0,0),3,1)
R1 <- matrix(list("r11",0,0,0,"r",0,0,0,"r"),3,3)
pi1 <- matrix(0,2,1); V1=diag(1,2)
model.list <- list(B=B1,U=U1,Q=Q1,Z=Z1,A=A1,R=R1,x0=pi1,V0=V1,tinitx=0)

The tinitx element tells MARSS whether the initial state for x is at t = 1 (tinitx=1) or t = 0 (tinitx=0).

MARSS has a number of text shortcuts for common parameter forms, such as "diagonal and unequal"; see below for the possible shortcuts.

Data

The data must be entered as a n × T matrix, or a ts() (or mts) object or vector (which will be converted to a n × T matrix).

Fit call

The call to fit the model is.

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

See ?MARSS for all other arguments. The common ones are method to change the fitting method from default of EM (slow). See below. control is used to pass in a list to control fitting, e.g. control=list(maxit=1000) to increase maximum iterations if the fit does not converge.

Example with simulated data:

library(MARSS)
set.seed(1234)
x <- rbind(arima.sim(n=50,list(ar=0.95), sd=0.4), 
           arima.sim(n=50,list(ar=0.95), sd=.02))
y <- Z1 %*% x + matrix(rnorm(3*50,0,0.1), 3, 50)
fit <- MARSS(y, model=model.list, silent=TRUE)
tidy(fit)

##    term     estimate    std.error      conf.low      conf.up
## 1  A.a1 0.0184598413 0.0269628326 -0.0343863395 0.0713060221
## 2 R.r11 0.0188070225 0.0062216241  0.0066128633 0.0310011817
## 3   R.r 0.0115645650 0.0024433098  0.0067757658 0.0163533641
## 4   B.b 0.8605716318 0.0634069752  0.7362962441 0.9848470195
## 5 Q.q11 0.1301323970 0.0288070155  0.0736716842 0.1865931098
## 6 Q.q12 0.0020987299 0.0029632765 -0.0037091853 0.0079066452
## 7 Q.q22 0.0001353807 0.0003058603 -0.0004640944 0.0007348558

Different fitting methods

The EM algorithm in the {MARSS} package is in R and on top of that EM algorithms are famously slow. You can try method="BFGS" and see if that is faster. For some models, it will be much faster and for others slower. The companion package {marssTMB} allows you to fit these models with TMB and will be the fastest, often much faster. Definitely if you are doing Dynamic Factor Analysis or working with large data sets, you will want to use {marssTMB}.

fit1 <- MARSS(y, model=model.list)
fit2 <- MARSS(y, model=model.list, method="BFGS")
fit3 <- MARSS(y, model=model.list, method="TMB")

method="BFGS" and method="TMB" are both using quasi-Newton methods to optimize and these can be sensitive to initial conditions. You can run EM a few iterations use that as initial conditions for BFGS or TMB, and this will guard against poor initial conditions issues.

fit1 <- MARSS(y, model=model.list, control = list(maxit=15))
fit2 <- MARSS(y, model=model.list, method="BFGS", inits = fit1)

Defaults for model list

Form of the model list is list(B=..., U=...) etc.

Tips

Use ggplot2::autoplot(fit) (or plot(fit)) to see a series of standard plots and diagnostics for your model. Use tidy(fit) for parameter estimates or coef(fit). Use fitted(fit) for model (y) estimates and tsSmooth(fit) for states (x) estimates. You can also use fit$states for the states.

Let’s say you specified your model with some text short-cuts, like Q="unconstrained", but you want the list matrix for for a next step. a <- summary(fit$model) returns that list (invisibly). Because the model argument of MARSS() will understand a list of list matrices, you can pass in model=a to specify the model. MARSSkfas(fit, return.kfas.model=TRUE) will return your model in {KFAS} format (class SSModel), thus you can use all the functions available in the {KFAS} package on your model.

Troubleshooting

Try MARSSinfo() if you get errors you don’t understand or fitting is taking a long time to converge. When fitting a model with MARSS(), pass in silent=2 to see what MARSS() is doing. This puts it in verbose mode. Use fit=FALSE to set up a model without fitting. Let’s say you do fit <- MARSS(..., fit=FALSE). Now you can do summary(fit$model) to see what MARSS() thinks you are trying to fit.

You can also try toLatex(fit$model) to make a LaTeX file and pdf version of your model (saved in the working directory). This loads the {Hmisc} package (and all its dependencies) and requires that you are able to process LaTeX files.

Z

Defaults to "identity". Can be a text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", or "onestate", or a length n vector of factors specifying which of the m hidden state time series correspond to which of the n observation time series. May be specified as a n × m list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric n × m matrix to use a custom fixed Z. "onestate" gives a n × 1 matrix of 1s. The text shortcuts all specify n × n matrices.

B

Defaults to "identity". Can be a text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "zero". Can also be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric m × m matrix to use custom fixed B, but in this case all the eigenvalues of B must fall in the unit circle.

U and x0

Defaults to "unequal". Can be a text string, "unconstrained", "equal", "unequal" or "zero". May be specified as a m × 1 list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric m × 1 matrix to use a custom fixed U or x₀.

A

Defaults to "scaling". Can be a text string, "scaling" ,"unconstrained", "equal", “unequal” or “zero”. May be specified as a n x 1 list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric n × 1 matrix to use a custom fixed A. Care must be taken so that the model is not under-constrained and unsolvable model. The default, "scaling", only applies to Z matrices that are design matrices (only 1s and 0s and all rows sum to 1). When a column in Z has multiple 1s, the first row in the A matrix associated with those Z rows is 0 and the other associated A rows have an estimated value. This treats A as an intercept where one intercept for each x (hidden state) is fixed at 0 and any other intercepts associated with that x have an estimated intercept. This ensures a solvable model when Z is a design matrix.

Q, R and V0

Can be a text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", “zero”. May be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric matrix to use a custom fixed matrix.

V0 defaults to "zero", which means that x₀ is treated as an estimated parameter.

D and C

Defaults to "zero" or if no covariates, defaults to "unconstrained". Can also be any of the options available for the variance matrices.

d and c

Defaults to "zero". Numeric matrix. No missing values allowed. Must have 1 column or the same number of columns as the data, The numbers of rows in must match the corresponding D or C matrix.

G and H

Defaults to "identity". Can be specified as a numeric matrix or array for time-varying cases. Size must match the corresponding Q or R matrix.

Covariates

Inputs, aka covariates, are in c and d. The are passed in via the model list and must be a numeric matrix (no missing values). C and D are the estimated parameters, aka covariate effects.

Let’s say you have temperature data and you want to include a linear effect of temperature that is different for each x time series:

temp <- matrix(rnorm(50, seq(0,1,1/50), 0.1),nrow=1)
C1 <- matrix(c("temp1","temp2"),2,1)
model.list$C <- C1
model.list$c <- temp

Fit as normal:

fit <- MARSS(y, model=model.list, method="BFGS")

A seasonal effect can be easily included via sine/cosine pairs. The fourier() function in the {forecast} package simplifies this. You will need to make your data into a ts/mts object.

yts <- ts(t(y), frequency = 12) # requires time down the rows
fcov <- forecast::fourier(yts,1) |> t()

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

If you want a factor effect, then you can use the seasonaldummy() function in {forecast}

yts <- ts(t(y), frequency = 12) # monthly data
mcov <- forecast::seasonaldummy(yts) |> t() # month factor

Linear constraints

Your model can have simple linear constraints within all the parameters except Q, R and Λ. For example 1 + 2a − 3b is a linear constraint. When entering this value for you matrix, you specify this as "1+2*a+-3*b". NOTE: +’s join parts so +-3*b to specify −3b. Anything after * is a parameter. So 1*1 has a parameter called "1". Example, let’s change the B and Q matrices in the previous model to: Q is fixed because Z is estimated and estimating both creates a statistically confounded model (both scale the variance of x).

This would be specified as (notice "1*z1+-1*z2" for z1-z2):

B1 <- matrix(list("-0.1+1*b",0,0,"0.1+1*b"),2,2)
Q1 <- matrix(list(1,0,0,1),2,2)
Z1 <- matrix(list("1*z1+-1*z2",0,"z2","2*z1","z1",0),3,2)
model.list <- list(B=B1,U=U1,Q=Q1,Z=Z1,A=A1,R=R1,x0=pi1,V0=V1,tinitx=0)

Fit as usual:

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

Time-varying parameters

The default model form allows you to pass in a 3-D array for a time-varying parameter (T is the number of time-steps in your data and is the 3rd dimension in the array):
Zeros are allowed on the diagonals of Q, R and Λ. NOTE(!!), the time indexing. Make sure you enter your arrays such that the correct parameter (or input) at time t lines up with x_t; e.g., it is common for state equations to have B_t − 1 lined up with x_t so you might need to enter the B array such that your B_t − 1 is entered at Bt[,,t] in your code.

The length of the 3rd dimension must be the same as your data. For example, say in your mean-reverting random walk model (the example on the first page) you wanted B(2, 2) to be one value before t = 20 and another value after but B(1, 1) to be time constant. You can pass in the following:

TT <- dim(y)[2]
B1 <- array(list(),dim=c(2,2,TT))
B1[,,1:20] <- matrix(list("b",0,0,"b_1"),2,2)
B1[,,21:TT] <- matrix(list("b",0,0,"b_2"),2,2)

Notice the specification is one-to-one to your B_t matrices on paper.