Reconcile Forecasts

This tutorial explains how use the ReconcileAdditiveForecasts class to create forecasts that satisfy inter-forecast additivity constraints.

The inputs are:

1. additive constraints to be satisfied

2. original (base) forecasts (timeseries)

3. actuals (timeseries)

The output is adjusted forecasts that satisfy the constraints.

Optimization Approach

The adjusted forecasts are computed as a linear transformation of the base forecasts. The linear transform is the solution to an optimization problem (Reconcile Forecasts).

In brief, the objective is to minimize the weighted sum of these error terms:

1. Training MSE: empirical MSE of the adjusted forecasts on the training set

2. Bias penalty: estimated squared bias of adjusted forecast errors

3. Variance penalty: estimated variance of adjusted forecast errors for an unbiased transformation, assuming base forecasts are unbiased (this underestimates the variance if the transformation is biased).

4. Adjustment penalty: regularization term that penalizes large adjustments

Subject to these constraints:

1. Adjusted forecasts satisfy inter-forecast additivity constraints (required)

2. Transform is unbiased (optional)

3. Transform matrix entries are between [lower, upper] bound (optional)

ReconcileAdditiveForecasts allows you to tune the optimization objective and constraints. It also exposes common methods as special cases of this optimization problem. The available methods are:

• "bottom_up" (bottom up)

• "ols" (OLS)

• "mint_sample" (MinT with sample covariance)

• "custom" (custom objective and constraints)

Note

"bottom_up" is applicable when the constraints can be represented as a tree. It produces reconciled forecasts by summing the leaf nodes. This is equivalent to the solution to the optimization that only penalizes adjustment to the leaf nodes’ forecasts.

"ols" and "mint_sample" include only the variance penalty and require that the transform be unbiased. The variance penalty depends on forecast error covariances. "ols" assumes base forecast errors are uncorrelated with equal variance. "mint_sample" uses sample covariance of the forecast errors.

Prepare Input Data

In this tutorial, we consider a 3-level tree with the parent-child relationships below.

        00        # level 0
      /   \
   10       11    # level 1
  / | \     /\
20 21 22   23 24  # level 2

We want the forecasts of parent nodes to equal the sum of the forecasts of their children.

First, we need to generate forecasts for each of the nodes. One approach is to generate the forecasts independently, using rolling window forecasting to get h-step ahead forecasts over time, for some constant h. This can be done with the benchmark class. (The variance penalty assumes the residuals have fixed covariance, and using constant h helps with that assumption.)

For this tutorial, we assume that forecasts have already been computed. Below, forecasts and actuals are pandas DataFrames in long format, where each column is a time series, and each row is a time step. The rows are sorted in ascending order.

import logging
import plotly
import warnings

import pandas as pd
import numpy as np

from greykite.algo.reconcile.convex.reconcile_forecasts import ReconcileAdditiveForecasts
from greykite.common.constants import TIME_COL
from greykite.common.data_loader import DataLoader
from greykite.common.viz.timeseries_plotting import plot_multivariate

logger = logging.getLogger()
logger.setLevel(logging.ERROR)  # reduces logging
warnings.simplefilter("ignore", category=UserWarning)  # ignores matplotlib warnings when rendering documentation

dl = DataLoader()
actuals = dl.load_data(data_name="daily_hierarchical_actuals")
forecasts = dl.load_data(data_name="daily_hierarchical_forecasts")
actuals.set_index(TIME_COL, inplace=True)
forecasts.set_index(TIME_COL, inplace=True)
forecasts.head().round(1)

	00	10	11	20	21	22	23	24
ts
2019-08-05	477.5	247.9	243.0	94.7	73.3	80.5	146.1	95.8
2019-08-06	473.0	250.4	233.3	76.2	88.1	86.6	130.8	102.6
2019-08-07	493.6	250.0	256.3	72.4	91.9	86.0	151.3	103.8
2019-08-08	500.1	248.7	265.1	75.1	87.8	84.8	163.3	101.1
2019-08-09	506.0	246.1	283.2	88.5	75.1	79.8	184.1	95.3

