Hierarchical Forecasting

Source Files

• twiga/hierarchy/structure.py - aggregate_panel, aggregate_levels, level_key

• twiga/hierarchy/alignment.py - align_to_level, alignment_multipliers

• twiga/hierarchy/reconcile.py - reconcile, Hierarchy, naive_residuals, shrunk_covariance

Products roll up into departments, stores and a total, and a plan is needed at every level. twiga.hierarchy offers two ways to make forecasts agree across levels: top-down alignment to one reference level, and reconciliation of forecasts made at every level.

Top-down alignment

A global model trained on thousands of products forecasts each one well, but errors it makes for every series at once (a missed trend, a level shift) add up when series are summed. On aggregated series (store, department, total) trend and level are far easier to see. Top-down alignment keeps the bottom model’s split across products and takes the level from a forecast made directly at a higher level:

\[ m_{g,t} = \frac{\hat{y}^{\text{level}}_{g,t}}{\sum_{i \in g} \hat{y}_{i,t}}, \qquad \tilde{y}_{i,t} = \hat{y}_{i,t}\,\bigl(1 + \lambda\,(m_{g,t} - 1)\bigr) \]

This is the adjustment the 2nd- and 5th-placed M5 Accuracy methods applied: N-BEATS forecasts at the top levels, or the recent store-department mean, as the reference. The damping \(\lambda\) (strength) plays the role of Anderer’s bias multiplier.

Workflow

from twiga.core.config import DataPipelineConfig
from twiga.forecaster.core import TwigaForecaster
from twiga.hierarchy import aggregate_panel, align_to_level

# 1. Bottom-level forecasts from a global model (panel mode).
bottom_preds, _ = bottom_forecaster.backtesting(panel)

# 2. Aggregate exactly the forecasted series to a level and forecast that small panel.
forecasted = panel[panel["id"].isin(bottom_preds["id"].unique())]
level = aggregate_panel(forecasted, ["store_id", "dept_id"], ["sales"], id_column="id", key_column="series")
level_config = DataPipelineConfig(
    target_feature="sales", period="1MS", lookback_window_size=24, forecast_horizon=6, id_column="series"
)
level_preds, _ = TwigaForecaster(level_config, level_model, cv_config).backtesting(level)

# 3. Rescale the bottom forecasts so each store-department sums to its level forecast.
aligned = align_to_level(
    bottom_preds,
    level_preds[["series", "timestamp", "forecast", "Folds"]],
    ["store_id", "dept_id"],
    id_column="id",
    tags=tags,  # one row per id with store_id, dept_id
    group_columns=["Folds"],
)

Things to get right:

• Same series on both sides. Aggregate only the series that have a bottom forecast. Series too short for a window, or launched later, would otherwise add to every level reference, and alignment would inflate the forecasted series to cover sales it never sees. When that set changes across backtest folds, build the level reference per fold.

• One reference per group and timestamp. level_forecasts must hold a single forecast per level group and timestamp (and fold); filter to one model first.

• Bounds. Multipliers are clipped to clip=(0.5, 2.0) by default to guard against near-zero bottom sums; strength < 1 damps the whole adjustment.

• Intervals and quantiles. Pass interval bounds or quantile columns in scale_columns to rescale them with the same factor.

Evaluate the result across levels with hierarchical_scores() (see Panel Metrics).

Example: M5 aggregated to months

All 30,490 series, 3 folds × 6 months; the bottom model is a global LightGBM, aligned to a level forecast made on the aggregated series. The panel forecasting tutorial runs the same steps on a 500-product sample:

Model	M5 score (12 levels)	± std over folds	Item-store WRMSSE
LightGBM aligned to Chronos-2 forecasts of store-department	0.992	0.076	1.016
LightGBM aligned to naive forecasts of store-department	1.056	0.016	1.022
Naive	1.077	0.021	1.124
LightGBM (no alignment)	1.178	0.349	1.030

Alignment fixed the global LightGBM’s shared under-forecast during rising sales. It made the model the best at every level and improved even the item levels. The alignment level and reference model were chosen on these same folds, so treat the best row as indicative.

Reconciliation

