自分のキャリアをあれこれ考えながら、Pythonで様々なデータを分析していくブログです

(その3-5) アップル引越しの需要予測を自己回帰移動平均(ARMA)モデルでやってみた

Data Analytics
Data Analytics

前回MAモデルの記事を書きました。

(その3-4) アップル引越しの需要予測を移動平均(MA)モデルでやってみた
前回ARモデルの記事を書きました。 結果は1年分の日次予測は難しいかなという結果でした。 今回はMAモデルを試してみます。 時系列モデル一覧 AR (自己回帰モデル) MA (移動平均モデル) ARMA (自己回帰移動平均モデル) ARIM...

結果はARモデルと同じく1年分の日次予測は難しいかなという結果でした。

今回はARMAモデルを試してみます。

時系列モデル一覧

  1. AR (自己回帰モデル)
  2. MA (移動平均モデル)
  3. ARMA (自己回帰移動平均モデル)
  4. ARIMA (自己回帰和分移動平均モデル)
  5. ARIMAX (自己回帰和分移動平均モデル with 外因性変数)
  6. Auto ARIMA (自動ARIMAモデル)
  7. SARIMA (季節ARIMAモデル)
  8. SARIMAX (季節ARIMAモデル with 外因性変数)
  9. ARCH (分散自己回帰モデル)
  10. GARCH (一般化分散自己回帰モデル)
  11. VAR (ベクトル自己回帰モデル)
  12. VARMA (ベクトル自己回帰移動平均モデル)

※ ARMAモデルについてはWikipedia参照

自己回帰移動平均モデル (…) は自己回帰モデルによる線形フィードバックと移動平均モデルによる線形フィードフォワードによりシステムを表現するモデルである。 (…) ARMAモデルは時系列データの将来値を予測するツールとして機能する。
引用: https://ja.wikipedia.org/wiki/自己回帰移動平均モデル

名前の通りARモデルとMAモデル両方を取り入れたモデルのようです。

スポンサーリンク

アップル引越しの引越し数をARMAモデルで予想する

モデル作成するパートまではARモデルと共通になります。

本記事では重複した内容になってしまうので、詳細な説明を省きます。

データの読み込み

# アップル引越しのデータセットを読み込む
import pandas as pd
from matplotlib import pyplot as plt
import pandas as pd
# 訓練データと予測付与用データの読み込み
df = pd.read_csv("/Users/hinomaruc/Desktop/blog/dataset/applehikkoshi/train.csv",parse_dates=['datetime'],index_col='datetime')
df_test = pd.read_csv("/Users/hinomaruc/Desktop/blog/dataset/applehikkoshi/test.csv",parse_dates=['datetime'],index_col='datetime')
# 描画設定
from IPython.display import HTML
import seaborn as sns
from matplotlib import ticker
import matplotlib.pyplot as plt
sns.set_style("whitegrid")
from matplotlib import rcParams
rcParams['font.family'] = 'Hiragino Sans' # Macの場合
#rcParams['font.family'] = 'Meiryo' # Windowsの場合
#rcParams['font.family'] = 'VL PGothic' # Linuxの場合
rcParams['xtick.labelsize'] = 12       # x軸のラベルのフォントサイズ
rcParams['ytick.labelsize'] = 12       # y軸のラベルのフォントサイズ
rcParams['axes.labelsize'] = 18        # ラベルのフォントとサイズ
rcParams['figure.figsize'] = 18,8      # 画像サイズの変更(inch)

日付の間隔をasfreqメソッドで変更する

asfreqメソッドはpandasのメソッドで日付の粒度を変更してくれたり、足りない日付も作成してくれる便利なメソッドです。

# https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.asfreq.html
# 日付の間隔を「日間隔」に変更 (元から日別のデータですが、、週や月間隔にも変更できます)
df = df.asfreq('d')
df.info()
Out[0]
DatetimeIndex: 2101 entries, 2010-07-01 to 2016-03-31
Freq: D
Data columns (total 5 columns):
 #   Column    Non-Null Count  Dtype
---  ------    --------------  -----
 0   y         2101 non-null   int64
 1   client    2101 non-null   int64
 2   close     2101 non-null   int64
 3   price_am  2101 non-null   int64
 4   price_pm  2101 non-null   int64
dtypes: int64(5)
memory usage: 98.5 KB

欠損値の確認

df.infoメソッドで確認した通り、欠損値は今回存在しませんのでスキップします。

引っ越し数(y)を描画

pandasのplotメソッドで描画できます。

df.y.plot()
Out[0]

yが正規分布に従うかどうかQQプロットを確認

統計手法を適用するにあたり、おおよそ正規分布を仮定する必要があります。

そのためデータが正規分布に従っているかどうかQQプロットという手法を使って確認します。

# 正規分布に従うかどうかをQQプロットを描画して確認
import scipy.stats as stats
import matplotlib.pyplot as plt

