add fevd method for VAR models

NAMESPACE

@@ -15,6 +15,7 @@ S3method(as_draws_rvars,mvgam)
 S3method(coef,mvgam)
 S3method(conditional_effects,mvgam)
 S3method(ensemble,mvgam_forecast)
+S3method(fevd,mvgam)
 S3method(find_predictors,mvgam)
 S3method(find_predictors,mvgam_prefit)
 S3method(fitted,mvgam)
@@ -41,6 +42,7 @@ S3method(nuts_params,mvgam)
 S3method(pairs,mvgam)
 S3method(plot,mvgam)
 S3method(plot,mvgam_conditional_effects)
+S3method(plot,mvgam_fevd)
 S3method(plot,mvgam_forecast)
 S3method(plot,mvgam_irf)
 S3method(plot,mvgam_lfo)
@@ -91,6 +93,7 @@ export(eval_mvgam)
 export(eval_smoothDothilbertDotsmooth)
 export(eval_smoothDotmodDotsmooth)
 export(eval_smoothDotmoiDotsmooth)
+export(fevd)
 export(forecast)
 export(get_mvgam_priors)
 export(gp)
@@ -180,6 +183,11 @@ importFrom(brms,set_prior)
 importFrom(brms,stancode)
 importFrom(brms,standata)
 importFrom(brms,student)
+importFrom(ggplot2,aes)
+importFrom(ggplot2,facet_wrap)
+importFrom(ggplot2,geom_bar)
+importFrom(ggplot2,ggplot)
+importFrom(ggplot2,labs)
 importFrom(ggplot2,scale_colour_discrete)
 importFrom(ggplot2,scale_fill_discrete)
 importFrom(ggplot2,theme_classic)

NEWS.md

