Facebook Prophetによる時系列予測 — 加法モデルで手軽に高精度予測

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from prophet import Prophet
import warnings
warnings.filterwarnings('ignore')

# AirPassengersデータの準備
from statsmodels.datasets import get_rdataset
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']

# Prophetの入力フォーマット（ds列とy列が必要）
df = pd.DataFrame({
    'ds': pd.date_range('1949-01-01', periods=len(data), freq='MS'),
    'y': data['passengers'].values.astype(float)
})

# 訓練/テスト分割
train = df[df['ds'] < '1959-01-01']
test = df[df['ds'] >= '1959-01-01']

# Prophetモデルの作成とフィッティング
model = Prophet(
    yearly_seasonality=True,     # 年次季節性
    weekly_seasonality=False,    # 週次季節性（月次データなので不要）
    daily_seasonality=False,     # 日次季節性（月次データなので不要）
    changepoint_prior_scale=0.05,  # トレンドの柔軟性
    seasonality_mode='multiplicative'  # 乗法的季節性
)
model.fit(train)

# 予測（テスト期間 + 24か月先）
future = model.make_future_dataframe(periods=len(test) + 24, freq='MS')
forecast = model.predict(future)

# プロット
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 予測結果
axes[0].plot(train['ds'], train['y'], 'o', markersize=3, label='Train')
axes[0].plot(test['ds'], test['y'], 'o', markersize=3, color='black', label='Test')
axes[0].plot(forecast['ds'], forecast['yhat'], linewidth=1.5,
             color='red', label='Forecast')
axes[0].fill_between(forecast['ds'], forecast['yhat_lower'], forecast['yhat_upper'],
                     color='red', alpha=0.1)
axes[0].set_title('Prophet Forecast — AirPassengers')
axes[0].legend()

# テスト期間の拡大
test_forecast = forecast[forecast['ds'].isin(test['ds'])]
axes[1].plot(test['ds'], test['y'], 'o-', markersize=4, label='Actual')
axes[1].plot(test_forecast['ds'], test_forecast['yhat'], 'o-',
             markersize=4, color='red', label='Predicted')
axes[1].fill_between(test_forecast['ds'],
                     test_forecast['yhat_lower'], test_forecast['yhat_upper'],
                     color='red', alpha=0.1)
axes[1].set_title('Test Period Detail')
axes[1].legend()

plt.tight_layout()
plt.show()

# 精度指標
y_true = test['y'].values
y_pred = test_forecast['yhat'].values
mae = np.mean(np.abs(y_true - y_pred))
mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
rmse = np.sqrt(np.mean((y_true - y_pred) ** 2))
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 prophet import Prophet
from statsmodels.datasets import get_rdataset
import warnings
warnings.filterwarnings('ignore')

# データ準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
df = pd.DataFrame({
    'ds': pd.date_range('1949-01-01', periods=len(data), freq='MS'),
    'y': data['passengers'].values.astype(float)
})

# モデル
model = Prophet(yearly_seasonality=True, weekly_seasonality=False,
                daily_seasonality=False, seasonality_mode='multiplicative')
model.fit(df)

future = model.make_future_dataframe(periods=24, freq='MS')
forecast = model.predict(future)

# 成分プロット
fig = model.plot_components(forecast, figsize=(12, 8))
plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from prophet import Prophet
from statsmodels.datasets import get_rdataset
import warnings
warnings.filterwarnings('ignore')

# データ準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
df = pd.DataFrame({
    'ds': pd.date_range('1949-01-01', periods=len(data), freq='MS'),
    'y': data['passengers'].values.astype(float)
})

scales = [0.001, 0.05, 0.5, 5.0]

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

for idx, scale in enumerate(scales):
    model = Prophet(changepoint_prior_scale=scale,
                    yearly_seasonality=True, weekly_seasonality=False,
                    daily_seasonality=False, seasonality_mode='multiplicative')
    model.fit(df)

    future = model.make_future_dataframe(periods=24, freq='MS')
    forecast = model.predict(future)

    axes[idx].plot(df['ds'], df['y'], 'o', markersize=2, alpha=0.5)
    axes[idx].plot(forecast['ds'], forecast['yhat'], linewidth=1.5, color='red')
    axes[idx].fill_between(forecast['ds'], forecast['yhat_lower'],
                           forecast['yhat_upper'], color='red', alpha=0.1)
    axes[idx].set_title(f'changepoint_prior_scale = {scale}')
    axes[idx].set_ylabel('Passengers')

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from prophet import Prophet
from prophet.diagnostics import cross_validation, performance_metrics
from statsmodels.datasets import get_rdataset
import warnings
warnings.filterwarnings('ignore')

