The forecast combination puzzle

This is a simulation of some of the ideas in a preprint by Claeskens et al, offering a new theoretical explanation of the forecast combination puzzle. The ideas were later published in the International Journal of Forecasting (full reference below).

The paper begins by explaining what the forecast combination puzzle is:

When several forecasts of the same event are available, it is natural to try and find a (linear) combination of the forecasts which is ‘best’ in some sense. Empirical evidence and extensive simulations show that a simple average often performs best– better than a theorectically derived optimum. This finding is known as the ‘forecast combination puzzle.’…

One important fact has, however, been overlooked in all (or almost all) previous research, namely the fact that the optimal weight derivation and its estimation are separated.

Setup

Consider the case of linearly combining two forecasts, y1 and y2, of an event:

\[ y_c = w y_1 + (1-w) y_2 \]

Most studies of the variance of yc, the combined forecast, assume the weight w is fixed. Claeskens et al. compare this with the more realistic case when w is random, thus creating additional variance. The simulations below focus on the specific case where (y1, y2, w) follows a trivariate normal distribution, although the paper explores other cases.

R helper functions

Data generation functions

A function to generate the trivariate normal distribution with a sample size n. (Fuller documentation of this function is in the appendix.)

library(MASS)
#' Generate Y1, Y2, weight from a trivariate normal distribution.
generateData <- function(n, mu=rep(0,3), sigma=diag(3)) {
  mat <- mvrnorm(n=n, mu, sigma)
  df <- data.frame(mat)
  colnames(df) <- c("Y1", "Y2", "weight")
  return(df)
}

A function for the variance of the traditional fixed-weight case– ignoring the random variable in the weight column.

fixedWeightVar <- function(data, fixed_weight) {
  y_combined <- fixed_weight * data$Y1 + (1-fixed_weight) * data$Y2
  df <- data.frame(fixed_weight, var(y_combined), "fixed")
  names(df) <- c("weight", "variance", "type")
  return(df)
}

A function for the variance of a combined forecast with random weights.

randomWeightVar <- function(data, true_mean_of_weight,
                            col_names=c("weight", "variance", "type")) {
  y_combined <- data$weight * data$Y1 + (1-data$weight) * data$Y2
  df <- data.frame(true_mean_of_weight, var(y_combined), "random")
  names(df) <- col_names
  return(df)
}

Simulation functions

simulate <- function(n, mean_rand_weight, fixed_weight, mu=rep(0,3),
                     sigma=diag(3)) {
  data <- generateData(n, mu=c(0 ,0, mean_rand_weight), sigma=sigma)
  results <- data.frame(randomWeightVar(data, mean_rand_weight))
  results<-rbind(results, fixedWeightVar(data, fixed_weight))
  return(results)
}
runSimulations <- function(num_simulations, n, weight_vector=seq(0,1,0.1),
                           mu=rep(0,3), sigma=diag(3)) {
  results <- data.frame()
  # I know for loops are inefficient, but everyone understands them.
  for (mean_rand_weight in weight_vector) {
    for (fixed_weight in weight_vector) {
      for (i in seq(num_simulations)) {
        results<-rbind(
          results,
          simulate(n, mean_rand_weight, fixed_weight, mu=mu, sigma=sigma))
      }
    }
  }
  return(results)
}

Random seed for replication

If you want to reproduce these results exactly, use the same random seed.

set.seed(1)

Running simulations and plotting results

Zero mean, equal variances, independent forecasts

Let’s run a simulation with simple size 10 and the default settings of zero means and no correlations.

library(ggplot2)
results <- runSimulations(10, 200)
ggplot(data=results, aes(x=weight, y=variance, colour=type)) + geom_point() +ylim(0,5)

The key lesson here is that the random weights cause the overall variance of the combined forecast to go up. A fixed weight combination would always have lower variance.

The minimum variance is at a weight of 0.5, meaning equal weights for the two forecasts will give the lowest combined forecast weight.

Zero mean, first has more variance, independent forecasts

Now let’s give a higher variance to the first forecast– double the variance of the 2nd forecast. Since this is a worse predictor, this higher-variance forecast would ideally be weighted less.

results <- runSimulations(10, 200, sigma=diag(c(2,1,1)))
ggplot(data=results, aes(x=weight, y=variance, colour=type)) + geom_point() + ylim(0,5)
## Warning: Removed 162 rows containing missing values (geom_point).

The overall variance has increased, as expected. And the optimal fixed weight has moved distinctly to the left, meaning a lower weight to the first forecast, as expected. (But it’s not as easy to see the optimum with the random weights, which always have higher variance.)

Let’s create orders of magnitue of difference: the first forecast has a variance of 2 while the second has a variance of 0.2

results <- runSimulations(10, 200, sigma=diag(c(2,0.2,1)))
ggplot(data=results, aes(x=weight, y=variance, colour=type)) + geom_point() + ylim(0,5)
## Warning: Removed 29 rows containing missing values (geom_point).

