Reconciliation Diagnostics

!pip install hierarchicalforecast statsforecast datasetsforecast

import numpy as np
import pandas as pd

from datasetsforecast.hierarchical import HierarchicalData, HierarchicalInfo
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, Naive

from hierarchicalforecast.core import HierarchicalReconciliation
from hierarchicalforecast.methods import BottomUp, TopDown, MinTrace

group_name = 'TourismSmall'
group = HierarchicalInfo.get_group(group_name)
Y_df, S_df, tags = HierarchicalData.load('./data', group_name)
S_df = S_df.reset_index(names="unique_id")
Y_df['ds'] = pd.to_datetime(Y_df['ds'])

# Train/test split
Y_test_df = Y_df.groupby('unique_id').tail(group.horizon)
Y_train_df = Y_df.drop(Y_test_df.index)

print(f"Hierarchy levels: {list(tags.keys())}")
print(f"Total series: {len(S_df)}")
print(f"Bottom series: {S_df.shape[1] - 1}")

Hierarchy levels: ['Country', 'Country/Purpose', 'Country/Purpose/State', 'Country/Purpose/State/CityNonCity']
Total series: 89
Bottom series: 56

fcst = StatsForecast(
    models=[AutoARIMA(season_length=group.seasonality), Naive()],
    freq="QE",
    n_jobs=-1
)
Y_hat_df = fcst.forecast(df=Y_train_df, h=group.horizon)
Y_hat_df.head()

# Run reconciliation with diagnostics
hrec = HierarchicalReconciliation(reconcilers=[BottomUp(), MinTrace(method='ols')])
Y_rec_df = hrec.reconcile(
    Y_hat_df=Y_hat_df,
    Y_df=Y_train_df,
    S_df=S_df,
    tags=tags,
    diagnostics=True  # Enable diagnostics
)

# View overall coherence verification
coherence_metrics = hrec.diagnostics.query(
    "level == 'Overall' and metric in "
    "['coherence_residual_mae_before', 'coherence_residual_mae_after', 'is_coherent', 'coherence_max_violation']"
)
coherence_metrics

# View coherence residuals by hierarchy level
residuals_by_level = hrec.diagnostics.query(
    "metric in ['coherence_residual_mae_before', 'coherence_residual_mae_after']"
).pivot(index='level', columns='metric')
residuals_by_level

# Run multiple reconciliation methods
hrec_compare = HierarchicalReconciliation(reconcilers=[
    BottomUp(),
    TopDown(method='forecast_proportions'),
    MinTrace(method='ols'),
    MinTrace(method='wls_struct'),
])
Y_rec_compare = hrec_compare.reconcile(
    Y_hat_df=Y_hat_df,
    Y_df=Y_train_df,
    S_df=S_df,
    tags=tags,
    diagnostics=True
)

# Compare adjustment magnitude across methods (Overall level)
adjustment_comparison = hrec_compare.diagnostics.query(
    "level == 'Overall' and metric in ['adjustment_mae', 'adjustment_rmse', 'adjustment_max']"
)
adjustment_comparison

# Compare adjustments by hierarchy level for AutoARIMA forecasts
adjustment_by_level = hrec_compare.diagnostics.query("metric == 'adjustment_mae'")

# Pivot for easier comparison
adjustment_pivot = adjustment_by_level.set_index('level').drop(columns=['metric'])
adjustment_pivot.columns = [c.replace('AutoARIMA/', '') for c in adjustment_pivot.columns]
adjustment_pivot = adjustment_pivot[[c for c in adjustment_pivot.columns if 'AutoARIMA' in c or 'Naive' not in c]]
adjustment_pivot

# Create forecasts with some negative values to demonstrate
Y_hat_with_negatives = Y_hat_df.copy()

# Introduce some negative base forecasts at random locations in the bottom level
bottom_ids = tags['Country/Purpose/State/CityNonCity']
mask = Y_hat_with_negatives['unique_id'].isin(bottom_ids[:10])
Y_hat_with_negatives.loc[mask, 'AutoARIMA'] -= 5000
Y_hat_with_negatives.loc[mask, 'Naive'] -= 5000

print(f"Negative forecasts introduced for AutoARIMA: {(Y_hat_with_negatives['AutoARIMA'] < 0).sum()}")
print(f"Negative forecasts introduced for Naive: {(Y_hat_with_negatives['Naive'] < 0).sum()}")

Negative forecasts introduced for AutoARIMA: 33
Negative forecasts introduced for Naive: 36

# Run reconciliation with diagnostics
hrec_neg = HierarchicalReconciliation(reconcilers=[
    BottomUp(),
    MinTrace(method='ols'),
    MinTrace(method='ols', nonnegative=True),  # Non-negative constraint
])
Y_rec_neg = hrec_neg.reconcile(
    Y_hat_df=Y_hat_with_negatives,
    Y_df=Y_train_df,
    S_df=S_df,
    tags=tags,
    diagnostics=True
)

# Check negative value metrics at Overall level
negative_metrics = hrec_neg.diagnostics.query(
    "level == 'Overall' and metric in "
    "['negative_count_before', 'negative_count_after', 'negative_introduced', 'negative_removed']"
)
negative_metrics

# Check which levels have negative value issues
negatives_by_level = hrec_neg.diagnostics.query(
    "metric in ['negative_count_before', 'negative_count_after']"
).pivot(index='level', columns='metric')
negatives_by_level

# Export full diagnostics report
# hrec.diagnostics.to_csv("reconciliation_diagnostics.csv", index=False)

# Or export a summary
summary = hrec.diagnostics.query("level == 'Overall'").copy()
summary