はじめに

弊社オークファンでは国内外のオークション・ショッピングサイトの過去の商品データを大量に保有しており、それらによる相場検索サービスを提供しています。
この記事では、この商品データの時系列分析をする想定で、次のような流れで実験をしてみます。

(1) 真の経時的な振る舞いが分かっている架空の商品データを生成
(2) 商品データを Prophet ライブラリで分析
(3) 真の振る舞いを正しく推測できているか確認

実行環境・バージョン

Windows
Python 3.9
Jupyter Notebook
Prophet 1.1.1
pandas 1.5.2
Matplotlib 3.6.2
NumPy 1.23.5

Prophet に関する背景知識

実験の前に、Prophet について簡単に説明します。

時系列分析と Prophet の関係

時系列分析における Prophet の位置づけを説明するために、時系列分析の手法を次の2つの切り口で分類してみます。

数理モデル ×（ツール・ライブラリ等の）実装

数理モデルの例としては以下があります。(参考)

古典的手法
- 指数平滑化法 etc.
ARIMA系
- ARIMA, SARIMA, SARIMAX etc.
状態空間モデル
機械学習系
- Prophet
- ノンパラメトリック
- ニューラルネットワーク
- etc...

実装方法の例としては以下があります。(参考)

クラウドサービス、製品等
- Amazon forecast, etc.
ライブラリ
- R言語
- Python
  - Darts
  - Prophet
  - Pycaret
  - Sktime
  - etc...

今回は数理モデルとして Prophet を扱います。
実装はMetaが公開しているライブラリの Prophet を利用します。これは名前の通り、数理モデル Prophet に特化した時系列分析ライブラリです。

数理モデル Prophet

数理モデルとしての Prophet について簡単に説明していきます。
詳しく知りたい方は、論文をご参照ください。（公式サイトからもリンクされているものです）

Prophet では、時系列データが次のような成分の和であるという仮定を置きます。

(時系列データ) = (トレンド成分) + (季節性成分) + (イベント成分) + (ノイズ成分)

イメージとしては、データを下の画像のように分解していきます。（公式ガイドより引用）

以下、各成分を簡単に見ていきます。

トレンド

折れ線グラフのように、各区間で線形に増減する成分です。
ハイパーパラメータとして、区間の分け目の数などをこちらで指定することもできます。

季節性

時系列データは１週間ごとや１年ごとの周期があることも多いため、それを表す成分です。
30日ごとなど、任意の周期を指定することもできます。

Prophetではフーリエ級数でモデリングしています。
フーリエ級数は、周期的なデータをsin波の足し算で近似するのに使われる手法です。（Wikipedia）

イベント

祝日など、周期性のない突発的な影響の成分です。

Prophet では、イベントの日にちを指定することで、イベント成分の大きさを推定することができます。

ノイズ

上記３成分の和と、実際の値との間の誤差です。

商品データとProphetの各成分の対応

前節の各成分のイメージが湧きやすいように、商品データに当てはめて考えてみます。

トレンド

例えば新商品の販売数は、次のようなトレンドとしてモデル化できそうです。

発売日までは販売数0件
発売日から販売数が増加
時間が経つと増加量が鈍化したり減少に転じたりする

季節性

例えばかき氷機は夏によく売れる、スノーボードは冬によく売れる、というように、売れ行きが季節により大きく異なる商品があることが知られています。
こうした増減が12ヶ月周期でデータに現れることになります。

イベント

ネットオークションは土日祝日に取引数が多くなる性質があります。
土日は季節性成分としても扱えますが、祝日は周期性がないためイベントとしてモデル化することになります。

サンプルデータの説明

ここから今回の実験の話に入ります。

まずは分析対象として、架空の商品のオークション落札数データを生成します。
条件は下にまとめています。

結果のグラフを見れば何となくわかると思うので、この節は必要に応じて読み返すと良いかもしれません。

条件の詳細

対象日

毎月1日

期間

2年
3年
5年
10年

落札件数

(落札件数) = (トレンド成分) + (季節性成分) + (ノイズ成分) とします。

トレンド

3年目の12カ月間 (25カ月目~36カ月目) に線形増加
変化量は下記のいずれか
- 500 件/月
- 300 件/月
- 100 件/月
- 0 件/月

季節性

夏に売れる商品を仮定しました。各月の落札件数は以下の通りです。

1月	2月	3月	4月	5月	6月	7月	8月	9月	10月	11月	12月
1000	1000	1000	1000	1000	2000	6000	8000	2000	1000	1000	1000

(イベント)

今回は時間の都合で割愛します。

ノイズ

今回は正規分布(連続分布、上限下限なし)の近似値を使うことにします。（参考: 正規分布とは）
- 近似の内容
  - 整数に丸める
  - 落札件数がマイナスにならないように切り上げる
- パラメータ
  - 平均 0
  - 標準偏差 100

実装

実装を始めます。
まずはライブラリのインストールからです。

インストール

pandas, Prophet をインストールします。

pip install pandas
pip install prophet

ライブラリの読み込み

NumPy, Matplotlib, pandas, Prophet を読み込みます。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from prophet import Prophet

データ生成

次のようなコードからデータを生成します。

n_year = (年数)
std = (ノイズの標準偏差)
trend_start = 24
trend_end = 36
trend_inc = (トレンドの増加量)
data_1y = (1年分データのNumPy配列)

data_all_year_no_noise = np.tile(data_1y, n_year)
data = add_noise(data_all_year_no_noise, std)
data = add_trend(data, trend_start, trend_end, trend_inc)

各成分は下記のように実装しています。

トレンド