# データ準備
data = get_rdataset("AirPassengers").data
data.columns = ['time', 'passengers']
df = pd.DataFrame({
    'ds': pd.date_range('1949-01-01', periods=len(data), freq='MS'),
    'y': data['passengers'].values.astype(float)
})

# モデルのフィッティング
model = Prophet(yearly_seasonality=True, weekly_seasonality=False,
                daily_seasonality=False, seasonality_mode='multiplicative')
model.fit(df)

# クロスバリデーション
# initial: 最初の訓練期間の長さ
# period: cutoffをずらす間隔
# horizon: 予測ホライズン
df_cv = cross_validation(model, initial='730 days', period='180 days',
                         horizon='365 days')
print("Cross-validation results:")
print(df_cv.head(10))

# 性能指標
df_metrics = performance_metrics(df_cv, rolling_window=1)
print("\nPerformance Metrics:")
print(df_metrics[['horizon', 'mape', 'rmse', 'mae']].to_string())

# 予測ホライズンごとのMAPE
fig, ax = plt.subplots(figsize=(10, 5))
df_cv['horizon_days'] = (df_cv['ds'] - df_cv['cutoff']).dt.days
df_cv['ape'] = np.abs((df_cv['y'] - df_cv['yhat']) / df_cv['y'])

# 月ごとのMAPEを集計
df_cv['horizon_month'] = df_cv['horizon_days'] // 30 + 1
monthly_mape = df_cv.groupby('horizon_month')['ape'].mean() * 100

ax.bar(monthly_mape.index, monthly_mape.values, color='steelblue', alpha=0.7)
ax.set_xlabel('Forecast Horizon (months)')
ax.set_ylabel('MAPE (%)')
ax.set_title('MAPE by Forecast Horizon (Cross-Validation)')
plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from prophet import Prophet
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.astype(float)
time_index = pd.date_range('1949-01-01', periods=len(passengers), freq='MS')

# 分割
split_date = '1958-01-01'
n_test = sum(time_index >= split_date)

# --- Prophet ---
df_prophet = pd.DataFrame({'ds': time_index, 'y': passengers})
train_prophet = df_prophet[df_prophet['ds'] < split_date]
test_prophet = df_prophet[df_prophet['ds'] >= split_date]

model_prophet = Prophet(yearly_seasonality=True, weekly_seasonality=False,
                        daily_seasonality=False,
                        seasonality_mode='multiplicative',
                        changepoint_prior_scale=0.05)
model_prophet.fit(train_prophet)
future = model_prophet.make_future_dataframe(periods=n_test, freq='MS')
forecast_prophet = model_prophet.predict(future)
pred_prophet = forecast_prophet[forecast_prophet['ds'] >= split_date]['yhat'].values

# --- SARIMA ---
ts = pd.Series(np.log(passengers), index=time_index)
train_sarima = ts[ts.index < split_date]

sarima_model = SARIMAX(train_sarima, order=(0, 1, 1),
                       seasonal_order=(0, 1, 1, 12)).fit(disp=False)
pred_sarima_log = sarima_model.forecast(n_test)
pred_sarima = np.exp(pred_sarima_log.values)

# --- 比較 ---
y_true = passengers[time_index >= split_date]

fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(time_index[time_index < split_date], passengers[time_index < split_date],
        linewidth=1, color='gray', label='Train')
ax.plot(time_index[time_index >= split_date], y_true,
        linewidth=2, color='black', label='Actual')
ax.plot(time_index[time_index >= split_date], pred_prophet,
        linewidth=2, linestyle='--', color='blue', label='Prophet')
ax.plot(time_index[time_index >= split_date], pred_sarima,
        linewidth=2, linestyle='--', color='red', label='SARIMA')
ax.axvline(x=pd.Timestamp(split_date), color='gray', linestyle=':', alpha=0.5)
ax.set_title('Prophet vs SARIMA — AirPassengers')
ax.legend()
plt.tight_layout()
plt.show()

# 精度比較
for name, preds in [('Prophet', pred_prophet), ('SARIMA', pred_sarima)]:
    mae = np.mean(np.abs(y_true - preds))
    mape = np.mean(np.abs((y_true - preds) / y_true)) * 100
    rmse = np.sqrt(np.mean((y_true - preds) ** 2))
    print(f"{name:10s}: MAE={mae:.1f}, MAPE={mape:.1f}%, RMSE={rmse:.1f}")