Note

To use the reconcile method, dataframe columns should contain only the forecasts or actuals timeseries. Time should not be its own column.

Above, we set time as the index using .set_index(). Index values are ignored by the reconcile method so you could also choose to drop the column.

The rows and columns in forecasts and actuals correspond to each other.

assert forecasts.index.equals(actuals.index)
assert forecasts.columns.equals(actuals.columns)

Next, we need to encode the constraints. In general, these can be defined by constraint_matrix. This is a c x m array encoding c constraints in m variables, where m is the number of timeseries. The columns in this matrix correspond to the columns in the forecasts/actuals dataframes below. The rows encode additive expressions that must equal 0.

constraint_matrix = np.array([
   # 00  10  11 20 21 22 23 24
    [-1,  1,  1, 0, 0, 0, 0, 0],  # 0 = -1*x_00 + 1*x_10 + 1*x_11
    [ 0, -1,  0, 1, 1, 1, 0, 0],  # 0 = -1*x_10 + 1*x_20 + 1*x_21 + 1*x_22
    [ 0,  0, -1, 0, 0, 0, 1, 1]   # 0 = -1*x_11 + 1*x_23 + 1*x_24
])

Alternatively, if the graph is a tree, you can use the levels parameter. This is a more concise way to specify additive tree # constraints, where forecasts of parent nodes must equal the sum of the forecasts of their children. It assumes the columns in forecasts and actuals are in the tree’s breadth first traversal order: i.e., starting from the root, scan left to right, top to bottom, as shown below for our example:

       0
    /     \
   1       2
 / | \    / \
3  4  5  6   7

Here is an equivalent specification using the levels parameter.

# The root has two children.
# Its children have 3 and 2 children, respectively.
levels = [[2], [3, 2]]
# Summarize non-leaf nodes by the number of children
# they have, and iterate in breadth first traversal order.
# Each level in the tree becomes a sublist of `levels`.
#
#          (2)     --> [2]
#         /   \
#      (3)     (2) --> [3, 2]

Note

More formally, levels specifies the number of children of each internal node in the tree. The ith inner list provides the number of children of each node in level i. Thus, the first sublist has one integer, the length of a sublist is the sum of the previous sublist, and all entries in levels are positive integers. All leaf nodes must have the same depth.

For illustration, we plot the inconsistency between forecasts of the root node, "00", and its children. Notice that the blue and orange lines do not perfectly overlap.

parent = "00"
children = ["10", "11"]
cols = {
    f"parent-{parent}": forecasts[parent],
    "sum(children)": sum(forecasts[child] for child in children)
}
cols.update({f"child-{child}": forecasts[child] for child in children})
cols[TIME_COL] = forecasts.index
parent_child_df = pd.DataFrame(cols)
fig = plot_multivariate(
    df=parent_child_df,
    x_col=TIME_COL,
    title=f"Forecasts of node '{parent}' and its children violate the constraint",
)
plotly.io.show(fig)

Forecast reconciliation

Training Evaluation

To reconcile these forecasts, we use the ReconcileAdditiveForecasts class.

raf = ReconcileAdditiveForecasts()

Fit

Call fit() to learn the linear transform. Available methods are "bottom_up", "ols", "mint_sample", "custom". Let’s start with the bottom up method.

_ = raf.fit(
    forecasts=forecasts,
    actuals=actuals,
    levels=levels,
    method="bottom_up",
)

Each row in the transform matrix shows how to compute the adjusted forecast as a linear combination of the base forecasts. For the “bottom up” transform, the matrix simply reflects the tree structure.

raf.transform_matrix

Out:

array([[0., 0., 0., 1., 1., 1., 1., 1.],
       [0., 0., 0., 1., 1., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 1.],
       [0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1.]])