plt.subplot(1,2,1)
sns.histplot(df["y"], kde=True)
plt.subplot(1,2,2)
stats.probplot(df["y"], dist="norm", plot=plt)
plt.show()
Out[0]

左側のグラフがヒストグラム、右側のグラフがQQプロットです。

青い点が赤い線に沿っていれば正規分布であるという表現方法になっています。

端の部分がずれているようですが、なんとなく正規分布に従っていると仮定します。

Mann-Kendall trend検定

H0: データにトレンドは存在しない
HA: データにトレンドが存在する

# https://pypi.org/project/pymannkendall/
import pymannkendall as mk
mk.original_test(df["y"])
Out[0]
Mann_Kendall_Test(trend='increasing', h=True, p=0.0, z=30.51302505819748, Tau=0.4440819564379774, s=979667.0, var_s=1030826418.3333334, slope=0.016624040920716114, intercept=14.544757033248079)

p <= 0.05なので帰無仮説は棄却され、データにトレンドは存在するという結果になりました。

そのため、拡張ディッキー・フラー検定ではトレンドがあるということを示すregression='ct'のオプションを追加しようと思います。

拡張ディッキー・フラー検定で時系列データが(弱)定常か非定常か確認する

帰無仮説と対立仮説は下記になります。

H0: 単位根がある
HA: 単位根がない

# https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.adfuller.html
from statsmodels.tsa.stattools import adfuller

# 検定結果を見やすく加工
result = adfuller(df["y"],regression='ct')
labels = ["検定統計量","P値","#ラグ数","#観測数"]
mod_result = pd.Series(result[0:4],index=labels)
for key,val in result[4].items():
    mod_result[f'臨界値 (有意水準:{key})'] = val;
mod_result
Out[0]
検定統計量               -3.139726
P値                   0.097095
#ラグ数                26.000000
#観測数              2074.000000
臨界値 (有意水準:1%)       -3.963142
臨界値 (有意水準:5%)       -3.412609
臨界値 (有意水準:10%)      -3.128298
dtype: float64

P値が0.05以下ではないため帰無仮説を棄却できず「単位根がないとは言えない」という結果になります。

つまりアップル引っ越しのデータセットは弱定常過程ではなく、非定常過程の可能性がありそうです。

非定常過程のデータを定常過程のデータに変換する (階差)

拡張ディッキー・フラー検定で単位根がないとは言えないという結果になったので、本データセットは非定常過程のデータであると思われます。従って階差によって定常過程のデータに変換してモデルを作成します。

# https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.diff.html
# 1階差分を取る
df_diff1 = df.diff(periods=1).dropna()
df_diff1.head()
Out[0]

y client close price_am price_pm
datetime
2010-07-02 1.0 0.0 0.0 0.0 0.0
2010-07-03 2.0 0.0 0.0 0.0 0.0
2010-07-04 0.0 0.0 0.0 0.0 0.0
2010-07-05 -6.0 0.0 0.0 0.0 0.0
2010-07-06 0.0 0.0 0.0 0.0 0.0
df_diff1.y.plot()
Out[0]

定常過程っぽいデータになりましたね

Mann-Kendall trend検定と拡張ディッキー・フラー検定を定常化したデータで確認する

下記記事と同じ内容になりますので、興味ある方はご確認ください。

(その3-3) アップル引越しの需要予測を自己回帰(AR)モデルでやってみた
やっと22年9月の最初の記事を投稿できます。 2週間くらい時系列モデルについて勉強していました。以前も取り組んだことはあるのですが、ARIMAモデルを動かした程度で終わってしまったので今回はもう少し時間をかけてみました。 定常性とかホワイト...

まとめると、
・トレンドはなさそう
・(弱)定常過程でありそう
という結果になります。

自己相関係数(ACF)と偏自己相関係数(PACF)の確認

ACFはMAのラグ数、PACFはARのラグ数の指針になるようです。

どちらも自己相関があると思われるラグの値を初期設定としてARMAモデルを構築して調整していくと良いようです。

# 自己相関係数の確認
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
import numpy as np

fig, ax = plt.subplots(figsize=(16,8))
plot_acf(df_diff1['y'], lags=50, zero=False, title="自己相関係数 (autocorrelation)", ax=ax, alpha=0.01)
plt.ylim([0,1])
plt.yticks(np.arange(0.1, 1.1, 0.1))
plt.xticks(np.arange(1, 51, 1))
plt.axhline(y=0.5, color="green")
plt.show()

fig2, ax2 = plt.subplots(figsize=(16,8))
plot_pacf(df_diff1['y'], lags=50, zero=False, title="偏自己相関係数 (partial autocorrelation)", method=('ywm'), ax=ax2, alpha=0.01)
plt.ylim([0,1])
plt.yticks(np.arange(0.1, 1.1, 0.1))
plt.xticks(np.arange(1, 51, 1))
plt.axhline(y=0.5, color="green")
plt.show()
Out[0]


