使用 THIEF 进行时间序列聚合

!pip install hierarchicalforecast statsforecast datasetsforecast

import numpy as np
import pandas as pd

from datasetsforecast.m3 import M3

m3_monthly, _, _ = M3.load(directory='data', group='Monthly')
m3_quarterly, _, _ = M3.load(directory='data', group='Quarterly')

m3_monthly = m3_monthly.groupby("unique_id", group_keys=False)\
                       .apply(lambda x: x.tail(len(x) //  12 * 12))\
                       .reset_index(drop=True)

m3_quarterly = m3_quarterly.groupby("unique_id", group_keys=False)\
                           .apply(lambda x: x.tail(len(x) //  4 * 4))\
                           .reset_index(drop=True)

horizon_monthly = 24
horizon_quarterly = 8

m3_monthly_test = m3_monthly.groupby("unique_id", as_index=False).tail(horizon_monthly)
m3_monthly_train = m3_monthly.drop(m3_monthly_test.index)

m3_quarterly_test = m3_quarterly.groupby("unique_id", as_index=False).tail(horizon_quarterly)
m3_quarterly_train = m3_quarterly.drop(m3_quarterly_test.index)

spec_temporal_monthly = {"yearly": 12, "semiannually": 6, "fourmonthly": 4, "quarterly": 3, "bimonthly": 2, "monthly": 1}
spec_temporal_quarterly = {"yearly": 4, "semiannually": 2, "quarterly": 1}

from hierarchicalforecast.utils import aggregate_temporal

# Monthly
Y_monthly_train, S_monthly_train, tags_monthly_train = aggregate_temporal(df=m3_monthly_train, spec=spec_temporal_monthly)
Y_monthly_test, S_monthly_test, tags_monthly_test = aggregate_temporal(df=m3_monthly_test, spec=spec_temporal_monthly)

# Quarterly
Y_quarterly_train, S_quarterly_train, tags_quarterly_train = aggregate_temporal(df=m3_quarterly_train, spec=spec_temporal_quarterly)
Y_quarterly_test, S_quarterly_test, tags_quarterly_test = aggregate_temporal(df=m3_quarterly_test,  spec=spec_temporal_quarterly)

S_monthly_train.iloc[:5, :5]

S_monthly_test.iloc[:5, :5]

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

c:\Users\ospra\miniconda3\envs\hierarchicalforecast\lib\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

Y_hats = []
id_cols = ["unique_id", "temporal_id", "ds", "y"]

# We loop over the monthly and quarterly data
for tags_train, tags_test, Y_train, Y_test in zip([tags_monthly_train, tags_quarterly_train], 
                                                  [tags_monthly_test, tags_quarterly_test],
                                                  [Y_monthly_train, Y_quarterly_train], 
                                                  [Y_monthly_test, Y_quarterly_test]):
    # We will train a model for each temporal level
    Y_hats_tags = []
    for level, temporal_ids_train in tags_train.items():
        # Filter the data for the level
        Y_level_train = Y_train.query("temporal_id in @temporal_ids_train")
        temporal_ids_test = tags_test[level]
        Y_level_test = Y_test.query("temporal_id in @temporal_ids_test")
        # For each temporal level we have a different frequency and forecast horizon. We use the timestamps of the first timeseries to automatically derive the frequency & horizon of the temporally aggregated series.
        unique_id = Y_level_train["unique_id"].iloc[0]
        freq_level = pd.infer_freq(Y_level_train.query("unique_id == @unique_id")["ds"])
        horizon_level = Y_level_test.query("unique_id == @unique_id")["ds"].nunique()
        # Train a model and create forecasts
        fcst = StatsForecast(models=[AutoARIMA()], freq=freq_level, n_jobs=-1)
        Y_hat_level = fcst.forecast(df=Y_level_train[["ds", "unique_id", "y"]], h=horizon_level)
        # Add the test set to the forecast
        Y_hat_level = pd.concat([Y_level_test.reset_index(drop=True), Y_hat_level.drop(columns=["unique_id", "ds"])], axis=1)
        # Put cols in the right order (for readability)
        Y_hat_cols = id_cols + [col for col in Y_hat_level.columns if col not in id_cols]
        Y_hat_level = Y_hat_level[Y_hat_cols]
        # Append the forecast to the list
        Y_hats_tags.append(Y_hat_level)

    Y_hat_tag = pd.concat(Y_hats_tags, ignore_index=True)
    Y_hats.append(Y_hat_tag)

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

reconcilers = [
    BottomUp(),
    MinTrace(method="wls_struct"),
]
hrec = HierarchicalReconciliation(reconcilers=reconcilers)
Y_recs = []
# We loop over the monthly and quarterly data
for Y_hat, S, tags in zip(Y_hats, 
                          [S_monthly_test, S_quarterly_test], 
                          [tags_monthly_test, tags_quarterly_test]):
    Y_rec = hrec.reconcile(Y_hat_df=Y_hat, S=S, tags=tags, temporal=True)
    Y_recs.append(Y_rec)

from hierarchicalforecast.evaluation import evaluate
from utilsforecast.losses import mae

Y_rec_monthly = Y_recs[0]
evaluation = evaluate(df = Y_rec_monthly.drop(columns = 'unique_id'),
                      tags = tags_monthly_test,
                      metrics = [mae],
                      id_col='temporal_id',
                      benchmark="AutoARIMA")

evaluation.columns = ['level', 'metric', 'Base', 'BottomUp', 'MinTrace(wls_struct)']
numeric_cols = evaluation.select_dtypes(include="number").columns
evaluation[numeric_cols] = evaluation[numeric_cols].map('{:.2f}'.format).astype(np.float64)

evaluation

Y_rec_quarterly = Y_recs[1]
evaluation = evaluate(df = Y_rec_quarterly.drop(columns = 'unique_id'),
                      tags = tags_quarterly_test,
                      metrics = [mae],
                      id_col='temporal_id',
                      benchmark="AutoARIMA")

evaluation.columns = ['level', 'metric', 'Base', 'BottomUp', 'MinTrace(wls_struct)']
numeric_cols = evaluation.select_dtypes(include="number").columns
evaluation[numeric_cols] = evaluation[numeric_cols].map('{:.2f}'.format).astype(np.float64)

evaluation