We can visualize this matrix to more easily see how forecasts are combined. The top row in this plot shows that the adjusted forecast for node “00” (tree root) is the sum of all the base forecasts of the leaf nodes. “10” and “11” are the sum of their children, and each leaf node keeps its original value.

fig = raf.plot_transform_matrix()
plotly.io.show(fig)

Transform

The transform() method applies the transform and returns the adjusted (consistent) forecasts. If we call it without arguments, it applies the transform to the training set.

adjusted_forecasts = raf.transform()
adjusted_forecasts.head().round(1)

	00	10	11	20	21	22	23	24
ts
2019-08-05	490.3	248.5	241.9	94.7	73.3	80.5	146.1	95.8
2019-08-06	484.3	250.9	233.3	76.2	88.1	86.6	130.8	102.6
2019-08-07	505.5	250.4	255.1	72.4	91.9	86.0	151.3	103.8
2019-08-08	512.1	247.7	264.4	75.1	87.8	84.8	163.3	101.1
2019-08-09	522.7	243.4	279.3	88.5	75.1	79.8	184.1	95.3

The adjusted forecasts on the training set are stored in the adjusted_forecasts attribute.

assert adjusted_forecasts.equals(raf.adjusted_forecasts)

Evaluate

Now that we have the actuals, forecasts, and adjusted forecasts, we can check how the adjustment affects forecast quality. Here, we do evaluation on the training set.

_ = raf.evaluate(
    is_train=True,         # evaluates on training set
    ipython_display=True,  # displays evaluation table
    plot=True,             # displays plots
    plot_num_cols=2,       # formats plots into two columns
)

Out:

    MAPE pp change  MedAPE pp change  RMSE % change  Base MAPE  Adjusted MAPE  Base MedAPE  Adjusted MedAPE  Base RMSE  Adjusted RMSE
00             0.4               0.2           15.4        2.9            3.3          2.6              2.8       17.8           20.5
10            -0.1              -0.0           -4.3        1.2            1.1          1.1              1.1        3.7            3.5
11             0.6               0.7            7.4        6.7            7.3          5.3              6.0       20.7           22.2
20             0.0               0.0            0.0        8.3            8.3          7.4              7.4        7.2            7.2
21             0.0               0.0            0.0        5.7            5.7          4.6              4.6        5.0            5.0
22             0.0               0.0            0.0        3.5            3.5          2.5              2.5        4.5            4.5
23             0.0               0.0            0.0       19.9           19.9         13.4             13.4       23.5           23.5
24             0.0               0.0            0.0        1.0            1.0          0.8              0.8        1.7            1.7

For better formatting in this documentation, let’s display the table again. evaluation_df contains the evaluation table for the training set. The errors for the leaf nodes are the same, as expected, because their forecasts have not changed. The error for nodes “00” and “11” have increased.

raf.evaluation_df.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	0.4	0.2	15.4	2.9	3.3	2.6	2.8	17.8	20.5
10	-0.1	-0.0	-4.3	1.2	1.1	1.1	1.1	3.7	3.5
11	0.6	0.7	7.4	6.7	7.3	5.3	6.0	20.7	22.2
20	0.0	0.0	0.0	8.3	8.3	7.4	7.4	7.2	7.2
21	0.0	0.0	0.0	5.7	5.7	4.6	4.6	5.0	5.0
22	0.0	0.0	0.0	3.5	3.5	2.5	2.5	4.5	4.5
23	0.0	0.0	0.0	19.9	19.9	13.4	13.4	23.5	23.5
24	0.0	0.0	0.0	1.0	1.0	0.8	0.8	1.7	1.7

Note

The ipython_display parameter controls whether to display the evaluation table.

• The “*change” columns show the change in error after adjustment.

• The “Base*” columns show evaluation metrics for the original base forecasts.

• The “Adjusted*” columns show evaluation metrics for the adjusted forecasts.

