はじめに

時系列予測とは、過去のデータから未来の値を予測するタスクです。売上予測、需要予測、株価予測など、ビジネスの様々な場面で活用されています。
本記事では、Pythonを用いてKaggleのStore Salesデータセットを分析し、時系列予測の基本的な手法を学びます。

実装例について
本記事では、KaggleのStore Salesコンペのデータを一例として使用します。

このコンペは、実際の小売店の売上データを用いて、将来の売上を予測するタスクです。多様な時系列予測手法を試すのに適したデータセットとなっています。

https://www.kaggle.com/competitions/store-sales-time-series-forecasting

前処理

データの読み込み

とりあえず、全データを読み込みます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import os
for dirname, _, filenames in os.walk('/kaggle/input'):
    for filename in filenames:
        print(os.path.join(dirname, filename))

from pathlib import Path
comp_dir = Path('/kaggle/input/store-sales-time-series-forecasting/')

oil = pd.read_csv(comp_dir / 'oil.csv')
holidays_events = pd.read_csv(comp_dir / 'holidays_events.csv')
stores = pd.read_csv(comp_dir / 'stores.csv')
train = pd.read_csv(comp_dir / 'train.csv')
test = pd.read_csv(comp_dir / 'test.csv')
transactions = pd.read_csv(comp_dir / 'transactions.csv')

各データの確認

各データに対しての情報を確認します。

trainデータ

train.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 3000888 entries, 0 to 3000887
Data columns (total 6 columns):
 #   Column       Dtype
---  ------       -----
 0   id           int64
 1   date         object
 2   store_nbr    int64
 3   family       object
 4   sales        float64
 5   onpromotion  int64
dtypes: float64(1), int64(3), object(2)
memory usage: 137.4+ MB

print(train['store_nbr'].unique())
print(train['family'].unique())

[ 1 10 11 12 13 14 15 16 17 18 19  2 20 21 22 23 24 25 26 27 28 29  3 30
 31 32 33 34 35 36 37 38 39  4 40 41 42 43 44 45 46 47 48 49  5 50 51 52
 53 54  6  7  8  9]
['AUTOMOTIVE' 'BABY CARE' 'BEAUTY' 'BEVERAGES' 'BOOKS' 'BREAD/BAKERY'
 'CELEBRATION' 'CLEANING' 'DAIRY' 'DELI' 'EGGS' 'FROZEN FOODS' 'GROCERY I'
 'GROCERY II' 'HARDWARE' 'HOME AND KITCHEN I' 'HOME AND KITCHEN II'
 'HOME APPLIANCES' 'HOME CARE' 'LADIESWEAR' 'LAWN AND GARDEN' 'LINGERIE'
 'LIQUOR,WINE,BEER' 'MAGAZINES' 'MEATS' 'PERSONAL CARE' 'PET SUPPLIES'
 'PLAYERS AND ELECTRONICS' 'POULTRY' 'PREPARED FOODS' 'PRODUCE'
 'SCHOOL AND OFFICE SUPPLIES' 'SEAFOOD']

oilデータ

oil.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1218 entries, 0 to 1217
Data columns (total 2 columns):
 #   Column      Non-Null Count  Dtype
---  ------      --------------  -----
 0   date        1218 non-null   object
 1   dcoilwtico  1175 non-null   float64
dtypes: float64(1), object(1)
memory usage: 19.2+ KB

データに欠損値があることがわかります。

fig = oil.plot()
fig.set_title('oil')
fig.set_xlabel('Date')

欠損値はありますが、石油価格は連続的に変動することを考慮し、前後の値で補完します。

holidays_eventsデータ

holidays_events.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 350 entries, 0 to 349
Data columns (total 6 columns):
 #   Column       Non-Null Count  Dtype
---  ------       --------------  -----
 0   date         350 non-null    object
 1   type         350 non-null    object
 2   locale       350 non-null    object
 3   locale_name  350 non-null    object
 4   description  350 non-null    object
 5   transferred  350 non-null    bool
dtypes: bool(1), object(5)
memory usage: 14.1+ KB

欠損値はありません。

しかし、祝日が移動している場合（transferred=True）は、祝日が移動していることを意味するため、その列は無視します。また、休日が重複している場合があるため、重複数を確認します。

display(holidays_events)

dup_days=holidays_events.duplicated(subset='date').sum()
print(f"重複数 {dup_days}")

fig, ax = plt.subplots(1,2, figsize=(15,5))
holidays_events['type'].value_counts().plot.bar(ax=ax[0])
ax[0].set_title('holidays_events type counts')
ax[0].set_xlabel('Type')
ax[0].set_ylabel('Counts')

holidays_events['locale'].value_counts().plot.bar(ax=ax[1])
ax[1].set_title('holidays_events locale counts')
ax[1].set_xlabel('Type')
ax[1].set_ylabel('Counts')