AR=[7,8,35,36]、MA=[7,14,21,28,35,42,49]の組み合わせを試して一番良いモデルのorderを算出した後、係数のP値を確認して妥当なモデルか確認していくという作業になっていくのかなと思います。

ARMAモデルの作成

まずは、ARMA(1,1)のモデルyt=c+ϕ1yt-11εt-1tを作成していこうと思います。

式の中身を分解すると下記になります。

・cはコンスタント値
・ytは現在の値
・yt-1は1ピリオド前の値
・εtは現在の残差(エラー)
・εt-1は1ピリオド前の残差(エラー)
・ϕは現在値を説明するのに過去の値がどれほど関連しているか (ARの部分)
・θは現在値を説明するのに過去の残差(エラー)がどれほど関連しているか (MAの部分)

ARモデルとMAモデルは対数尤度比検定(Log likelihood ratio test)でモデルん優劣を判断していましたが、簡略化のため、log likelihooodの値とAICの値を判断材料にしようと思います。log likelihooodの値は高いほどいいモデル、AICの値は低いほどいいモデルという判別になります。

それではARMAモデルの作成をしていきたいと思います。

**from statsmodels.tsa.arima_model import ARMA**というクラスは本記事で使用しているstatsmodels==0.13.2では実装から消されているので使えません。使おうとすると、NotImplementedErrorが発生するのでご注意ください。代わりに**statsmodels.tsa.arima.model.ARIMA**をお使いください。
from statsmodels.tsa.arima.model import ARIMA
# ARMA(1,1)モデルの作成
ar1ma1_model = ARIMA(df_diff1.y,order=(1,0,1))
ar1ma1_model_fit = ar1ma1_model.fit()
ar1ma1_model_fit.summary()
Out[0]

SARIMAX Results
Dep. Variable: y No. Observations: 2100
Model: ARIMA(1, 0, 1) Log Likelihood -7536.590
Date: Wed, 21 Sep 2022 AIC 15081.180
Time: 16:19:12 BIC 15103.779
Sample: 07-02-2010 HQIC 15089.458
- 03-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
const 0.0349 0.038 0.921 0.357 -0.039 0.109
ar.L1 0.6381 0.020 31.392 0.000 0.598 0.678
ma.L1 -0.9284 0.011 -87.638 0.000 -0.949 -0.908
sigma2 76.6691 1.784 42.974 0.000 73.172 80.166
Ljung-Box (L1) (Q): 0.28 Jarque-Bera (JB): 217.27
Prob(Q): 0.60 Prob(JB): 0.00
Heteroskedasticity (H): 2.53 Skew: -0.08
Prob(H) (two-sided): 0.00 Kurtosis: 4.57


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

ar.L1とma.L1のp値は0.0で意味のある係数のようです。AR係数の0.6381プラスの値なので現在と過去の値に正の傾向があるという意味で、MA係数の-0.9284はマイナスの値なので過去の値から離れるという意味のようです。

ARMA(2,2)モデルを作成してARMA(1,1)の結果と比較する

お試しでARMA(2,2)モデルを作成して、ARMA(1,1)とどちらのモデルが良いモデルか比べてみます。

# ARMA(2,2)モデルの作成
ar2ma2_model = ARIMA(df_diff1.y,order=(2,0,2))
ar2ma2_model_fit = ar2ma2_model.fit()
ar2ma2_model_fit.summary()
Out[0]

SARIMAX Results
Dep. Variable: y No. Observations: 2100
Model: ARIMA(2, 0, 2) Log Likelihood -7535.822
Date: Wed, 21 Sep 2022 AIC 15083.643
Time: 16:19:22 BIC 15117.541
Sample: 07-02-2010 HQIC 15096.059
- 03-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
const 0.0345 0.037 0.936 0.349 -0.038 0.107
ar.L1 0.2290 0.647 0.354 0.723 -1.039 1.497
ar.L2 0.2834 0.408 0.695 0.487 -0.516 1.083
ma.L1 -0.5335 0.650 -0.820 0.412 -1.808 0.741
ma.L2 -0.3729 0.602 -0.619 0.536 -1.554 0.808
sigma2 76.6121 1.781 43.007 0.000 73.121 80.104
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 219.50
Prob(Q): 0.88 Prob(JB): 0.00
Heteroskedasticity (H): 2.53 Skew: -0.08
Prob(H) (two-sided): 0.00 Kurtosis: 4.58


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
# log likelihoodとAICの値を確認
print(f"ARMA(1,1)\tlog likelihood: {ar1ma1_model_fit.llf}\tAIC: {ar1ma1_model_fit.aic}")
print(f"ARMA(2,2)\tlog likelihood: {ar2ma2_model_fit.llf}\tAIC: {ar2ma2_model_fit.aic}")
ARMA(1,1)   log likelihood: -7536.590178155023  AIC: 15081.180356310046
ARMA(2,2)   log likelihood: -7535.821547414062  AIC: 15083.643094828123

