Exogenous Variables

!pip install hierarchicalforecast statsforecast

import numpy as np
import pandas as pd

Y_df = pd.read_csv('https://raw.githubusercontent.com/Nixtla/transfer-learning-time-series/main/datasets/tourism.csv')
Y_df = Y_df.rename({'Trips': 'y', 'Quarter': 'ds'}, axis=1)
Y_df.insert(0, 'Country', 'Australia')
Y_df = Y_df[['Country', 'Region', 'State', 'Purpose', 'ds', 'y']]
Y_df['ds'] = Y_df['ds'].str.replace(r'(\d+) (Q\d)', r'\1-\2', regex=True)
Y_df['ds'] = pd.PeriodIndex(Y_df['ds'], freq='Q').to_timestamp()
Y_df.head()

np.random.seed(42)

# CPI: a national-level index, same across all bottom-level series per quarter
dates = Y_df['ds'].unique()
cpi_values = 100 + np.cumsum(np.random.normal(0.5, 0.3, len(dates)))
cpi_map = dict(zip(dates, cpi_values))
Y_df['cpi'] = Y_df['ds'].map(cpi_map)

# Marketing spend: varies by state, grows over time
state_spend_base = {s: np.random.uniform(50, 200) for s in Y_df['State'].unique()}
Y_df['marketing_spend'] = Y_df.apply(
    lambda r: state_spend_base[r['State']] * (1 + 0.02 * (r['ds'].year - 1998)) + np.random.normal(0, 5), 
    axis=1
).round(1)

Y_df.head()

spec = [
    ['Country'],
    ['Country', 'State'], 
    ['Country', 'Purpose'], 
    ['Country', 'State', 'Region'], 
    ['Country', 'State', 'Purpose'], 
    ['Country', 'State', 'Region', 'Purpose']
]

from hierarchicalforecast.utils import aggregate

Y_agg, S_df, tags = aggregate(
    Y_df, 
    spec, 
    exog_vars={
        'cpi': 'mean',              # rate variable -> average
        'marketing_spend': 'sum',   # additive variable -> sum
    }
)

# Top level: all of Australia
Y_agg[Y_agg['unique_id'] == 'Australia'].head()

# State level: Victoria
Y_agg[Y_agg['unique_id'] == 'Australia/Victoria'].head()

# Bottom level: a specific Region/Purpose combination
Y_agg[Y_agg['unique_id'] == 'Australia/Victoria/Melbourne/Holiday'].head()

Y_multi, _, _ = aggregate(
    Y_df, 
    spec, 
    exog_vars={
        'marketing_spend': ['sum', 'mean', 'std'],
    }
)

Y_multi[Y_multi['unique_id'] == 'Australia/Victoria'][
    ['unique_id', 'ds', 'y', 'marketing_spend_sum', 'marketing_spend_mean', 'marketing_spend_std']
].head()

Y_agg_no_exog, S_df_no_exog, tags_no_exog = aggregate(Y_df, spec)

horizon = 8

# With exogenous
Y_test_df = Y_agg.groupby('unique_id', as_index=False).tail(horizon)
Y_train_df = Y_agg.drop(Y_test_df.index)
X_test_df = Y_test_df[['unique_id', 'ds', 'cpi_mean', 'marketing_spend_sum']]

# Without exogenous
Y_test_no_exog = Y_agg_no_exog.groupby('unique_id', as_index=False).tail(horizon)
Y_train_no_exog = Y_agg_no_exog.drop(Y_test_no_exog.index)

print(f'Train: {Y_train_df.shape[0]:,} rows, Test: {Y_test_df.shape[0]:,} rows')

Train: 30,600 rows, Test: 3,400 rows

from statsforecast.models import AutoARIMA
from statsforecast.core import StatsForecast

/home/osprangers/Repositories/hierarchicalforecast/.venv/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

# With exogenous variables
fcst_exog = StatsForecast(models=[AutoARIMA(season_length=4)], freq='QS', n_jobs=-1)
Y_hat_exog = fcst_exog.forecast(df=Y_train_df, h=horizon, X_df=X_test_df, fitted=True)
Y_fitted_exog = fcst_exog.forecast_fitted_values()

Y_hat_exog.head()

# Without exogenous variables
fcst_no_exog = StatsForecast(models=[AutoARIMA(season_length=4)], freq='QS', n_jobs=-1)
Y_hat_no_exog = fcst_no_exog.forecast(df=Y_train_no_exog, h=horizon, fitted=True)
Y_fitted_no_exog = fcst_no_exog.forecast_fitted_values()

Y_hat_no_exog.head()

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

reconcilers = [
    BottomUp(),
    MinTrace(method='mint_shrink'),
]

# With exogenous
hrec_exog = HierarchicalReconciliation(reconcilers=reconcilers)
Y_rec_exog = hrec_exog.reconcile(
    Y_hat_df=Y_hat_exog, Y_df=Y_fitted_exog, S_df=S_df, tags=tags
)

Y_rec_exog.head()

# Without exogenous
hrec_no_exog = HierarchicalReconciliation(reconcilers=reconcilers)
Y_rec_no_exog = hrec_no_exog.reconcile(
    Y_hat_df=Y_hat_no_exog, Y_df=Y_fitted_no_exog, S_df=S_df_no_exog, tags=tags_no_exog
)

Y_rec_no_exog.head()

from hierarchicalforecast.evaluation import evaluate
from utilsforecast.losses import rmse, mase
from functools import partial

eval_tags = {
    'Total': tags['Country'],
    'Purpose': tags['Country/Purpose'],
    'State': tags['Country/State'],
    'Regions': tags['Country/State/Region'],
    'Bottom': tags['Country/State/Region/Purpose'],
}

# Evaluate with exogenous (drop exog columns before merge to avoid them being treated as models)
df_exog = Y_rec_exog.merge(Y_test_df[['unique_id', 'ds', 'y']], on=['unique_id', 'ds'])
eval_exog = evaluate(
    df=df_exog,
    tags=eval_tags,
    train_df=Y_train_df[['unique_id', 'ds', 'y']],
    metrics=[rmse, partial(mase, seasonality=4)],
)
eval_exog.columns = ['level', 'metric', 'Base (w/ exog)', 'BottomUp (w/ exog)', 'MinTrace (w/ exog)']

# Evaluate without exogenous
df_no_exog = Y_rec_no_exog.merge(Y_test_no_exog, on=['unique_id', 'ds'])
eval_no_exog = evaluate(
    df=df_no_exog,
    tags=eval_tags,
    train_df=Y_train_no_exog,
    metrics=[rmse, partial(mase, seasonality=4)],
)
eval_no_exog.columns = ['level', 'metric', 'Base (no exog)', 'BottomUp (no exog)', 'MinTrace (no exog)']

comparison = eval_exog.merge(eval_no_exog, on=['level', 'metric'])
comparison = comparison[[
    'level', 'metric', 
    'Base (no exog)', 'Base (w/ exog)', 
    'BottomUp (no exog)', 'BottomUp (w/ exog)',
    'MinTrace (no exog)', 'MinTrace (w/ exog)',
]]

numeric_cols = comparison.select_dtypes(include='number').columns
comparison[numeric_cols] = comparison[numeric_cols].map('{:.2f}'.format).astype(np.float64)

comparison.query('metric == "rmse"')

comparison.query('metric == "mase"')