時系列モデルを勉強していくぜ！～データの前処理編～（Signate引っ越し需要予測17位解法）

はじめに

時系列モデルってみなさん得意ですか？自分は食わず嫌いで避けてきましたが、（自称）データサイエンティストとしてレベルアップするために時系列データの勉強をしています。
勉強した成果をアウトプットとして記事にしておきます。全４回（予定）に渡ってお送りしますので、よろしければお付き合いください。

第１回（データの前処理）
第２回（LightGBMモデルの構築、Pycaretによるアンサンブルモデルの構築）
第３回（Prophetによる時系列モデルの構築、アンサンブル学習）

題材

SignateのSOTA(State-of-the-Art) Challengeから「アップル 引越し需要予測」を選びました。
テーマとしては分かりやすいですし、予想も立てやすい（3～4月が繁忙期だろう等）からです。
また、データも訓練データとテストデータのみなので、とっつきやすいと思います。

SOTA(State-of-the-Art) Challengeとは？

・過去に終了済みのコンペに再び参加できる
・精度向上を極限まで挑戦し、知見を共有することを目的とした終わりのないコンペ

結論

先にモデルの構築の流れと、結果をお伝えしておきます。

モデル構築

①前処理
②LightGBMによる回帰モデルの構築と特徴量の選択
③特徴量選択したデータをPycaretを用いて欠損値補完と回帰モデル構築（Extra Trees Regressor、Gradient Boosting Regressor、Random Forest Regressorのアンサンブル）
④欠損値補完したデータを用いてProphetによる時系列モデルの構築
⑤Pycaretで構築したアンサンブル回帰モデルとProphetで構築した時系列モデルをアンサンブル

結果

現状、17位/775人です！銀メダルは確実ですが、金メダルまでもう少し・・・！！

EDA

まずはデータを眺めることから始めます。
訓練データは2010-07-01から2016-03-31までのデータ、テストデータは2016-04-01から2017-03-31までのデータとなっています。説明変数はすべて数値型で、目的変数'y'（引っ越し数）はなぜかobject型でしたので、後ほどint型に変換することとします。

import pandas as pd
import numpy as np
from IPython.display import display
import matplotlib.pyplot as plt
import sweetviz as sv

train_df = pd.read_csv('**/train.csv')
test_df = pd.read_csv('**/test.csv')
display(train_df.head(), test_df.head())

実行結果：

datetime：日時（YYYY-MM-DD）
y：引越し数
client：法人が絡む特殊な引越し日フラグ
close：休業日
price_am：午前の料金区分（-1は欠損を表す。5が最も料金が高い）
price_pm：午後の料金区分（-1は欠損を表す。5が最も料金が高い）

省略しますが、以下のような基本情報も確認します。

display(train_df.shape, test_df.shape)

display(train_df.info(), test_df.info())

display(train_df.describe(), test_df.describe())

もう少し具体的にデータを見るため、SweetvizというEDAライブラリを使用します。Sweetvizでは訓練データとテストデータの比較ができるので両者の間でデータの偏りがないか確認できるのでとても便利です。
Sweetviz詳細は以下記事を参照ください。

それでは、Sweetvizでデータを見ていきます。sv.config_parser.read_string("[General]\nuse_cjk_font=1")を記述しておくことで、日本語が化けて豆腐にならずに済みます。今回はカラム名もデータも日本語は入っていないので必要ありませんが、紹介しておきます。

sv.config_parser.read_string("[General]\nuse_cjk_font=1")
report = sv.compare(train_df, test_df, target_feat='y')
report.show_html('hikkoshi_eda.html')