Even with orders of magnitude of difference in variance, the optimal fixed weight is not very far from 0.5 (equal weights), and a weight of 0.5 doesn’t perform much worse than the optimum. Fixed weights only begin to do as poorly as random weights for a very wrong weight, above 0.75.

Zero mean, exploring correlation structures

So far, we have assumed zero correlation everywhere.

Let’s add a correlation between the forecasts (which we are combining) but zero correlation between each forecast and the random weights.

matrix <- rbind(cbind(1, 0.5, 0), cbind(0.5, 1, 0), cbind(0, 0, 1))
matrix
##      [,1] [,2] [,3]
## [1,]  1.0  0.5    0
## [2,]  0.5  1.0    0
## [3,]  0.0  0.0    1
results <- runSimulations(10, 200, sigma=matrix)
ggplot(data=results, aes(x=weight, y=variance, colour=type)) + geom_point() + ylim(0,5)

This raised the variance in the fixed weight combination a bit and lowered variance in the random weight combination a lot, so now random and fixed are much closer.

Let’s try the reverse– each predictor has some correlation with the weight but none with each other.

matrix <- rbind(cbind(1, 0, 0.5), cbind(0, 1, 0.5), cbind(0.5, 0.5, 1))
matrix
##      [,1] [,2] [,3]
## [1,]  1.0  0.0  0.5
## [2,]  0.0  1.0  0.5
## [3,]  0.5  0.5  1.0
results <- runSimulations(10, 200, sigma=matrix)
ggplot(data=results, aes(x=weight, y=variance, colour=type)) + geom_point() + ylim(0,5)
## Warning: Removed 1 rows containing missing values (geom_point).

This looks similar to the original no-correlation case.

Conclusion: if you have highly correlated predictors, then combining them with random weights is almost as good as using fixed weights. But why not just use fixed weights?

References

Claeskens, Gerda & Magnus, Jan R. & Vasnev, Andrey L. & Wang, Wendun, 2016. “The forecast combination puzzle: A simple theoretical explanation,” International Journal of Forecasting, Elsevier, vol. 32(3), pages 754-762.

Appendix: Documentation for data-generating functions

As the documentation took up a lot of rows of data, I have moved it here.

library(MASS)
#' Generate Y1, Y2, weight from a trivariate normal distribution.
#'
#' @param n Sample size to return-- number of rows in the data.frame.
#' @param mu Optional means of Y1, Y2, and weight.  Should be a 3-element
#'           vector.  By default all means are zero.
#' @param sigma Optional covariance matrix.  Should be a 3x3 matrix.  By
#'          default forecasts are independent with variance of 1.
#' @return An n x 3 data.frame with columns Y1, Y2, and weight.
#' @examples
#' set.seed(1) # To reproduce the exact same results as below.
#' generateData(2)
#' #           Y1         Y2     weight
#' # 1  0.3295078 -0.8356286 -0.6264538
#' # 2 -0.8204684  1.5952808  0.1836433
generateData <- function(n, mu=rep(0,3), sigma=diag(3)) {
  mat <- mvrnorm(n=n, mu, sigma)
  df <- data.frame(mat)
  colnames(df) <- c("Y1", "Y2", "weight")
  return(df)
}

#' Calculates variance of a combined forecast with a fixed weight.
#'
#' This uses a combined forecast of
#'    fixed_weight * data$Y1 + (1-fixed_weight) * data$Y2
#'
#' @param data A data frame with columns Y1 and Y2
#' @param fixed_weight The weight to apply to Y1.  (1-fixed_weight) will be
#'             applied to Y2.
#' @return A one-row data.frame with columns weight, variance, and type.
#'            The weight echoes the fixed_weight.  The variance is the
#'            calculated variance of the forecast combination.  The type
#'            is 'fixed', meaning a fixed weight.
fixedWeightVar <- function(data, fixed_weight) {
  y_combined <- fixed_weight * data$Y1 + (1-fixed_weight) * data$Y2
  df <- data.frame(fixed_weight, var(y_combined), "fixed")
  names(df) <- c("weight", "variance", "type")
  return(df)
}

#' Calculates variance of a combined forecast with random weights.
#'
#' This uses a combined forecast of
#'    data$weight * data$Y1 + (1-data$weight) * data$Y2
#'
#' @param data A data frame with columns Y1, Y2, and weight.
#' @param true_mean_of_weight The true mean that generated the weight column.
#'            This is only used to echo back the return data.
#' @return A one-row data.frame with columns weight, variance, and type.
#'            The weight echoes the true_mean_of_weight.  The variance is the
#'            calculated variance of the forecast combination.  The type
#'            is 'random', meaning a random weight.
randomWeightVar <- function(data, true_mean_of_weight,
                            col_names=c("weight", "variance", "type")) {
  y_combined <- data$weight * data$Y1 + (1-data$weight) * data$Y2
  df <- data.frame(true_mean_of_weight, var(y_combined), "random")
  names(df) <- col_names
  return(df)
}