時系列データの特徴把握（ACF・PACF・季節分解・単位根検定）

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
plt.rcParams['font.family'] = 'MS Gothic'

# データの読み込み
file_path = '../00_data/01_kaggle_Air_Passengers/AirPassengers.csv'
df = pd.read_csv(file_path)

# カラム名変更
df = df.rename(columns={'#Passengers':'乗客数'})

# 年月データを型変換
df['Month'] = pd.DatetimeIndex(df['Month'])

# 差分取得（前月比±何%といった形式で算出）
df['前月比%'] = df['乗客数']/df['乗客数'].shift(1) -1

# 年月をインデックスへ
index_cols = ['Month']
df = df.set_index(index_cols)

# 元データから可視化用に必要箇所取得
# 乗客数
use_cols = ['乗客数']
df_rawdata = df[use_cols]
display(df_rawdata.head(2))

# 前月比%
use_cols = ['前月比%']
df_zengetsuhi = df[use_cols]
# 前月比を算出できない月は削除
df_zengetsuhi = df_zengetsuhi[df_zengetsuhi['前月比%'].notnull()]
display(df_zengetsuhi.head(2))

# 可視化
fig, axes = plt.subplots(3, 2, figsize=(16, 16))

# 元データ分析の可視化
axes[0][0].plot(df_rawdata)
axes[0][0].set_title("米国の月間航空機搭乗者数推移（1949-1960年）")
axes[0][0].set_xlabel("年月")
axes[0][0].set_ylabel("乗客数（万人）")

fig = sm.graphics.tsa.plot_acf(df_rawdata, lags=40, ax=axes[1][0])
axes[1][0].set_title("元データの自己相関係数（ACF）")
axes[1][0].set_xlabel("ラグ（月数）")
axes[1][0].set_ylabel("ACF")

fig = sm.graphics.tsa.plot_pacf(df_rawdata, lags=40, ax=axes[2][0])
axes[2][0].set_title("元データの偏自己相関係数（PACF）")
axes[2][0].set_xlabel("ラグ（月数）")
axes[2][0].set_ylabel("PACF")


# 差分化モデルの可視化
axes[0][1].plot(df_zengetsuhi.dropna(), alpha=1)
axes[0][1].set_title("差分化モデル（前月比差分%）")
axes[0][1].set_xlabel("年月")
axes[0][1].set_ylabel("乗客数（万人）")

fig = sm.graphics.tsa.plot_acf(df_zengetsuhi, lags=40, ax=axes[1][1])
axes[1][1].set_title("前月比差分の自己相関係数（ACF）")
axes[1][1].set_xlabel("ラグ（月数）")
axes[1][1].set_ylabel("ACF")

fig = sm.graphics.tsa.plot_pacf(df_zengetsuhi, lags=40, ax=axes[2][1])
axes[2][1].set_title("前月比差分の偏自己相関係数（PACF）")
axes[2][1].set_xlabel("ラグ（月数）")
axes[2][1].set_ylabel("PACF")

# 画像出力
out_path = '../02_work/03_時系列予測_ARIMA/01_可視化/01_可視化_搭乗者数推移_ACF_PAXF_前月比差分.png'
fig.savefig(out_path, dpi=150, bbox_inches='tight')

# 成分分解
# コレログラムでの結果をふまえて季節周期は12か月として成分分解
period = 12
dcmp = sm.tsa.seasonal_decompose(df_rawdata, model='additive', period=period)

# 可視化
# 元データ、トレンド成分、季節成分、残差（不規則変動）の順に縦に並べて可視化
titles = ['元データ', 'トレンド成分', '季節成分', '残差成分']
datas = [dcmp.observed, dcmp.trend, dcmp.seasonal, dcmp.resid]

fig, axes = plt.subplots(4, 1, figsize=(12, 8), sharex=True)

for i, data in enumerate(datas):
    axes[i].set_title(titles[i])
    axes[i].plot(data)

# 画像出力
out_path = '../02_work/03_時系列予測_ARIMA/01_可視化/01_可視化_搭乗者数推移_季節分解.png'
fig.savefig(out_path, dpi=150, bbox_inches='tight')

plt.show()

# 単位根検定（ADF検定）実施
# 元データ
adf_rawdata = sm.tsa.stattools.adfuller(df_rawdata)
print('ADF検定結果（元データ）：', adf_rawdata)
# 差分化後データ
adf_zengetsuhi = sm.tsa.stattools.adfuller(df_zengetsuhi)
print('ADF検定結果（差分化データ）：', adf_zengetsuhi)

# adfullerの出力は順に、テスト統計値、p値、ラグ数、使用したデータ数、検定統計量の臨界値(1%，5%，10%水準)、AIC
# ADF検定結果（元データ）： (0.8153688792060597, 0.9918802434376411, 13, 130, {'1%': -3.4816817173418295, '5%': -2.8840418343195267, '10%': -2.578770059171598}, 996.692930839019)
# ADF検定結果（差分化データ）： (-2.926108907545377, 0.04239503972524038, 14, 128, {'1%': -3.4825006939887997, '5%': -2.884397984161377, '10%': -2.578960197753906}, -425.9255544019345)