はい、便利～(ﾟ∀ﾟ)
例として以下のことが読み取れます。
・'client'（法人フラグ）はテストデータに偏り
・'client'（法人フラグ）が1だと、引っ越し数多い
・'close'（休業日）はもちろん引っ越し数は0、訓練/テストデータ間で偏りなし
・'price_am'、'price_pm'は両社とも高い方が引っ越し数多い（ビジネスなので、需要に比例して価格を上げるのは必然）
・'price_am'、'price_pm'の訓練/テストデータ間で'-1'（欠損値）に偏り有り。訓練データに欠損値が多い。
・相関係数を見ると、引っ越し数と価格の相関が高い

ここからはもう少し時系列らしいEDAをしていきます。以下の記事が参考になります。

成分分解をする前に、データに定常性があるか（P値（有意水準）が0.05（5%）以下かどうか）を確認します。

from statsmodels.tsa.stattools import adfuller

def adf_test(series):
    adf_df = pd.DataFrame(
        [
            adfuller(series)[1]
        ],
        columns=['P値']
    )
    adf_df['P値'] = adf_df['P値'].round(decimals=3).astype(str)
    print(adf_df)

adf_test(series=train_df['y'])

実行結果：
P値
0.575

うん、定常性はないですね。ちなみに差分を取ることで定常性を確保できましたが、今回は定常性を前提としたモデル（ARIMAやSARIMA）は使用しないので、差分は取らずに進めます（差分で進めると、予測結果も差分になって、スケールを戻す作業が発生してめんどくさい）。

まずは訓練データとテストデータをマージして、日付をindexにしてしまいましょう。
欠損値があるとこの後の分析ができなくなるので、テストデータの'y'はとりあえず0で埋めておきます。

train_df['is_train'] = 1
test_df['is_train'] = 0

train_cols = train_df.columns.tolist()
test_cols = test_df.columns.tolist()
common_cols = list(set(train_cols) & set(test_cols))
train_only = list(set(train_cols) - set(common_cols))
test_only = list(set(test_cols) - set(common_cols))

for col in train_only:
    test_df[col] = None

for col in test_only:
    train_df[col] = None

merged_df = pd.concat([train_df, test_df], ignore_index=True, sort=False)
merged_df.set_index('datetime', inplace=True)
merged_df.index = pd.to_datetime(merged_df.index)
merged_df['y'] = merged_df['y'].fillna(0).astype(int)
merged_df.head()

実行結果：

目的変数のトレンドや季節性を確認します。

from statsmodels.tsa.seasonal import seasonal_decompose

price_am_seasonal = seasonal_decompose(merged_df['y'], model='additive', period=365)

# 分解結果のプロット
fig = price_am_seasonal.plot()
plt.show()

実行結果：

2016年4月以降はテストデータなので、無視するとして、トレンドを見るとなんとなく引っ越し数は年々上昇傾向にあるみたいですね。定住する人（家を買う人）が減っているのですかね。
一年周期の季節変動も見て取れます。

自己相関も確認しておきましょう。

from statsmodels.graphics.tsaplots import plot_acf

plot_acf(train_df['y'], lags=730)
plt.show()

実行結果：

一年周期の季節変動があるので、1年前との自己相関は高いですね。1年前のラグ特徴量を追加するのも良さそうです（そもそもテストデータが1年間のデータなので、目的変数のラグ特徴量は1年以上前からしか取れません）。

他の特徴量も同様に成分分解と自己相関を確認しました。
・'price_am'、'price_pm'はおおよそ7日前までと1年前の自己相関係数高め
・'price_am'、'price_pm'にトレンドは無いが、季節性の周期はある
・'close'は2日前までと1年前の自己相関係数高め

また、'price_am'と'price_pm'には2010-07-01から2011-01-03までがすべて欠損値（'-1'表記）となっていました。目的変数の予測に重要な特徴量ですので、欠損値というのはよろしくないですね。しかし、連続して半年以上の欠損値なので、補完は困難と判断し、削除しました。

merged_df = merged_df.loc[(merged_df.index< pd.Timestamp('2010-07-01')) |
            (merged_df.index >= pd.Timestamp('2011-01-01'))]

