時系列分析による時系列データの解析と未来予測（ARIMA, SARIMA）

この記事は個人ブログと同じ内容です

www.ritolab.com

時系列分析の基礎を確認しつつ、「データの確認・理解」「定常データへの変換」「モデル構築（ARIMA, SARIMA）」と一連の時系列分析の流れを実施し、時系列データの未来予測を行っていきます。

時系列分析とは

時系列分析は、時間的な順序で取られたデータ（＝時系列データ）の特性やパターンを分析し、過去の振る舞いから将来の振る舞いを予測するための統計的手法です。

時系列データとは

時系列データは、一定の時間間隔（日次、週次、月次など）で観測されたデータポイントから構成されるデータです。

このような時系列データは、経済・株価・気象・トラフィックデータなど、多くの実世界の現象を表現するために使用されます。

  • ある地域の毎日の気温
  • ある店舗の日次の売上

通常の時系列データでは観測者によって観測の時間間隔が設定されます。

時系列データの特徴

  • 時系列データは一度しか観測されない
  • 観測値から平均や分散などを推定することはできない
    • 時間の非独立性
      • 時系列データの観測値は時間的に依存しており、過去の値が現在の値に影響を与える可能性がある。したがって、単純な平均や分散の推定では、データの時間的なパターンや相関関係を考慮することができない。
    • 季節性や周期性
      • 時系列データには季節性や周期性が存在する場合がある。これらの要素は平均や分散に影響を与える可能性があり、単純な統計的手法では適切にモデリングできない。
    • トレンド
      • 時系列データには長期的なトレンドが存在することがある。トレンドは平均値の変動を引き起こし、推定結果に影響を与える可能性がある。

時系列データとは区別するべきデータ

地震観測データや為替取引の Tick データは、「点過程（ポイントプロセス）データ」または「マーク付き点過程データ」と呼ばれます。これらのデータは、発生時間と発生間隔に意味があり、観測者が発生間隔を設定できません。

これらも時系列分析の対象となり得ますが、通常の時系列データとは異なる特性を持つため、それに応じた解析手法やモデルの適用が必要となります。

時系列分析で出来ること・わかること

時系列分析では、以下について知ることができます。

  • トレンド
    • 時系列データに現れる長期的な変化や傾向（トレンド）を把握できます。トレンドの有無や方向、変動のパターンを特定することができます。
  • 季節性
    • 時系列データに現れる季節的な変動を把握できます。特定の季節パターンや周期性を検出し、季節要素の影響を理解することができます。
  • 周期性
    • 時系列データに周期的な変動がある場合、その周期や周期の長さを特定することができます。サイクルの長さや振幅の変動を分析することができます。
  • 予測
    • 過去の時系列データを基に将来の値を予測することができます。予測モデルを構築し、将来の傾向や変動を推定することができます。
  • 異常値検出
    • 時系列データから異常値や外れ値を検出することができます。異常な振る舞いや予測モデルからの逸脱を特定し、異常値の原因や特徴を分析することができます。
  • 時系列データ間の相関関係
    • 他の時系列データや外部要因との関係を分析し、相互の影響や連動性を理解することができます。

これらの分析を通じて、時系列データの特徴や変動要因を理解し、将来の予測や意思決定に活用することができます。

時系列分析の手順

  1. データの特性の理解
    • 時系列データの基本的な特性を調査し、傾向（トレンド）、季節性、周期性、ランダム性などのパターンを特定。
  2. モデリング
    1. 時系列データにモデルを適用し、データの生成プロセスを表現するための数学的なモデルを構築。
    2. 代表的なモデル
      • ARIMA（Autoregressive Integrated Moving Average）
      • SARIMA（Seasonal ARIMA）
      • VAR（Vector AutoRegression）など
  3. 予測
    • 構築したモデルを使用し将来のデータポイントを予測。
  4. モデルの診断と改善
    • 構築したモデルの適合度や残差の診断を行い、モデルの改善や修正を行う。モデルの信頼性と予測精度を向上させていく。

サンプルデータ

今回は、航空機の乗客数データを使って未来の予測や季節性などの理解を行っていきます。

その際に、ARIMA モデル、及び SARIMA モデルを作成します。

R 言語の標準データセットとして提供されている「AirPassengers」を利用しますが、今回は Python で進めるため、Kaggle で AirPassengers データセットをダウンロードしておきます。

Air Passengers - Kaggle

データの読み込みと調整

ライブラリの読み込みを行い、データを読み込みます。

import pandas as pd
import numpy as np
import statsmodels
import statsmodels.api as sm

from matplotlib import pylab as plt
import seaborn as sns
%matplotlib inline
sns.set()

from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 16, 8

# データ読み込み
df = pd.read_csv('AirPassengers.csv')

読み込んだデータを確認します。

# df

    Month   #Passengers
0   1949-01 112
1   1949-02 118
2   1949-03 132
3   1949-04 129
4   1949-05 121
... ... ...
139 1960-08 606
140 1960-09 508
141 1960-10 461
142 1960-11 390
143 1960-12 432


# df.dtypes

Month          object
#Passengers     int64
dtype: object

データを使いやすくするためにカラム名をリネームし、Month をインデックスにしておきます。

## カラム名リネーム
df = df.rename(columns={'#Passengers': 'Passengers'})
## datetime 化
df['Month'] = pd.to_datetime(df['Month'] )

# Month をインデックスにする
df.set_index('Month', inplace=True)

再度データを確認します。

# df

      Passengers
Month   
1949-01-01  112
1949-02-01  118
1949-03-01  132
1949-04-01  129
1949-05-01  121
... ...
1960-08-01  606
1960-09-01  508
1960-10-01  461
1960-11-01  390
1960-12-01  432
144 rows × 1 columns


# df.dtypes

Passengers    int64
dtype: object

Month をインデックスに変換したのは、整列・集計・プロットなどを行いやすくするためです。

データの確認

データの読み込みと整理が済んだので改めてデータを見ていくと、1949 年 1 月から、1960 年 12 月までの乗客数データが収録されています。 プロットして見てみます。

