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.

 86 import logging
 87 import plotly
 88 import warnings
 89
 90 import pandas as pd
 91 import numpy as np
 92
 93 from greykite.algo.reconcile.convex.reconcile_forecasts import ReconcileAdditiveForecasts
 94 from greykite.common.constants import TIME_COL
 95 from greykite.common.data_loader import DataLoader
 96 from greykite.common.viz.timeseries_plotting import plot_multivariate
 97
 98 logger = logging.getLogger()
 99 logger.setLevel(logging.ERROR)  # reduces logging
100 warnings.simplefilter("ignore", category=UserWarning)  # ignores matplotlib warnings when rendering documentation
101
102 dl = DataLoader()
103 actuals = dl.load_data(data_name="daily_hierarchical_actuals")
104 forecasts = dl.load_data(data_name="daily_hierarchical_forecasts")
105 actuals.set_index(TIME_COL, inplace=True)
106 forecasts.set_index(TIME_COL, inplace=True)
107 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.

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.

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

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.

157 # The root has two children.
158 # Its children have 3 and 2 children, respectively.
159 levels = [[2], [3, 2]]
160 # Summarize non-leaf nodes by the number of children
161 # they have, and iterate in breadth first traversal order.
162 # Each level in the tree becomes a sublist of `levels`.
163 #
164 #          (2)     --> [2]
165 #         /   \
166 #      (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.

182 parent = "00"
183 children = ["10", "11"]
184 cols = {
185     f"parent-{parent}": forecasts[parent],
186     "sum(children)": sum(forecasts[child] for child in children)
187 }
188 cols.update({f"child-{child}": forecasts[child] for child in children})
189 cols[TIME_COL] = forecasts.index
190 parent_child_df = pd.DataFrame(cols)
191 fig = plot_multivariate(
192     df=parent_child_df,
193     x_col=TIME_COL,
194     title=f"Forecasts of node '{parent}' and its children violate the constraint",
195 )
196 plotly.io.show(fig)

Forecast reconciliation

Training Evaluation

To reconcile these forecasts, we use the ReconcileAdditiveForecasts class.

206 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.

215 _ = raf.fit(
216     forecasts=forecasts,
217     actuals=actuals,
218     levels=levels,
219     method="bottom_up",
220 )

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.

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.

233 fig = raf.plot_transform_matrix()
234 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.

241 adjusted_forecasts = raf.transform()
242 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.

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.

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

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.

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.)

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

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.)

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.

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.

365 raf = ReconcileAdditiveForecasts()
366 raf.fit_transform(  # fits and transforms the training data
367     forecasts=forecasts_train,
368     actuals=actuals_train,
369     levels=levels,
370     method="mint_sample"
371 )
372 assert raf.transform_matrix is not None    # train fit result, set by fit
373 assert raf.adjusted_forecasts is not None  # train transform result, set by transform
374 fig = raf.plot_transform_matrix()
375 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.

381 raf.evaluate(is_train=True)
382 assert raf.evaluation_df is not None         # train evaluation result, set by evaluate
383 assert raf.figures is not None               # train evaluation figures, set by evaluate
384 assert raf.constraint_violation is not None  # train constraint violation, set by evaluate
385 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.

393 raf.transform_evaluate(  # transform and evaluates on test data
394     forecasts_test=forecasts_test,
395     actuals_test=actuals_test,
396     ipython_display=False,
397     plot=False,
398 )
399 assert raf.adjusted_forecasts_test is not None    # test transform result, set by transform
400 assert raf.evaluation_df_test is not None         # test evaluation result, set by evaluate
401 assert raf.figures_test is not None               # test evaluation figures, set by evaluate
402 assert raf.constraint_violation_test is not None  # test constraint violation, set by evaluate
403 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.

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

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

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”.

473 raf.transform_evaluate(
474     forecasts_test=forecasts_test,
475     actuals_test=actuals_test,
476     ipython_display=False,
477     plot=False
478 )
479 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.

492 # the order of `weights` corresponds to `forecasts.columns`
493 weight = np.array([1, 1, 1, 1, 1, 1, 1, 5])  # weight is 5x higher for node "24"
494 raf = ReconcileAdditiveForecasts()
495 _ = raf.fit_transform_evaluate(
496     forecasts=forecasts_train,
497     actuals=actuals_train,
498     fit_kwargs=dict(
499         levels=levels,
500         method="custom",
501         lower_bound=None,
502         upper_bound=None,
503         unbiased=True,
504         lam_adj=1.0,
505         lam_bias=1.0,
506         lam_train=1.0,
507         lam_var=1.0,
508         covariance="sample",
509         weight_adj=weight,    # apply the weights to adjustment penalty
510         weight_bias=None,
511         weight_train=None,
512         weight_var=None,
513     )
514 )

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.

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.

539 fig = raf.plot_transform_matrix()
540 plotly.io.show(fig)
543 plotly.io.show(raf.figures["adj_size"])

The test error looks better than before.

548 # Transform and evaluate on the test set.
549 raf.transform_evaluate(
550     forecasts_test=forecasts_test,
551     actuals_test=actuals_test,
552     ipython_display=False,
553     plot=True,
554     plot_num_cols=2,
555 )
556 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


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:

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

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).

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.]])
  1. The constraint violation should be 0 for the actuals.

Out:

{'adjusted': 1.874446823677435e-15, 'forecast': 0.04015962097753542, 'actual': 1.9465566690613992e-16}
  1. 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.”).

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.

648 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.

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.

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]])}
  1. Check the convex optimization problem.

663 print(type(raf.prob))
664 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)))])

Gallery generated by Sphinx-Gallery