特徴量の作成

精度を見ながら特徴量を追加・削減しました。以下のすべて使ったわけではないですが、参考までに載せておきます。

ラグ特徴量、移動平均

# ラグ特徴量の追加
merged_df['y_lag_365'] = merged_df['y'].shift(365)

for lag in [1, 7, 90, 365]:
    merged_df[f'price_am_lag_{lag}'] = merged_df['price_am'].shift(lag)

for lag in [1, 7, 90, 365]:
    merged_df[f'price_pm_lag_{lag}'] = merged_df['price_pm'].shift(lag)

for lag in [1, 365]:
    merged_df[f'close_lag{lag}'] = merged_df['close'].shift(lag)

# 移動平均の追加
windows = [7, 30, 90]
for window in windows:
    merged_df[f'price_am_moving_avg_{window}'] = merged_df['price_am'].rolling(window=window).mean()

windows = [7, 30, 90]
for window in windows:
    merged_df[f'price_pm_moving_avg_{window}'] = merged_df['price_pm'].rolling(window=window).mean()

年月日

merged_df['year'] = list(pd.Series(merged_df.index).apply(lambda x: x.year))
merged_df['month'] = list(pd.Series(merged_df.index).apply(lambda x: x.month))
merged_df['day'] = list(pd.Series(merged_df.index).apply(lambda x: x.day))

曜日フラグ、週末フラグ等

思いつく限り追加します。

