How to make out-of-sample forecasts with Dynamic Linear Model of MARSS package?

I'm trying to understand how to use Dynamic Linear Modeling for forecasting. I found an example of the DLM functionality of the MARSS package in R being used for forecasting. Below is all the code in the example, starting with loading the data and ending with creating the in-sample forecasts. What I don't understand is how I would make an out-of-sample forecast? The code below generates "in-sample" forecasts, where it uses already-known information to generate predictions about already-existing data. Say I want to forecast the Salmon Survival tomorrow rather than throughout the last several weeks. How would I do that?

Any help would be appreciated.

# load the data
data(SalmonSurvCUI, package = "MARSS")

# get time indices
years <- SalmonSurvCUI[, 1]

# number of years of data
TT <- length(years)

# get response variable: logit(survival)
dat <- matrix(SalmonSurvCUI[, 2], nrow = 1)

# get predictor variable
CUI <- SalmonSurvCUI[, 3]

## z-score the CUI
CUI.z <- matrix((CUI - mean(CUI))/sqrt(var(CUI)), nrow = 1)

# number of regr params (slope + intercept)
m <- dim(CUI.z)[1] + 1

# for process eqn
B <- diag(m)  ## 2x2; Identity
U <- matrix(0, nrow = m, ncol = 1)  ## 2x1; both elements = 0
Q <- matrix(list(0), m, m)  ## 2x2; all 0 for now
diag(Q) <- c("q.alpha", "q.beta")  ## 2x2; diag = (q1,q2)

# for observation eqn
Z <- array(NA, c(1, m, TT))  ## NxMxT; empty for now
Z[1, 1, ] <- rep(1, TT)  ## Nx1; 1's for intercept
Z[1, 2, ] <- CUI.z  ## Nx1; predictor variable
A <- matrix(0)  ## 1x1; scalar = 0
R <- matrix("r")  ## 1x1; scalar = r

# only need starting values for regr parameters
inits.list <- list(x0 = matrix(c(0, 0), nrow = m))

# list of model matrices & vectors
mod.list <- list(B = B, U = U, Q = Q, Z = Z, A = A, R = R)

# fit univariate DLM
dlm1 <- MARSS(dat, inits = inits.list, model = mod.list)

# get list of Kalman filter output
kf.out <- MARSSkfss(dlm1)
## forecasts of regr parameters; 2xT matrix
eta <- kf.out$xtt1

## ts of E(forecasts)
fore.mean <- vector()
for (t in 1:TT) {
    fore.mean[t] <- Z[, , t] %*% eta[, t, drop = FALSE]
}
# variance of regr parameters; 1x2xT array
Phi <- kf.out$Vtt1
## obs variance; 1x1 matrix
R.est <- coef(dlm1, type = "matrix")$R
## ts of Var(forecasts)
fore.var <- vector()
for (t in 1:TT) {
    tZ <- matrix(Z[, , t], m, 1)  ## transpose of Z
    fore.var[t] <- Z[, , t] %*% Phi[, , t] %*% tZ + R.est
}

Solution

The model of the beta and alpha is a random walk without drift so the prediction of beta(TT+k) and alpha(TT+k) will just be beta(TT) and alpha(TT) where TT is the last time step in the data (in this case CUI.z).

So your prediction is

logit.survival(TT+k) = alpha(TT) + beta(TT)*CUI.z(TT+k)

alpha(TT) and beta(TT) would be output via kf.out$xtT[,TT], i.e. last state estimate. You will need to provide a CUI.z at t=TT+k.

MARSS version 3.11.0 will have predict function and will output these predictions along with the prediction intervals. But release date is sometime late summer 2020. The functionality is in the GitHub development site (under the resids_update branch) but final testing is still being done.