このコンペでは、店舗store_nbrごと・商品カテゴリfamilyごとに売上salesを予測する必要があるということが、trainとtestデータの差分からわかります。
また、その期間は、trainデータが2013年1月1日から2017年8月15日まで、testデータが2017年8月16日から2017年8月31日までとなっています。

データの結合と前処理

分析を行うために、必要なデータを結合してしまいましょう。
このとき、date列は日付型に変換しておくと便利です。（読み込み時に変換しても良い）
ただし、dateindexに設定すると、結合が面倒になる場合があるため、date列はそのままにしておきます。

今回は、trainデータを中心に、storesデータ、oilデータ、holidays_eventsデータを結合します。
ただし、holidays_eventsデータは重複があるため、休日かどうかのみを判定して結合します。
transferred=Trueの場合は、祝日が移動していることを意味するため、その列は無視します。
また、oilデータの欠損値は前後の値で補完します。これは、石油価格が連続的に変動することを考慮しています。

# 休日フラグの作成
holidays_events = holidays_events[holidays_events['transferred'] == False]
holidays_events['is_holiday'] = 1
holidays_events = holidays_events[['date', 'is_holiday']].drop_duplicates(subset=['date'])

# データの結合
train_merged = train.merge(stores, on='store_nbr', how='left')
train_merged = train_merged.merge(oil, on='date', how='left')
train_merged['dcoilwtico'] = train_merged['dcoilwtico'].ffill().bfill()
train_merged = train_merged.merge(holidays_events, on='date', how='left')
train_merged['is_holiday'] = train_merged['is_holiday'].fillna(0)

# データ型の変換とソート
train_merged = train_merged.sort_values(['store_nbr', 'family', 'date'])
train_merged = train_merged.astype({
    'store_nbr': 'str',
    'family': 'str',
    'city': 'str',
    'state': 'str',
    'type': 'str',
    'cluster': 'str',
    'is_holiday': 'uint8',
    'dcoilwtico': 'float32',
})
train_merged['date'] = pd.to_datetime(train_merged['date'])

# 日付関連の特徴量作成
train_merged['year'] = train_merged.date.dt.year
train_merged['month'] = train_merged.date.dt.month
train_merged['dayofmonth'] = train_merged.date.dt.day
train_merged['dayofweek'] = train_merged.date.dt.dayofweek
train_merged['dayname'] = train_merged.date.dt.strftime('%A')

train_merged.head()

読み込み時に使用するものを確定している場合は、データの前処理として、日時データはparse_datesで日付型に変換し、index_colでインデックスに設定する。
また、各列のデータ型をdtypesで指定して読み込みます。

import numpy as np
import pandas as pd

df = pd.read_csv(
    'path/to/your/data.csv',
    use_cols=['store_nbr', 'family', 'date', 'sales', 'onpromotion'],
    dtype={
        'store_nbr': 'str',
        'family': 'str',
        'sales': 'float32',
        'onpromotion': 'uint32',
    },
    parse_dates=['date_column'],
    index_col='date_column'
)

後から、date型に変換や、インデックス設定を行う場合は以下のようにします。

df['date'] = pd.to_datetime(df['date'])
# inplace=Trueで元のDataFrameを直接変更
df.set_index('date', inplace=True)

# datetimeIndexに変換
df.index = pd.DatetimeIndex(df.index)

# 型を変換
df['store_nbr'] = df['store_nbr'].astype('str')

データ分析（EDA）

モデルの作成の前に、データの可視化・分析を行い、データの理解を深めます。

EDA（探索的データ分析：Exploratory Data Analysis）とは、データ分析の初期段階で行われ、データに触れて特徴やパターンを掴み、仮説を立てるためのプロセスのことです。

時系列データの可視化・分析

時系列データの理解を深めるために、まずはデータの可視化と基本的な分析を行います。

前処理として、日次データと月次データにリサンプリングします。

daily_df = train_merged.copy()
daily_df.set_index('date', inplace=True)
daily_df.index = pd.DatetimeIndex(daily_df.index)
daily_df = daily_df['sales'].resample('D').mean().fillna(0)

month_df = train_merged
month_df.set_index('date', inplace=True)
month_df.index = pd.DatetimeIndex(month_df.index)
month_df = month_df['sales'].resample('ME').mean().fillna(0)

STL分解

STL分解（Seasonal-Trend decomposition using Loess）は、時系列データを以下の3つの成分に分解する手法です。

Trend（トレンド）: 長期的な傾向を示す成分
Seasonal（季節性）: 定期的なパターンを示す成分
Residual（残差）: トレンドと季節性を除いた後のランダムな成分

import statsmodels.api as sm

stl_month = sm.tsa.STL(month_df).fit()

fig = stl_month.plot()
fig.set_size_inches(15,10)

