売上高の時系列分析（２）

美容室経営研究所Refineの井上です。

今回は、前回（売上高の時系列分析（１））に続き、売上の時系列データに対して、ボックス・ジェンキンス法（Box–Jenkins method）を適用して分析してみようと思います。

データは私が自作したもので、前回は24か月分でしたが、今回は36か月分になっています！

ボックス・ジェンキンス法について

ボックス・ジェンキンス法は古典的な手法であり、ボックスさんとジェンキンスさんの2名のお名前からとった分析手法です。

簡単に言うと、自己回帰移動平均（ARMA）モデルや自己回帰和分移動平均（ARIMA）モデルを適用して、データに最も適合するものを求めるものです。

この方法は、主に次のステップで構成されています

①同定（Identification）

時系列データを分析して、データの特性を把握します。

この段階で、データが定常性（時間の経過による統計的特性の変化がないこと）を持っているかどうかを評価し、必要に応じて変換（差分取りなど）を行います。（ADF検定にて行います。）

拡張ディッキー–フラー検定（augmented Dickey–Fuller test, ADF test）は以下のようなものです。

統計学と計量経済学において、時系列標本が単位根を持つかどうかの仮説検定である。

これは大きくより複雑な時系列モデルに対するディッキー–フラー検定の拡張版となっている。

検定で用いられる拡張ディッキー–フラー検定統計量は負の値を取る。より大きな負の値を取れば取るほど、ある有意水準の下で単位根が存在するという仮説を棄却する可能性が強くなる。

Wikipedia

ADF（Augmented Dickey-Fuller）検定の式：

$$\Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \sum_{i=1}^{p} \phi_i \Delta y_{t-i} + \varepsilon_t$$

ただし、

• \(\Delta{y_t}\)：時刻\(t\)での時系列データの1階差分。これは\(y_{t} – y_{t−1}\)を意味します。
• \(\alpha\)：切片（またはドリフト）項。時系列データに線形トレンドがない場合は0になります。
• \(\beta_{t}\)：時間の経過に対する係数と時間\(t\)の積。時系列データに線形トレンドがある場合に使用されます。
• \(\gamma\)：ラグ1の時系列データの係数。これは、単位根過程（ユニットルートプロセス）の存在を検定するためのものです。
• \( \sum_{i=1}^{p} \phi_i \Delta y_{t-i} \)：過去の差分の合計。ここで\(p\)はラグの数で、検定の精度を高めるために含められます。
• \(\phi_i\)：ラグi の差分の係数。
• \(\epsilon_t\)：時刻tでの誤差項。これは、ホワイトノイズと仮定されています。

②推定（Estimation）

データに適合するモデルのパラメータを推定します。ARIMA（自己回帰積分移動平均）モデルが一般的に使用されます。

このモデルは、自己回帰（AR）項、差分（I）項、移動平均（MA）項を組み合わせて構成されています。

ARIMAモデルの式：

$$(1 – \sum_{i=1}^{p} \phi_i L^i) (1 – L)^d y_t = (1 + \sum_{j=1}^{q} \theta_j L^j) \varepsilon_t$$

ただし、

•  \(Yt\)は時刻 \(t\)での時系列データです。
•  \(L\)はバックシフト（ラグ）演算子で、 \(LY_t = Y_{t-1}\)です。
• \((1-L)^d\)は差分演算子で、時系列データの非定常性（トレンドや季節性など）を除去するために使用されます。
• \(\sum_{i=1}^{p} \phi_i B^i\)はAR項を表し、時系列の過去の値が現在の値にどのように影響するかを示します。
• \(\sum_{j=1}^{q} \theta_j B^j\)はMA項を表し、時系列の過去の誤差が現在の値にどのように影響するかを示します。
• \(\varepsilon_t\)は時刻 \(t\)での誤差項（ホワイトノイズ）です。

参考に自己回帰モデル（autoregressive model：AR model）と移動平均モデル（moving average model：MA model）の式も載せておきます。

AR(p)の式：

$$X_t = c + \sum_{i=1}^{p} \phi_i X_{t-i} + \varepsilon_t$$

ただし、

\(\phi_1, \ldots, \phi_p\)はモデルのパラメーターであり、\(c\)は定数項、\(\epsilon_t\) はホワイトノイズ。

MA(q)の式：

$$y_t = \mu + \sum_{k=0}^{q} \theta_k \varepsilon_{t-k} = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \ldots + \theta_q \varepsilon_{t-q}$$

ただし、

 \(\mu\)は定数、\(\theta_k\)はパラメータ (\(\theta_0\))、\(\epsilon_t\) は時刻\(t\)におけるホワイトノイズ。

③診断チェック（Diagnostic Checking）

推定されたモデルがデータに適しているかを評価します。