• MAPE/MedAPE = mean/median absolute percentage error, RMSE = root mean squared error, pp = percentage point.

We can check the diagnostic plots for more information. The “Base vs Adjusted” and “Adjustment Size” plots show that the forecasts for “00” and “11” are higher after adjustment. The “Forecast Error” plot shows that this increased the forecast error. (Plots are automatically shown when plot=True. To make plots appear inline in this tutorial, we need to explicitly show the figures.)

plotly.io.show(raf.figures["base_adj"])

plotly.io.show(raf.figures["adj_size"])

plotly.io.show(raf.figures["error"])

Note

The plot parameter controls whether to display diagnostic plots to adjusted to base forecasts.

• “Base vs Adjusted Forecast” shows base forecast (blue) vs adjusted forecast (orange)

• “Adjustment Size (%)” shows the size of the adjustment.

• “Forecast Error (%)” shows the % error before (blue) and after (orange) adjustment. Closer to 0 is better.

• Note that the y-axes are independent.

For completeness, we can verify that the actuals and adjusted forecasts satisfy the constraints. constraint_violation shows constraint violation on the training set, defined as root mean squared violation (averaged across time points and constraints), divided by root mean squared actual value. It should be close to 0 for “adjusted” and “actual”. (This is not necessary to check, because a warning is printed during fitting if actuals do not satisfy the constraints or if there is no solution to the optimization problem.)

raf.constraint_violation

Out:

{'adjusted': 1.2097724692944673e-16, 'forecast': 0.03741571732039168, 'actual': 2.0074667673074238e-16}

Test Set Evaluation

Evaluation on the training set is sufficient for the "bottom_up" and "ols" methods, because they do not use the forecasts or actuals to learn the transform matrix. The transform depends only on the constraints.

The "mint_sample" and "custom" methods use forecasts and actuals in addition to the constraints, so we should evaluate accuracy on an out-of-sample test set.

Information used by each method

	constraints	forecasts	actuals
bottom_up	X
ols	X
mint_sample	X	X	X
custom	X	X	X

"custom" always uses the constraints. Whether it uses forecasts and actuals depends on the optimization terms:

• forecasts: used for adjustment penalty, train penalty, variance penalty with “sample” covariance, preset weight options (“MedAPE”, “InverseMedAPE”).

• actuals: used for bias penalty, train penalty, variance penalty with “sample” covariance, preset weight options (“MedAPE”, “InverseMedAPE”).

Train

We’ll fit to the first half of the data and evaluate accuracy on the second half.

train_size = forecasts.shape[0]//2
forecasts_train = forecasts.iloc[:train_size,:]
actuals_train = actuals.iloc[:train_size,:]
forecasts_test = forecasts.iloc[train_size:,:]
actuals_test = actuals.iloc[train_size:,:]

Let’s try the "mint_sample" method. First, fit the transform and apply it on the training set. The transform matrix is more complex than before.

raf = ReconcileAdditiveForecasts()
raf.fit_transform(  # fits and transforms the training data
    forecasts=forecasts_train,
    actuals=actuals_train,
    levels=levels,
    method="mint_sample"
)
assert raf.transform_matrix is not None    # train fit result, set by fit
assert raf.adjusted_forecasts is not None  # train transform result, set by transform
fig = raf.plot_transform_matrix()
plotly.io.show(fig)

Now, evaluate accuracy on the training set. In our example, all the reconciled forecasts have lower error than the base forecasts on the training set.

