R/forecast.R

In dsge: Dynamic Stochastic General Equilibrium Models

# Dynamic forecasting for DSGE models

#' Forecast from a Fitted DSGE Model
#'
#' Produces dynamic multi-step forecasts from a fitted DSGE model.
#' Forecasts are generated by iterating the state-space solution forward
#' from the last filtered state, with forecast-error variance computed
#' analytically by iterating the state covariance.
#'
#' @param object A `dsge_fit` object.
#' @param horizon Integer. Number of periods to forecast ahead. Default is 12.
#' @param ... Additional arguments (currently unused).
#'
#' @return An object of class `"dsge_forecast"` containing:
#'   \describe{
#'     \item{forecasts}{Data frame with columns `period`, `variable`,
#'       `value`, and `sd` (one-step forecast standard deviation at each
#'       horizon).  Use `value +/- qnorm(0.5+level/2)*sd` to construct
#'       a confidence band at any level.}
#'     \item{horizon}{The forecast horizon.}
#'     \item{states}{Matrix of forecasted state vectors.}
#'     \item{obs_matrix}{Forecast point estimates in matrix form.}
#'     \item{obs_sd}{Matrix of forecast standard deviations (same shape
#'       as `obs_matrix`).}
#'     \item{history}{Matrix of in-sample observed data (un-demeaned) for
#'       use in plotting forecasts alongside history.}
#'   }
#'
#' @export
forecast.dsge_fit <- function(object, horizon = 12L, ...) {
  sol <- object$solution

  if (!sol$stable) {
    stop("Cannot forecast from an unstable model.", call. = FALSE)
  }

  G <- sol$G
  H <- sol$H
  M <- sol$M
  D <- sol$D

  # Get last filtered state and its covariance
  n_T   <- nrow(object$kalman$filtered_states)
  x_last <- object$kalman$filtered_states[n_T, ]
  P_last <- if (!is.null(object$kalman$filtered_P))
              object$kalman$filtered_P[[n_T]]
            else matrix(0, length(x_last), length(x_last))

  obs_vars <- object$model$variables$observed
  n_obs <- length(obs_vars)
  Z <- D %*% G
  Q <- M %*% t(M)

  # Iterate forward (no shocks in forecast period)
  forecasted_obs    <- matrix(0, horizon, n_obs)
  forecasted_states <- matrix(0, horizon, length(x_last))
  forecast_obs_sd   <- matrix(0, horizon, n_obs)

  x_current <- x_last
  P_current <- P_last

  for (h in seq_len(horizon)) {
    # Mean
    x_current <- as.numeric(H %*% x_current)
    y_current <- as.numeric(Z %*% x_current)
    # Covariance
    P_current <- H %*% P_current %*% t(H) + Q
    P_current <- (P_current + t(P_current)) / 2
    Sigma_y   <- Z %*% P_current %*% t(Z)
    Sigma_y   <- (Sigma_y + t(Sigma_y)) / 2

    forecasted_states[h, ] <- x_current
    forecasted_obs[h, ]    <- y_current
    forecast_obs_sd[h, ]   <- sqrt(pmax(diag(Sigma_y), 0))
  }

  # Add back data means (point forecast levels)
  forecasted_obs <- sweep(forecasted_obs, 2, object$data_means)

  colnames(forecasted_obs)  <- obs_vars
  colnames(forecast_obs_sd) <- obs_vars

  # Build tidy data frame (period, variable, value, sd)
  fc_list <- list()
  for (j in seq_along(obs_vars)) {
    fc_list[[j]] <- data.frame(
      period   = seq_len(horizon),
      variable = obs_vars[j],
      value    = forecasted_obs[, j],
      sd       = forecast_obs_sd[, j],
      stringsAsFactors = FALSE
    )
  }
  fc_df <- do.call(rbind, fc_list)

  # In-sample data (un-demean) for plotting alongside the forecast
  hist_y <- object$data
  if (!is.null(hist_y) && !is.null(object$data_means)) {
    hist_y <- sweep(hist_y, 2, object$data_means, FUN = "+")
  }
  if (!is.null(hist_y)) colnames(hist_y) <- obs_vars

  structure(
    list(
      forecasts  = fc_df,
      horizon    = horizon,
      states     = forecasted_states,
      obs_matrix = forecasted_obs,
      obs_sd     = forecast_obs_sd,
      history    = hist_y
    ),
    class = "dsge_forecast"
  )
}