# fig, axes = plt.subplots(4, 1, figsize=(15,10))
# month_df.plot(ax=axes[0])
# stl.trend.plot(ax=axes[1])
# stl.seasonal.plot(ax=axes[2])
# stl.resid.plot(ax=axes[3])

一番上に表示されているグラフは元のデータです。

結論

Trend（トレンド）: グラフを見ると右肩上がりです。つまり、店全体の売上は年々成長していることがわかります。
Seasonal（季節性）: 定期的な波が見られます。特に1年単位での強いパターンがあることが示唆されます。
Residual（残差）: トレンドと季節性を除いた不規則な動きです。ここが大きい場合、突発的なイベント（地震や特異なセールなど）の影響がある可能性があります。
トレンドがあることからデータは非定常であり、予測モデルを構築する際にはトレンド成分を考慮する必要があります。

自己相関の確認

自己相関とは、時系列データが過去の値とどの程度関連しているかを示す指標です。自己相関を確認することで、データの周期性やトレンドを把握できます。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# コレログラムを描画
fig, ax = plt.subplots(2, 1, figsize=(15,7))
plot_acf (daily_df, lags=40, ax=ax[0]) # 自己相関係数
plot_pacf(daily_df, lags=27, ax=ax[1]) # 偏自己相関係数
plt.xlabel('Lag [day]')
plt.tight_layout()
plt.show()

ACF（上段のグラフ）は、過去の値との相関を示しています。

PACF（下段のグラフ）は、間の影響を取り除いて、「そのラグの値」が「現在の値」に直接与える影響力を示します（例：昨日の影響を除いた、一昨日の純粋な影響）。

結論

ACFの値がなかなか下がらないことから、データにトレンド成分が含まれている（非定常）であることを示唆しています。（これは、STL分解で確認したトレンド成分と一致します）
Lag 1は高い値を示し、Lag 2で急激に減少していることから、直近の売上が翌日の売上に強く影響していることがわかります。
Lag 7でも高い自己相関が見られ、これは週次の周期性を示唆しています。

月次データの場合

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# コレログラムを描画
fig, ax = plt.subplots(2, 1, figsize=(15,7))
plot_acf (month_df, lags=40, ax=ax[0]) # 自己相関係数
plot_pacf(month_df, lags=27, ax=ax[1]) # 偏自己相関係数
plt.xlabel('Lag [month]')
plt.tight_layout()
plt.show()

結論

なかなかACFの値が下がらないことから、データにトレンド成分が含まれている（非定常）であることを示唆しています。（これは、STL分解で確認したトレンド成分と一致します）
PACFでは、Lag 1が高い値を示し、それ以外のラグでは急激に減少していることから、直近の売上が翌月の売上に強く影響していることがわかります。

ラグプロット

実際にラグごとのプロットを確認してみます。

from learntools.time_series.utils import plot_lags

plot_lags(daily_df, lags=8, nrows=2);

結論

1日のラグ（Lag 1）と7日のラグ（Lag 7）で強い相関が見られます。これは、先ほどの自己相関の結果と一致しています。

月次データの場合

plot_lags(month_df, lags=16, nrows=2);

結論

1ヶ月のラグ（Lag 1）で強い相関が見られます。これは、先ほどの自己相関の結果と一致しています。

周期性の確認

どのくらいの周期でデータが繰り返されているかを確認するために、周期プロットとパワースペクトル密度（Periodogram）を使用します。

import datetime

from learntools.time_series.style import *  # plot style settings
from learntools.time_series.utils import plot_periodogram, seasonal_plot

df = daily_df.to_frame(name='sales')

df['week'] = df.index.isocalendar().week
df['day'] = df.index.dayofweek
seasonal_plot(df, y='sales', period='week', freq='day');

結論

週次データの場合、曜日ごとの売上の違いが明確に見られます。
例えば、金曜日と土曜日の売上が他の曜日に比べて高い傾向があります。

月次データの場合

import datetime

from learntools.time_series.style import *  # plot style settings
from learntools.time_series.utils import plot_periodogram, seasonal_plot

df = month_df.to_frame(name='sales')

df['year'] = df.index.year
df['month'] = df.index.month
seasonal_plot(df, y='sales', period='year', freq='month');

結論

12月に売上がピークになる傾向が見られます。これは、年末のホリデーシーズンに関連している可能性があります。

パワースペクトル密度（Periodogram）

どの周期が支配的かを周波数解析しています。

from learntools.time_series.utils import plot_periodogram

fig, ax = plt.subplots(1,1, figsize=(15,5))
plot_periodogram(daily_df.squeeze(), ax=ax);

結論

Weekly (週次)が最も支配的な周期であることがわかります。これは、先ほどの周期プロットの結果と一致しています。

date列以外の特徴量分析

