SARIMAモデルによる季節時系列予測 — 季節性を捉える拡張ARIMAの理論と実装

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset

# AirPassengersデータの読み込み
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
n = len(passengers)
time_index = pd.date_range('1949-01', periods=n, freq='MS')

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

# 元データ
axes[0].plot(time_index, passengers, linewidth=1.5)
axes[0].set_title('Original Data (AirPassengers)')
axes[0].set_ylabel('Passengers')

# 対数変換（分散安定化）
log_passengers = np.log(passengers)
axes[1].plot(time_index, log_passengers, linewidth=1.5, color='darkorange')
axes[1].set_title('Log Transformed')
axes[1].set_ylabel('log(Passengers)')

# 季節差分
seasonal_diff = log_passengers[12:] - log_passengers[:-12]
axes[2].plot(time_index[12:], seasonal_diff, linewidth=1.5, color='green')
axes[2].set_title('Seasonal Difference (d=0, D=1, s=12)')
axes[2].set_ylabel('Δ₁₂ log(Passengers)')
axes[2].axhline(y=0, color='gray', linestyle='--', alpha=0.5)

# 季節差分 + 1次差分
double_diff = seasonal_diff[1:] - seasonal_diff[:-1]
axes[3].plot(time_index[13:], double_diff, linewidth=1.5, color='purple')
axes[3].set_title('Seasonal + First Difference (d=1, D=1, s=12)')
axes[3].set_ylabel('Δ₁Δ₁₂ log(Passengers)')
axes[3].axhline(y=0, color='gray', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# データの準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')

# 対数変換 + 季節差分 + 1次差分
log_pass = np.log(passengers)
seasonal_diff = log_pass[12:] - log_pass[:-12]
double_diff = seasonal_diff[1:] - seasonal_diff[:-1]

# ADF検定
adf_result = adfuller(double_diff)
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.6f}")
print(f"=> {'Stationary' if adf_result[1] < 0.05 else 'Non-stationary'}")

# ACF/PACF
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

plot_acf(double_diff, lags=36, ax=axes[0])
axes[0].set_title('ACF of Δ₁Δ₁₂ log(Passengers)')

plot_pacf(double_diff, lags=36, ax=axes[1])
axes[1].set_title('PACF of Δ₁Δ₁₂ log(Passengers)')

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import warnings
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.statespace.sarimax import SARIMAX
import itertools

warnings.filterwarnings('ignore')

# データ準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
log_pass = np.log(passengers)
ts = pd.Series(log_pass, index=pd.date_range('1949-01', periods=len(passengers), freq='MS'))

# パラメータの候補
p_range = range(0, 3)
q_range = range(0, 3)
P_range = range(0, 2)
Q_range = range(0, 2)

# d=1, D=1, s=12 は固定
d, D, s = 1, 1, 12

results = []
for p, q, P, Q in itertools.product(p_range, q_range, P_range, Q_range):
    try:
        model = SARIMAX(ts, order=(p, d, q),
                        seasonal_order=(P, D, Q, s),
                        enforce_stationarity=False,
                        enforce_invertibility=False)
        fit = model.fit(disp=False, maxiter=200)
        results.append({
            'order': f'({p},{d},{q})({P},{D},{Q})_{s}',
            'p': p, 'q': q, 'P': P, 'Q': Q,
            'AIC': fit.aic,
            'BIC': fit.bic
        })
    except Exception:
        continue

results_df = pd.DataFrame(results).sort_values('AIC').reset_index(drop=True)
print("Top 10 Models by AIC:")
print(results_df.head(10).to_string(index=False))

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
warnings.filterwarnings('ignore')

# データ準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
log_pass = np.log(passengers)
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')
ts = pd.Series(log_pass, index=time_index)

# 訓練/テスト分割（最後の24か月をテスト）
train = ts[:'1958-12']
test = ts['1959-01':]

# SARIMA(0,1,1)(0,1,1)_12 のフィッティング
model = SARIMAX(train, order=(0, 1, 1),
                seasonal_order=(0, 1, 1, 12),
                enforce_stationarity=False,
                enforce_invertibility=False)
