R Forecasting with covariates MARSS package

I've written a model with the MARSS package for R.

The main idea behind the model is to forecast the observable vector for at least 10 quarters, however I can't seem to do it using the MARSSsimulate function (I believe it is because of the inclussion the exogenous vector, named season, in the estimation process). I would really appreciate your help.

Thanks in advance!

The dataset can be downloaded here

I've used the following code

info        <- read.table("series_kalman2.txt",header=T,dec=".")
dat_est_spa <- t(info[,3:6])

Sigma <-  sqrt(apply(dat_est_spa, 1, var, na.rm=TRUE))
y.bar <-  apply(dat_est_spa, 1, mean, na.rm=TRUE)
dat.z <-  (dat_est_spa - y.bar) * (1/Sigma)
rownames(dat.z) = rownames(dat_est_spa)

N.ts <- dim(dat_est_spa)[1]

season <- rbind(rep(c(1,0,0,0),ceiling(dim(dat_est_spa)[2]/4)), 
                rep(c(0,1,0,0),ceiling(dim(dat_est_spa)[2]/4)), 
                rep(c(0,0,1,0),ceiling(dim(dat_est_spa)[2]/4)),
                rep(c(0,0,0,1),ceiling(dim(dat_est_spa)[2]/4)))

rownames(season) <- c("Q1","Q2","Q3","Q4")

season <- season[,-((dim(dat_est_spa)[2]+1):dim(season)[2])]

### Model

cntl.list = list(minit=200, maxit=60000, allow.degen=FALSE)
mod_est_spa   <- list(A="zero", R="diagonal and equal", m=3)
estim_est_spa <- MARSS(dat.z, model=mod_est_spa, control=cntl.list, 
                       form="dfa", covariates=season)

### Forecast

MARSSsimulate(estim_est_spa, tSteps = 10)

Solution

I cannot download the data file from the linked site as it says "permission denied".

Anyway, you are correct in that the inclusion of the season covariate will preclude the use of MARSSsimulate(), but you do have another option that comes with a large caveat.

You are fitting a DFA model with 3 latent trends, which themselves are merely unbiased random walks. Therefore, you could easily simulate the random walks by drawing innovations from a multivariate normal with mean vector

mu = matrix(0, m, 1)

and variance-covariance matrix

Sigma = coef(estim_est_spa, "matrix")$Q.

You can get the estimated states for the last time step T, which would be the starting place for your forecasted states by

X_T = estim_est_spa$states[,dim(dat_est_spa)[2]].

The loadings are

Z = coef(estim_est_spa, "matrix")$Q,

which would need to be rotated. See the DFA example in the MARSS User's Guide for the matrix math.

Note, however, that forecasting with a DFA model is not likely to be very promising because the latent trends are random walks, which generally make very bad forecast models.