merged_df['day_of_week'] = merged_df.index.dayofweek
merged_df['week_of_year'] = merged_df.index.weekofyear
merged_df['week_of_month'] = ((merged_df['day'] - 1) // 7) + 1
merged_df["is_wknd"] = (merged_df['day_of_week'] // 4)
merged_df['is_month_start'] = merged_df.index.is_month_start.astype(int)
merged_df['is_month_end'] = merged_df.index.is_month_end.astype(int)
merged_df['quarter'] = merged_df.index.quarter
merged_df['is_quarter_start'] = merged_df.index.is_quarter_start.astype(int)
merged_df['is_quarter_end'] = merged_df.index.is_quarter_end.astype(int)
merged_df['is_year_start'] = merged_df.index.is_year_start.astype(int)
merged_df['is_year_end'] = merged_df.index.is_year_end.astype(int)

祝日フラグ

import holidays

jp_holidays = holidays.Japan()

def label_holidays(row, holiday_list):
    if row.name in holiday_list:
        return 1  # 祝日の場合は1を返す
    return 0  # 祝日でなければ0を返す

# 祝日の列をデータフレームに追加
merged_df['is_holiday'] = merged_df.apply(label_holidays, holiday_list=jp_holidays, axis=1)

季節フラグ

def get_season(month):
    if month in [3, 4, 5]:
        return 'Spring'
    elif month in [6, 7, 8]:
        return 'Summer'
    elif month in [9, 10, 11]:
        return 'Autumn'
    else:
        return 'Winter'

merged_df['season'] = merged_df['month'].apply(get_season)

六曜フラグ

下記githubを参照。

from qreki import Kyureki

def get_rokuyou(date):
    kyureki = Kyureki.from_date(date)
    return kyureki.rokuyou

# 六曜フラグを追加
merged_df['rokuyo'] = merged_df.index.map(get_rokuyou)

大型連休フラグ

def is_golden_week(date):
    return (date.month == 4 and date.day >= 29) or (date.month == 5 and date.day <= 5)

merged_df['golden_week'] = merged_df.index.map(is_golden_week).astype(int)

def is_obon(date):
    return date.month == 8 and date.day >= 13 and date.day <= 16

merged_df['obon'] = merged_df.index.map(is_obon).astype(int)

def is_new_year(date):
    return (date.month == 12 and date.day >= 29) or (date.month == 1 and date.day <= 3)

merged_df['new_year'] = merged_df.index.map(is_new_year).astype(int)

merged_df['golden_week'] = merged_df.index.map(is_golden_week).astype(int)
merged_df['obon'] = merged_df.index.map(is_obon).astype(int)
merged_df['new_year'] = merged_df.index.map(is_new_year).astype(int)

経済指標等

なんとなく経済指標も効いてきそうですよね。FREDから取得するのが簡単だと思います。
IDが良く変更されるので、FREDのサイトできちんと確認したほうが良いです。
こういう指標は月単位だったり4半期単位、年単位が多いので.ffil().bfill()
で前方補完と後方補完をしています。

import pandas_datareader.data as web

start_date = merged_df.index.min()
end_date = merged_df.index.max()

# 日本のGDP成長率 (Quarterly)
gdp_growth_japan = web.DataReader('JPNRGDPEXP', 'fred', start_date, '2017-04-01').resample('D').interpolate(method='time').ffill().bfill()

# 日本の失業率 (Monthly)
unemployment_rate_japan = web.DataReader('LRUN64TTJPM156S', 'fred', start_date, end_date).resample('D').interpolate(method='time').ffill().bfill()

# 日本の消費者物価指数 (CPI) (Monthly)
cpi_japan = web.DataReader('JPNCPIALLMINMEI', 'fred', start_date, end_date).resample('D').interpolate(method='time').ffill().bfill()

# 日本の不動産価格指数 (Quarterly)
real_estate_price_index_japan = web.DataReader('QJPN628BIS', 'fred', start_date, end_date).resample('D').interpolate(method='time').ffill().bfill()

# 賃金
mortage_lone_japan = web.DataReader('LCEAPR03JPM661S', 'fred', start_date, end_date).resample('D').interpolate(method='time').ffill().bfill()

# 消費者態度指数
cpi_houce = web.DataReader('JPNCPIHOUMINMEI', 'fred', start_date, end_date).resample('D').interpolate(method='time').ffill().bfill()

# 中央銀行金利
central_bank_rates = web.DataReader('IRSTCB01JPM156N', 'fred', start_date, end_date).resample('D').interpolate(method='time').ffill().bfill()


merged_df = merged_df.merge(gdp_growth_japan, left_index=True, right_index=True, how='left', suffixes=('', '_gdp_growth_japan'))
merged_df = merged_df.merge(unemployment_rate_japan, left_index=True, right_index=True, how='left', suffixes=('', '_unemployment_rate_japan'))
merged_df = merged_df.merge(cpi_japan, left_index=True, right_index=True, how='left', suffixes=('', '_cpi_japan'))
merged_df = merged_df.merge(real_estate_price_index_japan, left_index=True, right_index=True, how='left', suffixes=('', 'real_estate_price_index_japan'))
merged_df = merged_df.merge(mortage_lone_japan, left_index=True, right_index=True, how='left', suffixes=('', 'mortage_lone_japan')) 
merged_df = merged_df.merge(cpi_houce, left_index=True, right_index=True, how='left', suffixes=('', 'cpi_houce'))
merged_df = merged_df.merge(central_bank_rates, left_index=True, right_index=True, how='left', suffixes=('', 'central_bank_rates'))

訓練データとテストデータの分離

train_data = merged_df[merged_df['is_train'] == 1].copy()
test_data = merged_df[merged_df['is_train'] == 0].copy()

train_data.drop(['is_train'], axis=1, inplace=True)
test_data.drop(['y','is_train'], axis=1, inplace=True)

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

はい、データの準備ができました。

おわりに

今回はデータの前処理まで行いました。次回はLightGBMによる単純な回帰モデルを構築していきます。
まだまだ勉強中なので、もっとこんなデータみろよ！とか、こんな特徴量追加しろよ！とかございましたらコメントよろしくお願いします。

自分のような時系列初心者の参考に少しでもなれば幸いです。

それでは次の記事でお会いしましょう～('ω')ノ