@@ -1,6 +1,7 @@
 # mvgam 1.1.4 (development version; not yet on CRAN)
 ## New functionalities
 * Added a `stability.mvgam` method to compute stability metrics from models fit with Vector Autoregressive dynamics (#21)
+* Added a `fevd.mvgam` method to compute forecast error variance decompositions from models fit with Vector Autoregressive dynamics
 
 ## Bug fixes
 * Fixed a minor bug in the way `trend_map` recognises levels of the `series` factor

R/fevd.mvgam.R

@@ -0,0 +1,153 @@
+#' Calculate latent VAR forecast error variance decompositions
+#'
+#' Compute forecast error variance decompositions from
+#' \code{mvgam} models with Vector Autoregressive dynamics
+#'
+#'@name fevd.mvgam
+#'@param object \code{list} object of class \code{mvgam} resulting from a call to [mvgam()]
+#'that used a Vector Autoregressive latent process model (either as `VAR(cor = FALSE)` or
+#'`VAR(cor = TRUE)`)
+#'@param h Positive \code{integer} specifying the forecast horizon over which to calculate
+#'the IRF
+#'@param ... ignored
+#'@details A forecast error variance decomposition is useful for quantifying the amount
+#'of information each series that in a Vector Autoregression contributes to the forecast
+#'distributions of the other series in the autoregression. This function calculates
+#'the forecast error variance decomposition using the
+#'orthogonalised impulse response coefficient matrices \eqn{\Psi_h}, which can be used to
+#'quantify the contribution of series \eqn{j} to the
+#'h-step forecast error variance of series \eqn{k}:
+#'\deqn{
+#'\sigma_k^2(h) = \sum_{j=1}^K(\psi_{kj, 0}^2 + \ldots + \psi_{kj,
+#' h-1}^2) \quad
+#' }
+#' If the orthogonalised impulse reponses \eqn{(\psi_{kj, 0}^2 + \ldots + \psi_{kj, h-1}^2)}
+#'are divided by the variance of the forecast error \eqn{\sigma_k^2(h)},
+#'this yields an interpretable percentage representing how much of the
+#'forecast error variance for \eqn{k} can be explained by an exogenous shock to \eqn{j}.
+#'@return An object of class \code{mvgam_fevd} containing the posterior forecast error
+#'variance decompositions. This
+#'object can be used with the supplied S3 functions \code{plot}
+#'@author Nicholas J Clark
+#'@references Lütkepohl, H (2006).
+#'New Introduction to Multiple Time Series Analysis. Springer, New York.
+#'@seealso \code{\link{VAR}}, \code{\link{plot.mvgam_irf}}, \code{\link{stability}}
+#' @examples
+#' \donttest{
+#' # Simulate some time series that follow a latent VAR(1) process
+#' simdat <- sim_mvgam(family = gaussian(),
+#'                     n_series = 4,
+#'                     trend_model = VAR(cor = TRUE),
+#'                     prop_trend = 1)
+#' plot_mvgam_series(data = simdat$data_train, series = 'all')
+#'
+#' # Fit a model that uses a latent VAR(1)
+#' mod <- mvgam(y ~ -1,
+#'              trend_formula = ~ 1,
+#'              trend_model = VAR(cor = TRUE),
+#'              family = gaussian(),
+#'              data = simdat$data_train,
+#'              silent = 2)
+#'
+#' # Calulate forecast error variance decompositions for each series
+#' fevds <- fevd(mod, h = 12)
+#'
+#' # Plot median contributions to forecast error variance
+#' plot(fevds)
+#' }
+#'@export
+fevd <- function(object, ...){
+  UseMethod("fevd", object)
+}
+
+#'@rdname fevd.mvgam
+#'@method fevd mvgam
+#'@export
+fevd.mvgam <- function(object,
+                      h = 1,
+                      ...){
+  validate_pos_integer(h)
+  trend_model <- attr(object$model_data, 'trend_model')
+  if(!trend_model %in% c('VAR', 'VARcor', 'VAR1', 'VAR1cor')){
+    stop('Only VAR(1) models currently supported for calculating IRFs',
+         call. = FALSE)
+  }
+  beta_vars <- mcmc_chains(object$model_output, 'A')
+  sigmas <- mcmc_chains(object$model_output, 'Sigma')
+  n_series <- object$n_lv
+
+  if(is.null(n_series)){
+    n_series <- nlevels(object$obs_data$series)
+  }
+
+  all_fevds <- lapply(seq_len(NROW(beta_vars)), function(draw){
+
+    # Get necessary VAR parameters into a simple list format
+    x <- list(K = n_series,
+              A = matrix(beta_vars[draw,],
+                         nrow = n_series,
+                         ncol = n_series,
+                         byrow = TRUE),
+              Sigma = matrix(sigmas[draw,],
+                             nrow = n_series,
+                             ncol = n_series,
+                             byrow = TRUE),
+              p = 1)
+
+    # Calculate the IRF
+    gen_fevd(x, h = h)
+  })
+  class(all_fevds) <- 'mvgam_fevd'
+  return(all_fevds)
+}
+
+#### Functions to compute forecast error variance decompositions
+# Much of this code is modified from R code generously provided in the vars
+# package https://github.com/cran/vars ####
+#' Forecast error variance decomposition
+#' @noRd
+gen_fevd <- function(x, h = 6, ...){
+
+    K <- x$K
+    ynames <- paste0('process_', 1:K)
+    msey <- var_fecov(x, h = h)
+    Phi <- var_phi(x, h = h)
+    mse <- matrix(NA, nrow = h, ncol = K)
+    Omega <- array(0, dim = c(h, K, K))
+    for(i in 1 : h){
+      mse[i, ] <- diag(msey[, , i])
+      temp <- matrix(0, K, K)
+      for(j in 1 : i){
+        temp <- temp + Phi[ , , j]^2
+      }
+      temp <- temp / mse[i, ]
+      for(j in 1 : K){
+        Omega[i, ,j] <- temp[j, ]
+      }
+    }
+    result <- list()
+    for(i in 1 : K){
+      result[[i]] <- matrix(Omega[, , i], nrow = h, ncol = K)
+      colnames(result[[i]]) <- ynames
+    }
+    names(result) <- ynames
+    return(result)
+  }
+
+#' Forecast error covariance matrix
+#' @noRd
+var_fecov = function(x, h) {
+  sigma_yh <- array(NA, dim = c(x$K, x$K, h))
+  Phi <- var_phi(x, h = h)
+  sigma_yh[, , 1] <- Phi[, , 1] %*% t(Phi[, , 1])
+  if (h > 1) {
+    for (i in 2:h) {
+      temp <- matrix(0, nrow = x$K, ncol = x$K)
+      for (j in 2:i) {
+        temp <- temp + Phi[, , j] %*% t(Phi[, , j])
+      }
+      sigma_yh[, , i] <- temp + sigma_yh[, , 1]
+    }
+  }
+  return(sigma_yh)
+}

R/irf.mvgam.R

@@ -217,22 +217,3 @@ var_psi = function(x, h=10){
   }
   return(Psi)
 }