log likelihoodの値は ARMA(2.2)の方が高いのですが、AICの値も高くなっているのでARMA(1,1)より良いモデルとは言えないという結果になりました。

確認方法がつかめたのでARMAモデルでARとMAそれぞれの次数の値を探っていこうと思います。

最適な次数のMAモデルを求める

# Log likelihood ratio testの関数作成
"""
Log likelihood ratio test
H0: 2つのモデルの比に有意差はない
HA: 2つのモデルの比に有意差がある
引数1: モデル1のLog Likelihood
引数2: モデル2のLog Likelihood
引数3: 自由度
リターン: p値
"""
def llr_test(llf1,llf2,dof):
    LL1 = llf1
    LL2 = llf2
    LLR = 2 * (LL2-LL1)

    from scipy.stats.distributions import chi2
    p = chi2.sf(LLR,dof).round(3)

    return p
# ACFとPACFのグラフを確認した結果のとおり、
# まずはAR=[7,8,35,36]、MA=[7,14,21,28,35,42,49]の組み合わせを試して良さそうなスタートラインを決める

from statsmodels.tsa.arima.model import ARIMA
import datetime

best_arma_model_fit = None

for ar in [7,8,35,36]:
    for ma in [7,14,21,28,35,42,49]:   
        # ARMAモデル作成
        arma_model = ARIMA(df_diff1["y"],order=(ar,0,ma))
        arma_model_fit = arma_model.fit()

        # ARMA(7,7)は自動的にベストモデルになる
        if (ar == 7) and  (ma==7):
            best_arma_model_fit = arma_model_fit
        else:
        # ARMA(7,7) 以降の処理
            llf1 = best_arma_model_fit.llf
            aic1 = best_arma_model_fit.aic
            llf2 = arma_model_fit.llf
            aic2 = arma_model_fit.aic

            # llfは高いほど良く、AICは低いほど良い
            if llf2 > llf1 and aic2 < aic1:
                #条件に合致すれば、ベストモデルを更新
                best_arma_model_fit = arma_model_fit
                print(f"{datetime.datetime.now()}: best model is ARMA({ar,ma}) -> Log Likelihood:{best_arma_model_fit.llf}、AIC:{best_arma_model_fit.aic}")
                print("-----")

途中warningが出てしいるのは問題ないようです。
ちなみにConvergenceWarning: Maximum Likelihood optimization failed to convergeはモデルの当てはまりが悪いなどでパラメーターの探索がうまくいってない場合に発生するようです。

This message suggests that the optimizer is having a hard time determining good values for all of the parameters. There are a number of possible reasons why this could happen, but one reason is if the model is not a good fit for your data.
引用: https://github.com/statsmodels/statsmodels/issues/8314

Out[0]

※ warningはアウトプットから削っています。かなり出るのでwarningsをoffにすればよかったかも

2022-09-21 17:08:52.092158: best model is ARMA((7, 14)) -> Log Likelihood:-7326.532247414697、AIC:14699.064494829394
2022-09-21 17:09:37.148874: best model is ARMA((7, 21)) -> Log Likelihood:-7310.3116106352245、AIC:14680.623221270449
2022-09-21 17:11:00.768425: best model is ARMA((7, 28)) -> Log Likelihood:-7302.741788220932、AIC:14679.483576441864
2022-09-21 17:13:16.580240: best model is ARMA((7, 35)) -> Log Likelihood:-7294.690771513018、AIC:14677.381543026037
2022-09-21 17:17:43.116703: best model is ARMA((7, 42)) -> Log Likelihood:-7282.795994795316、AIC:14667.591989590632
2022-09-21 18:06:28.780062: best model is ARMA((35, 14)) -> Log Likelihood:-7272.621157888031、AIC:14647.242315776062
2022-09-21 18:10:09.148433: best model is ARMA((35, 21)) -> Log Likelihood:-7244.432971703251、AIC:14604.865943406501
2022-09-21 18:14:13.904902: best model is ARMA((35, 28)) -> Log Likelihood:-7225.9263071699925、AIC:14581.852614339985

ARMA(35, 28)が選択されました。

# ベストモデルのサマリ情報の確認
best_arma_model_fit.summary()
Out[0]