def add_trend(a, i_start, i_end, increase):
    increases = np.zeros(a.shape)
    increases[i_start:i_end] = np.arange(1, i_end - i_start + 1) * increase
    increases[i_end:] = increases[i_end-1]
    return a + increases

季節性

data_1y = np.array([1000, 1000, 1000,
                    1000, 1000, 2000,
                    6000, 8000, 2000,
                    1000, 1000, 1000,])

ノイズ

def generate_noise(limit, std, ndim):
    noise = np.random.normal(loc=0, scale=std, size=ndim)
    noise[noise < -limit] = -limit
    return np.around(noise)

def add_noise(a, std):
    return a + generate_noise(limit=np.min(a), std=std, ndim=a.shape)

予測と可視化

前節のデータ生成から、データを Prophet モデルに当てはめて、予測値をグラフにするまでを次の関数で実装しました。

def predict_and_plot(data_1y, n_year, std, trend_start=None, trend_end=None, trend_inc=None):
    # データ生成
    data_all_year_no_noise = np.tile(data_1y, n_year)
    data = add_noise(data_all_year_no_noise, std)

    if trend_start and trend_end and trend_inc:
        data = add_trend(data, trend_start, trend_end, trend_inc)

    # データの変換
    df = pd.DataFrame({"ds": make_ds(n_year=n_year),
                       "y": data})

    # モデルの適用と予測
    m = Prophet()
    m.fit(df)
    predict_months = 24
    future = m.make_future_dataframe(periods=predict_months, freq="MS")
    forecast = m.predict(future)

    # 可視化
    fig1 = m.plot(forecast)
    fig2 = m.plot_components(forecast)
    # 毎月1日のみの季節性
    plt.figure()
    fig, ax = plt.subplots()
    ax.plot(np.arange(1, 13), forecast["yearly"][0:12:])
    ax.set_xlabel("month")
    ax.set_ylabel("yearly seasonality of month")

    return forecast

以下、各パートの解説です。（データの生成は前節参照）

データの変換

Prophet が受け付けるのは次の2列を持った pandas の DataFrame となります。

ds (datestamp): YYYY-MM-DD 又は YYYY-MM-DD HH:MM:SS
y : 予測対象となる値

先ほどのデータからこの DataFrame を用意します。

def make_ds(n_year):
    mm_dd = [f"{i+1:02}" + "-01" for i in range(12)]
    yyyy = [str(i+2010) for i in range(n_year)]
    return [y + "-" + md for y in yyyy for md in mm_dd]

df = pd.DataFrame({"ds": make_ds(n_year=n_year),
                   "y": data})

モデルの適用と予測

Prophetに DataFrame を渡して、学習を実行します。

m = Prophet()
m.fit(df)
predict_months = 24
future = m.make_future_dataframe(periods=predict_months, freq="MS")
forecast = m.predict(future)

今回は学習データが毎月1日なので、予測も毎月1日のみとするために make_future_dataframe() に freq="MS" を渡しています。

可視化

学習ができたので、予測値とトレンド、季節性成分を可視化してみます。
Prophetに可視化用のメソッドがすでに用意されています。（便利ですね）

fig1 = m.plot(forecast)
fig2 = m.plot_components(forecast)

ただし、季節性成分はデフォルトだと毎月1日以外のデータも予測を出すため、今回の要件に合いません。
毎月1日の季節性成分は下記のように可視化します。

plt.figure()
fig, ax = plt.subplots()
ax.plot(np.arange(1, 13), forecast["yearly"][0:12:])
ax.set_xlabel("month")
ax.set_ylabel("yearly seasonality of month")

結果

ここから各条件でシミュレーションしていきます。

シンプルな条件の場合

まずは5年分、トレンドなしというシンプルな条件で結果を見てみます。

predict_and_plot(data_1y=data_1y, n_year=5, std=100)

3枚のグラフが出力されるので、順にみていきます。

点が実データ、線が過去の推定値(ノイズなしの値)と将来の予測値です。
薄い青色で塗りつぶされているのは、推定や予測の誤差の大きさの目安です。（デフォルトはMCMCで得た推定値サンプルの+-80パーセンタイルのようです）
大まかには予測できていそうです。

トレンドは上下の幅が300件程度(正解は0件)になり、少しダウントレンドの予測になってしまいました。

一方で、季節性成分は年間でかなり乱れているように見えますが、これは毎月1日以外のデータも予測を出しているためです。

毎月1日のみのグラフを見てみます。

トレンド(約2000件)を足すとほぼ真の値(サンプルデータの説明の節を参照)と一致していることが分かります。

過去データの年数を変化させた場合

データの年数を変化させて、予測への影響を見てみます。

predict_and_plot(data_1y=data_1y, n_year=10, std=100)
predict_and_plot(data_1y=data_1y, n_year=3, std=100)
predict_and_plot(data_1y=data_1y, n_year=2, std=100)

10年分のデータのグラフ

3年分のデータのグラフ

2年分のデータのグラフ

2年分のデータだと、 周期性のない大きなノイズの乗ったデータである という推定結果になってしまいました。
今回の条件では3年分以上のデータが必要になるようです。

トレンドを変化させた場合

トレンドを追加してみます。

predict_and_plot(data_1y=data_1y, n_year=5, std=100, trend_start=24, trend_end=36, trend_inc=500)
predict_and_plot(data_1y=data_1y, n_year=5, std=100, trend_start=24, trend_end=36, trend_inc=300)
predict_and_plot(data_1y=data_1y, n_year=5, std=100, trend_start=24, trend_end=36, trend_inc=100)