Temporal hierarchy over a count time series
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)
Figure 1: Carpart - monthly car part sales.
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) <- levelsThe 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)
}
Figure 2: The aggregated time series of the temporal hierarchy.
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 summing matrix \(\mathbf{S}\), 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)
# Summing matrix
S <- recon.matrices$S
# A 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{S}\) 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(
S,
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
