Modify predict function.

Output from stan is a 3d array of shape (p, p + 1, horizon). At each steps ahead h = 1, ..., horizon,
result[:, 1, h] is the point prediction for h steps ahead
result[:, 2:(p+1), h] is the covariance matrix for h steps ahead

Logic:

• No change to existing code to get forecasts for 1 step ahead. This gives us 4 things:
  • a = E(alpha_{n + 1} | Y_n)
  • P = Cov(alpha_{n + 1} | Y_n)
  • d + Z*a = E(y_{n + 1} | Y_n)
  • F = Cov(y_{n + 1} | Y_n) = ZPZ' + H = to_symmetric_matrix(quad_form(P, Z') + H)
• if horizon > 1, for loop from h = 2 to horizon.
  • E[alpha_{n + h} | Y_n] = E[c + T \alpha_{n + h - 1} + R eta_{t + h - 1} | Y_n] = c + T a_{n + h - 1}. Update a = c + T a
  • Cov(alpha_{n + h} | Y_n) = Cov(c + T \alpha_{n + h - 1} + R eta_{t + h - 1} | Y_n) = TPT' + RQR'
  • E(y_{n + h} | Y_n) = E(d + Z alpha_{n + h} + epsilon_{n + h}) = d + Z*a
  • Cov(y_{n + h} | Y_n) = Cov(d + Z alpha_{n + h} + epsilon_{n + h}) = ZPZ' + H = to_symmetric_matrix(quad_form(P, Z') + H)

The intercept parameter will be a single real number. To handle this in the filtering, set up a column vector c of length m with first entry equal to the intercept parameter and all others 0.

in the model section of the stan file, rather than calculating the likelihood, calculate the log predictive density h steps after the last time point. Probably want to write a helper function to do this calculation.

Function:
lssm_constant_forecast_log_score (or similar)

Inputs:

• LSSM matrices
• observed data y
• horizon h

Logic:

• initialize return value = 0
• for t in 1:(n - h):
  • update filtering up through time t
  • obtain mean and covariance of forecast distribution at time t + h
  • return value += log predictive density at time t + h evaluated at y[t + h]

Then in the model section of stan file, call
target += lssm_constant_forecast_log_score(...)

SARIMA.stan, SARIMA_predict.stan, functions.stan to stan from stan_models