初めに：
このブログはAidemy Premiumのカリキュラムの一環で、受講修了条件を満たすために公開しています。

実行環境：
Python3.10.12
Windows10
Google Colaboratory

目的：
時系列データ解析の実践
SARIMAの実践と予測の精度を高める

概要：
まずは、SARIMAモデルに使用するパラメーターsを決定するために、時系列解析を実践して、データの周期性やトレンドを確認する。
原系列に外れ値があるため、原系列の外れ値を欠損させ、それをスプライン補完し、自然な系列に近づけたデータを作成し、原系列と補完後データのSARIMAを比較してみる。

工程：

必要なモジュールのインポートとデータを読み込み
原系列の時系列解析
2-1. 自己相関係数と偏自己相関係数グラフで周期性の確認
2-2. トレンド、季節変動、不規則変動(残差)に分け可視化
原系列の外れ値を欠損値にし、スプライン補完
スプライン補完後データの時系列解析
4-1. 自己相関係数と偏自己相関係数グラフで周期性の確認
4-2. トレンド、季節変動、不規則変動(残差)に分け可視化
原系列とスプライン補完後のSARIMAを比較

Appendix： SARIMAの最小のBICの特定と、すべてのパラメーターのリストアップの自動化

1. 必要なモジュールのインポートとデータを読み込み

#必要なモジュールのインポート
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import numpy as np
import warnings
import itertools
from sklearn.neighbors import LocalOutlierFactor
from scipy import interpolate

#Google Driveをマウント
from google.colab import drive
drive.mount('/content/drive')

#データを読み込み、日付の列をindex化する
df = pd.read_excel('/content/drive/MyDrive/Colab Notebooks/DATA/製品A時系列解析データ.xlsx',sheet_name='製品A時系列解析データ')
df.index=df.CreatedDate
df=df[['TotalQuantity']]
print(df)
df.plot()

元データはこのように月別製品個数の36か月分のデータ。例えば2021年12月が5000個台となっており、四分位範囲上は外れ値のあるデータである。

2. 原系列の時系列解析
2-1 自己相関係数と偏自己相関係数グラフで周期性の確認

「パラメーター s に関しては事前に時系列データや次に説明する偏自己相関の可視化を行うことによって調べておく」
とのことで、SARIMAで使用するパラメーターsを特定するために周期性を確認する。

# グラフの枠をセット
fig=plt.figure(figsize=(12, 8))

# 自己相関係数のグラフを出力します
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df, lags=35, ax=ax1)

# 偏自己相関係数のグラフを出力します
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df, lags=17, ax=ax2)
plt.show()

原系列のコレログラム：

コレログラムにおいて、水色の領域を超える自己相関の値は、偶然の範囲を超えてデータ間に相関があることを示しており、この範囲を超えない自己相関の値は、偶然によるものと考えられる。
偏自己相関の方で3と6が比較的数値が高いが、すべて水色の範囲に収まるので、偶然によるものと考えられ、自己相関はなさそう。周期性がはっきりとは特定できないデータ。

2-2 トレンド、季節変動、不規則変動(残差)に分け可視化

#tsa.seasonal_decompose()を用い、原系列をトレンド、季節変動、残差に分けて出力
fig = sm.tsa.seasonal_decompose(df, period=8).plot()
plt.show()

period=8とした。周期性がわかりにくいため自己相関からは特定できず、結局パラメーターを変えながらSARIMAに当てはめて可視化を繰り返して、過去との類似度が高い8を選択した。再現性上精度が高そうな期間だったため。

出力結果：

原系列＝トレンド + 季節性 + 残差
Trend：なにか傾向があれば右肩上がりなどの傾向を示すとのこと。どちらがというとアップダウンがあり、直近は低下傾向。
Seasonal：2～3か月おきに低下を繰り返すパターンということだろうか。period=8でよいのか疑問。
Resid：これが定常過程とのこと。2022年以降は比較的0のあたりに集中しており何等かの傾向がありそう。グラフ内のゼロはResidualの平均だそう。±5000のグラフでバラツキが大きいと感じる。

