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

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

## ----finTShier, fig.cap="**Figure 1**: financial time series hierarchy.", out.width = '80%', 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 <- rbind(buis$upper_rec_samples, buis$bottom_rec_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 <- rbind(buis$upper_rec_samples, buis$bottom_rec_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 <- rbind(buis$upper_rec_samples, buis$bottom_rec_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 <- rbind(buis$upper_rec_samples, buis$bottom_rec_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))