date以外の情報についても分析をしてみます。
この分析によって、どの特徴量を使用するべきかなどを考えます。

全データの相関分析

!pip install dython

from dython.nominal import associations

result = associations(train_merged, figsize=(10, 10), cmap='coolwarm')
corr_matrix = result['corr']
print(corr_matrix)

dythonをインストールするとstatsmodelsのバージョンが上がり、既存のコードが動かなくなる可能性があります。

結論

salesとfamilyに強い相関があります。つまり、商品カテゴリを考慮することが売上予測に重要である可能性があります。
salesとonpromotionにも正の相関があります。つまり、販促活動が売上にプラスの影響を与えている可能性があります。

onpromotion

onpromotion（販促活動を行っているかどうか）の影響を確認します。

試したコード例

plt.figure()
train_merged.plot.scatter(x='onpromotion',y='sales')

train_merged.groupby('store_nbr')[['onpromotion','sales']].sum().plot.scatter(x='onpromotion',y='sales')

train_merged.groupby('family')[['onpromotion','sales']].sum().plot.scatter(x='onpromotion',y='sales')

train_merged.groupby('family')[['onpromotion','sales']].sum().plot.scatter(x='onpromotion',y='sales', logy=True)

結論

onpromotionとsalesには正の相関が見られます。つまり、販促活動が売上にプラスの影響を与えている可能性があります。
salesはlog scaleの方が良さそうです。

cluster(stores)

import matplotlib.cm as cm

fig, ax = plt.subplots(1,1, figsize=(15,5))

clusters = train_merged['cluster'].unique()

colors = cm.tab20(np.linspace(0, 1, len(clusters)))
for i, cluster in enumerate(clusters):
    df = cluster_daily.query('cluster == @cluster')[['date', 'sales']]
    df.set_index('date',inplace=True)

    plot_periodogram(df.squeeze(), ax=ax)
    ax.lines[-1].set_color(colors[i])
    ax.lines[-1].set_label(f'Cluster {cluster}')

結論

クラスターごとの周期性はそれほど変わらないようです。

is_holiday

sns.stripplot(x="is_holiday", y="sales", data=train_merged)

モデルの構築

データ準備

まずモデルを構築する前に、データの前処理を train と test の両方に対して行います。

ここで、train と test を結合して一括で前処理を行うと、同じ前処理が両方に適用されるため、コードが簡潔になります。

データ結合・前処理の関数

def create_all_df(train, test, stores, oil, is_holiday_df):
    train_withflg = train.copy()
    train_withflg['is_train'] = True

    test_withflg = test.copy()
    test_withflg['is_train'] = False

    all_df = pd.concat([train_withflg, test_withflg], axis=0, ignore_index=True).sort_values('date')

    all_df = all_df.merge(stores, on='store_nbr', how='left')

    all_df = all_df.merge(oil, on='date', how='left')
    all_df['dcoilwtico'] = all_df['dcoilwtico'].ffill().bfill()
    all_df['oil_ma7'] = all_df['dcoilwtico'].rolling(7).mean()

    all_df = all_df.merge(is_holiday_df, on='date', how='left')
    all_df['is_holiday'] = all_df['is_holiday'].fillna(0)

    all_df = all_df.sort_values(['store_nbr', 'family', 'date'])
    all_df = all_df.astype({
        'store_nbr': 'uint8',
        'family': 'str',
        'city': 'str',
        'state': 'str',
        'type': 'str',
        'cluster': 'str',
        'is_holiday': 'uint8',
        'dcoilwtico': 'float32',
    })
    all_df['date'] = pd.to_datetime(all_df['date'])

    return all_df

データを結合した後、モデル予測に用いる特徴量を追加します。

追加する特徴量の例として、以下のようなものがあります。

deterministic Process特徴量

この特徴量は、statsmodelsのDeterministicProcessを使用して作成します。

トレンド成分（線形トレンドなど）
季節成分（曜日ダミー、Fourier項など）

ラグ特徴量

過去の売上データを用いて、以下のようなラグ特徴量を作成します。

どのラグを使用するかは、先ほどの自己相関分析の結果を参考にします。

lag_1: 1日前の売上
lag_7: 7日前の売上
lag_14: 14日前の売上
roll_mean_7: 過去7日間の移動平均

その他の特徴量

年、月、日、曜日、週末フラグなどの日時特徴量
カテゴリ変数のエンコード（Label Encodingなど）

特徴量追加の関数

from sklearn.preprocessing import LabelEncoder
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess

def create_lags(df):
    # lag特徴量の作成
    df['lag_1'] = df.groupby(['store_nbr', 'family'])['sales'].shift(1)
    df['lag_7'] = df.groupby(['store_nbr', 'family'])['sales'].shift(7)
    df['lag_14'] = df.groupby(['store_nbr', 'family'])['sales'].shift(14)
    df['roll_mean_7'] = df.groupby(['store_nbr', 'family'])['lag_1'].transform(lambda x: x.rolling(7).mean())
    return df