SARIMAX Results
Dep. Variable: y No. Observations: 2100
Model: ARIMA(35, 0, 28) Log Likelihood -7225.926
Date: Wed, 21 Sep 2022 AIC 14581.853
Time: 21:34:58 BIC 14949.083
Sample: 07-02-2010 HQIC 14716.359
- 03-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
const 0.0198 0.010 2.034 0.042 0.001 0.039
ar.L1 -0.7257 0.242 -2.999 0.003 -1.200 -0.251
ar.L2 -0.3908 0.184 -2.123 0.034 -0.752 -0.030
ar.L3 0.1043 0.139 0.752 0.452 -0.168 0.376
ar.L4 0.3029 0.131 2.308 0.021 0.046 0.560
ar.L5 0.0886 0.146 0.607 0.544 -0.198 0.375
ar.L6 -0.2113 0.128 -1.653 0.098 -0.462 0.039
ar.L7 -0.3062 0.132 -2.315 0.021 -0.565 -0.047
ar.L8 -0.0558 0.109 -0.514 0.607 -0.268 0.157
ar.L9 0.0316 0.103 0.307 0.759 -0.170 0.233
ar.L10 0.1369 0.103 1.335 0.182 -0.064 0.338
ar.L11 -0.0662 0.088 -0.755 0.450 -0.238 0.106
ar.L12 -0.1526 0.087 -1.746 0.081 -0.324 0.019
ar.L13 -0.2797 0.097 -2.891 0.004 -0.469 -0.090
ar.L14 -0.0338 0.101 -0.335 0.738 -0.232 0.164
ar.L15 0.2545 0.103 2.460 0.014 0.052 0.457
ar.L16 0.0302 0.099 0.306 0.759 -0.163 0.223
ar.L17 -0.0639 0.096 -0.668 0.504 -0.251 0.124
ar.L18 -0.2908 0.087 -3.334 0.001 -0.462 -0.120
ar.L19 -0.2082 0.121 -1.718 0.086 -0.446 0.029
ar.L20 -0.0355 0.108 -0.328 0.743 -0.248 0.177
ar.L21 0.3376 0.102 3.298 0.001 0.137 0.538
ar.L22 0.2314 0.122 1.902 0.057 -0.007 0.470
ar.L23 -0.1042 0.110 -0.943 0.345 -0.321 0.112
ar.L24 -0.5207 0.120 -4.340 0.000 -0.756 -0.286
ar.L25 -0.4876 0.134 -3.641 0.000 -0.750 -0.225
ar.L26 -0.0731 0.125 -0.584 0.559 -0.318 0.172
ar.L27 0.3077 0.108 2.844 0.004 0.096 0.520
ar.L28 0.6300 0.138 4.565 0.000 0.360 0.901
ar.L29 0.1812 0.050 3.656 0.000 0.084 0.278
ar.L30 0.1631 0.040 4.109 0.000 0.085 0.241
ar.L31 0.1870 0.045 4.143 0.000 0.099 0.275
ar.L32 0.2667 0.053 5.038 0.000 0.163 0.370
ar.L33 0.1455 0.054 2.691 0.007 0.039 0.251
ar.L34 0.0781 0.042 1.864 0.062 -0.004 0.160
ar.L35 0.0217 0.036 0.597 0.551 -0.049 0.093
ma.L1 0.2575 0.242 1.065 0.287 -0.217 0.732
ma.L2 -0.0341 0.212 -0.161 0.872 -0.450 0.382
ma.L3 -0.3813 0.187 -2.043 0.041 -0.747 -0.015
ma.L4 -0.4205 0.162 -2.588 0.010 -0.739 -0.102
ma.L5 -0.0920 0.170 -0.543 0.587 -0.424 0.240
ma.L6 0.2480 0.155 1.603 0.109 -0.055 0.551
ma.L7 0.3390 0.155 2.190 0.029 0.036 0.642
ma.L8 -0.0253 0.127 -0.199 0.842 -0.274 0.224
ma.L9 -0.1316 0.124 -1.058 0.290 -0.375 0.112
ma.L10 -0.2298 0.130 -1.774 0.076 -0.484 0.024
ma.L11 0.0072 0.106 0.068 0.946 -0.201 0.215
ma.L12 0.0486 0.108 0.450 0.652 -0.163 0.260
ma.L13 0.2452 0.112 2.194 0.028 0.026 0.464
ma.L14 -0.0268 0.115 -0.233 0.816 -0.252 0.198
ma.L15 -0.3096 0.121 -2.557 0.011 -0.547 -0.072
ma.L16 0.0282 0.113 0.249 0.803 -0.194 0.250
ma.L17 -0.0001 0.113 -0.001 0.999 -0.221 0.221
ma.L18 0.2253 0.108 2.088 0.037 0.014 0.437
ma.L19 0.1103 0.139 0.795 0.427 -0.162 0.382
ma.L20 -0.0456 0.130 -0.352 0.725 -0.300 0.209
ma.L21 -0.3466 0.120 -2.885 0.004 -0.582 -0.111
ma.L22 -0.1843 0.133 -1.389 0.165 -0.444 0.076
ma.L23 0.1885 0.121 1.561 0.119 -0.048 0.425
ma.L24 0.4590 0.137 3.359 0.001 0.191 0.727
ma.L25 0.2895 0.133 2.184 0.029 0.030 0.549
ma.L26 -0.1433 0.125 -1.146 0.252 -0.388 0.102
ma.L27 -0.3896 0.128 -3.041 0.002 -0.641 -0.139
ma.L28 -0.6192 0.155 -3.996 0.000 -0.923 -0.315
sigma2 61.5569 1.735 35.476 0.000 58.156 64.958
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 324.95
Prob(Q): 0.92 Prob(JB): 0.00
Heteroskedasticity (H): 2.29 Skew: -0.23
Prob(H) (two-sided): 0.00 Kurtosis: 4.87


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