残差（モデルからの予測値と実際のデータの差）を分析し、モデルがデータの情報を適切に捉えているかをチェックします。

これには下記のLjung-Box検定やJarque-Bera検定を使用します。

Ljung-Box検定（リュング・ボックスけんてい、英: Ljung-Box test）は、ある時系列の自己相関の集まりが０と異なるかに関する統計的検定の一種である。

個々の時間差（ラグ）に関する無作為性を検定するのではなく、複数期間に対する時間差に基づく全般的な無作為性を検定することから、かばん検定である。

Wikipedia

ジャック＝ベラ検定（ジャック＝ベラけんてい、英: Jarque–Bera test）とは、統計学において標本データが正規分布に従う尖度と歪度を有しているかどうかを調べる適合度検定である。検定名はCarlos JarqueとAnil K. Beraにちなんで名づけられた。

Wikipedia

Ljung-Box検定の式：

$$Q = n(n+2) \sum_{k=1}^{m} \frac{\hat{\rho}_k^2}{n-k}$$

ただし、

\(\hat{\rho}_k\)は時間差\(k\) における標本自己相関、\(m\)は検定する時間差の数

Jarque-Bera検定の式：

$$JB = \frac{n}{6} \left( S^2 + \frac{1}{4} (K – 3)^2 \right)$$

ただし、

$$S = \frac{\frac{1}{n} \sum_{i=1}^{n} (X_i – \bar{X})^3}{\left(\frac{1}{n} \sum_{i=1}^{n} (X_i – \bar{X})^2\right)^{3/2}}$$

$$K = \frac{\frac{1}{n} \sum_{i=1}^{n} (X_i – \bar{X})^4}{\left(\frac{1}{n} \sum_{i=1}^{n} (X_i – \bar{X})^2\right)^2}$$

④予測（Forecasting）

モデルを用いて未来のデータ点を予測します。

では、見ていきましょう♪

実践

読み込み

まずはデータセットを読み込んでいこうと思います。

今回は、statsmodelsのtsa（時系列分析）から、holtwintersのExponentialSmoothingを使用しますので先にインポートしておきます。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

#読み込み
data = pd.read_excel('revenue36.xlsx')

data.head()

今回のデータは上記のようにいたってシンプルです。

美容室の月間売上高が36カ月分入っています。

可視化

データの様子を見ていこうと思います。

# 日付をインデックスに設定
data['月'] = pd.to_datetime(data['月'])
data.set_index('月', inplace=True)

plt.figure(figsize=(12, 6))
plt.plot(data['売上'])
plt.title('Monthly Revenue')
plt.xlabel('Month')
plt.ylabel('Revenue')
plt.grid(True)
plt.show()

どうでしょうか？

緩やかな上昇トレンドと、周期性が見られますね。

前回のデータに12か月分追加されて36か月分になっています。

周期は美容業界でよく繁忙月と言われる7月、12月、3月になるようにデータを作成しました。

この周期性をとらえた分析ができればいいなと思います。

コレログラム