raf.evaluate(is_train=True)
assert raf.evaluation_df is not None         # train evaluation result, set by evaluate
assert raf.figures is not None               # train evaluation figures, set by evaluate
assert raf.constraint_violation is not None  # train constraint violation, set by evaluate
raf.evaluation_df.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	-0.0	-0.1	-4.0	2.8	2.8	2.4	2.3	17.0	16.4
10	-0.2	-0.1	-14.2	1.2	1.0	1.0	0.9	3.6	3.1
11	-1.1	-0.2	-17.1	7.3	6.2	5.3	5.2	21.6	17.9
20	-1.2	-1.1	-19.4	7.9	6.8	6.9	5.8	7.9	6.4
21	-0.0	0.1	-0.8	5.7	5.7	5.0	5.1	5.0	5.0
22	-0.7	0.3	-24.8	3.9	3.3	2.8	3.1	4.9	3.7
23	-4.0	-0.6	-21.6	17.0	13.0	12.0	11.5	24.0	18.8
24	-0.5	-0.6	-33.9	1.4	0.9	1.3	0.7	1.9	1.3

Test

Next, apply the transform to the test set and evaluate accuracy. Not all forecasts have improved on the test set. This demonstrates the importance of test set evaluation.

raf.transform_evaluate(  # transform and evaluates on test data
    forecasts_test=forecasts_test,
    actuals_test=actuals_test,
    ipython_display=False,
    plot=False,
)
assert raf.adjusted_forecasts_test is not None    # test transform result, set by transform
assert raf.evaluation_df_test is not None         # test evaluation result, set by evaluate
assert raf.figures_test is not None               # test evaluation figures, set by evaluate
assert raf.constraint_violation_test is not None  # test constraint violation, set by evaluate
raf.evaluation_df_test.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	-0.0	-0.4	1.9	2.9	2.9	2.8	2.4	18.4	18.8
10	-0.2	-0.1	-17.3	1.3	1.0	1.1	1.0	3.8	3.1
11	-0.0	-0.4	2.7	6.1	6.1	5.3	4.9	19.7	20.2
20	-1.1	-2.3	3.5	8.7	7.6	7.6	5.4	6.5	6.7
21	-0.2	-0.3	-1.3	5.6	5.4	4.5	4.2	5.0	5.0
22	0.0	0.0	8.8	3.0	3.1	2.4	2.4	4.1	4.5
23	-3.2	-3.0	-7.9	22.9	19.6	16.0	13.0	22.9	21.1
24	-0.0	-0.1	0.9	0.7	0.7	0.6	0.6	1.4	1.4

Note

The results for the test set are in the corresponding attributes ending with "_test".

As a summary, here are some key attributes containing the results:

transform_matrix :          transform learned from train set
adjusted_forecasts :        adjusted forecasts on train set
adjusted_forecasts_test :   adjusted forecasts on test set
evaluation_df :             evaluation result on train set
evaluation_df_test :        evaluation result on test set
constraint_violation :      normalized constraint violations on train set
constraint_violation_test : normalized constraint violations on test set
figures :                   evaluation plots on train set
figures_test :              evaluation plots on test set

For full attribute details, see ReconcileAdditiveForecasts.

Model Tuning

Now that you understand the basic usage, we’ll introduce some tuning parameters. If you have enough holdout data, you can use the out of sample evaluation to tune the model.

First, try the presets for the method parameter: "bottom_up", "ols", "mint_sample", "custom".

If you’d like to tune further, use the "custom" method to tune the optimization objective and constraints. The tuning parameters and their default values are shown below. See ReconcileAdditiveForecasts for details.

raf = ReconcileAdditiveForecasts()
_ = raf.fit_transform_evaluate(  # fits, transforms, and evaluates on training data
    forecasts=forecasts_train,
    actuals=actuals_train,
    fit_kwargs=dict(  # additional parameters passed to fit()
        levels=levels,
        method="custom",
        # tuning parameters, with their default values for the custom method
        lower_bound=None,     # Lower bound on each entry of ``transform_matrix``.
        upper_bound=None,     # Upper bound on each entry of ``transform_matrix``.
        unbiased=True,        # Whether the resulting transformation must be unbiased.
        lam_adj=1.0,          # Weight for the adjustment penalty (adj forecast - forecast)
        lam_bias=1.0,         # Weight for the bias penalty (adj actual - actual).
        lam_train=1.0,        # Weight for the training MSE penalty (adj forecast - actual)
        lam_var=1.0,          # Weight for the variance penalty (variance of adjusted forecast errors for an unbiased transformation, assuming base forecasts are unbiased)
        covariance="sample",  # Variance-covariance matrix of base forecast errors, used to compute the variance penalty ("sample", "identity" or numpy array)
        weight_adj=None,      # Weight for the adjustment penalty to put a different weight per-timeseries.
        weight_bias=None,     # Weight for the bias penalty to put a different weight per-timeseries.
        weight_train=None,    # Weight for the train MSE penalty to put a different weight per-timeseries.
        weight_var=None,      # Weight for the variance penalty to put a different weight per-timeseries.
    ),
    evaluate_kwargs=dict()  # additional parameters passed to evaluate()
)