Forecasts made separately at each level rarely add up: the store forecasts do not sum to the total forecast. Reconciliation takes base forecasts \(\hat{y}\) at every node of the hierarchy and maps them to coherent forecasts, those that are sums of bottom-level forecasts through the summing matrix \(S\) (nodes × bottom series).

Theory

The linear reconcilers are projections. MinT (Wickramasuriya, Athanasopoulos and Hyndman, 2019) chooses the one that minimises the trace of the reconciled forecast-error covariance, given the covariance \(W\) of the base forecast errors:

\[ \tilde{y} = S\,(S^\top W^{-1} S)^{-1} S^\top W^{-1}\,\hat{y}. \]

It is unbiased whenever the base forecasts are, and it is the generalised least-squares fit of \(\hat{y} = S b + \varepsilon\) with \(\operatorname{Cov}(\varepsilon) = W\). With \(C = [\,I \mid -S_{\text{agg}}\,]\), one row per aggregate node so that \(C y = 0\) exactly when \(y\) is coherent, the same projection is

\[ \tilde{y} = \hat{y} - W C^\top (C W C^\top)^{-1} C\,\hat{y}, \]

which twiga solves. It needs one solve over the \(m\) aggregate nodes instead of the \(n_b\) bottom series (\(m = 12{,}350\) against \(n_b = 30{,}490\) on M5), and \(C W C^\top\) stays sparse when \(W\) is diagonal. The methods differ only in \(W\):

method	\(W\)	Reference	Needs residuals
bottom_up	none: every node is the sum of its bottom base forecasts		no
ols	\(I\)	Hyndman et al. (2011)	no
wls_struct	\(\operatorname{diag}(S\mathbf{1})\), the number of bottom series under each node	Athanasopoulos, Hyndman, Kourentzes and Petropoulos (2017)	no
wls_var	\(\operatorname{diag}(\hat{W}_1)\), each node’s mean squared residual	Hyndman, Lee and Wang (2016)	yes
mint_shrink	\(\hat{W}_1\) with its correlations shrunk toward zero, intensity \(\hat\lambda\) estimated from the data	Wickramasuriya, Athanasopoulos and Hyndman (2019); Schäfer and Strimmer (2005)	yes

\(\hat{W}_1\) is the uncentred covariance of one-step in-sample residuals; the \(h\)-step covariance is taken proportional to it, and a scalar factor does not change the projection. MinT-shrink keeps each node’s variance and shrinks its correlations:

\[ W = D^{1/2}\bigl(\hat\lambda\, I + (1-\hat\lambda)\,\hat{R}\bigr)D^{1/2}, \qquad \hat\lambda = \frac{\sum_{i\ne j}\widehat{\operatorname{Var}}(\hat r_{ij})}{\sum_{i\ne j}\hat r_{ij}^2}, \qquad \widehat{\operatorname{Var}}(\hat r_{ij}) = \frac{\overline{w^2}_{ij} - \hat r_{ij}^2}{T-1}, \]

where \(w_{tij} = z_{ti} z_{tj}\) for standardised residuals \(z\). This matches hts::shrink.estim, which the tests reproduce to machine precision.

Computation

The shrunk covariance is a diagonal plus a matrix of rank \(T\), the number of residual periods:

\[ W = \operatorname{diag}(a) + U U^\top, \qquad U = \sqrt{(1-\hat\lambda)/T}\; D^{1/2} Z^\top \quad (n \times T). \]

So \(C W C^\top = A + B B^\top\) with \(A = C \operatorname{diag}(a) C^\top\) sparse and \(B = C U\) only \(T\) columns wide. Twiga factorises \(A\) once (sparse LU) and applies the Woodbury identity

\[ (A + B B^\top)^{-1} = A^{-1} - A^{-1}B\,(I_T + B^\top A^{-1} B)^{-1} B^\top A^{-1}, \]

which costs a \(T \times T\) solve. The shrinkage intensity uses the same structure: both sums over pairs \(i \ne j\) are computed from \(T \times T\) Gram matrices (\(\sum_{ij} \hat r_{ij}^2 = \lVert Z Z^\top \rVert_F^2 / T^2\)), so no \(n \times n\) matrix is ever formed.