import statsmodels.api as sm
#  自己相関のグラフ
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(data, lags=35, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(data, lags=12, ax=ax2)

ADF検定

次にADF検定で定常性の検定をしていきます。

from statsmodels.tsa.stattools import adfuller

# Perform Augmented Dickey-Fuller test
adf_test_result = adfuller(data['売上'])

# Extracting the p-value and test statistics
adf_statistic, adf_p_value = adf_test_result[0], adf_test_result[1]

adf_statistic, adf_p_value

p値0.88で棄却できませんでした。

1階の差分系列で再度検定してみましょう。

# データの１階差分をとる
data_diff = data['売上'].diff().dropna()

# Perform Augmented Dickey-Fuller test on the differenced data
adf_test_diff = adfuller(data_diff)

# Extracting the p-value and test statistics for the differenced data
adf_statistic_diff, adf_p_value_diff = adf_test_diff[0], adf_test_diff[1]

adf_statistic_diff, adf_p_value_diff

p値0.03で棄却できました！

というわけでARIMA(p,1,q)であるとして、最適なｐとｑを求めていきます。

from statsmodels.tsa.arima.model import ARIMA
import numpy as np
p_values = range(0, 5)
#ADF検定より
d_value = 1
q_values = range(0, 5)

#繰り返し演算
def calculate_aic_corrected(data, p_values, d_value, q_values):
    aic_values = []
    for p in p_values:
        for q in q_values:
            try:
                model = ARIMA(data, order=(p, d_value, q))
                model_fit = model.fit()
                aic_values.append((p, q, model_fit.aic))
            except Exception as e:
                aic_values.append((p, q, None))
    return aic_values

# AIC再計算
aic_results_corrected = calculate_aic_corrected(data['売上'], p_values, d_value, q_values)

# AICがNoneの場合
aic_results_filtered = [result for result in aic_results_corrected if result[2] is not None]

# AIC最小の組み合わせ
aic_results_sorted_corrected = sorted(aic_results_filtered, key=lambda x: x[2])
#トップ５出力
aic_results_sorted_corrected[:5] 

AICが最小になるモデルはARIMA(1,1,0)と分かりました。

ARIMAモデル学習

先ほど求めたモデルを学習させます。

# 学習
model = ARIMA(data['売上'], order=(1, 1, 0))
model_fit = model.fit()

# 結果表示
model_fit_summary = model_fit.summary()
model_fit_summary

結果、特にジャックベラ検定とリュングボックス検定についてみてみようと思います。

Ljung-Box検定の結果：
Q=3.20
p値：0.07
Ljung-Box検定のp値が0.07であるため、5%の有意水準で帰無仮説（残差に自己相関がない）を棄却するにはわずかに不足しています。

これは、残差に自己相関がない、つまりモデルがデータの自己相関を適切に捉えている可能性があることを示唆しています。ただし、このp値は10%の有意水準では帰無仮説を棄却する可能性があります。

Jarque-Bera検定の結果：
JB=2.38
p値：0.30
Jarque-Bera検定のp値が0.30であることから、帰無仮説（残差が正規分布に従う）を棄却できません。

これは、残差が正規分布に近いことを示しており、モデルの残差が正規分布の仮定に合致していると考えられます。

＜注意＞

多くの統計的検定では、研究者が立証しようとしている仮説を対立仮説（H1）として設定します。

しかし、Jarque-Bera検定のような正規性の検定では、この一般的なアプローチが少し異なり、帰無仮説が棄却されない（残差が正規分布に従う）ことが「望ましい」結果とされることが多いです。

また、Ljung-Box検定においても帰無仮説が棄却されないことが好ましい結果とされます。

Ljung-Box検定は時系列モデルの残差に自己相関がないかどうかを検定しますので、自己相関がないということは、モデルが時系列データの情報を十分に取り込んでいることを意味し、残差が単なるホワイトノイズであることを示唆します。

これは、モデルの残差が統計的手法の標準的な前提に合致していることを意味します。このように、検定の種類によっては帰無仮説の内容とその解釈が異なるため、検定の目的とコンテキストを理解することが重要です。

予測

このモデルを用いて予測を行っていきます。

# 学習したARIMA(1,1,0)で未来12か月を予測
forecast = model_fit.forecast(steps=12)

# DataFrameにして表示
forecast_dates = pd.date_range(start=data.index[-1], periods=13, freq='MS')[1:]
forecast_df = pd.DataFrame({'Forecast': forecast}, index=forecast_dates)
forecast_df

このような結果になりました。

プロットしてみようと思います。

#原系列と予測のプロット
plt.figure(figsize=(14, 7))
plt.plot(data['売上'], label='Revenue')
plt.plot(forecast_df['Forecast'], label='Forecasted Revenue', color='orange')

plt.title('Revenue and Revenue Forecast')
plt.xlabel('Date')
plt.ylabel('Revenue')
plt.legend()
plt.grid(True)
plt.show()

95％信頼区間も表示してみます。

# 予測と信頼区間を得る
forecast_results = model_fit.get_forecast(steps=12)
forecast_mean = forecast_results.predicted_mean
forecast_conf_int = forecast_results.conf_int(alpha=0.05)

# DataFrameにまとめる
forecast_dates = pd.date_range(start=data.index[-1], periods=13, freq='MS')[1:]
forecast_df = pd.DataFrame({'Forecast': forecast_mean}, index=forecast_dates)
forecast_df['Lower CI'] = forecast_conf_int['lower 売上']
forecast_df['Upper CI'] = forecast_conf_int['upper 売上']

# プロット
plt.figure(figsize=(14, 7))
plt.plot(data['売上'], label='Revenue')
plt.plot(forecast_df['Forecast'], label='Forecasted Revenue', color='orange')
plt.fill_between(forecast_df.index,
                 forecast_df['Lower CI'],
                 forecast_df['Upper CI'],
                 color='orange', alpha=0.3)

plt.title('Revenue Forecast with 95% Confidence Intervals')
plt.xlabel('Date')
plt.ylabel('Revenue')
plt.legend()
plt.grid(True)
plt.show()

図示してみると、平坦な予測で、95％区間も広く、あまりいい結果とは言えないですね…

データ数が少ないこと、季節性を加味していないことなどが原因と考えられます。

とりあえず今回は、ボックス・ジェンキンス法をやってみるという試みでしたので、この辺にしておこうと思います。

実務で用いるにはデータについてもっと詳細に調べる必要がありそうです。

では♪

参考書籍

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装 

馬場 真哉 (著)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) 

沖本 竜義 (著)