Using "custom" with default settings, we find good training set performance overall.

raf.evaluation_df.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	-0.0	0.0	-3.3	2.8	2.8	2.4	2.4	17.0	16.5
10	-0.0	0.1	-1.9	1.2	1.2	1.0	1.1	3.6	3.5
11	-1.1	-0.2	-15.8	7.3	6.2	5.3	5.1	21.6	18.2
20	-1.2	-1.9	-17.7	7.9	6.7	6.9	5.0	7.9	6.5
21	-0.0	0.3	-0.2	5.7	5.7	5.0	5.3	5.0	5.0
22	-0.6	-0.0	-18.8	3.9	3.4	2.8	2.8	4.9	4.0
23	-3.4	-1.4	-19.5	17.0	13.6	12.0	10.7	24.0	19.3
24	-0.1	-0.1	-7.3	1.4	1.3	1.3	1.2	1.9	1.8

Test set performance is also good, except for node “24”.

raf.transform_evaluate(
    forecasts_test=forecasts_test,
    actuals_test=actuals_test,
    ipython_display=False,
    plot=False
)
raf.evaluation_df_test.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	-0.1	-0.4	-1.8	2.9	2.8	2.8	2.4	18.4	18.1
10	-0.3	-0.4	-20.2	1.3	1.0	1.1	0.8	3.8	3.0
11	-0.1	-0.4	1.1	6.1	6.0	5.3	4.9	19.7	19.9
20	-0.9	-2.0	1.1	8.7	7.8	7.6	5.7	6.5	6.6
21	-0.0	0.1	-0.1	5.6	5.6	4.5	4.6	5.0	5.0
22	-0.1	-0.2	1.1	3.0	3.0	2.4	2.2	4.1	4.1
23	-3.1	-3.2	-8.0	22.9	19.7	16.0	12.8	22.9	21.1
24	0.3	0.4	39.8	0.7	1.0	0.6	1.0	1.4	1.9

Notice from the tables that node “24” had the most accurate base forecast of all nodes. Therefore, we don’t want its adjusted forecast to change much. It’s possible that the above transform was overfitting this node.

We can increase the adjustment penalty for node “24” so that its adjusted forecast will be closer to the original one. This should allow us to get good forecasts overall and for node “24” specifically.

# the order of `weights` corresponds to `forecasts.columns`
weight = np.array([1, 1, 1, 1, 1, 1, 1, 5])  # weight is 5x higher for node "24"
raf = ReconcileAdditiveForecasts()
_ = raf.fit_transform_evaluate(
    forecasts=forecasts_train,
    actuals=actuals_train,
    fit_kwargs=dict(
        levels=levels,
        method="custom",
        lower_bound=None,
        upper_bound=None,
        unbiased=True,
        lam_adj=1.0,
        lam_bias=1.0,
        lam_train=1.0,
        lam_var=1.0,
        covariance="sample",
        weight_adj=weight,    # apply the weights to adjustment penalty
        weight_bias=None,
        weight_train=None,
        weight_var=None,
    )
)

Note