Step	Dense MinT	Twiga
Memory for \(W\)	\(O(n^2)\): 14.7 GB for M5’s 42,812 nodes	\(O(nT)\): 8 MB with \(T = 24\)
Shrinkage intensity	\(O(n^2 T)\)	\(O(n T^2)\)
Projection	\(O(m^3)\) dense solve	sparse LU of \(A\), plus \(O(m T^2)\)

All twelve M5 levels, 42,812 nodes, six months reconciled at once on a laptop CPU:

method	Time	Most negative forecast
bottom_up	0.3 s	0
ols	0.7 s	−24,332
wls_struct	0.7 s	−65
wls_var	0.7 s	−265
mint_shrink (\(\hat\lambda = 0.86\))	0.7 s	−340
wls_var, nonnegative="clip"	0.7 s	0
wls_var, nonnegative="optimal"	46 s	0

Groups that cover the same series, such as several models in one fold, share one hierarchy. For OLS, the incoherence of the total is spread equally over 30,000 products, which is why it produces large negative forecasts; weighted methods put it where the errors are.

Non-negative demand

A projection can turn small product forecasts negative.

• nonnegative="clip" (or True): sets negative bottom forecasts to zero and sums up again. The result is coherent and fast, but not the best non-negative solution.

• nonnegative="optimal": solves the constrained problem \(\min_{b \ge 0} (\hat{y} - S b)^\top W^{-1} (\hat{y} - S b)\) (Wickramasuriya, Turlach and Hyndman, 2020) as a bounded sparse least-squares problem per timestamp. It needs a diagonal \(W\); with mint_shrink it falls back to clipping, with a warning.

Usage

from twiga.hierarchy import aggregate_levels, naive_residuals, reconcile

LEVELS = {"total": [], "store": ["store_id"], "store-dept": ["store_id", "dept_id"], "item-store": ["id"]}

# Base forecasts at every level, one row per level, node and timestamp (and fold).
# Forecast each aggregated panel from aggregate_levels(history, LEVELS, ["sales"], id_column="id")
# with any twiga model, then stack the predictions with a "level" column.
residuals = naive_residuals(history_before_origin, tags, LEVELS, "sales", id_column="id", window=24)
reconciled = reconcile(
    base,
    tags,
    LEVELS,
    id_column="id",
    method="mint_shrink",
    residuals=residuals,
    group_columns=["Folds"],
    nonnegative="clip",
)
reconciled.attrs["diagnostics"]  # MinT shrinkage intensity per group

Things to get right:

• Per-group hierarchies. The hierarchy is built per group (e.g. per fold) from the bottom series that have a base forecast there, so the set of series may grow across folds. The aggregate base forecasts must describe that same set of series, as for alignment.

• Residuals. wls_var and mint_shrink should read the in-sample residuals of the base models themselves. When those are not at hand, naive_residuals gives one-step (seasonal) naive errors, the scale RMSSE uses: a proxy that ranks nodes by how hard they are to forecast. Pass only history before the forecast origin.

• Unobserved periods. A product before its launch has no error, not a zero error: variances use only the periods a node was observed. Nodes with no residuals at all borrow the median variance of their level, with a warning.

• Enough periods for MinT. The shrinkage needs at least 3 residual periods; with few periods relative to the number of nodes, \(\hat\lambda\) is close to 1 and MinT-shrink approaches wls_var.

• Point forecasts only. Reconciling quantiles is not a linear projection; reconcile sample paths, or the median as a point forecast.

• Alignment or reconciliation. Alignment trusts one level completely and keeps the bottom model’s split; reconciliation weighs every level by how reliable it is. Compare both with hierarchical_scores on backtest folds.

References

• Athanasopoulos, G., Hyndman, R. J., Kourentzes, N. and Petropoulos, F. (2017). Forecasting with temporal hierarchies. European Journal of Operational Research, 262(1), 60–74.

• Hyndman, R. J., Ahmed, R. A., Athanasopoulos, G. and Shang, H. L. (2011). Optimal combination forecasts for hierarchical time series. Computational Statistics & Data Analysis, 55(9), 2579–2589.