def add_feature(df):
    df['year'] = df.date.dt.year
    df['month'] = df.date.dt.month
    df['dayofmonth'] = df.date.dt.day
    df['dayofweek'] = df.date.dt.dayofweek
    df['weekend'] = (df['dayofweek'] >= 5).astype(int)


    le = LabelEncoder()

    categorical_cols = ['family', 'city', 'state', 'type']
    for col in categorical_cols:
        df[col] = df[col].astype(str) # 欠損埋めのために文字列化
        df[col + '_enc'] = le.fit_transform(df[col])

    # store_nbrとclusterは元々数値ですが、念のためエンコード扱いの列も作っておく
    df['store_enc'] = le.fit_transform(df['store_nbr'])
    df['cluster_enc'] = le.fit_transform(df['cluster'])


    # 時間特徴量の作成
    unique_dates = np.sort(pd.to_datetime(df['date'].unique()))
    period_index = pd.date_range(
        start=unique_dates.min(),
        end=unique_dates.max(),
        freq='D'
    )

    # fourier_weekly = CalendarFourier(freq="W", order=4)
    fourier_yearly = CalendarFourier(freq="A", order=10)
    dp = DeterministicProcess(
        index=period_index,
        constant=True,       # 切片
        order=1,             # 線形トレンド
        seasonal=True,       # 曜日ダミー (Weeklyの強いスパイク対策)
        additional_terms=[fourier_yearly],
        drop=True
    )

    X_time = dp.in_sample()
    X_time['date'] = X_time.index

    dp_features = [c for c in X_time.columns if c != 'date']

    df = df.merge(X_time, on='date', how='left')

    df = create_lags(df)

    return df, dp_features

all_df = create_all_df(train, test, stores, oil, is_holiday_df)

all_df_with_feature, dp_features = add_feature(all_df)

xgb_features = [
    'oil_ma7', 'onpromotion',
    'family_enc', 'store_enc', 'cluster_enc', 'type_enc',
    'lag_1', 'lag_7', 'roll_mean_7',
    'is_holiday', 'weekend',
    'month', 'dayofweek', 'dayofmonth'
]

train_mask = (all_df_with_feature['date'] >= '2016-01-01') & (all_df_with_feature['is_train'] == True)
train_data = all_df_with_feature[train_mask].copy()
train_data = train_data.dropna(subset=['lag_1', 'lag_7', 'roll_mean_7'])

test_data = all_df_with_feature[all_df_with_feature['is_train'] == False]

display(train_data.head())
display(test_data.head())

モデル選択と構築

この章は、以下の記事の章を参考にしています。

https://medium.com/@alan.vourch/forecasting-a-practical-guide-6173f421c1ed

時系列予測にはさまざまなアプローチがありますが、次の一時点の予測だけではなく、複数の将来の時点を同時に予測する（multi-step forecasting）ことが多いです。

実際には以下のような図となっています。

横軸が時間とした時の、青色が予測に使用できるデータ、赤色が予測したいデータです。

origin: 予測に使用できる最新のデータが存在する時点
lead time: 予測を開始するまでの時間
forecast horizon: 予測したい将来の時点の数

lead timeは、予測に使用できる最新のデータが存在する時点から、予測を開始するまでの時間を指します。

例えば、ある日のデータを使って翌日の売上を予測する場合、lead timeは1日となります。

マルチステップ予測の戦略

マルチステップ予測には、いくつかの異なる戦略があります。代表的なものを以下に示します。

戦略	説明	メリット	デメリット
Multioutputモデル	一つのモデルが複数の将来の時点を同時に予測します。例えば、1日後、2日後、3日後の売上を一度に予測します。	- 計算効率が良い - 異なる時点間の関係性を学習 - 予測値の一貫性が保たれやすい	- 複雑なモデル構造が必要 - 長期予測には向かない - 出力が増えるため訓練が困難
Recursiveモデル	一つのモデルが次の時点を予測し、その予測値を使ってさらに次の時点を予測します。例えば、1日後の売上を予測し、その予測値を使って2日後の売上を予測します。	- モデルが単純 - 長期予測に対応可能 - 訓練データ量が少なくて済む	- 予測誤差が累積 - 計算量が多い - 長期予測の精度低下
Directモデル	各将来の時点に対して個別のモデルを構築します。例えば、1日後の売上を予測するモデル、2日後の売上を予測するモデル、3日後の売上を予測するモデルをそれぞれ作成します。	- 各時点に最適化可能 - 誤差が累積しない - 並列計算が可能	- 訓練と管理が複雑 - 計算量が多い - 時点間の関係性を考慮しにくい
Hybridモデル	上記の戦略を組み合わせたアプローチです。例えば、トレンド成分を線形モデルで、残差をXGBoostで予測するなど。	- 複数の手法の長所を活かせる - 予測精度が向上する可能性 - 柔軟な設計が可能	- モデル設計が複雑 - 手動調整が必要 - 解釈性が低下

