`bayesRecon`

This vignette shows how to perform *probabilistic
reconciliation* with the `bayesRecon`

package. We provide
three examples:

*Temporal hierarchy for a count time series*: we build a temporal hierarchy over a count time series, produce the base forecasts using`glarma`

and reconcile them via Bottom-Up Importance Sampling (BUIS).*Temporal hierarchy for a smooth time series*: we build a temporal hierarchy over a smooth time series, compute the base forecasts using`ets`

and we reconcile them in closed form using Gaussian reconciliation. The covariance matrix is diagonal.*Hierarchical of smooth time series*: this is an example of a cross-sectional hierarchy. We generate the base forecasts using`ets`

and we reconcile them via Gaussian reconciliation. The covariance matrix is full and estimated via shrinkage.

The package, available at this CRAN page, can be installed and loaded with the usual commands:

Load the package:

We select a monthly time series of counts from the *carparts*
dataset, available from the expsmooth package (R.
J. Hyndman 2015). The data set contains time series of sales of
cars part from Jan. 1998 to Mar. 2002. For this example we select time
series #2655, which we make available as
`bayesRecon::carparts_example`

.

This time series has a skewed distribution of values.

```
layout(mat = matrix(c(1, 2), nrow = 1, ncol = 2), widths = c(2, 1))
plot(carparts_example, xlab = "Time", ylab = "Car part sales", main = NULL)
hist(carparts_example, xlab = "Car part sales", main = NULL)
```

We divide the time series into train and test; the test set
contains the last 12 months.

```
train <- window(carparts_example, end = c(2001, 3))
test <- window(carparts_example, start = c(2001, 4))
```

We build the temporal hierarchy using the
`temporal aggregation`

function. We specify the aggregation
levels using the `agg_levels`

argument; in this case they are
*2-Monthly*, *Quarterly*, *4-Monthly*,
*Biannual*, and *Annual*.

```
train.agg <- bayesRecon::temporal_aggregation(train, agg_levels = c(2, 3, 4, 6, 12))
levels <- c("Annual", "Biannual", "4-Monthly", "Quarterly", "2-Monthly", "Monthly")
names(train.agg) <- levels
```

The function returns a list of aggregated time series, ordered from the most aggregated (top of the hierarchy) to the most disagreggated (bottom of the hierarchy). We plot them below.

```
par(mfrow = c(2, 3), mai = c(0.6, 0.6, 0.5, 0.5))
for (l in levels) {
plot(train.agg[[l]], xlab = "Time", ylab = "Car part sales", main = l)
}
```

We compute the *base forecasts* using the package `glarma`

,
a package specific for forecasting count time series. We forecast 12
steps ahead at the monthly level, 4 steps ahead at the quarterly level,
etc. by iterating over the levels of the hierarchy, At each level, we
fit a `glarma`

model with Poisson predictive distribution and
we forecast; each forecast is constituted by 20000 integer samples.

Eventually we collect the samples of the 28 predictive distributions
(one at the *Annual* level, two at the *Biannual* level,
etc.) in a list. The code below takes about 30 seconds on a standard
computer.

```
# install.packages("glarma", dependencies = TRUE)
#library(glarma)
fc.samples <- list()
D <- 20000
fc.count <- 1
# iterating over the temporal aggregation levels
for (l in seq_along(train.agg)) {
f.level <- frequency(train.agg[[l]])
print(paste("Forecasting at ", levels[l], "...", sep = ""))
# fit an independent model for each aggregation level
model <- glarma::glarma(
train.agg[[l]],
phiLags = if (f.level == 1) 1 else 1:(min(6, f.level - 1)),
thetaLags = if (f.level == 1) NULL else f.level,
X = cbind(intercept = rep(1, length(train.agg[[l]]))),
offset = cbind(intercept = rep(0, length(train.agg[[l]]))),
type = "Poi"
)
# forecast 1 year ahead
h <- f.level
tmp <- matrix(data = NA, nrow = h, ncol = D)
for (s in (1:D)) {
# each call to 'forecast.glarma' returns a simulation path
tmp[, s] <- glarma::forecast(
model,
n.ahead = h,
newdata = cbind(intercept = rep(1, h)),
newoffset = rep(0, h)
)$Y
}
# collect the forecasted samples
for (i in 1:h) {
fc.samples[[fc.count]] <- tmp[i, ]
fc.count <- fc.count + 1
}
}
#> [1] "Forecasting at Annual..."
#> [1] "Forecasting at Biannual..."
#> [1] "Forecasting at 4-Monthly..."
#> [1] "Forecasting at Quarterly..."
#> [1] "Forecasting at 2-Monthly..."
#> [1] "Forecasting at Monthly..."
```