• Hyndman, R. J., Lee, A. J. and Wang, E. (2016). Fast computation of reconciled forecasts for hierarchical and grouped time series. Computational Statistics & Data Analysis, 97, 16–32.

• Schäfer, J. and Strimmer, K. (2005). A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics. Statistical Applications in Genetics and Molecular Biology, 4(1).

• Wickramasuriya, S. L., Athanasopoulos, G. and Hyndman, R. J. (2019). Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. Journal of the American Statistical Association, 114(526), 804–819.

• Wickramasuriya, S. L., Turlach, B. A. and Hyndman, R. J. (2020). Optimal non-negative forecast reconciliation. Statistics and Computing, 30(5), 1167–1182.

API Reference

Hierarchy structure: aggregate a long-format bottom-level panel to one or every level of the hierarchy.

twiga.hierarchy.structure.aggregate_levels(data, levels, value_columns, *, id_column, date_column='timestamp', tags=None, key_column='series', level_column='level', group_columns=())

Sum a bottom-level panel to every level of a hierarchy and stack the results.

The stacked frame holds one row per level, group and timestamp, which is the layout reconcile() reads base forecasts and residuals in.

Parameters:

• data (DataFrame) – Bottom-level long frame (history or forecasts).

• levels (Mapping[str, Sequence[str]]) – Level name to the columns that define it, e.g. {"total": [], "store": ["store_id"], "item-store": ["id"]}; the bottom level is the one defined by [id_column].

• value_columns (Sequence[str]) – Summable columns to aggregate.

• id_column (str) – Bottom-level series identifier.

• date_column (str) – Timestamp column.

• tags (DataFrame | None) – One row per bottom series with the hierarchy columns, when data lacks them.

• key_column (str) – Name of the group key column in the result.

• level_column (str) – Name of the level-name column in the result.

• group_columns (Sequence[str]) – Extra columns kept separate while summing (e.g. Model, Folds).

Return type:

DataFrame

Returns:

[level_column, key_column, *group_columns, date_column, *value_columns].

twiga.hierarchy.structure.aggregate_panel(data, level_columns, value_columns, *, id_column, date_column='timestamp', tags=None, key_column='series', group_columns=())

Sum a bottom-level panel to one level of the hierarchy.

The result is itself a long panel, one row per level group and timestamp, so any twiga model can forecast it (DataPipelineConfig(id_column=key_column)).

Parameters:

• data (DataFrame) – Bottom-level long frame (history or forecasts).

• level_columns (Sequence[str]) – Columns that define the level, e.g. ["store_id", "dept_id"]; empty for the grand total.

• value_columns (Sequence[str]) – Summable columns to aggregate (sales, forecasts, revenue, …).

• id_column (str) – Bottom-level series identifier.

• date_column (str) – Timestamp column.

• tags (DataFrame | None) – One row per bottom series with the hierarchy columns; only needed when data does not already carry the level columns.

• key_column (str) – Name of the level-group key column in the result.

• group_columns (Sequence[str]) – Extra columns kept separate while summing (e.g. Model, Folds).

Return type:

DataFrame

Returns:

[*group_columns, *level_columns, date_column, *value_columns, key_column].

twiga.hierarchy.structure.level_key(frame, level_columns)

One string key per row naming its group at a level ("CA_1/FOODS"; "Total" for the top).

Parameters:

• frame (DataFrame) – Frame holding the level columns.

• level_columns (Sequence[str]) – Columns that define the level; empty for the grand total.

Return type:

Series

Returns:

Key per row, aligned with frame.

Top-down alignment of bottom-level forecasts to forecasts made at a higher level.

A global model trained on thousands of products forecasts each one well, but errors it makes for every series at once (a missed trend, a level shift) add up when series are summed. Trend and level are far easier to forecast on aggregated series. Alignment keeps the bottom model’s split across products and takes the level from the aggregate forecast:

multiplier[g, t] = level_forecast[g, t] / sum_{i in g} bottom_forecast[i, t] aligned[i, t] = bottom_forecast[i, t] * (1 + strength * (multiplier[g, t] - 1))

This is the adjustment the 2nd- and 5th-placed M5 Accuracy methods applied (N-BEATS forecasts at the top levels, or the recent store-department mean, as the reference). strength < 1 damps the adjustment, like Anderer’s bias multiplier.