Hybrid Multioutputモデル

時系列予測には、さまざまな手法がありますが、まずは代表的なliner回帰モデル+XGBoostを用いたHybrid Multioutputモデルを使用してみます。

STL分解において、トレンド成分が存在することがわかっているため、トレンド成分を線形回帰モデルで予測し、残差成分をXGBoostで予測するハイブリッドアプローチを採用します。

class BoostedHybrid:
    def __init__(self, model_1, model_2):
        self.model_1 = model_1
        self.model_2 = model_2
        self.y_columns = None

    def fit(self, X_1, X_2, y):
        self.model_1.fit(X_1, y)
        y_fit = self.model_1.predict(X_1)

        y_resid = y - y_fit

        self.model_2.fit(X_2, y_resid)

        self.y_fit_ = y_fit
        self.y_resid_ = y_resid

    def predict(self, X_1, X_2):
        y_pred = self.model_1.predict(X_1)
        y_pred += self.model_2.predict(X_2)
        return y_pred

モデルのoutputには、recursive戦略を採用してみます。

実際のコード例

from tqdm import tqdm

def recursive_predict(model, history_df, target_df, features_1, features_2):

    predictions = []

    target_dates = target_df['date'].unique()
    target_dates = np.sort(target_dates)

    print(f"Starting recursive prediction for {len(target_dates)} days...")

    work_df = pd.concat([history_df, target_df], ignore_index=True).sort_values(['store_nbr', 'family', 'date'])

    mask_future = work_df['date'].isin(target_dates)
    work_df.loc[mask_future, 'sales'] = np.nan

    for target_date in tqdm(target_dates):
        # 1. その時点でのデータでラグを作成
        # (直前の予測値が埋まった状態のwork_dfを使う)
        work_df = create_lags(work_df)

        # 2. 当日データの抽出
        current_data = work_df[work_df['date'] == target_date].copy()

        # 3. 予測実行
        # NaNが含まれるとエラーになる場合があるので、fillna等が必要な場合はここで行う
        # current_data = current_data.fillna(0) # 必要に応じて

        pred_log = model.predict(current_data[features_1], current_data[features_2])
        pred_sales = np.expm1(pred_log)
        pred_sales = np.clip(pred_sales, 0, None) # 負の値は0に

        # 4. 予測値を保存
        current_data['pred_sales'] = pred_sales
        predictions.append(current_data[['id', 'date', 'store_nbr', 'family', 'pred_sales']])

        # 5. 【重要】予測値をwork_dfのsales列に書き戻す（翌日のラグ計算用）
        work_df.loc[work_df['date'] == target_date, 'sales'] = pred_sales

    # 結果を結合
    pred_result = pd.concat(predictions)
    return pred_result

from sklearn.linear_model import LinearRegression
from xgboost import XGBRegressor

def run_model_pipeline1(train_df, target_df, features_1, features_2):

    y_train = np.log1p(train_df['sales'])
    X_train_1 = train_df[features_1]
    X_train_2 = train_df[features_2]

    model_linear = LinearRegression()
    model_xgb = XGBRegressor(
        n_estimators=1000,
        learning_rate=0.05,
        max_depth=8,
        subsample=0.8,
        colsample_bytree=0.8,
        random_state=42,
        n_jobs=-1,
        # early_stopping_rounds=50 if mode == 'validation' else None
    )

    model = BoostedHybrid(model_linear, model_xgb)

    model.fit(X_train_1, X_train_2, y_train)

    history_start_date = train_df['date'].max() - pd.Timedelta(days=380)
    history_df = train_df[train_df['date'] >= history_start_date].copy()

    pred_result = recursive_predict(model, history_df, target_df, features_1, features_2)

    return model, pred_result

学習データをtrainとvalidationに分割して、モデルの学習と評価を行ってみます。

プロット用の関数

def plot_pred(train_df, pred_df, actual_df=None, lookback_days=180):
    last_train_date = train_df['date'].max()
    start_plot_date = last_train_date - pd.Timedelta(days=lookback_days)
    train_subset_plot = train_df[train_df['date'] >= start_plot_date].copy()

    daily_train_actual = train_subset_plot.groupby('date')['sales'].sum()
    if actual_df is not None:
        daily_target_actual = actual_df.groupby('date')['sales'].sum()
    daily_target_pred = pred_df.groupby('date')['pred_sales'].sum()

    plt.figure(figsize=(15, 6))

    plt.plot(daily_train_actual.index, daily_train_actual, label='Train Actual', color='gray', alpha=0.5)
    if actual_df is not None:
        plt.plot(daily_target_actual.index, daily_target_actual, label='Valid Actual', color='tab:blue')
    plt.plot(daily_target_pred.index, daily_target_pred, label='Forecast', color='red', linestyle='--')

    plt.title(f"Train History & Forecast Check (RMSLE: {rmsle:.4f})")
    plt.xlabel("Date")
    plt.ylabel("Sales")
    plt.legend()
    plt.grid(True, alpha=0.5)

    plt.axvline(x=last_train_date, color='black', linestyle=':', alpha=0.5, label='Split Date')