Reconciliation requires the aggregation matrix \(\mathbf{A}\), which we obtain using the
function `get_reconc_matrices`

. It requires:

- the aggregation factors of the hierarchy, which in this example are \(\{2, 3, 4, 6, 12\}\);
- the length of the forecasting horizon at the bottom level, which is 12 in this example.

```
recon.matrices <- bayesRecon::get_reconc_matrices(agg_levels = c(2, 3, 4, 6, 12), h = 12)
# Aggregation matrix
A <- recon.matrices$A
```

To reconcile using Bottom-Up Important Sampling (BUIS) we we use the
function `reconc_BUIS`

, passing to it the \(\mathbf{A}\) matrix, the *base
forecasts*, the type of the base forecasts
(`in_type`

=“samples”) and whether the samples are discrete or
integer (`distr`

= “discrete”).

```
recon.res <- bayesRecon::reconc_BUIS(
A,
base_forecasts = fc.samples,
in_type = "samples",
distr = "discrete",
seed = 42
)
```

Here we obtain samples from the reconciled forecast distribution.

We now compute the Mean Absolute Error (MAE) and the Continuous
Ranked Probability Score (CRPS) for the bottom (i.e., *monthly*)
time series. For computing CRPS, we use the package `scoringRules`

.

```
# install.packages("scoringRules", dependencies = TRUE)
library(scoringRules)
ae.fc <- list()
ae.reconc <- list()
crps.fc <- list()
crps.reconc <- list()
for (h in 1:length(test)) {
y.hat_ <- median(fc.samples[[nrow(A) + h]])
y.reconc_ <- median(recon.res$bottom_reconciled_samples[, h])
# Compute Absolute Errors
ae.fc[[h]] <- abs(test[h] - y.hat_)
ae.reconc[[h]] <- abs(test[h] - y.reconc_)
# Compute Continuous Ranked Probability Score (CRPS)
crps.fc[[h]] <-
scoringRules::crps_sample(y = test[h], dat = fc.samples[[nrow(A) + h]])
crps.reconc[[h]] <-
scoringRules::crps_sample(y = test[h], dat = recon.res$bottom_reconciled_samples[, h])
}
mae.fc <- mean(unlist(ae.fc))
mae.reconc <- mean(unlist(ae.reconc))
crps.fc <- mean(unlist(crps.fc))
crps.reconc <- mean(unlist(crps.reconc))
metrics <- data.frame(
row.names = c("MAE", "CRPS"),
base.forecasts = c(mae.fc, crps.fc),
reconciled.forecasts = c(mae.reconc, crps.reconc)
)
metrics
#> base.forecasts reconciled.forecasts
#> MAE 1.2500000 1.0416667
#> CRPS 0.7093284 0.6840278
```

In this second example, we select a smooth monthly time series
(N1485) from the M3 forecasting competition (Makridakis and Hibon 2000). The data set is
available in the Mcomp package (R. Hyndman
2018) and we make available this time series as
`bayesRecon::m3_example`

.

We build the temporal hierarchy using the
`temporal_aggregation`

function.

Without specifying `agg_levels`

, the function generates by
default all the feasible aggregation: 2-Monthly, Quarterly, 4-Monthly,
Biannual, and Annual.

```
train.agg <- bayesRecon::temporal_aggregation(M3_example$train)
levels <- c("Annual", "Biannual", "4-Monthly", "Quarterly", "2-Monthly", "Monthly")
names(train.agg) <- levels
```

We generate the base forecasts using `ets`

from the forecast package
(R. J. Hyndman and Khandakar 2008). At
every level we predict 18 months ahead. All the predictive distributions
are Gaussian.