plt.plot(df['Passengers'])

視覚的には、どことなく似たような周期のアップダウンを繰り返しながら上昇しているように見えます。

この時系列データが非定常性を持つかを数値的に判断するために、ディッキーフラー検定を実施してみます。

Dickey-Fuller 検定（ディッキーフラー検定）

ディッキーフラー検定は、自己回帰モデルにおける単位根の有無の検定です。

  • 帰無仮説：「データ系列に単位根が存在する」
  • 対立仮説：「データ系列は定常性を有す」

時系列データの平均、分散、自己相関係数などの計算をはじめとした時系列データの分析をする際には、事前に検定を行い、分析対象の時系列データが定常性を有するか確かめることが必要です。チャートを見ると非定常性が明らかな感じがしますが、見た目だけ判断せず検定を行って確認をしておきます。

単位根

単位根（unit root） は、時系列データの性質を表す統計的な概念です。単位根を持つ時系列データは、長期的なトレンドや構造的な変化が存在し、平均値が一定でない（恒久的にランダムウォークをするような特性を持っている）ことを示します。単位根は統計的仮説検定によって検出され、データの非定常性を考慮する必要があります。

つまり、Dickey-Fuller 検定によって帰無仮説が棄却されない場合、単位根が存在すると判断し、その時系列データは非定常性である。と結論づけられます。

非定常性

非定常性（non-stationary） とは、時系列データの統計的性質が時間に依存し、一定のパターンや特性を持たない状態を指します。非定常性を持つ時系列データは、平均や分散が時間の経過とともに変動する傾向 があります。

非定常性を示す時系列データには、以下のような特徴があります:

  • トレンド（Trend）
    • データが長期的に上昇または下降する傾向がある場合、トレンドが存在します。トレンドは、統計的に見て平均値が時間とともに変化することを示します。
  • 季節性（Seasonality）
    • データに周期的なパターンや季節的な変動がある場合、季節性が存在します。季節性は、統計的に見て周期的な変動があることを示します。
  • 周期性（Cyclicity）
    • データに定期的な周期性があり、季節性とは異なる場合、周期性が存在します。周期性は、統計的に見て一定の期間で変動が生じることを示します。
  • 自己相関（Autocorrelation）
    • データの過去の値との相関関係があり、現在の値が過去の値に依存する場合、自己相関が存在します。自己相関は、統計的に見てデータが時間的な依存関係を持つことを示します。

非定常性を持つ時系列データは、定常性の仮定を満たさないため、統計モデリングや予測の精度を下げる可能性があります。そのため、非定常性のデータを扱う場合は、定常性を復元するための前処理や変換が必要となる場合があります。例えば、トレンドの除去や差分化、季節性の調整などが一般的なアプローチとして使用されます。これにより、データの統計的性質が一定である定常な状態に変換され、より正確な予測やモデリングが可能になります。

ディッキー・フラー検定の実施

Python では statsmodels の adfuller 関数で実施します。

statsmodels.tsa.stattools.adfuller

# ディッキー・フラー検定の実施
result = sm.tsa.stattools.adfuller(df['Passengers'])