from sklearn.metrics import mean_squared_log_error

split_date = '2017-08-01'

train_cv = train_data[train_data['date'] < split_date].copy()
valid_cv = train_data[train_data['date'] >= split_date].copy()

print(f"Train: ~ {train_cv.date.max()}, Valid: {valid_cv.date.min()} ~")

val_model, val_pred_result = run_model_pipeline1(
    train_df=train_cv,
    target_df=valid_cv,
    features_1=dp_features,
    features_2=xgb_features,
)

val_merged = valid_cv[['id', 'date', 'sales']].merge(
    val_pred_result[['id', 'pred_sales']],
    on='id',
    how='left'
)
rmsle = np.sqrt(mean_squared_log_error(val_merged['sales'], val_merged['pred_sales']))

plot_pred(train_cv, val_pred_result, valid_cv)

青色がvalidationデータの実績値、赤色が予測値です。
結果の図を見ると、予測値が実績値に比べて非常に高い値を示してしまっていることがわかります。
これは、trend成分を過剰評価してしまっている可能性があります。

STL分解の結果を再度確認すると、trend成分は右肩上がりだが、年が経つにれて成長率が鈍化しているように見えます。
また、recursive戦略を採用しているため、予測誤差が累積してしまっている可能性もあります。

STL分解の結果

そのため、学習に使うデータを、2017年以降に変更なども検討してみましたが、あまり改善しませんでした。

train_mask = (all_df_with_feature['date'] >= '2017-01-01') & (all_df_with_feature['is_train'] == True)

XGBoost

先ほどの反省から、trend成分をlinear回帰モデルで予測するのではなく、XGBoost単体で予測を行ってみます。
また、マルチステップ予測の戦略も、recursiveではなくDirRec戦略を採用してみます。

関数のコード例

def create_multistep_dataset(df, horizon=16):
    df_wide = df.copy()
    df_wide = df_wide.sort_values(['store_nbr', 'family', 'date'])
    target_cols = []

    # 1日後〜16日後のsalesを列として追加
    for i in range(1, horizon + 1):
        col_name = f'target_{i}'
        df_wide[col_name] = df_wide.groupby(['store_nbr', 'family'])['sales'].shift(-i)
        target_cols.append(col_name)

    # 未来の正解がない行（直近16日分）は学習に使えないので削除
    df_wide = df_wide.dropna(subset=target_cols)
    return df_wide, target_cols

from sklearn.multioutput import RegressorChain

def run_model_pipeline2(X_train, y_train_log):
    model = RegressorChain(
        base_estimator=XGBRegressor(
            n_estimators=500,
            learning_rate=0.05,
            max_depth=6,
            subsample=0.8,
            colsample_bytree=0.8,
            n_jobs=-1,
            random_state=42
        ),
        order=None # Noneなら target_1 -> target_2 -> ... の順で連鎖
    )
    model.fit(X_train, y_train_log)

    return model

def plot_pred2(train_data, split_date, pred_matrix, y_true_matrix=None, history_days=90):
    plt.figure(figsize=(15, 6))

    last_train_date = pd.to_datetime(split_date)
    history_mask = (train_data['date'] < split_date) & (train_data['date'] >= last_train_date - pd.Timedelta(days=history_days))
    history_data = train_data[history_mask].groupby('date')['sales'].sum()

    start_pred_date = last_train_date + pd.Timedelta(days=1)
    valid_dates = pd.date_range(start=start_pred_date, periods=16, freq='D')
    daily_pred = pd.Series(pred_matrix.sum(axis=0), index=valid_dates)

    plt.plot(history_data.index, history_data, label='History (Train)', color='gray', alpha=0.5)
    plt.plot(daily_pred.index, daily_pred, label='Forecast', color='red', linestyle='--')

    if y_true_matrix is not None:
        daily_true = pd.Series(y_true_matrix.sum(axis=0), index=valid_dates)
        plt.plot(daily_true.index, daily_true, label='Actual', color='tab:blue')

    plt.title(f"Pure XGBoost with RegressorChain (RMSLE: {rmsle:.4f})")
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

同様に、学習データをtrainとvalidationに分割して、モデルの学習と評価を行ってみます。