```
# install.packages("forecast", dependencies = TRUE)
library(forecast)
#> Registered S3 method overwritten by 'quantmod':
#> method from
#> as.zoo.data.frame zoo
H <- length(M3_example$test)
H
#> [1] 18
fc <- list()
level.idx <- 1
fc.idx <- 1
for (level in train.agg) {
level.name <- names(train.agg)[level.idx]
# fit an ETS model for each temporal level
model <- ets(level)
# generate forecasts for each level within 18 months
h <- floor(H / (12 / frequency(level)))
print(paste("Forecasting at ", level.name, ", h=", h, "...", sep = ""))
level.fc <- forecast(model, h = h)
# save mean and sd of the gaussian predictive distribution
for (i in 1:h) {
fc[[fc.idx]] <- list(mean = level.fc$mean[[i]],
sd = (level.fc$upper[, "95%"][[i]] - level.fc$mean[[i]]) / qnorm(0.975))
fc.idx <- fc.idx + 1
}
level.idx <- level.idx + 1
}
#> [1] "Forecasting at Annual, h=1..."
#> [1] "Forecasting at Biannual, h=3..."
#> [1] "Forecasting at 4-Monthly, h=4..."
#> [1] "Forecasting at Quarterly, h=6..."
#> [1] "Forecasting at 2-Monthly, h=9..."
#> [1] "Forecasting at Monthly, h=18..."
```

Using the function `get_reconc_matrices`

, we get matrix
\(\mathbf{A}\).

```
rmat <- get_reconc_matrices(agg_levels = c(2, 3, 4, 6, 12), h = 18)
par(mai = c(1,1,0.5,0.5))
image(1:ncol(rmat$A), 1:nrow(rmat$A),
t(apply(t(rmat$A),1,rev)),
xaxt='n', yaxt='n', ylab = "", xlab = levels[6])
axis(1, at=1:ncol(rmat$A), label=1:ncol(rmat$A), las=1)
axis(2, at=c(23,22,19,15,9), label=levels[1:5], las=2)
```

The function `reconc_gaussian`

implements the exact
Gaussian reconciliation. We also run `reconc_BUIS`

, to check
the consistency between the two approaches.

```
recon.gauss <- bayesRecon::reconc_gaussian(
A = rmat$A,
base_forecasts.mu = sapply(fc, "[[", 1),
base_forecasts.Sigma = diag(sapply(fc, "[[", 2)) ^ 2
)
reconc.buis <- bayesRecon::reconc_BUIS(
A = rmat$A,
base_forecasts = fc,
in_type = "params",
distr = "gaussian",
num_samples = 20000,
seed = 42
)
# check that the algorithms return consistent results
round(rbind(
c(rmat$S %*% recon.gauss$bottom_reconciled_mean),
rowMeans(reconc.buis$reconciled_samples)
))
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
#> [1,] 74977 35913 39063 41491 23520 25136 26321 27174 17464 18450 19251 19812
#> [2,] 74946 35897 39049 41470 23532 25091 26323 27134 17462 18435 19231 19818
#> [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
#> [1,] 20412 21079 11527 11993 12393 12743 13047 13274 13531 13643 14317 5694
#> [2,] 20425 21046 11547 11985 12365 12726 13041 13282 13482 13652 14336 5694
#> [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
#> [1,] 5833 5936 6056 6138 6255 6324 6419 6508 6538 6619 6655 6759
#> [2,] 5853 5914 6070 6127 6237 6308 6418 6505 6536 6612 6670 6729
#> [,37] [,38] [,39] [,40] [,41]
#> [1,] 6772 6881 6762 7160 7157
#> [2,] 6753 6943 6709 7152 7184
```

We now compare *base forecasts* and *reconciled
forecasts*:

```
yhat.mu <- tail(sapply(fc, "[[", 1), 18)
yhat.sigma <- tail(sapply(fc, "[[", 2), 18)
yhat.hi95 <- qnorm(0.975, mean = yhat.mu, sd = yhat.sigma)
yhat.lo95 <- qnorm(0.025, mean = yhat.mu, sd = yhat.sigma)
yreconc.mu <- rowMeans(reconc.buis$bottom_reconciled_samples)
yreconc.hi95 <- apply(reconc.buis$bottom_reconciled_samples, 1,
function(x) quantile(x, 0.975))
yreconc.lo95 <- apply(reconc.buis$bottom_reconciled_samples, 1,
function(x) quantile(x, 0.025))
ylim_ <- c(min(M3_example$train, M3_example$test, yhat.lo95, yreconc.lo95) - 1,
max(M3_example$train, M3_example$test, yhat.hi95, yreconc.hi95) + 1)
plot(M3_example$train, xlim = c(1990, 1995.6), ylim = ylim_,
ylab = "y", main = "N1485 Forecasts")
lines(M3_example$test, lty = "dashed")
lines(ts(yhat.mu, start = start(M3_example$test), frequency = 12),
col = "coral", lwd = 2)
lines(ts(yreconc.mu, start = start(M3_example$test), frequency = 12),
col = "blue2", lwd = 2)
xtest <- time(M3_example$test)
polygon(c(xtest, rev(xtest)), c(yhat.mu, rev(yhat.hi95)),
col = "#FF7F5066", border = "#FF7F5066")
polygon(c(xtest, rev(xtest)), c(yhat.mu, rev(yhat.lo95)),
col = "#FF7F5066", border = "#FF7F5066")
polygon(c(xtest, rev(xtest)), c(yreconc.mu, rev(yreconc.hi95)),
col = "#0000EE4D", border = "#0000EE4D")
polygon(c(xtest, rev(xtest)), c(yreconc.mu, rev(yreconc.lo95)),
col = "#0000EE4D", border = "#0000EE4D")
```

In this example, we consider the hierarchical time series
*infantgts*, which is available from the `hts`

package
(R. Hyndman et al. 2021) We make it
available also in our package as
`bayesRecon::infantMortality`

.

It contains counts of infant mortality (deaths) in Australia, disaggregated by state and sex (male and female).

We forecast one year ahead using `auto.arima`

from the `forecast`

package. We collect the residuals, which we will later use to compute
the covariance matrix.

```
# install.packages("forecast", dependencies = TRUE)
library(forecast)
fc <- list()
residuals <- matrix(NA,
nrow = length(infantMortality$total),
ncol = length(infantMortality))
fc.idx <- 1
for (s in infantMortality) {
s.name <- names(infantMortality)[fc.idx]
print(paste("Forecasting at ", s.name, "...", sep = ""))
# fit an auto.arima model and forecast with h=1
model <- auto.arima(s)
s.fc <- forecast(model, h = 1)
# save mean and sd of the gaussian predictive distribution
fc[[s.name]] <- c(s.fc$mean,
(s.fc$upper[, "95%"][[1]] - s.fc$mean) / qnorm(0.975))
residuals[, fc.idx] <- s.fc$residuals
fc.idx <- fc.idx + 1
}
#> [1] "Forecasting at total..."
#> [1] "Forecasting at NSW..."
#> [1] "Forecasting at VIC..."
#> [1] "Forecasting at QLD..."
#> [1] "Forecasting at SA..."
#> [1] "Forecasting at WA..."
#> [1] "Forecasting at NT..."
#> [1] "Forecasting at ACT..."
#> [1] "Forecasting at TAS..."
#> [1] "Forecasting at male..."
#> [1] "Forecasting at female..."
#> [1] "Forecasting at NSW male..."
#> [1] "Forecasting at NSW female..."
#> [1] "Forecasting at VIC male..."
#> [1] "Forecasting at VIC female..."
#> [1] "Forecasting at QLD male..."
#> [1] "Forecasting at QLD female..."
#> [1] "Forecasting at SA male..."
#> [1] "Forecasting at SA female..."
#> [1] "Forecasting at WA male..."
#> [1] "Forecasting at WA female..."
#> [1] "Forecasting at NT male..."
#> [1] "Forecasting at NT female..."
#> [1] "Forecasting at ACT male..."
#> [1] "Forecasting at ACT female..."
#> [1] "Forecasting at TAS male..."
#> [1] "Forecasting at TAS female..."
```

Now we build the \(\mathbf{A}\) matrix.

```
# we have 16 bottom time series, and 11 upper time series
A <- matrix(data = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,
1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,
0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1), byrow=TRUE, ncol = 16)
# plot of A
par(mai = c(1.5,1,0.5,0.5))
image(1:ncol(A), 1:nrow(A),
t(apply(t(A),1,rev)),
xaxt='n', yaxt='n', ann=FALSE)
axis(1, at=1:ncol(A), label=names(infantMortality)[12:27], las=2)
axis(2, at=c(1:11), label=rev(names(infantMortality)[1:11]), las=2)
```