inst/doc/reconciliation_properties.R

In bayesRecon: Probabilistic Reconciliation via Conditioning

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(bayesRecon)

## ----out.width = '50%', echo = FALSE------------------------------------------
knitr::include_graphics("img/finTS_hier.jpg")

## -----------------------------------------------------------------------------
# Hierarchy composed by 6 time series: 5 bottom and 1 upper
n_b <- 5
n_u <- 1
n <- n_b + n_u   

A <- matrix(1, ncol = n_b, nrow = n_u)  # aggregation matrix

# Actual values:
actuals <- data.frame(extr_mkt_events)  # convert to data frame
# Base forecasts:
base_fc <- extr_mkt_events_basefc

N <- nrow(actuals)  # number of days (3508)

# # If you want to run only N reconciliations (instead of 3508):
# N <- 200
# actuals <- actuals[1:N,]
# base_fc$mu <- base_fc$mu[1:N,]

## -----------------------------------------------------------------------------
# We need to save the mean and median of the reconciled distribution
# in order to compute the squared error and the absolute error:
rec_means   <- matrix(NA, ncol = n, nrow = N)
rec_medians <- matrix(NA, ncol = n, nrow = N)

# We need to save the lower and upper quantiles of the reconciled distribution
# in order to compute the interval score:
rec_L <- matrix(NA, ncol = n, nrow = N)
rec_U <- matrix(NA, ncol = n, nrow = N)
int_cov = 0.9   # use 90% interval         
q1 <- (1 - int_cov) / 2
q2 <- (1 + int_cov) / 2

# Set the number of samples to draw from the reconciled distribution:
N_samples <- 1e4

start <- Sys.time()
for (j in 1:N) {
  # Base forecasts:
  base_fc_j <- c()   
  for (i in 1:n) {
    base_fc_j[[i]] <- list(size = base_fc$size[[i]], mu = base_fc$mu[[j,i]])
  }
  
  # Reconcile via importance sampling:
  buis <- reconc_BUIS(A, base_fc_j, "params", "nbinom", 
                      num_samples = N_samples, seed = 42)
  samples_y <- buis$reconciled_samples
  
  # Save mean, median, and lower and upper quantiles:
  rec_means[j,]   <- rowMeans(samples_y)
  rec_medians[j,] <- apply(samples_y, 1, median)
  rec_L[j,]       <- apply(samples_y, 1, quantile, q1)
  rec_U[j,]       <- apply(samples_y, 1, quantile, q2)
}
stop <- Sys.time()
cat("Computational time for ", N, " reconciliations: ",     
    round(difftime(stop, start, units = "secs"), 2), "s")

## -----------------------------------------------------------------------------
base_means   <- base_fc$mu
base_medians <- matrix(NA, ncol = n, nrow = N)
base_L       <- matrix(NA, ncol = n, nrow = N)
base_U       <- matrix(NA, ncol = n, nrow = N)

for (i in 1:n) {
  base_medians[,i] <- sapply(base_means[,i], 
                       function(mu) qnbinom(p=0.5, size=base_fc$size[[i]], mu=mu))
  base_L[,i]       <- sapply(base_means[,i], 
                       function(mu) qnbinom(p=q1,  size=base_fc$size[[i]], mu=mu))
  base_U[,i]       <- sapply(base_means[,i], 
                       function(mu) qnbinom(p=q2,  size=base_fc$size[[i]], mu=mu))
}

## -----------------------------------------------------------------------------
# Compute the squared errors
SE_base <- (base_means - actuals)^2
SE_rec  <- (rec_means - actuals)^2

# Compute the absolute errors
AE_base <- abs(base_medians - actuals)
AE_rec  <- abs(rec_medians - actuals)

# Define the function for the interval score
int_score <- function(l, u, actual, int_cov = 0.9) {
  is <- (u - l) +
  2 / (1-int_cov) * (actual - u) * (actual>u) +
  2 / (1-int_cov) * (l - actual) * (l>actual) 
  return(is)
}