train_wide, target_cols = create_multistep_dataset(train_data, horizon=16)

features = xgb_features + dp_features

split_date = '2017-07-30'

X_train = train_wide[train_wide['date'] < split_date][features]
y_train = train_wide[train_wide['date'] < split_date][target_cols]
y_train_log = np.log1p(y_train)

X_valid = train_wide[train_wide['date'] == split_date][features]
y_valid_true = train_wide[train_wide['date'] == split_date][target_cols]

val_model2 = run_model_pipeline2(X_train, y_train_log)

pred_log = val_model2.predict(X_valid)
pred_val_sales2 = np.clip(np.expm1(pred_log), 0, None)

y_true_flat = y_valid_true.values.flatten()
y_pred_flat = pred_val_sales2.flatten()

rmsle = np.sqrt(mean_squared_log_error(y_true_flat, y_pred_flat))

plot_pred2(train_data, split_date, pred_val_sales2, y_valid_true.values)

先ほどに比べて、予測値が実績値に近づいていることがわかります。

まだ、ノイズが大きいので、さらなる改善が必要ですが、XGBoost単体で予測を行う方が良さそうです。

ノイズを抑えるために、seed値を変えて学習して平均を取ってみました。

関数のコード例

from sklearn.multioutput import RegressorChain

def run_model_pipeline2(X_train, y_train_log, X_test_features=None, n_seeds=3):
    models = []
    preds_log_list = []

    seeds = list(range(1, n_seeds + 1))

    print(f"Training {n_seeds} models for ensembling...")

    for seed in seeds:
        print(f"  - Training with seed {seed}...")
        model = RegressorChain(
            base_estimator=XGBRegressor(
                n_estimators=500,
                learning_rate=0.05,
                max_depth=6,
                subsample=0.8,       # データを間引いて多様性を出す
                colsample_bytree=0.8,
                n_jobs=-1,
                random_state=seed    # ここを変える！
            ),
            order=None
        )
        model.fit(X_train, y_train_log)
        models.append(model)

        # もしテストデータが渡されていたら予測しておく
        if X_test_features is not None:
            pred = model.predict(X_test_features)
            preds_log_list.append(pred)

    # テスト予測がある場合は平均をとって返す
    if X_test_features is not None:
        # 平均をとる (Ensemble)
        avg_pred_log = np.mean(preds_log_list, axis=0)
        return models, avg_pred_log

    return model, None

先ほどよりもノイズが抑えられ、予測値が実績値に近づいていることがわかります。（RMSLEの値が改善しています。）

最後のモデルで作成した予測値を用いて、提出ファイルを作成します。

提出ファイル作成のコード例

# --- 3. 提出用フォーマットへの変換 (Reshape) ---
# 予測結果は (店舗x商品数, 16日間) の行列になっています。
# これを「縦持ち」に変換して、submissionファイルのIDと紐付けます。

# A. 行列をDataFrame化（列名を日付にする）
# 予測期間: 2017-08-16 〜 2017-08-31
predict_dates = pd.date_range(start='2017-08-16', periods=16, freq='D')
pred_df = pd.DataFrame(pred_sales2, columns=predict_dates)

# B. キー情報（store_nbr, family）を付与
pred_df['store_nbr'] = X_test['store_nbr'].values
pred_df['family'] = X_test['family'].values

# C. 横持ちから縦持ちへ変換 (Melt)
pred_melt = pred_df.melt(
    id_vars=['store_nbr', 'family'],
    var_name='date',
    value_name='sales'
)
pred_melt['date'] = pd.to_datetime(pred_melt['date'])

# D. 提出用ファイル(test.csv)の読み込みと結合
# kaggleのtest.csvには 'id' があるので、それをdate, store, familyで紐付けます
# test変数が読み込まれていない場合は、pd.read_csv('test.csv') などをしてください
test_data_base = test.copy() # 元のtestデータ
test_data_base['date'] = pd.to_datetime(test_data_base['date'])

submission = test_data_base.merge(
    pred_melt,
    on=['date', 'store_nbr', 'family'],
    how='left'
)

# --- 4. 保存 ---
submission = submission[['id', 'sales']]
submission.to_csv('submission2.csv', index=False)

print("Submission file saved: submission2.csv")
print(submission.head())

結果の値は、0.44572となりました。まずまずの結果ですね。

まとめ

時系列データ予測において、分析を行いながら特徴量を設計してモデルも作成してみました。
より精度を上げるには、以下のことをやってみると良いかもしれません。

ハイパラサーチ
NNモデルの検討（LSTM, Transformerなど）
外部データの活用（災害情報など、特別な時期は除いて学習するなど）
familyごとにモデルを分ける
その他特徴量エンジニアリング（例: lag特徴量の拡充、相互作用特徴量など）

参考にさせていただいた記事・資料