result = model.fit(disp=False)

print(result.summary())

# 予測（テスト期間 + 12か月先）
n_forecast = len(test) + 12
forecast = result.get_forecast(steps=n_forecast)
forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int(alpha=0.05)

# 元のスケール（指数変換）に戻す
train_orig = np.exp(train)
test_orig = np.exp(test)
forecast_orig = np.exp(forecast_mean)
ci_lower = np.exp(forecast_ci.iloc[:, 0])
ci_upper = np.exp(forecast_ci.iloc[:, 1])

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 対数スケールでのプロット
axes[0].plot(train.index, train, label='Train', linewidth=1.5)
axes[0].plot(test.index, test, label='Test', linewidth=1.5, color='darkorange')
axes[0].plot(forecast_mean.index, forecast_mean, label='Forecast',
             linewidth=1.5, color='red', linestyle='--')
axes[0].fill_between(forecast_ci.index,
                     forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1],
                     color='red', alpha=0.1)
axes[0].set_title('SARIMA(0,1,1)(0,1,1)₁₂ — Log Scale')
axes[0].set_ylabel('log(Passengers)')
axes[0].legend()

# 元スケールでのプロット
axes[1].plot(train_orig.index, train_orig, label='Train', linewidth=1.5)
axes[1].plot(test_orig.index, test_orig, label='Test', linewidth=1.5, color='darkorange')
axes[1].plot(forecast_orig.index, forecast_orig, label='Forecast',
             linewidth=1.5, color='red', linestyle='--')
axes[1].fill_between(ci_lower.index, ci_lower, ci_upper,
                     color='red', alpha=0.1)
axes[1].set_title('SARIMA(0,1,1)(0,1,1)₁₂ — Original Scale')
axes[1].set_ylabel('Passengers')
axes[1].legend()

plt.tight_layout()
plt.show()

# 予測精度（テスト期間のみ）
test_forecast = forecast_orig[:len(test)]
mae = np.mean(np.abs(test_orig.values - test_forecast.values))
mape = np.mean(np.abs((test_orig.values - test_forecast.values) / test_orig.values)) * 100
rmse = np.sqrt(np.mean((test_orig.values - test_forecast.values) ** 2))

print(f"\n予測精度（テスト期間: 24か月）:")
print(f"MAE  = {mae:.2f}")
print(f"MAPE = {mape:.2f}%")
print(f"RMSE = {rmse:.2f}")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# データ準備（訓練データのみ）
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
passengers = data['passengers'].values
log_pass = np.log(passengers)
time_index = pd.date_range('1949-01', periods=len(passengers), freq='MS')
ts = pd.Series(log_pass, index=time_index)
train = ts[:'1958-12']

# フィッティング
model = SARIMAX(train, order=(0, 1, 1),
                seasonal_order=(0, 1, 1, 12),
                enforce_stationarity=False,
                enforce_invertibility=False)
result = model.fit(disp=False)

# 残差診断
residuals = result.resid

fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 残差の時系列
axes[0, 0].plot(residuals, linewidth=0.8)
axes[0, 0].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[0, 0].set_title('Residuals')

# 残差のACF
plot_acf(residuals, lags=36, ax=axes[0, 1])
axes[0, 1].set_title('ACF of Residuals')

# Q-Qプロット
stats.probplot(residuals, dist="norm", plot=axes[1, 0])
axes[1, 0].set_title('Q-Q Plot')

# 残差のヒストグラム
axes[1, 1].hist(residuals, bins=25, density=True, alpha=0.7, color='steelblue')
x_range = np.linspace(residuals.min(), residuals.max(), 100)
axes[1, 1].plot(x_range,
                stats.norm.pdf(x_range, residuals.mean(), residuals.std()),
                'r-', linewidth=2)
axes[1, 1].set_title('Distribution of Residuals')

plt.tight_layout()
plt.show()

# Ljung-Box検定
lb = acorr_ljungbox(residuals, lags=[12, 24, 36], return_df=True)
print("Ljung-Box Test:")
print(lb)