# Compute the interval scores
IS_base <- mapply(int_score, base_L, base_U, data.matrix(actuals))
IS_base <- matrix(IS_base, nrow = N)
IS_rec <- mapply(int_score, rec_L, rec_U, data.matrix(actuals))
IS_rec <- matrix(IS_rec, nrow = N)

## -----------------------------------------------------------------------------
SS_AE <- (AE_base - AE_rec) / (AE_base + AE_rec) * 2
SS_SE <- (SE_base - SE_rec) / (SE_base + SE_rec) * 2
SS_IS <- (IS_base - IS_rec) / (IS_base + IS_rec) * 2
SS_AE[is.na(SS_AE)] <- 0
SS_SE[is.na(SS_SE)] <- 0
SS_IS[is.na(SS_IS)] <- 0

mean_skill_scores <- c(round(colMeans(SS_IS), 2),
                       round(colMeans(SS_SE), 2),
                       round(colMeans(SS_AE), 2))
mean_skill_scores <- data.frame(t(matrix(mean_skill_scores, nrow = n)))
colnames(mean_skill_scores) <- names(actuals)
rownames(mean_skill_scores) <- c("Interval score", "Squared error", "Absolute error")
knitr::kable(mean_skill_scores)

## -----------------------------------------------------------------------------
# Now we draw a larger number of samples:
N_samples <- 1e5

### Example of concordant-shift effect
j <- 124
base_fc_j <- c()   
for (i in 1:n) base_fc_j[[i]] <- list(size = extr_mkt_events_basefc$size[[i]],
                                      mu = extr_mkt_events_basefc$mu[[j,i]])
# Reconcile
buis <- reconc_BUIS(A, base_fc_j, "params", "nbinom", 
                    num_samples = N_samples, seed = 42)
samples_y <- buis$reconciled_samples

# The reconciled mean is lower than both the base and bottom-up means:
means <- c(round(extr_mkt_events_basefc$mu[[j,1]], 2),
           round(sum(extr_mkt_events_basefc$mu[j,2:n]), 2),
           round(mean(samples_y[1,]), 2))
col_names <- c("Base upper mean", "Bottom-up upper mean", "Reconciled upper mean")
knitr::kable(matrix(means, nrow=1), col.names = col_names)

### Example of combination effect 
j <- 1700
base_fc_j <- c()   
for (i in 1:n) base_fc_j[[i]] <- list(size = extr_mkt_events_basefc$size[[i]],
                                      mu = extr_mkt_events_basefc$mu[[j,i]])
# Reconcile 
buis <- reconc_BUIS(A, base_fc_j, "params", "nbinom", 
                    num_samples = N_samples, seed = 42)
samples_y <- buis$reconciled_samples

# The reconciled mean is between the base and the bottom-up mean:
means <- c(round(extr_mkt_events_basefc$mu[[j,1]], 2),
           round(sum(extr_mkt_events_basefc$mu[j,2:n]), 2),
           round(mean(samples_y[1,]), 2))
col_names <- c("Base upper mean", "Bottom-up upper mean", "Reconciled upper mean")
knitr::kable(matrix(means, nrow=1), col.names = col_names)

## -----------------------------------------------------------------------------
j <- 2308
base_fc_j <- c()   
for (i in 1:n) base_fc_j[[i]] <- list(size = extr_mkt_events_basefc$size[[i]],
                                      mu = extr_mkt_events_basefc$mu[[j,i]])
# Reconcile
buis <- reconc_BUIS(A, base_fc_j, "params", "nbinom", num_samples = N_samples, seed = 42)
samples_y <- buis$reconciled_samples

# Compute the variance of the base and reconciled bottom forecasts
base_bottom_var <- mapply(function(mu, size) var(rnbinom(n = 1e5, size = size, mu = mu)),
                          extr_mkt_events_basefc$mu[j,2:n], 
                          extr_mkt_events_basefc$size[2:n])
rec_bottom_var <- apply(samples_y[2:n,], MARGIN = 1, var) 

# The reconciled variance is greater than the base variance:
bottom_var <- rbind(base_bottom_var, rec_bottom_var)
rownames(bottom_var) <- c("var base", "var reconc")
knitr::kable(round(bottom_var, 2))