# 結果の表示
print('ADF統計量（ADF Statistic）:', result[0])
print('p値（p-value）:', result[1])
print('臨界値（Critical Values）:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))

結果出力

ADF統計量（ADF Statistic）: 0.8153688792060421
p値（p-value）: 0.9918802434376409
臨界値（Critical Values）:
    1%: -3.482
    5%: -2.884
    10%: -2.579

ディッキー・フラー検定の結果を解釈するためには、ADF 統計量と p 値を考慮し、臨界値と比較して帰無仮説を棄却するかどうかを確認します。

ADF 統計量は 0.8153688792060421 です。この値を臨界値と比較することで、時系列データの非定常性を評価します。臨界値は 1%、5%、および 10% の有意水準に対する値であり、以下のようになります

  • 1%: -3.482
  • 5%: -2.884
  • 10%: -2.579

結果を解釈するためには、ADF 統計量を臨界値と比較します。ADF 統計量が臨界値よりも小さい場合は定常性の存在が示唆され、帰無仮説を棄却できます。逆に、ADF 統計量が臨界値よりも大きい場合、定常性の存在が示されず、帰無仮説を棄却できません。

今回の場合、ADF 統計量の値（0.8153688792060421）は臨界値（-3.482、-2.884、-2.579）よりも大きいため、定常性の存在が示されません。また、p 値が非常に高い（0.9918802434376409）ため、帰無仮説を支持し、定常性の存在を示す証拠はほとんどありません。

したがって、与えられた結果からは「帰無仮説を棄却できない」ということになり、「データ系列に単位根が存在する」つまり「定常性を有するとは言えず、非定常性を有する可能性が高い」と結論付けられます。

非定常性・定常性を持つデータであるか

乗客数の原系列データ（未加工のもとのデータ）を見れば、年々数値が上昇（トレンドがある）していたり、そのトレンドの中でも似たようなアップダウン（周期性・季節性）があることがなんとなくわかると思います。

時系列分析におけるポイントとして、予測のためのモデル構築を行う際には、こういった「非定常データ」を「定常データ」という、トレンドや周期性・季節性を排除したデータへ変換してからモデルを構築する必要がある。という背景があります。

非定常データから定常データへの変換を行っていく過程で、変換後のデータが定常性を持つデータになっているかを確認するために、ディッキー・フラー検定は有効です。

時系列データの理解

続いて、今回の時系列データがどのような性質をもっているのかを確認していきます。

自己相関と偏自己相関を視覚化し、データの構造やパターンを理解します。

自己相関（Autocorrelation）

自己相関は、時系列データ内の異なる時間ラグ間の相互関係を測る統計量です。具体的には、ある時刻のデータと一定の時間ラグだけずれた時刻のデータとの間の相関を計算します。自己相関を調べることで、時系列データの周期性やパターンの特徴を把握することができます。

自己相関のグラフは、時系列データの過去の値と現在の値との間の相関関係を示します。このグラフを通じて、データがどれくらい自己相関を持っているか、周期性や季節性のパターンが存在するか、および他のタイプの相関関係があるかを視覚的に確認することができます。自己相関が高い場合、データが過去の値に依存しており、トレンドや周期性のパターンが存在する可能性があります。

また、自己相関、及び次に紹介する偏自己相関のグラフは一般に コレログラム（Correlogram） と呼ばれることもあります。

# 自己相関
fig = sm.graphics.tsa.plot_acf(df['Passengers'], lags=40)

x=0 の地点にある値が基点となるデータで、x=1 の地点にある値がラグ 1 となるデータです。つまり今回の乗客数データでは、x=0 が基点月で、x=1 が前月、x=1 が二ヶ月前で... という読み方ができます。

このコレログラムを見ると、 基点月と前月は強い正の自己相関を示しています。

また、全体的に一定のパターンで推移していそうなことも読み取れます。

グラフの青色の範囲は 95%信頼区間（=ACF{自己相関関数}の推定値が統計的に有意でないと考えられる範囲）を示してしますが、範囲を抜けているラグが 1 年くらいというのもあること、チャートの推移から1 年周期でのパターン（12の周期）があることが読み取れます。

偏自己相関（Partial Autocorrelation）

偏自己相関は、ある時刻と別の時刻の間で、他の時刻の影響を取り除いた相関を測る統計量です。つまり、一連の中間時刻を通じて制約されない2つの時刻間の直接的な相関を計算します。偏自己相関を調べることで、直接的な相関関係を評価することができます。

偏自己相関のグラフは、時系列データの過去の値と現在の値との間の直接的な相関関係を示します。このグラフを通じて、データがどの時点で自己相関を持つのか、他の時点の影響を排除した直接の相関関係が存在するかを視覚的に確認することができます。

# 偏自己相関
fig = sm.graphics.tsa.plot_pacf(df['Passengers'], lags=40)

0 地点のデータは 1 地点前のデータと強い正の自己相関にあり、前月の乗客数が多ければ当月も多くなることがわかります。

また、12 ヶ月周期の相関が見られるので、周期性ないし季節性があるとわかります。

和分課程〜非定常データから定常データへ

これらの時系列データから予測を立てられるようにするには、先ほど確認した時系列的な変動パターンがあると予測が立てられないため、周期変動を除去し、非定常から定常性のデータに変換していく必要があります。

時系列分析において、非定常データを定常データに変換することにはいくつかの重要な意味があります。

  • 定常性の仮定
    • 時系列データが定常である場合、データの統計的特性が時間に依存しなくなります。この特性は、時系列分析において重要な仮定です。定常データは、平均や分散が一定であり、自己共分散（自己相関）が時間に依存しないという特徴を持ちます。定常データである場合、統計的手法がより正確で信頼性の高い結果を提供しやすくなります。
  • モデルの安定性
    • 非定常データをそのまま使った場合、データの統計的特性が時間に依存しているため、モデルのパラメータが時間経過とともに変化してしまう可能性があります。これによって、予測の精度が低下したり、モデルの安定性が損なわれたりすることがあります。定常データに変換することで、モデルのパラメータを安定させることができます。
  • データ解釈と比較
    • 定常データは、時間に依存しない統計的特性を持つため、データの解釈や異なる時点のデータとの比較が容易になります。また、非定常データではトレンドや季節性の影響がデータ全体に影響を及ぼすため、異なる時点のデータを比較することが困難になることがあります。

非定常データを定常データに変換するために、和分過程（差分化）や季節調整などの方法が利用されます。これらの手法を用いることで、データの非定常性を取り除き、より信頼性のある分析や予測モデルの構築が可能となります。

モデルを構築する前に、これら時系列データの非定常性を除去していく様子を見ていこうと思います。

ちなみに、データの差分を取ることで非定常性（単位根）を除去した後の過程を 和分過程（Integrated Process） といいます。

差分系列（階差系列）

差分系列（階差系列, Difference series）は、時系列データの、連続する観測値間の差分を計算した系列（1時点離れたデータとの差をとったデータ）です。差分系列を計算することで、トレンドや季節性の影響を除去し、定常性を持ったデータに近づけることができます。

差分系列は数式で表すと以下で表現できます。

\displaystyle \Delta y_t = y_t - y_{t-1}

y_t は時点 t における観測値を表し、y_{t-1} はその直前の時点 t-1 における観測値を示します。差分系列 \Delta y_t は、現在の観測値 y_t から直前の観測値 y_{t-1} を引いた値として計算されます。

Passengers（乗客数）の値から階差系列を出力します。

# 差分系列
df['difference'] = df['Passengers'].diff()

# df

    Passengers  difference
Month       
1949-01-01  112 NaN
1949-02-01  118 6.0
1949-03-01  132 14.0
1949-04-01  129 -3.0
1949-05-01  121 -8.0
... ... ...
1960-08-01  606 -16.0
1960-09-01  508 -98.0
1960-10-01  461 -47.0
1960-11-01  390 -71.0
1960-12-01  432 42.0
144 rows × 2 columns

一番はじめ、1949-01-01 の difference が NaN なのは、1 地点前のデータが存在しないためです。

差分系列をプロットしてみます。

plt.plot(df['difference'])

トレンドは除去されたようです。

しかし振れ幅がだんだん大きくなっていっているので分散はまだ一定にはなっていないことがわかります。

したがって、この時点では完全に定常性を持ったデータにはなっていないことがわかります。

対数系列（Log series）

対数系列は、時系列データの各観測値に対して対数変換を適用した系列です。対数変換により非線形な変動を緩和し、データの特性を正規分布に近づけることができます。

数式では以下で表現します。

\displaystyle \log y_t = log(y_t)

y_t は個々の時点 tにおける観測値や測定値を示し、 \log y_t は y_t の自然対数を表します。この数式では、各時点の y_t に対して自然対数を適用し、その結果を \log y_t として表現しています。

定常性を持ったデータに近づけるため、こちらも試してみます。

# 対数系列（底が10となる対数を作成する）
df['log'] = np.log10(df['Passengers'])

対数系列データを確認します。

# df

Passengers  difference  log
Month           
1949-01-01  112 NaN 2.049218
1949-02-01  118 6.0 2.071882
1949-03-01  132 14.0    2.120574
1949-04-01  129 -3.0    2.110590
1949-05-01  121 -8.0    2.082785
... ... ... ...
1960-08-01  606 -16.0   2.782473
1960-09-01  508 -98.0   2.705864
1960-10-01  461 -47.0   2.663701
1960-11-01  390 -71.0   2.591065
1960-12-01  432 42.0    2.635484
144 rows × 3 columns

グラフに描画します。

# 対数系列をプロット
plt.plot(df['log'])

トレンドはあるものの、振れ幅（分散）はあまり変わっていないように見えます。

対数系列を取ると、分散がある程度一定に抑えられることがわかりました。

対数差分系列（Log Difference series）

階差系列への変換ではトレンドが除去でき、対数系列への変換では分散をある程度一定に抑えられることがわかったので、対数差分系列への変換を行ってみます。

対数差分系列は、時系列データの対数変換と差分の組み合わせです。具体的には、各時点の観測値に対して対数変換を適用し、その後に差分を計算します。

対数差分系列は、対数変換によって非線形性を緩和し、変動の幅を縮小させる効果があります。また、差分を取ることでトレンドや季節性の影響を除去することができます。対数変換は正規分布に近い性質を持つデータに対して効果的であり、対数差分系列は定常性の要件を満たすことが多いです。

対数差分系列の数式は次のように表されます

\displaystyle \Delta \log y_t = \log y_t - \log y_{t-1}

ここで、 y_t は時点 t における観測値を表し、 \log y_t はその対数変換を示します。差分系列 \Delta \log y_t は、現在の対数変換後の観測値 \log y_t から直前の対数変換後の観測値 \log y_{t-1}を引いた値として計算されます。

# 対数差分系列
df['log_difference'] = np.log10(df['Passengers']).diff(periods=1)

対数差分系列データを確認します。

# df

    Passengers  difference  log log_difference
Month               
1949-01-01  112 NaN 2.049218    NaN
1949-02-01  118 6.0 2.071882    0.022664
1949-03-01  132 14.0    2.120574    0.048692
1949-04-01  129 -3.0    2.110590    -0.009984
1949-05-01  121 -8.0    2.082785    -0.027804
... ... ... ... ...
1960-08-01  606 -16.0   2.782473    -0.011318
1960-09-01  508 -98.0   2.705864    -0.076609
1960-10-01  461 -47.0   2.663701    -0.042163
1960-11-01  390 -71.0   2.591065    -0.072636
1960-12-01  432 42.0    2.635484    0.044419
144 rows × 4 columns

グラフに描画します。

# 対数差分系列をプロット
plt.plot(df['log_difference'])

トレンドが除去でき、分散もおおよそ一定に近づいたではないでしょうか。

一方で、周期性は残っているように見えます。年単位で周期があるように見受けられるため、周期性というより季節性かもしれません。

一旦、ここまで非定常データを定常に近づけるために変換処理を行ったデータに対して、再度ディッキー・フラー検定を実施し、自己相関・偏自己相関も見てみます。

# 対数差分系列の先頭行は NaN のため除去
list = df['log_difference'][1:]

# ディッキー・フラー検定の実施
result = sm.tsa.stattools.adfuller(list)

# 結果の表示
print('ADF統計量（ADF Statistic）:', result[0])
print('p値（p-value）:', result[1])
print('臨界値（Critical Values）:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))

# 出力された結果

ADF統計量（ADF Statistic）: -2.7171305983880982
p値（p-value）: 0.07112054815086455
臨界値（Critical Values）:
    1%: -3.483
    5%: -2.884
    10%: -2.579

ADF統計量が -2.7171305983880982 であり、p値が 0.07112054815086455 です。この場合、有意水準 5% での臨界値（-2.884）を下回っていませんが、有意水準 10% での臨界値（-2.579）を下回っています。

つまり、p 値が 0.05 より大きいため、5% の有意水準で帰無仮説を棄却できませんが、10% の有意水準で帰無仮説を棄却できる可能性があります。

したがって、データが非定常かどうかについてはやや曖昧な結果となりますが、有意水準 10% での臨界値を下回っていることから、データには定常性がある可能性も出てきました。変換処理に効果が出ているということですね。

自己相関と偏自己相関も見てみます。

# 自己相関をプロット
fig = sm.graphics.tsa.plot_acf(list, lags=40)


# 偏自己相関をプロット
fig = sm.graphics.tsa.plot_pacf(list, lags=40)

まだ季節性が残っていることが見受けられます。

よって完全なる定常データへ変換できたか？という問いに対しては「まだ季節性が残っているため十分ではない」とし、さらなる変換が必要ということになります。

現時点で全ての非定常性要素を除去できたわけではなりませんが、このように、通常であれば非定常な時系列データを定常データに変換していき、それを使って予測モデル（ARIMA, SARIMA etc...）を構築していく流れになります。

ARIMA モデル

本記事では割愛していますが、ARMA モデル（自己回帰移動平均モデル）というものがあります。

この ARMA モデルは、定常データに対しては説明能力は良いが非定常データには使えません。そのため、差分系列をとって定常過程に変換してから、ARMA モデルを適用することを考えます。これを ARIMA（アリマ, 自己回帰和分移動平均モデル、Autoregressive Integrated Moving Average）モデルといいます。

差分系列へ ARMA モデルを適用する場合、d 次和分過程を I(d) と書くので、真ん中に入れて “ARIMA” と呼ばれます。

ARIMA = AR 過程（Auto Regressive process, 自己回帰過程）＋和分課程（I）＋MA 過程（Moving Average process, 平均移動過程）

ということで、データは原系列データではなく、差分系列データを用いてモデルを作成していきます。

データの分割

データを学習用とテスト用に分けます。

今回のデータでは、1949-01 〜 1957-12-01 までのデータを学習用として使い、ARIMA モデルを構築します。

1958-01 〜 1960-12 までのデータはテスト用とし、作成したモデルで同期間の予測を行った後、このテスト用データと付け合わせて実際の予測がどれだけの精度かを確認します。

# データを学習用とテスト用に分ける
df_train_arima_diff = df['difference'][:'1957-12-01'].dropna()
df_test_arima_diff = df['difference']['1958-01-01':]

次数の決定

ARIMAモデルは、3 つの主要なパラメータ（p、d、q）によって定義されます。これらのパラメータは、モデルの自己回帰（AR）成分、積分（I）成分、および移動平均（MA）成分を制御します。

  • p（自己回帰次数）
    • 自己回帰成分の次数を示します。直前の時刻の値にどれだけの過去の値を使用するかを示します。大きな値は、過去の多くの時刻の値が予測に影響を与えることを意味します。
    • AR 過程の次数
  • d（積分次数）
    • 積分成分の次数を示します。これは、元の時系列データが非定常過程（トレンドや季節性の影響を受けて変動する）であるかどうかを示します。dの値が 0 であれば、元の時系列データは定常過程と見なされます。d の値が 1 以上であれば、差分を取ることによってデータを定常化します。
  • q（移動平均次数）
    • 移動平均成分の次数を示します。誤差項に過去の誤差値をどれだけ使用するかを示します。大きな値は、過去の誤差値が予測に影響を与えることを意味します。
    • MA 過程の次数

これらの次数を最適値を見つけるために、statsmodels の arma_order_select_ic 関数を利用します。

statsmodels.tsa.stattools.arma_order_select_ic

この関数は、情報基準（Information Criterion）を使用して異なる次数の組み合わせを比較し、最適な次数を見つけるのに役立ちます。

一般的に、AIC（Akaike Information Criterion）やBIC（Bayesian Information Criterion）などの情報基準が使用されます。

sm.tsa.arma_order_select_ic(df_train_arima_diff, ic='aic', trend='n')

結果出力

{'aic':              0           1           2
 0  1001.530812  990.101618  987.950157
 1   994.820617  987.280756  982.138924
 2   990.473898  981.180360  983.831761
 3   991.560168  983.089715  978.733996
 4   982.579395  984.165016  978.372978,
 'aic_min_order': (4, 2)}

AIC が最も低いものが最も良いモデルとされ、p=4, q=2 が　AIC=978.372978 と最も低く最適な次数であるという結果が出たので、これを使っていきます。

モデル作成

ARIMA モデルは statsmodels の ARIMA() で作成できます。

statsmodels.tsa.arima.model.ARIMA

from statsmodels.tsa.arima.model import ARIMA

# ARIMA モデル作成
model = ARIMA(data_diff, order=(4, 1, 2))

result = model.fit()

推定されたパラメータを確認してみます。

# result.params

ar.L1      -0.408674
ar.L2       0.041754
ar.L3      -0.208824
ar.L4      -0.333281
ma.L1      -0.179896
ma.L2      -0.817283
sigma2    853.547001
dtype: float64

モデルの適合度と仮定の検証

作成したモデルがデータにどれだけ適合しているかを、残差を使って確認します。

ARIMAモデルは、一定の仮定に基づいています。例えば、残差は平均ゼロ、定常性を持ち、自己相関を持たないという仮定があります。残差のプロットや統計テストを通じて、これらの仮定が満たされているかどうかを確認できます。

残差の正規性を見たいのでヒストグラムをプロットします。

# 残差
res = result.resid

# ヒストグラムをプロット
plt.hist(res, bins=16)

左に寄っています。正規分布に従うとは言えないでしょう。

続いて自己相関や偏自己相関を見て、周期性ないし季節性が無いことを確認します。

fig = plt.figure(figsize=(12,8))

# 自己相関
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(res.values.squeeze(), lags=40, ax=ax1)

# 偏自己相関
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(res, lags=40, ax=ax2)

周期性が残ってしまっています。

季節性成分に対応する

モデルの検証を行った結果、残念ながらこれでは予測に使えないという結論になりました。

ARIMAモデルは、データの定常性や自己相関などの統計的な特性に基づいて設計されているため、周期性が残ってしまうとモデルがデータに適合しづらくなり、予測精度が低下する可能性があります。

特に、季節性が強いデータでは、ARIMAモデルだけでは十分な適合が難しいことがあります。

サンプルデータは、前項で行った変換処理を通じて季節性のあるデータであることはわかっていましたがその上で ARIMA モデルを作成したので当然の結果ではあります。

ちなみに参考までに。

この「使えない」モデルで予測を行うと以下のようになります。青い線が正解データで、赤い線が作成したモデルを使って予測したものです。

このように、周期性が残っているデータに対しては ARIMA モデルでは限界があるため、次に季節変動ありの ARIMA である SARIMA モデルを作成していきます。

SARIMAモデル

SARIMA（Seasonal ARIMA）モデルは、ARIMA モデルに季節性成分を追加したモデルです。季節性のあるデータに適用され、季節性成分をモデル化するためのパラメータを追加します。SARIMAモデルは、季節的な変動を捉えるため、季節性パターンを特定するのに有用です。

ARIMA と SARIMA は過去の値や誤差項を用いて現在の値を予測します。ARIMA と SARIMAは一般的に季節性のないデータに適用されますが、SARIMA は季節性を持つデータにも対応できます。

データの分割

データを学習用とテスト用に分けます。分ける範囲は ARIMA のときと同じです。

  • 学習用：1949-01 〜 1957-12-01
  • テスト用：1958-01 〜 1960-12

学習用データを使って SARIMA モデルを作成後、テスト用範囲を予測、結果をテスト用データと突合させて精度を確認する。という流れです。

今回も差分系列データを使っていきます。

# データを学習用とテスト用に分ける
df_train_diff = df['difference'][:'1957-12-01'].dropna()
df_test_diff = df['difference']['1958-01-01':]

SARIMAX モジュール

SARIMAX モデルは statsmodels の SARIMAX 関数で作成できます。

statsmodels.tsa.statespace.sarimax.SARIMAX

from statsmodels.tsa.statespace.sarimax import SARIMAX

次数の決定

SARIMA モデルの次数は、元のARIMAモデルの次数（p、d、q）と季節性成分の次数（P、D、Q、s）から構成されます。

  • P（季節自己回帰次数）
    • 季節性の自己回帰成分の次数です。季節性成分のパターンをモデル化します。
  • D（季節積分次数）
    • 季節性成分の積分次数です。季節性の差分を取る回数を示します。
  • Q（季節移動平均次数）
    • 季節性の移動平均成分の次数です。誤差項に季節性の過去の誤差値をどれだけ使用するかを示します。
  • s（季節周期）
    • 季節性の周期を示します。月次データであれば 12（1年の周期）、四半期データであれば 4（1年の四半期）などです。

これらの次数を見つけるために、総当りで SARIMA モデルの次数を探索し、AIC が最も低いモデルを見つけます。

（ちなみに総当りは手段の 1 つです。データの量が多い場合や計算時間が制約されている場合には他の方法も検討できます。）

処理が多いので関数として定義します。

def find_model_with_lowest_aic(df_train):
    """
    総当りで SARIMA モデルを作成し、最も低い AIC 値を持つモデルのパラメータと AIC を出力します。

    Parameters:
    df_train (DataFrame): テスト用データ

    Returns:
    string: 最も低い AIC 値を持つモデル

    """
    
    # 各パラメータの候補値リスト
    ## p, d, q
    p_list = [1, 2]
    d_list = [1]
    q_list = [1, 2]
    
    ## 季節項（P, D, Q）
    sp_list = [1, 2]
    sd_list = [1]
    sq_list = [1, 2]

    ## 何ヶ月か（s）
    m = 12

    parameter = []
    results   = []
    
    # 総当りで SARIMA モデルを作成
    for p in p_list:
        for d in d_list:
            for q in q_list:
                for sp in sp_list:
                    for sd in sd_list:
                        for sq in sq_list:
                            # パラメータを格納
                            parameter.append([p, d, q, sp, sd, sq])
                            # モデル作成
                            model = SARIMAX(df_train, order=(p, d, q), seasonal_order=(sp, sd, sq, m))
                            aic = model.fit(disp=0).aic
                            results.append({'parmeter': [p, d, q, sp, sd, sq], 'aic': aic})
                            print('parmeter', [p, d, q, sp, sd, sq], ', AIC=', aic)

    # 比較用: 最小 AIC 値の初期値を設定
    min_aic = float('inf')
    best_result = None

    # 最小 AIC 値を持つモデルを見つける
    for result in results:
        aic = result['aic']
        if aic < min_aic:
            min_aic = aic
            best_result = result
    
    print("最も低い AIC 値を持つモデル:", best_result)

各パラメータの候補値リスト（p_list, d_list, q_list, sp_list, sd_list, sq_list`）に定義されている数字は、それぞれ SARIMA モデルの次数に対する候補値です。これらの値を組み合わせてモデルを構築し、AIC を計算して最適なモデルの次数を決定します。

具体的な値はデータやドメイン知識に基づいて選択する必要がある。として、決め打ちです。とはいえ過学習のリスクや計算量コストを考えると、基本的には以下に沿うことになるかなと思っています。

  • p_list（自己回帰次数の候補値）
    • SARIMA モデルの自己回帰次数（AR次数）に関する候補値を指定
    • 通常、小さな整数値（1から3程度）を選択
  • d_list（積分次数の候補値）
    • SARIMA モデルの積分次数（差分の取る回数）に関する候補値を指定
    • データの定常性に応じて、0 または 1 を選択
  • q_list（移動平均次数の候補値）
    • SARIMA モデルの移動平均次数（MA次数）に関する候補値を指定
    • 通常、小さな整数値（1から3程度）を選択
  • sp_list（季節自己回帰次数の候補値）
    • 季節性の自己回帰次数（季節AR次数）に関する候補値を指定
    • 通常、小さな整数値（1から3程度）を選択
  • sd_list（季節積分次数の候補値）
    • 季節性の積分次数（季節差分の取る回数）に関する候補値を指定
    • データの季節性に応じて、0 または 1 を選択
  • sq_list（季節移動平均次数の候補値）
    • 季節性の移動平均次数（季節MA次数）に関する候補値を指定
    • 通常、小さな整数値（1から3程度）を選択

それではこの関数を使って次数を決定していきます。

# warning 多ければ抑制
# import warnings
# warnings.filterwarnings('ignore')

# 関数実行
find_model_with_lowest_aic(df_train_diff)

出力された結果は以下

parmeter [1, 1, 1, 1, 1, 1] , AIC= 706.6661313300135
parmeter [1, 1, 1, 1, 1, 2] , AIC= 700.7887100270527
parmeter [1, 1, 1, 2, 1, 1] , AIC= 702.2413978900132
parmeter [1, 1, 1, 2, 1, 2] , AIC= 701.9189059594827
parmeter [1, 1, 2, 1, 1, 1] , AIC= 707.7131926442661
parmeter [1, 1, 2, 1, 1, 2] , AIC= 702.5540894596875
parmeter [1, 1, 2, 2, 1, 1] , AIC= 703.8802139331958
parmeter [1, 1, 2, 2, 1, 2] , AIC= 703.5023730470815
parmeter [2, 1, 1, 1, 1, 1] , AIC= 708.2181279257733
parmeter [2, 1, 1, 1, 1, 2] , AIC= 702.7085048503061
parmeter [2, 1, 1, 2, 1, 1] , AIC= 704.100556663309
parmeter [2, 1, 1, 2, 1, 2] , AIC= 703.745258100657
parmeter [2, 1, 2, 1, 1, 1] , AIC= 709.7110131616942
parmeter [2, 1, 2, 1, 1, 2] , AIC= 700.8798954969151
parmeter [2, 1, 2, 2, 1, 1] , AIC= 702.3974872249754
parmeter [2, 1, 2, 2, 1, 2] , AIC= 702.2549954185089

最も低い AIC 値を持つモデル: {'parameter': [1, 1, 1, 1, 1, 2], 'aic': 700.7887100270527}

最も低いAIC値を持つモデルは、パラメータ p=1, d=1, q=1, P=1, D=1, Q=2 のものとなりました。 これらを当てはめて再度モデルを作成します。

# SARIMA モデル作成
#        SARIMAX(df_train_diff, order=(p,d,q), seasonal_order=(P,D,Q,m)).fit()
r_diff = SARIMAX(df_train_diff, order=(1,1,1), seasonal_order=(1,1,2,12)).fit()

# r_diff.summary()

SARIMAX Results

Dep. Variable:  difference  No. Observations:   107
Model:  SARIMAX(1, 1, 1)x(1, 1, [1, 2], 12) Log Likelihood  -347.466
Date:   Fri, 18 Aug 2023    AIC 706.931
Time:   08:31:16    BIC 722.191
Sample: 02-01-1949  HQIC    713.095
      - 12-01-1957      
Covariance Type:    opg     

          coef      std err z   P>|z|    [0.025  0.975]
ar.L1     -0.2249   0.087   -2.599  0.009   -0.395  -0.055
ma.L1     -0.9964   0.427   -2.332  0.020   -1.834  -0.159
ar.S.L12  0.5576    0.727   0.768   0.443   -0.866  1.982
ma.S.L12  -0.8143   0.772   -1.055  0.291   -2.327  0.698
ma.S.L24  0.2828    0.213   1.330   0.184   -0.134  0.700
sigma2  77.6646 32.494  2.390   0.017   13.977  141.352

Ljung-Box (L1) (Q): 0.04    
Jarque-Bera (JB):   3.46
Prob(Q):    0.85    
Prob(JB):   0.18
Heteroskedasticity (H): 1.36    
Skew:   0.44
Prob(H) (two-sided):    0.40    
Kurtosis:   2.69

モデルの適合度と仮定の検証

作成したモデルがデータにどれだけ適合しているかを、残差を使って確認します。

ARIMA のときはヒストグラムやコレログラムなどを個別に見ていきましたが、statsmodels の plot_diagnostics 関数を使うと複数の残差の診断プロットを一回で生成できます。

statsmodels.tsa.arima.model.ARIMAResults.plot_diagnostics

# 残差の診断プロット
r_diff.plot_diagnostics(lags=20);

plot_diagnostics 関数は、ARIMA および SARIMA モデルの残差の診断を行うために使用されます。生成される 4つのチャートは、モデルの残差が一定の基準を満たしているかどうかを確認するためのものです。

Standardized Residuals Plot (標準化残差プロット)
このプロットは、モデルの残差を標準化したものを時系列で表示します。残差がランダムにばらついていることが期待されます。もし残差にパターンやトレンドが見られる場合、モデルがデータに適合していない可能性があります。
Histogram Plus Estimated Density Plot (ヒストグラムと推定密度プロット)
このプロットは、残差のヒストグラムとカーネル密度推定を表示します。残差が正規分布に近いかどうかを確認するのに役立ちます。理想的には、ヒストグラムと密度推定の形状が正規分布に近い形になることが望ましいです。
Normal Q-Q (Quantile-Quantile) Plot (正規Q-Qプロット)
このプロットは、残差の分位数を正規分布の分位数と比較したものです。正規分布に従う場合、プロットされた点は対角線に近い位置に分布します。プロットされた点が対角線から外れている場合、残差が正規分布から逸脱している可能性があります。
Correlogram (自己相関プロット)
このプロットは、残差の自己相関係数をタイムラグに対してプロットします。残差が白色雑音（無相関性）である場合、自己相関プロットはほとんどのラグでゼロに近くなるはずです。自己相関プロットにおけるラグがゼロでない値を示す場合、残差に時系列的なパターンが残っている可能性があります。

パターンやトレンドはなく、まあまあ正規分布に近く、QQ プロットはおおむね対角線に乗っている。コレログラムでも、周期性・季節性は見られない。ということで、おおそよ良いのではないかと思われます。

では、作成したモデルを使って実際に予測を行い、テスト用データと突合させてみます。

pred_diff = r_diff.predict('1958-01-01', '1960-12-01')

予測結果

1958-01-01    10.862560
1958-02-01   -12.957117
1958-03-01    51.532015
1958-04-01    -6.048657
1958-05-01     6.919694
1958-06-01    68.024899
1958-07-01    44.618030
1958-08-01    -2.338867
1958-09-01   -60.685910
1958-10-01   -57.188398
1958-11-01   -41.271487
1958-12-01    34.165328
1959-01-01    10.742535
1959-02-01   -13.216177
1959-03-01    53.727288
1959-04-01    -6.715341
1959-05-01     7.652743
1959-06-01    71.907824
1959-07-01    44.975607
1959-08-01    -0.446448
1959-09-01   -63.921007
1959-10-01   -59.853134
1959-11-01   -42.304585
1959-12-01    33.934304
1960-01-01    11.246330
1960-02-01   -13.367209
1960-03-01    55.074747
1960-04-01    -6.993050
1960-05-01     8.162161
1960-06-01    74.172291
1960-07-01    45.274499
1960-08-01     0.708273
1960-09-01   -65.625626
1960-10-01   -61.239695
1960-11-01   -42.781264
1960-12-01    33.904902

ここまで差分系列データで進めてきたため、結果も差分系列データのままです。実際の乗客数データに変換して戻します。

def inverse_difference(initial_value, diff_series):
    """
    差分系列から原系列に戻す

    Parameters:
    initial_value (int): 予測開始した月の前月の乗客数
    diff_series (int): 差分系列データ

    Returns:
    diff_series: 原系列データと同じ次元の予測データ

    """
    cum_sum = diff_series.cumsum()
    original_series = cum_sum + initial_value
    return original_series

# 1957-12-01 の乗客数
initial_value = df['Passengers']['1957-12-01']
# 差分系列から原系列に戻す
pred_original = inverse_difference(initial_value, pred_diff)

変換結果

1958-01-01    346.862560
1958-02-01    333.905443
1958-03-01    385.437459
1958-04-01    379.388802
1958-05-01    386.308495
1958-06-01    454.333395
1958-07-01    498.951424
1958-08-01    496.612558
1958-09-01    435.926648
1958-10-01    378.738250
1958-11-01    337.466764
1958-12-01    371.632091
1959-01-01    382.374626
1959-02-01    369.158449
1959-03-01    422.885737
1959-04-01    416.170397
1959-05-01    423.823140
1959-06-01    495.730964
1959-07-01    540.706571
1959-08-01    540.260123
1959-09-01    476.339117
1959-10-01    416.485983
1959-11-01    374.181398
1959-12-01    408.115702
1960-01-01    419.362032
1960-02-01    405.994823
1960-03-01    461.069570
1960-04-01    454.076520
1960-05-01    462.238681
1960-06-01    536.410972
1960-07-01    581.685471
1960-08-01    582.393744
1960-09-01    516.768118
1960-10-01    455.528423
1960-11-01    412.747158
1960-12-01    446.652061

予測した結果をプロットしてみます。

# 正解データと予測結果をプロット
plt.plot(df['Passengers'])
plt.plot(pred_original, "r")

青い線が正解データ。赤い線が実際に予測した部分です。

多少上ぶれている部分はあるものの、そこそこトレンドや季節性を掴めていることが確認できます。

精度検証

予測が出来たので、精度を検証してみます。

時系列データに使える精度指標としては以下になります。

  • RMSE
  • MAPE

RMSE（Root Mean Square Error, 二乗平均平方根誤差）

RMSE は、予測モデルの予測と実際の観測値との間での誤差を評価するための指標です。RMSE が小さいほど予測が実際のデータに近いことを示します。

\displaystyle RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}

n はデータポイントの数、 y_i は実際の値、 \hat{y}_i は予測値です。

scikit-learn の mean_squared_error 関数を利用します。

sklearn.metrics.mean_squared_error

from sklearn.metrics import mean_squared_error

test_original = df['Passengers']['1958-01-01':]

np.sqrt(mean_squared_error(test_original, pred_original))

# => RMSE: 21.32544528858937

RMSE は 21.3 でした。これは、予測値と実際の値との平均誤差がおおよそ 21.3 であることを意味します。

300 〜 600 くらいの予測なのでこの平均誤差が多いか少ないかは要件によるところもありますが、まだ改善はできそうです。これが一桁の平均誤差ならかなり精度の高い予測ができていると言えそうです。

MAPE（Mean Absolute Percentage Error, 平均絶対パーセント誤差）

MAPE は、予測モデルの予測値と実際の観測値との間での誤差を評価するための指標です。予測の正確さをパーセントで示します。MAPEの値が小さいほど予測が正確であることを示し、大きいほど予測の誤差が大きいことを示します。

\displaystyle MAPE = \frac{1}{n} \sum_{i=1}^{n} \left| \frac{y_i - \hat{y}_i}{y_i} \right| \times 100

n はデータポイントの数、 y_i は実際の値、 \hat{y}_i は予測値です。

np.mean(np.abs((pred_original - test_original) / test_original))

# => MAPE: 0.045705066956362256

MAPE 0.0457 でした。これは予測モデルの予測が実際の観測値と比べて、平均して約 4.57 ％の誤差があることを示しています。つまり、予測値が実際の値から約 4.57 ％程度離れていることが平均的な誤差として示されています。

値としては小さいため、モデルの予測が相対的に実際の観測値に近いと判断できそうです。

モデルの改善

ここまでで、「データの確認・理解」「定常データへの変換」「モデル構築」と一連の時系列分析の流れを実施しました。

あとは作成したモデルの精度を上げてく工程がありますが、本記事ではここまでとします。

「モデルのチューニング・パラメータ調整」「データの変換」など、アプローチは色々とあるため、次の機会に実施したいと思います。

トレンドや季節性成分の抽出

最後に、原系列データからトレンドや季節性成分を抽出してみます。

python の Pmdarima というパッケージを利用します。

alkaline-ml/pmdarima

Pmdarima（Pyramid ARIMA）は、Pythonで時系列データの予測モデリングを行うためのツールキットです。 ARIMA, SARIMA モデルのパラメーター推定と予測モデリングのプロセスを簡素化し、ユーザーが容易に時系列データを分析できるように支援してくれます。

Pmdarima パッケージを利用すると、モデル作成のためのパラメータ（次数）を自動で選択してくれたり、モデルの適合度を評価・診断プロットを生成する機能があったりと、今回実施した手順を簡略化してモデル構築が行えたりします。

Pmdarima パッケージを利用したモデル構築は次の機会に試すとして、ここではトレンドや季節性成分の抽出を行ってみます。

# データセット読み込み
column_names = ['Month', 'Passengers']
df=pd.read_csv('AirPassengers.csv', index_col='Month', parse_dates=True, names=column_names, header=None, skiprows=1)

# df.head()
Month   Passengers
1949-01-01  112
1949-02-01  118
1949-03-01  132
1949-04-01  129
1949-05-01  121

この原系列データからトレンドや季節性成分を抽出します。

decompose 関数で成分の分解を行い、decomposed_plot 関数でそれぞれの結果をプロットします。

from pmdarima import utils
from pmdarima import arima

data = df['Passengers'].values

utils.decomposed_plot(
    arima.decompose(data,'additive',m=12), figure_kwargs = {'figsize': (16, 12)} 
)

  • data
    • 元のデータ
  • trend
    • トレンド
  • seasonal
    • 季節性成分
  • random
    • 残差成分
    • 残差成分はトレンドと季節性を除いた残りの部分を表現しています。通常、残差はランダムであるべきで、特定のパターンや構造がないことが望まれます。

トレンドと季節性成分を抽出できました。

Pmdarima パッケージも便利そうなので次の機会に深掘ってみたいと思います。

現在 back check 開発チームでは一緒に働く仲間を募集中です。 herp.careers