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} | Yn] = E[c + T \alpha{n + h - 1} + R eta_{t + h - 1} | Yn] = c + T a{n + h - 1}. Update a = c + T a
Cov(alpha_{n + h} | Yn) = Cov(c + T \alpha{n + h - 1} + R eta_{t + h - 1} | Y_n) = TPT' + RQR'
E(y_{n + h} | Yn) = E(d + Z alpha{n + h} + epsilon_{n + h}) = d + Z*a
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:
to_symmetric_matrix(quad_form(P, Z') + H)a = c + T ato_symmetric_matrix(quad_form(P, Z') + H)