twiga.hierarchy.alignment.align_to_level(bottom, level_forecasts, level_columns, *, id_column, date_column='timestamp', forecast_column='forecast', level_key_column='series', level_forecast_column='forecast', tags=None, group_columns=(), strength=1.0, clip=(0.5, 2.0), scale_columns=())

Rescale bottom-level forecasts so each level group sums to its level forecast.

Each product keeps its share of the group; the group’s level (and so its trend) comes from the level forecast. Groups or timestamps without a level forecast are left unchanged.

The level forecasts must describe the same set of series as bottom: aggregate only the series that have a bottom forecast. Series without one (too short for a window, launched later) would otherwise add to every level reference, and alignment would inflate the forecasted series to cover sales it never sees.

Parameters:

• bottom (DataFrame) – Bottom-level forecasts, e.g. the predictions frame from backtesting().

• level_forecasts (DataFrame) – Forecasts of the aggregated series (see alignment_multipliers()).

• level_columns (Sequence[str]) – Columns that define the level, e.g. ["store_id", "dept_id"].

• id_column (str) – Bottom-level series identifier.

• date_column (str) – Timestamp column.

• forecast_column (str) – Bottom forecast column to align.

• level_key_column (str) – Level-group key column in level_forecasts.

• level_forecast_column (str) – Forecast column in level_forecasts.

• tags (DataFrame | None) – Hierarchy columns per bottom id, when bottom lacks the level columns.

• group_columns (Sequence[str]) – Columns matched between the frames besides key and timestamp (e.g. Folds).

• strength (float) – 1 applies the full multiplier, 0 none; values in between damp the adjustment.

• clip (tuple[float, float] | None) – Bounds on the multiplier. None disables.

• scale_columns (Sequence[str]) – Further columns rescaled by the same factor (interval bounds, quantiles).

Return type:

DataFrame

Returns:

bottom with forecast_column (and scale_columns) aligned and an alignment column holding the factor applied to each row.

twiga.hierarchy.alignment.alignment_multipliers(bottom, level_forecasts, level_columns, *, id_column, date_column='timestamp', forecast_column='forecast', level_key_column='series', level_forecast_column='forecast', tags=None, group_columns=(), clip=(0.5, 2.0))

Per level group and timestamp, the factor that makes bottom forecasts sum to the level forecast.

Parameters:

• bottom (DataFrame) – Bottom-level forecasts (long frame).

• level_forecasts (DataFrame) – Forecasts of the aggregated series, keyed by level_key_column (the key aggregate_panel() produces).

• level_columns (Sequence[str]) – Columns that define the level.

• id_column (str) – Bottom-level series identifier.

• date_column (str) – Timestamp column (in both frames).

• forecast_column (str) – Bottom forecast column.

• level_key_column (str) – Level-group key column in level_forecasts.

• level_forecast_column (str) – Forecast column in level_forecasts.

• tags (DataFrame | None) – Hierarchy columns per bottom id, when bottom lacks the level columns.

• group_columns (Sequence[str]) – Columns matched between the frames besides the key and timestamp (e.g. Folds); a bottom Model column is aligned per model automatically.

• clip (tuple[float, float] | None) – Bounds on the multiplier; guards against a near-zero bottom sum. None disables.

Return type:

DataFrame

Returns:

[*group, key, date, bottom_sum, level_forecast, multiplier] per level group and timestamp.

Forecast reconciliation: make forecasts at every level of a hierarchy add up.

Forecasts made separately at each level (total, store, department, product) rarely sum: the store forecasts do not add to the total forecast. Reconciliation projects the vector of base forecasts at every node, \(\hat{y}\), onto the coherent subspace spanned by the summing matrix \(S\) (which maps bottom series to every node):

\[\tilde{y} = S\,(S^\top W^{-1} S)^{-1} S^\top W^{-1}\,\hat{y} = \hat{y} - W C^\top (C W C^\top)^{-1} C\,\hat{y},\]