3. 原系列の外れ値を欠損値にし、スプライン補完

#　補完後用のdf2と列を作成
df2 = df.copy()
df2['TotalQuantity_NaN']=df2['TotalQuantity']

#　Scikit-learnによるLOFの実装
lof = LocalOutlierFactor(n_neighbors=7)
outliers = lof.fit_predict(df2[['TotalQuantity_NaN']])
df2['Outlier'] = outliers

# 外れ値は-1とし欠損値にする
df2.loc[df2['Outlier'] == -1, 'TotalQuantity_NaN'] = np.nan

# スプライン補完
df3 = df2.dropna()
tck = interpolate.splrep(df3.index, df3['TotalQuantity_NaN'], s=0)
df2['InterpolatedQty'] = interpolate.splev(df2.index, tck, der=0)

print(df2)
df2['TotalQuantity'].plot()
df2['InterpolatedQty'].plot()

出力結果：
このように外れ値は-1と表示され、その月はNaNとなり、そこが補完されたものがInterpolatedQty列。

元データとスプライン補完後データのグラフ：
青色が原系列で、外れ値となった月のデータをなだらかに補完してくれている。

4. スプライン補完後データの時系列解析
4-1 自己相関係数と偏自己相関係数グラフで周期性の確認

#SARIMA用に補完データInterpolatedQty列以外をdrop
df2.drop(columns=['TotalQuantity', 'TotalQuantity_NaN', 'Outlier'], inplace=True)

# 自己相関用グラフの枠のセット
fig=plt.figure(figsize=(12, 8))

# 自己相関係数のコレログラフを出力
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df2, lags=35, ax=ax1)

# 偏自己相関係数のコレログラフを出力
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df2, lags=13, ax=ax2)
plt.show()

スプライン補完後のコレログラム：

補完後でも周期性がわかりにくくsは特定できず、こちらもパラメーターを変えながらSARIMAに当てはめて可視化を繰り返して、精度が高そうな期間のperiod=10を使用することにした。偏自己相関係コレログラフで、1の次に増加しているのが9番目という点も参考にした。

4-2 トレンド、季節変動、不規則変動(残差)に分け可視化

#tsa.seasonal_decompose()を用い、補完後データをトレンド、季節変動、残差に分けて出力
fig = sm.tsa.seasonal_decompose(df2, period=10).plot()
plt.show()

出力結果：

Trendはさきほどよりはっきりと低下傾向に。
Seasonalも10カ月に一回のバンプがある周期性があるようにみえる。
Residは原型列のときよりバラつき大きくなってしまったかと思ったが、y軸の幅が±5000から±1000になっておりバラつき低下。また、residが2022-05あたりで増加しているのは思い当たるフシがある。

補完したデータ方が予測に向いていそうなデータにみえる。

5. 原系列とスプライン補完後のSARIMAを比較
SARIMAモデルのパラメーター、(p, d, q) (sp, sd, sq, s)をBIC(ベイズ情報量基準) を使って最適化するための関数selectparameterの定義

def selectparameter(DATA, s):
    p = d = q = range(0, 2)
    pdq = list(itertools.product(p, d, q))
    seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
    parameters = []
    BICs = np.array([])
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(DATA,
                                            order=param,
                                            seasonal_order=param_seasonal)
                results = mod.fit()
                parameters.append([param, param_seasonal, results.bic])
                BICs = np.append(BICs,results.bic)
            except:
                continue
    return parameters[np.argmin(BICs)]

#さきほど決定したsで、原系列、補完後ともに関数に当てはめパラメーター(p, d, q) (sp, sd, sq, s)を特定
selectparameter(df, 8)
selectparameter(df2, 10)

# SARIMAに当てはめ
SARIMA_df = sm.tsa.statespace.SARIMAX(df,order=(0, 1, 1),seasonal_order=(0, 1, 1, 8)).fit()
SARIMA_df2 = sm.tsa.statespace.SARIMAX(df2,order=(1, 1, 0),seasonal_order=(0, 1, 1, 10)).fit()