・AR係数のar.L01~ ar.L35のうち18個が P>|z| が 0.05以上で意味をなさない係数
・MA係数のma.L1 ~ ma.L28のうち17個が P>|z| が 0.05以上で意味をなさない係数
でした。

これはもっとシンプルなモデルにできることを示唆しているようですが、今回はこのまま進めていこうと思います。

作成したARMAモデルの残差分析

残差の正規性はQQプロットで確認することにします。残差の間に相関がない(独立している)ことはLjung-Box検定で確認します。

resid = best_arma_model_fit.resid
# ARMA(35,28)の残差データを確認
resid.plot()
Out[0]

残差の平均と分散を確認

print(f"mean:{resid.mean()}、variance:{resid.var()}")
Out[0]
mean:0.0003919578747833204、variance:56.7044965544352

残差の正規性の確認

# 正規分布に従うかどうかをQQプロットを描画して確認
import scipy.stats as stats
import matplotlib.pyplot as plt

plt.subplot(1,2,1)
sns.histplot(resid, kde=True)
plt.subplot(1,2,2)
stats.probplot(resid, dist="norm", plot=plt)
plt.show()
Out[0]

概ね正規分布に従っていそうです。

自己相関係数の確認

# 自己相関係数の確認
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
import numpy as np

fig, ax = plt.subplots(figsize=(16,8))
plot_acf(resid, lags=50, zero=False, title="自己相関係数 (autocorrelation)", ax=ax, alpha=0.01)
plt.ylim([0,1])
plt.yticks(np.arange(0.1, 1.1, 0.1))
plt.xticks(np.arange(1, 51, 1))
plt.show()
Out[0]

全てのラグが青いエリア内なので問題なさそうです。

Ljung-Box検定で残差の独立性の確認

H0: 残差は独立分布している
HA: 残差は独立分布ではない

残差が独立分布していることが望ましいので帰無仮説を棄却しない(P値 >= 0.05)結果になることを期待したいと思います。

# https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.acorr_ljungbox.html

import statsmodels.api as sm
sm.stats.acorr_ljungbox(resid,lags=50)
Out[0]

lb_stat lb_pvalue
1 0.011686 0.913915
2 0.143040 0.930978
3 0.150921 0.985094
4 0.417733 0.980999
5 0.693918 0.983301
6 0.702819 0.994430
7 1.164911 0.991712
8 2.749300 0.949092
9 3.577558 0.936958
10 5.307432 0.869718
11 5.581831 0.899760
12 6.452172 0.891597
13 7.408150 0.879912
14 8.243099 0.876280
15 8.316958 0.910448
16 9.234596 0.903439
17 10.899463 0.861742
18 10.910356 0.898122
19 10.973101 0.924740
20 11.404043 0.935064
21 11.439417 0.953532
22 11.537977 0.966201
23 14.138049 0.922850
24 14.435594 0.936193
25 14.504890 0.952202
26 14.655830 0.963308
27 14.657846 0.973962
28 15.025384 0.978166
29 15.704724 0.978705
30 15.739981 0.984838
31 16.161825 0.986957
32 16.590876 0.988737
33 16.630544 0.992082
34 16.779439 0.994146
35 18.892458 0.987967
36 19.098703 0.990629
37 19.106045 0.993429
38 19.252047 0.995096
39 19.294917 0.996561
40 19.645050 0.997157
41 20.069758 0.997551
42 20.390804 0.998010
43 21.062967 0.998028
44 21.097238 0.998639
45 23.134654 0.997179
46 23.134671 0.998044
47 26.216152 0.993923
48 27.987032 0.990720
49 28.023933 0.993058
50 28.571382 0.993617

1~50のラグで検定してみました。全てp値が0.05以上なので帰無仮説を棄却できず残差は独立分布ではないとは言えないという結果になりましたので問題なさそうです。

残差分析の結果

残差分析の結果は、ARMA(35,28)の残差が平均0の正規分布に従い、かつ残差の間には相関がありませんでしたので問題ありませんでした。

ARMA(35,28)モデルを使って引っ越し数を予測する

ARモデルとMAモデルの時は長期間の予測は難しそうでしたが、ARMAモデルだとどうなるのでしょうか?楽しみです。

# 2016-04-01 ~ 2017-03-31までの階差を予測
best_arma_model_fit.predict(start = df_test.index[0],end = df_test.index[-1]).plot()
Out[0]

ARモデルのときやMAモデルのときよりは改善されました。

階差の予測値なので引っ越し数の予測値に変換します

# 学習データの一番最後の引っ越し数を確認
df.iloc[-1]
Out[0]
y           105
client        1
close         0
price_am      5
price_pm      4
Name: 2016-03-31 00:00:00, dtype: int64