-
-#' Forecast error decomposition
-#' @noRd
-var_fecov = function(x, h) {
-  sigma_u <- x$Sigma
-  sigma_yh <- array(NA, dim = c(x$K, x$K, h))
-  sigma_yh[, , 1] <- sigma_u
-  Phi <- var_phi(x, h = h)
-  if (h > 1) {
-    for (i in 2:h) {
-      temp <- matrix(0, nrow = x$K, ncol = x$K)
-      for (j in 2:i) {
-        temp <- temp + Phi[, , j] %*% sigma_u %*% t(Phi[, , j])
-      }
-      sigma_yh[, , i] <- temp + sigma_yh[, , 1]
-    }
-  }
-  return(sigma_yh)
-}

R/mvgam_fevd-class.R

@@ -0,0 +1,62 @@
+#' `mvgam_fevd` object description
+#'
+#' A \code{mvgam_fevd} object returned by function \code{\link{fevd}}.
+#' Run `methods(class = "mvgam_fevd")` to see an overview of available methods.
+#' @details A `mvgam_fevd` object contains a list of posterior forecast
+#' error variance decompositions, each stored as
+#' its own list
+#' @seealso [mvgam], [VAR]
+#' @author Nicholas J Clark
+#' @name mvgam_fevd-class
+NULL
+
+#'Plot forecast error variance decompositions from an mvgam_fevd object
+
+#'This function takes an \code{mvgam_fevd} object and produces
+#'a plot of the posterior median contributions to forecast variance for each series
+#'in the fitted Vector Autoregression
+#'
+#'@importFrom ggplot2 ggplot aes geom_bar facet_wrap labs
+#'@param x \code{list} object of class \code{mvgam_fevd}. See [fevd()]
+#'@param ... ignored
+#'@return A \code{\link[ggplot2]{ggplot}} object,
+#'which can be further customized using the \pkg{ggplot2} package
+#'@author Nicholas J Clark
+#'@export
+plot.mvgam_fevd = function(x, ...){
+
+  ynames <- names(x[[1]])
+  fevd_df = function(x){
+    do.call(rbind, lapply(seq_len(length(x)), function(process){
+      data.frame(horizon = 1:NROW(x[[process]]),
+                 evd = as.vector(x[[process]]),
+                 Series = sort(rep(colnames(x[[process]]),
+                                   NROW(x[[process]]))),
+                 target = ynames[process])
+    }))
+  }
+
+  # Calculate posterior median error variance contributions
+  fevd_dfs <- do.call(rbind, lapply(seq_len(length(x)),
+                                    function(draw){
+                                      fevd_df(x[[draw]])
+                                    })) %>%
+    dplyr::group_by(horizon, target, Series) %>%
+    dplyr::summarise(mean_evd = median(evd)) %>%
+    dplyr::ungroup() %>%
+    dplyr::group_by(horizon, target) %>%
+    dplyr::mutate(total_evd = sum(mean_evd),
+                  mean_evd = mean_evd / total_evd) %>%
+    dplyr::ungroup() -> mean_evds
+
+  # Plot as a ggplot object
+  ggplot2::ggplot(mean_evds,
+                  ggplot2::aes(fill = Series,
+                               y = mean_evd,
+                               x = horizon)) +
+    ggplot2::geom_bar(position = "stack", stat = "identity") +
+    ggplot2::facet_wrap(~target) +
+    ggplot2::labs(x = 'Forecast horizon',
+                  y = 'Median contribution to forecast variance')
+
+}