The default weight=None puts equal weight on all nodes. Weight can also be "MedAPE" (proportional to MedAPE of base forecasts), "InverseMedAPE" (proportional to 1/MedAPE of base forecasts), or a numpy array that specifies the weight for each node.

Note

When the transform is unbiased (unbiased=True), the bias penalty is zero, so lam_bias and weight_bias have no effect.

The training error looks good.

raf.evaluation_df.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	-0.0	-0.1	-4.0	2.8	2.8	2.4	2.3	17.0	16.4
10	-0.2	-0.1	-13.2	1.2	1.0	1.0	0.9	3.6	3.1
11	-1.1	-0.0	-17.0	7.3	6.2	5.3	5.3	21.6	17.9
20	-1.2	-1.2	-19.2	7.9	6.7	6.9	5.7	7.9	6.4
21	-0.0	0.2	-0.8	5.7	5.7	5.0	5.2	5.0	5.0
22	-0.7	0.2	-24.2	3.9	3.2	2.8	3.0	4.9	3.7
23	-3.9	-0.7	-21.2	17.0	13.1	12.0	11.4	24.0	18.9
24	-0.2	-0.2	-14.3	1.4	1.2	1.3	1.1	1.9	1.6

Plots of the transform matrix and adjustment size show that node “24“‘s adjusted forecast is almost the same as its base forecast.

fig = raf.plot_transform_matrix()
plotly.io.show(fig)

plotly.io.show(raf.figures["adj_size"])

The test error looks better than before.

# Transform and evaluate on the test set.
raf.transform_evaluate(
    forecasts_test=forecasts_test,
    actuals_test=actuals_test,
    ipython_display=False,
    plot=True,
    plot_num_cols=2,
)
raf.evaluation_df_test.round(1)

	MAPE pp change	MedAPE pp change	RMSE % change	Base MAPE	Adjusted MAPE	Base MedAPE	Adjusted MedAPE	Base RMSE	Adjusted RMSE
00	-0.1	-0.3	0.5	2.9	2.9	2.8	2.5	18.4	18.5
10	-0.3	-0.3	-20.6	1.3	1.0	1.1	0.8	3.8	3.0
11	-0.1	-0.3	2.2	6.1	6.0	5.3	5.0	19.7	20.1
20	-1.0	-2.3	2.5	8.7	7.7	7.6	5.4	6.5	6.7
21	-0.1	-0.2	-1.0	5.6	5.5	4.5	4.3	5.0	5.0
22	-0.0	-0.2	5.8	3.0	3.0	2.4	2.2	4.1	4.3
23	-3.3	-3.4	-8.1	22.9	19.6	16.0	12.6	22.9	21.1
24	-0.0	-0.0	0.4	0.7	0.7	0.6	0.6	1.4	1.4

plotly.io.show(raf.figures_test["base_adj"])

plotly.io.show(raf.figures_test["adj_size"])

plotly.io.show(raf.figures_test["error"])

Tuning Tips

If you have enough data, you can use cross validation with multiple test sets for a better estimate of test error. You can use test error to select the parameters.

To tune the parameters,

1. Try all four methods.

2. Tune the lambdas and the weights for the custom method.

For example, start with these lambda settings to see which penalties are useful:

lambdas = [
    # lam_adj, lam_bias, lam_train, lam_var
    (0, 0, 0, 1),  # the same as "mint_sample" if other params are set to default values.
    (0, 0, 1, 1),
    (1, 0, 0, 1),
    (1, 0, 1, 1),  # the same as "custom" if other params are set to default values.
    (1, 1, 1, 1),  # try this one with unbiased=False
]

Tips:

• var penalty is usually helpful

• train, adj, bias penalties are sometimes helpful

• You can increase the lambda for penalties that are more helpful.

To try a biased transform, set (unbiased=False, lam_bias>0). Avoid (unbiased=False, lam_bias=0), because that can result in high bias.

Choose weights that fit your needs. For example, you may care about the accuracy of some forecasts more than others.