yの値を使います。

predict_diff = best_arma_model_fit.predict(start = df_test.index[0],end = df_test.index[-1])
predict_diff
Out[0]
2016-04-01    -6.041302
2016-04-02    -2.132627
2016-04-03    -8.123813
2016-04-04   -10.215175
2016-04-05    -4.719134
                ...    
2017-03-27    -3.798514
2017-03-28    -5.038416
2017-03-29    -2.744863
2017-03-30     1.081610
2017-03-31     3.390011
Freq: D, Name: predicted_mean, Length: 365, dtype: float64
# 残差予測から引っ越し数の予測値を作成
y_pred = df.y.iloc[-1]
y_pred_list = []
for idx,val in enumerate(predict_diff):
    y_pred_list.append((y_pred + val))

# テストデータに予測結果を追加
df_test["y_pred"] = y_pred_list
# 学習データの引っ越し数の値(青い線)と予測した引っ越し数の値(オレンジの線)を描画
plt.figure(figsize=(20, 3))

ytrue = df.y.iloc[-31:]
ypred = df_test.y_pred

plt.plot(ytrue, label="Training Data")
plt.plot(ypred, label="Forecasts")

plt.title("apple-hikkoshi prediction")
_ = plt.legend()
Out[0]

39日目以降は横一直線になってしまい予測できていません。またしてもこれでは予測精度を確認しても意味がないので本記事ではここまでの確認とします。

(オプション) 残差の予測ではなく、引っ越し数をそのままARMAモデルに投入してみる

定常データにするために引っ越し数に対して1階差分をとっていましたが、そのまま引っ越し数で予測したらどうなるか試してみます。

とりあえず気にせず同じorderであるARMA(35,28)で作成してみます。