# predに予測データを代入する
SARIMA_pred = SARIMA_df.predict("2023-3-1", "2024-11-1")
SARIMA_pred2 = SARIMA_df2.predict("2023-3-1", "2024-11-1")

# predデータと元の系列データの可視化
plt.plot(df)
plt.plot(SARIMA_pred, color="g")
plt.plot(SARIMA_pred2, color="r")
plt.show()

出力結果：

青＝原系列
赤＝スプライン補完後のSARIMA
緑＝原系列のSARIMA

原系列と期間の重なりのある部分は、再現度が補完後SARIMAの方があるようにみえたのだが、その後の期間についてはここまで低下するのだろうか...という印象。周期やトレンドを忠実に反映しているようにみえるものの。原系列SARIMAは再現度は低く、逆に増加傾向になっている。

Appendix： SARIMAの最小のBICの特定と、すべてのパラメーターのリストアップを自動化
BICが小さければ小さいほど、sの予測精度が高いとのことで、一番小さいBICの特定、すべてのパラメーター抽出、各BICのリストアップを一気に出すコード。

param_list = []
for s in range(2, 25):
    param = selectparameter(df, s)[2]
    param_list.append(param)
    print(selectparameter(df, s))
print(min(param_list))

上記を走らせると、時間はかかるが、一番小さいBICとs=2～24の各パラメーターを自動的にリストアップし、warningを除外すると下記が残る。

出力結果：
429.1211688
[(0, 1, 0), (0, 1, 1, 11), 473.98611471143175]
[(0, 1, 1), (0, 1, 1, 12), 429.1211687894908]
[(0, 1, 1), (0, 1, 1, 4), 610.2703531661796]
[(0, 1, 1), (0, 1, 1, 5), 586.7788574114014]
[(0, 1, 1), (0, 1, 1, 6), 573.015323245988]
...
周期性がわかりづらいデータだと自己相関からsを特定しづらく、ひとつひとつパラメーターを走らせては可視化を繰り返さなければならないので、リストアップしておくと楽だった。

解決できたことやわかったこと：

スプライン補完後のデータの方が分散が低下し、予測に使いやすそう。時系列解析時のデータ前処理に取り入れていきたい。
BICは小さければ小さい程よいと思っていたが、dfのBICはS=23のBIC=4が一番小さかったが、実際SARIMAに当てはめると非現実的なグラフに。元データが36か月と少ないためsが大きい値は適さない様子で、必ずしもBICが小さければ最適というわけではなかった。
周期性、コレログラムの理解が少し進んだ。

考察・今後の課題：

予測方法について

「SARIMAモデルとはARIMAモデルをさらに季節周期を持つ時系列データにも拡張できるようにしたモデル」なので、周期性の少ないデータはSARIMAでなくARIMAの方が適切であれば、今後試してみたい。
データ数が少ないことも今回は特定しづらい要因だった可能性があるそうだが、そういうケースの機械的予測方法があるのか。
売上や受注個数予測は、気象や火山などの何十年にも渡る周期性を持っていそうなデータとは性質が違い、周期性があったとしても、その年の外的要因やビジネスプランなどにより一定ではないケースが多そう。調査における実験法と観測法の分類のように、気象などは観察法で周期性を重視するが、ビジネスは実験法なので、明らかにわかっている説明変数は取り入れて予測する方が精度が上がるかもしれない。
SARIMAは、グラフで、原系列と予測の重なった部分の再現度をみながら精度の高さを人間が判断するので、それが精度が高いかどうかを機械的数値的に判断する方法。

周期性について

tsa.seasonal_decomposeもperiodの設定で出力結果が結構変わる。周期性の少ないデータの周期の特定の仕方。
BIC=tsa.seasonal_decomposeの周期性は、最適なものであれば、基本的には一致するのか。
periodは小数点の使用も可能なのか。

補完方法について

スプライン補完以外にも補完方法があるようなので調べる。

外れ値の定義について

LocalOutlierFactorの外れ値の定義をn_neighbors=7の選択について。
他にも外れ値の定義の仕方があるようなので、他も調べて実践してみたい。