library(bayesRecon)
# Consider a simple hierarchy with two bottom and one upper
A <- matrix(c(1, 1), nrow = 1)
# The bottom forecasts are Poisson with lambda=15
lambda <- 15
n_tot <- 60
base_fc_bottom <- list()
base_fc_bottom[[1]] <- apply(matrix(seq(0, n_tot)), MARGIN = 1,
FUN = \(x) dpois(x, lambda = lambda))
base_fc_bottom[[2]] <- apply(matrix(seq(0, n_tot)), MARGIN = 1,
FUN = \(x) dpois(x, lambda = lambda))
# The upper forecast is a Normal with mean 40 and std 5
base_fc_upper <- list(mean = 40, cov = matrix(5^2))
# Reconcile with reconc_MixCond
res.mixCond <- reconc_MixCond(A, base_fc_bottom, base_fc_upper)
# Note that the bottom distributions are slightly shifted to the right
PMF_summary(res.mixCond$bottom_rec_pmf[[1]])
PMF_summary(base_fc_bottom[[1]])
PMF_summary(res.mixCond$bottom_rec_pmf[[2]])
PMF_summary(base_fc_bottom[[2]])
# The upper distribution is slightly shifted to the left
PMF_summary(res.mixCond$upper_rec_pmf[[1]])
PMF_get_var(res.mixCond$upper_rec_pmf[[1]])
library(bayesRecon)
# Consider a simple hierarchy with two bottom and one upper
A <- matrix(c(1, 1), nrow = 1)
# The bottom forecasts are Poisson with lambda=15
lambda <- 15
n_tot <- 60
base_fc_bottom <- list()
base_fc_bottom[[1]] <- apply(matrix(seq(0, n_tot)), MARGIN = 1,
FUN = \(x) dpois(x, lambda = lambda))
base_fc_bottom[[2]] <- apply(matrix(seq(0, n_tot)), MARGIN = 1,
FUN = \(x) dpois(x, lambda = lambda))
# The upper forecast is a Normal with mean 40 and std 5
base_fc_upper <- list(mean = 40, cov = matrix(c(5^2)))
# Reconcile with reconc_TDcond
res.TDcond <- reconc_TDcond(A, base_fc_bottom, base_fc_upper)
# Note that the bottom distributions are shifted to the right
PMF_summary(res.TDcond$bottom_rec_pmf[[1]])
PMF_summary(base_fc_bottom[[1]])
PMF_summary(res.TDcond$bottom_rec_pmf[[2]])
PMF_summary(base_fc_bottom[[2]])
# The upper distribution remains similar
PMF_summary(res.TDcond$upper_rec_pmf[[1]])
PMF_get_var(res.TDcond$upper_rec_pmf[[1]])
## Example 2: reconciliation with unbalanced hierarchy
# We consider the example in Fig. 9 of Zambon et al. (2024).
# The hierarchy has 5 bottoms and 3 uppers
A <- matrix(c(
1, 1, 1, 1, 1,
1, 1, 0, 0, 0,
0, 0, 1, 1, 0
), nrow = 3, byrow = TRUE)
# Note that the 5th bottom only appears in the highest level, this is an unbalanced hierarchy.
n_upper <- nrow(A)
n_bottom <- ncol(A)
# The bottom forecasts are Poisson with lambda=15
lambda <- 15
n_tot <- 60
base_fc_bottom <- list()
for (i in seq(n_bottom)) {
base_fc_bottom[[i]] <- apply(matrix(seq(0, n_tot)), MARGIN = 1,
FUN = \(x) dpois(x, lambda = lambda))
}
# The upper forecasts are a multivariate Gaussian
mean <- c(75, 30, 30)
cov <- matrix(c(
5^2, 5, 5,
5, 10, 0,
5, 0, 10
), nrow = 3, byrow = TRUE)
base_fc_upper <- list(mean = mean, cov = cov)
if (FALSE) {
# If we reconcile with reconc_TDcond it won't work (unbalanced hierarchy)
res.TDcond <- reconc_TDcond(A, base_fc_bottom, base_fc_upper)
}
# We can balance the hierarchy by duplicating the node b5:
# i) consider the time series observations for b5 as the upper u4,
# ii) fit the multivariate ts model for u1, u2, u3, u4.
# In this example we simply assume that the forecast for u1-u4 is
# Gaussian with the mean and variance of u4 given by the parameters in b5.
mean_b5 <- lambda
var_b5 <- lambda
mean <- c(75, 30, 30, mean_b5)
cov <- matrix(c(
5^2, 5, 5, 5,
5, 10, 0, 0,
5, 0, 10, 0,
5, 0, 0, var_b5
), nrow = 4, byrow = TRUE)
base_fc_upper <- list(mean = mean, cov = cov)
# We also need to update the aggregation matrix
A <- matrix(c(
1, 1, 1, 1, 1,
1, 1, 0, 0, 0,
0, 0, 1, 1, 0,
0, 0, 0, 0, 1
), nrow = 4, byrow = TRUE)
# We can now reconcile with TDcond
res.TDcond <- reconc_TDcond(A, base_fc_bottom, base_fc_upper)
# Note that the reconciled distribution of b5 and u4 are identical,
# keep this in mind when using the results of your reconciliation!
max(abs(res.TDcond$bottom_rec_pmf[[5]] - res.TDcond$upper_rec_pmf[[4]]))
Run the code above in your browser using DataLab