# 引っ越し数を変換かけずに投入し、ARMA(35,28)モデルを作成する
model = ARIMA(df.y,order=(35,0,28))
model_fit = model.fit()
model_fit.summary()
Out[0]

    /Users/hinomaruc/Desktop/blog/my-venv/lib/python3.8/site-packages/statsmodels/tsa/statespace/sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
      warn('Non-stationary starting autoregressive parameters'
    /Users/hinomaruc/Desktop/blog/my-venv/lib/python3.8/site-packages/statsmodels/tsa/statespace/sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
      warn('Non-invertible starting MA parameters found.'
    /Users/hinomaruc/Desktop/blog/my-venv/lib/python3.8/site-packages/statsmodels/base/model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
      warnings.warn("Maximum Likelihood optimization failed to "

SARIMAX Results
Dep. Variable: y No. Observations: 2101
Model: ARIMA(35, 0, 28) Log Likelihood -7253.589
Date: Thu, 22 Sep 2022 AIC 14637.178
Time: 09:11:05 BIC 15004.439
Sample: 07-01-2010 HQIC 14771.692
- 03-31-2016
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
const 34.0967 636.834 0.054 0.957 -1214.076 1282.269
ar.L1 -0.6410 1.682 -0.381 0.703 -3.939 2.657
ar.L2 -0.3892 0.816 -0.477 0.633 -1.988 1.209
ar.L3 0.0576 0.804 0.072 0.943 -1.518 1.633
ar.L4 0.2821 0.421 0.671 0.502 -0.542 1.106
ar.L5 0.1041 0.475 0.219 0.826 -0.826 1.034
ar.L6 -0.1403 0.389 -0.361 0.718 -0.902 0.622
ar.L7 -0.3071 0.325 -0.944 0.345 -0.945 0.331
ar.L8 -0.0925 0.622 -0.149 0.882 -1.311 1.126
ar.L9 0.0288 0.219 0.132 0.895 -0.400 0.458
ar.L10 0.2444 0.227 1.079 0.281 -0.200 0.688
ar.L11 0.0341 0.438 0.078 0.938 -0.825 0.893
ar.L12 -0.0736 0.170 -0.433 0.665 -0.406 0.259
ar.L13 -0.1077 0.177 -0.608 0.543 -0.455 0.239
ar.L14 0.1369 0.247 0.554 0.579 -0.347 0.621
ar.L15 0.3513 0.312 1.128 0.259 -0.259 0.962
ar.L16 0.0453 0.577 0.079 0.937 -1.086 1.176
ar.L17 -0.0137 0.181 -0.076 0.940 -0.368 0.340
ar.L18 -0.2680 0.179 -1.496 0.135 -0.619 0.083
ar.L19 -0.1476 0.510 -0.290 0.772 -1.147 0.852
ar.L20 0.1024 0.228 0.449 0.653 -0.344 0.549
ar.L21 0.4672 0.193 2.426 0.015 0.090 0.845
ar.L22 0.2978 0.819 0.364 0.716 -1.307 1.903
ar.L23 -0.0693 0.398 -0.174 0.862 -0.849 0.710
ar.L24 -0.4068 0.241 -1.686 0.092 -0.880 0.066
ar.L25 -0.3653 0.655 -0.558 0.577 -1.649 0.919
ar.L26 -0.0432 0.552 -0.078 0.938 -1.126 1.039
ar.L27 0.1642 0.222 0.739 0.460 -0.271 0.600
ar.L28 0.4254 0.308 1.381 0.167 -0.178 1.029
ar.L29 0.2692 0.759 0.355 0.723 -1.219 1.758
ar.L30 0.2674 0.345 0.774 0.439 -0.409 0.944
ar.L31 0.2097 0.397 0.529 0.597 -0.568 0.987
ar.L32 0.2526 0.291 0.868 0.385 -0.318 0.823
ar.L33 0.1426 0.401 0.356 0.722 -0.643 0.928
ar.L34 0.1248 0.192 0.649 0.517 -0.252 0.502
ar.L35 0.0533 0.207 0.258 0.796 -0.352 0.458
ma.L1 1.1709 1.683 0.696 0.487 -2.128 4.470
ma.L2 1.1645 1.638 0.711 0.477 -2.046 4.375
ma.L3 0.8395 1.867 0.450 0.653 -2.819 4.498
ma.L4 0.4640 1.151 0.403 0.687 -1.792 2.720
ma.L5 0.4160 0.729 0.571 0.568 -1.013 1.845
ma.L6 0.6375 0.543 1.173 0.241 -0.428 1.703
ma.L7 0.9822 0.908 1.082 0.279 -0.797 2.761
ma.L8 0.9301 1.557 0.597 0.550 -2.121 3.982
ma.L9 0.7559 1.349 0.560 0.575 -1.889 3.401
ma.L10 0.4198 1.046 0.401 0.688 -1.631 2.470
ma.L11 0.4221 0.542 0.779 0.436 -0.640 1.484
ma.L12 0.4562 0.647 0.705 0.481 -0.812 1.724
ma.L13 0.5480 0.619 0.885 0.376 -0.666 1.762
ma.L14 0.3911 0.835 0.468 0.640 -1.246 2.028
ma.L15 0.0503 0.542 0.093 0.926 -1.012 1.113
ma.L16 0.1468 0.154 0.953 0.341 -0.155 0.449
ma.L17 0.1667 0.286 0.582 0.560 -0.394 0.728
ma.L18 0.4283 0.277 1.548 0.122 -0.114 0.970
ma.L19 0.4695 0.727 0.646 0.519 -0.956 1.895
ma.L20 0.2890 0.657 0.440 0.660 -0.999 1.577
ma.L21 -0.1022 0.478 -0.214 0.830 -1.038 0.834
ma.L22 -0.2135 0.285 -0.750 0.453 -0.771 0.344
ma.L23 0.0909 0.333 0.273 0.785 -0.562 0.744
ma.L24 0.5241 0.269 1.949 0.051 -0.003 1.051
ma.L25 0.7622 0.831 0.917 0.359 -0.867 2.392
ma.L26 0.6583 1.161 0.567 0.571 -1.617 2.933
ma.L27 0.4631 1.037 0.447 0.655 -1.570 2.496
ma.L28 0.0751 0.723 0.104 0.917 -1.343 1.493
sigma2 61.6912 1.731 35.638 0.000 58.298 65.084
Ljung-Box (L1) (Q): 0.73 Jarque-Bera (JB): 302.75
Prob(Q): 0.39 Prob(JB): 0.00
Heteroskedasticity (H): 2.28 Skew: -0.22
Prob(H) (two-sided): 0.00 Kurtosis: 4.81


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

p値を確認するとほとんど0.05以上なので採用してはいけないモデルかも知れませんが、お試しなのでこのままやります。

# 学習データの引っ越し数の値(青い線)と予測した引っ越し数の値(オレンジの線)を描画
plt.figure(figsize=(20, 3))

ytrue = df.y.iloc[-31:]
ypred = model_fit.predict(start = df_test.index[0],end = df_test.index[-1])

plt.plot(ytrue, label="Training Data")
plt.plot(ypred, label="Forecasts")

plt.title("apple-hikkoshi prediction")
_ = plt.legend()
Out[0]

結果は対して変わらずでした。。

まとめ

ARMAモデルでも1年分の日別予測は難しかったです。

次はAuto ARIMAモデルを試してみたいと思います。

スポンサーリンク

本記事で利用したライブラリのバージョン

import pandas as pd
import numpy as np
import statsmodels as sm
import matplotlib as mpl
import pymannkendall as mk
import scipy as scp
print('pandas',pd.__version__)
print('numpy',np.__version__)
print('statsmodels',sm.__version__)
print('matplotlib',mpl.__version__)
print('pymannkendall',mk.__version__)
print('scipy',scp.__version__)
Out[0]
pandas 1.4.3
numpy 1.23.2
statsmodels 0.13.2
matplotlib 3.5.3
pymannkendall 1.4.2
scipy 1.9.1
タイトルとURLをコピーしました