Setting weight_adj to "InverseMedAPE" is a convenient way to penalize adjustment to base forecasts that are already accurate.

Setting weight_bias, weight_train, or weight_var to "MedAPE" is a convenient way to improve the error on base forecasts that start with high error.

Debugging

Some tips if you need to debug:

1. Make sure the constraints are properly encoded (for the bottom up method, another way is to check the transform matrix).

raf.constraint_matrix

Out:

array([[-1.,  0.,  0.,  1.,  1.,  1.,  1.,  1.],
       [ 0., -1.,  0.,  1.,  1.,  1.,  0.,  0.],
       [ 0.,  0., -1.,  0.,  0.,  0.,  1.,  1.]])

2. The constraint violation should be 0 for the actuals.

raf.constraint_violation
raf.constraint_violation_test

Out:

{'adjusted': 1.874446823677435e-15, 'forecast': 0.04015962097753542, 'actual': 1.9465566690613992e-16}

3. Check the transform matrix to understand predictions.

fig = raf.plot_transform_matrix()
plotly.io.show(fig)

4. For all methods besides “bottom_up”, check if a solution was found to the optimization problem. If False, then the transform_matrix may be set to a fallback option (bottom up transform, if available). A warning is printed when this happens (“Failed to find a solution. Falling back to bottom-up method.”).

raf.is_optimization_solution

Out:

True

5. Check prob.status for details about cvxpy solver status and look for printed warnings for any issues. You can pass solver options to the fit method. See ReconcileAdditiveForecasts for details.

raf.prob.status

Out:

'optimal'

6. Inspect the objective function value at the identified solution and its breakdown into components. This shows the terms in the objective after multiplication by the lambdas/weights.

raf.objective_fn_val

Out:

{'adj': 0.0003175260797858461, 'bias': 1.9020021907581692e-30, 'train': 0.0037819550216529265, 'var': 0.0037819550216529143, 'total': 0.007881436123091686}

7. Check objective function weights, to make sure covariance, etc., match expectations.

raf.objective_weights

Out:

{'weight_adj': array([[0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0.5, 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0.5, 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0. , 2.5]]), 'weight_bias': array([[1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1.]]), 'weight_train': array([[1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1.]]), 'weight_var': array([[1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1.]]), 'covariance': array([[ 0.00846011, -0.00054813,  0.00880111,  0.00252624, -0.00150595,
        -0.00141349,  0.00878735, -0.00038106],
       [-0.00054813,  0.00036918, -0.0013657 , -0.00042228,  0.00038015,
         0.0003675 , -0.0016771 ,  0.00015219],
       [ 0.00880111, -0.0013657 ,  0.0135727 ,  0.00320355, -0.00189314,
        -0.00254459,  0.01486313, -0.00096638],
       [ 0.00252624, -0.00042228,  0.00320355,  0.00182758, -0.00092728,
        -0.00099977,  0.00344105, -0.00023999],
       [-0.00150595,  0.00038015, -0.00189314, -0.00092728,  0.00073333,
         0.00052669, -0.00206424,  0.00015938],
       [-0.00141349,  0.0003675 , -0.00254459, -0.00099977,  0.00052669,
         0.00070826, -0.0028976 ,  0.00022745],
       [ 0.00878735, -0.0016771 ,  0.01486313,  0.00344105, -0.00206424,
        -0.0028976 ,  0.0167878 , -0.00114732],
       [-0.00038106,  0.00015219, -0.00096638, -0.00023999,  0.00015938,
         0.00022745, -0.00114732,  0.00010755]])}

8. Check the convex optimization problem.

print(type(raf.prob))
raf.prob

Out:

<class 'cvxpy.problems.problem.Problem'>

Problem(Minimize(Expression(CONVEX, NONNEGATIVE, ())), [Equality(Expression(AFFINE, UNKNOWN, (3, 8)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, UNKNOWN, (8, 5)), Constant(CONSTANT, NONNEGATIVE, (8, 5)))])