where \(C = [\,I \mid -S_{\text{agg}}\,]\) holds one aggregation constraint per aggregate node. The second form needs a solve over the aggregate nodes only and keeps \(W\) sparse when it is diagonal, so it scales to tens of thousands of series. The methods differ in \(W\):

• bottom_up: no projection; every node is the sum of its bottom base forecasts.

• ols: \(W = I\) (Hyndman et al., 2011).

• wls_struct: \(W = \operatorname{diag}(S\mathbf{1})\), the number of bottom series under each node (Athanasopoulos et al., 2017); needs no residuals.

• wls_var: \(W = \operatorname{diag}(\hat{\sigma}^2)\), each node’s residual variance.

• mint_shrink: \(W\) is the residual covariance with its correlations shrunk toward zero (Wickramasuriya, Athanasopoulos and Hyndman, 2019; Schäfer-Strimmer shrinkage).

The shrunk covariance is a diagonal plus a rank-\(T\) matrix (\(T\) residual periods), so MinT-shrink is solved with a sparse LU and a \(T \times T\) Woodbury correction, and its shrinkage intensity from \(T \times T\) Gram matrices: no \(n \times n\) matrix is ever formed, and all methods scale to every node of M5.

nonnegative="optimal" replaces the projection by the non-negative weighted least-squares solution (Wickramasuriya, Turlach and Hyndman, 2020) for the diagonal-weight methods.

Residuals are in-sample errors per node, in the layout aggregate_levels() produces; naive_residuals() builds one-step naive errors from history when the base models’ own residuals are not at hand.

class twiga.hierarchy.reconcile.Hierarchy(nodes, bottom, summing, n_aggregate)

Bases: object

The nodes of a hierarchy over a set of bottom series, and its summing matrix.

Variables:

• nodes – One row per node, aggregates first then the bottom series, with columns level and series (the key level_key() gives).

• bottom – Bottom series ids, in the column order of summing.

• summing – Sparse summing matrix S (nodes x bottom series).

• n_aggregate – Number of aggregate nodes (the first rows of nodes).

aggregate(bottom_values)

Every node’s value from the bottom values, S @ bottom_values.

Parameters:

bottom_values (ndarray) – (n_bottom,) or (n_bottom, k) values in bottom order.

Return type:

ndarray

Returns:

(n_nodes,) or (n_nodes, k) values in nodes order.

bottom: Index

classmethod from_tags(tags, levels, *, id_column, ids=None)

Build the hierarchy of the bottom series in tags (optionally only ids).

Parameters:

• tags (DataFrame) – One row per bottom series with id_column and every level column.

• levels (Mapping[str, Sequence[str]]) – Level name to the columns that define it; exactly one level must be [id_column] (the bottom level).

• id_column (str) – Bottom-level series identifier.

• ids (Sequence | None) – Restrict the hierarchy to these bottom series, e.g. the ones forecast in a fold.

Return type:

Hierarchy

Returns:

The Hierarchy.

n_aggregate: int

property n_nodes: int

Number of nodes (aggregates plus bottom series).

node_index()

(level, series) index of the nodes, in nodes order.

Return type:

MultiIndex

nodes: DataFrame

summing: csr_matrix

twiga.hierarchy.reconcile.naive_residuals(history, tags, levels, target_column, *, id_column, date_column='timestamp', seasonality=1, window=None, key_column='series', level_column='level', residual_column='residual')

In-sample (seasonal) naive one-step errors at every node, for wls_var and mint_shrink.

The error of the naive forecast is the scale RMSSE uses, so these weights treat each node by how hard it is to forecast. Pass only the history before a fold’s origin.

Parameters:

• history (DataFrame) – Bottom-level long history.

• tags (DataFrame) – One row per bottom series with the level columns (when history lacks them).

• levels (Mapping[str, Sequence[str]]) – Level name to defining columns.

• target_column (str) – Target column.

• id_column (str) – Bottom-level series identifier.

• date_column (str) – Timestamp column.

• seasonality (int) – Naive lag (1 = last value, 12 = same month last year on monthly data).

• window (int | None) – Keep the last window timestamps of errors; None keeps all.

• key_column (str) – Node key column in the result.

• level_column (str) – Level-name column in the result.

• residual_column (str) – Residual column in the result.

Return type:

DataFrame

Returns:

[level_column, key_column, date_column, residual_column].

twiga.hierarchy.reconcile.reconcile(base, tags, levels, *, id_column, method='wls_struct', residuals=None, date_column='timestamp', forecast_column='forecast', residual_column='residual', level_column='level', key_column='series', group_columns=(), nonnegative=False)

Reconcile base forecasts made at every level of a hierarchy.

The hierarchy is built per group (e.g. per backtest fold) from the bottom series that have a base forecast in that group, so the set of series may change across folds; base forecasts of the aggregate nodes must describe that same set of series.

Parameters:

• base (DataFrame) – Base forecasts, one row per level, node, timestamp (and group), with level_column, key_column, date_column and forecast_column; node keys as aggregate_levels() names them. bottom_up only reads the bottom level.

• tags (DataFrame) – One row per bottom series with id_column and every level column.

• levels (Mapping[str, Sequence[str]]) – Level name to defining columns; one level must be [id_column].

• id_column (str) – Bottom-level series identifier (column of tags).

• method (Literal['bottom_up', 'ols', 'wls_struct', 'wls_var', 'mint_shrink']) – "bottom_up", "ols", "wls_struct", "wls_var" or "mint_shrink".

• residuals (DataFrame | None) – In-sample residuals per node (level_column, key_column, date_column, residual_column, and optionally group_columns), for wls_var and mint_shrink; see naive_residuals().

• date_column (str) – Timestamp column.

• forecast_column (str) – Base forecast column; the result holds the reconciled forecast here.

• residual_column (str) – Residual column of residuals.

• level_column (str) – Level-name column.

• key_column (str) – Node key column.

• group_columns (Sequence[str]) – Columns reconciled separately (e.g. Folds, Model).

• nonnegative (Union[bool, Literal['clip', 'optimal']]) – Keep the reconciled forecasts non-negative (demand, sales). "clip" (or True) sets negative bottom forecasts to zero and sums up again: coherent and fast, but not the best non-negative solution. "optimal" solves the non-negative weighted least-squares problem instead (diagonal-weight methods only; mint_shrink falls back to "clip").

Return type:

DataFrame

Returns:

[*group_columns, level_column, key_column, date_column, forecast_column, "base_forecast"] for every node; base_forecast is NaN where a node had none (always the case for aggregates under bottom_up). Diagnostics such as the MinT shrinkage intensity per group are in result.attrs["diagnostics"].

twiga.hierarchy.reconcile.shrinkage_intensity(standardised)

Schäfer-Strimmer shrinkage intensity of the residual correlations toward zero.

\[\hat\lambda = \frac{\sum_{i \ne j} \widehat{\operatorname{Var}}(\hat r_{ij})}{\sum_{i \ne j} \hat r_{ij}^2}, \qquad \widehat{\operatorname{Var}}(\hat r_{ij}) = \frac{1}{T - 1}\Bigl(\overline{w^2}_{ij} - \hat r_{ij}^2\Bigr), \quad w_{tij} = z_{ti} z_{tj},\]

as in Schäfer and Strimmer (2005) and the MinT-shrink estimator of Wickramasuriya, Athanasopoulos and Hyndman (2019) (hts::shrink.estim). Both sums are computed through T x T Gram matrices, O(n T^2) time and O(n T) memory, so the n x n correlation matrix is never formed.

Parameters:

standardised (ndarray) – (T, n) residuals divided by their root mean square (see _scaled_residuals()), zero where unobserved.

Return type:

float

Returns:

The intensity in [0, 1]; 1 when there are no off-diagonal correlations.

twiga.hierarchy.reconcile.shrunk_covariance(residuals)

Residual covariance with its correlations shrunk toward zero (Schäfer and Strimmer, 2005).

The dense form of the MinT-shrink weight matrix, for inspection and small hierarchies; reconcile() uses its diagonal-plus-low-rank form and never builds it.

Parameters:

residuals (ndarray) – (T, n) residuals, one column per node; NaN where a node is unobserved (e.g. a product before its launch).

Return type:

tuple[ndarray, float]

Returns:

(covariance, shrinkage) – the (n, n) shrunk covariance and the intensity in [0, 1].