演習：予測モデルを構築・比較しよう

import pandas as pd
import numpy as np
from pmdarima import auto_arima
from prophet import Prophet
from sklearn.metrics import mean_absolute_error, mean_squared_error
import matplotlib.pyplot as plt

# --- データ準備 ---
train = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/train.csv',
                     parse_dates=['date'])
grocery = train[(train['store_nbr'] == 1) & (train['family'] == 'GROCERY I')]
grocery = grocery.set_index('date')['sales']

train_data = grocery[:'2017-07-31']
test_data = grocery['2017-08-01':'2017-08-15']

# --- SARIMAモデル ---
sarima_model = auto_arima(
    train_data,
    seasonal=True, m=7,
    start_p=0, max_p=3,
    start_q=0, max_q=3,
    start_P=0, max_P=2,
    start_Q=0, max_Q=2,
    information_criterion='aic',
    trace=True,
    error_action='ignore',
    suppress_warnings=True,
    stepwise=True
)
sarima_forecast = sarima_model.predict(n_periods=15)
sarima_forecast = pd.Series(sarima_forecast, index=test_data.index)

# --- Prophetモデル ---
holidays_ec = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/holidays_events.csv',
                           parse_dates=['date'])
national_holidays = holidays_ec[holidays_ec['locale'] == 'National'][['date', 'description']]
national_holidays.columns = ['ds', 'holiday']

prophet_df = train_data.reset_index()
prophet_df.columns = ['ds', 'y']

prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=True,
    holidays=national_holidays,
    changepoint_prior_scale=0.05
)
prophet_model.fit(prophet_df)

future = prophet_model.make_future_dataframe(periods=15)
prophet_pred = prophet_model.predict(future)
prophet_forecast = prophet_pred.set_index('ds')['yhat']['2017-08-01':'2017-08-15']
prophet_forecast = prophet_forecast.clip(lower=0)

# --- 評価 ---
def evaluate(actual, predicted, model_name):
    mae = mean_absolute_error(actual, predicted)
    rmse = np.sqrt(mean_squared_error(actual, predicted))
    mape = np.mean(np.abs((actual - predicted) / actual)) * 100
    print(f"{model_name}: MAE={mae:.2f}, RMSE={rmse:.2f}, MAPE={mape:.2f}%")
    return {'model': model_name, 'MAE': mae, 'RMSE': rmse, 'MAPE': mape}

results = []
results.append(evaluate(test_data, sarima_forecast, 'SARIMA'))
results.append(evaluate(test_data, prophet_forecast, 'Prophet'))

# --- 可視化 ---
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(test_data.index, test_data.values, 'ko-', label='実績')
ax.plot(test_data.index, sarima_forecast.values, 'b--', label='SARIMA')
ax.plot(test_data.index, prophet_forecast.values, 'r--', label='Prophet')
ax.legend()
ax.set_title('SARIMA vs Prophet: 15日間予測比較')
ax.set_ylabel('売上')
plt.tight_layout()
plt.show()

import lightgbm as lgb
from sklearn.metrics import mean_absolute_error
import numpy as np

# --- データ準備（全店舗・GROCERYカテゴリ） ---
grocery_all = train[train['family'] == 'GROCERY I'].copy()
oil = pd.read_csv('/kaggle/input/store-sales-time-series-forecasting/oil.csv',
                   parse_dates=['date'])
oil = oil.set_index('date').reindex(grocery_all['date'].unique()).interpolate()

grocery_all = grocery_all.merge(oil.reset_index(), on='date', how='left')
grocery_all['dcoilwtico'] = grocery_all['dcoilwtico'].interpolate()

# --- 基本特徴量 ---
grocery_all['dayofweek'] = grocery_all['date'].dt.dayofweek
grocery_all['month'] = grocery_all['date'].dt.month
grocery_all['year'] = grocery_all['date'].dt.year
grocery_all['day'] = grocery_all['date'].dt.day

base_features = ['store_nbr', 'dayofweek', 'month', 'year', 'day']

# --- 拡張特徴量 ---
grocery_all = grocery_all.sort_values(['store_nbr', 'date'])
for lag in [7, 14, 28]:
    grocery_all[f'lag_{lag}'] = grocery_all.groupby('store_nbr')['sales'].shift(lag)

grocery_all['rolling_mean_7'] = grocery_all.groupby('store_nbr')['sales'].transform(
    lambda x: x.shift(1).rolling(7).mean()
)
grocery_all['rolling_std_7'] = grocery_all.groupby('store_nbr')['sales'].transform(
    lambda x: x.shift(1).rolling(7).std()
)
grocery_all['rolling_mean_28'] = grocery_all.groupby('store_nbr')['sales'].transform(
    lambda x: x.shift(1).rolling(28).mean()
)

# 給与日フラグ（15日と末日）
grocery_all['is_payday'] = grocery_all['day'].isin([15, 28, 29, 30, 31]).astype(int)

# 祝日フラグ
holiday_dates = set(holidays_ec[holidays_ec['locale'] == 'National']['date'])
grocery_all['is_holiday'] = grocery_all['date'].isin(holiday_dates).astype(int)

extended_features = base_features + [
    'lag_7', 'lag_14', 'lag_28',
    'rolling_mean_7', 'rolling_std_7', 'rolling_mean_28',
    'dcoilwtico', 'is_payday', 'is_holiday'
]

# --- 訓練/テスト分割 ---
train_mask = grocery_all['date'] <= '2017-07-31'
test_mask = (grocery_all['date'] >= '2017-08-01') & (grocery_all['date'] <= '2017-08-15')

# --- ベースラインモデル ---
X_train_base = grocery_all[train_mask][base_features]
X_test_base = grocery_all[test_mask][base_features]
y_train = grocery_all[train_mask]['sales']
y_test = grocery_all[test_mask]['sales']

lgb_base = lgb.LGBMRegressor(n_estimators=500, learning_rate=0.05, num_leaves=31)
lgb_base.fit(X_train_base, y_train)
pred_base = lgb_base.predict(X_test_base).clip(min=0)
mape_base = np.mean(np.abs((y_test - pred_base) / y_test.replace(0, 1))) * 100

# --- 拡張モデル ---
grocery_clean = grocery_all.dropna(subset=extended_features)
train_mask_c = grocery_clean['date'] <= '2017-07-31'
test_mask_c = (grocery_clean['date'] >= '2017-08-01') & (grocery_clean['date'] <= '2017-08-15')

X_train_ext = grocery_clean[train_mask_c][extended_features]
X_test_ext = grocery_clean[test_mask_c][extended_features]
y_train_ext = grocery_clean[train_mask_c]['sales']
y_test_ext = grocery_clean[test_mask_c]['sales']

lgb_ext = lgb.LGBMRegressor(n_estimators=500, learning_rate=0.05, num_leaves=31)
lgb_ext.fit(X_train_ext, y_train_ext)
pred_ext = lgb_ext.predict(X_test_ext).clip(min=0)
mape_ext = np.mean(np.abs((y_test_ext - pred_ext) / y_test_ext.replace(0, 1))) * 100

print(f"ベースライン MAPE: {mape_base:.2f}%")
print(f"拡張モデル MAPE: {mape_ext:.2f}%")
print(f"改善率: {(mape_base - mape_ext) / mape_base * 100:.1f}%")

# --- 特徴量重要度 ---
importance = pd.DataFrame({
    'feature': extended_features,
    'importance': lgb_ext.feature_importances_
}).sort_values('importance', ascending=False)

fig, ax = plt.subplots(figsize=(10, 6))
importance.head(10).plot.barh(x='feature', y='importance', ax=ax)
ax.set_title('特徴量重要度 Top 10')
ax.invert_yaxis()
plt.tight_layout()
plt.show()

from sklearn.linear_model import Ridge
from itertools import product

# --- 各モデルの予測値を揃える（store_nbr=1のテスト期間） ---
# sarima_forecast, prophet_forecast は Mission 1 で取得済み
# LightGBM の store_nbr=1 予測
lgb_store1 = pred_ext[grocery_clean[test_mask_c]['store_nbr'] == 1]
lgb_forecast = pd.Series(lgb_store1, index=test_data.index)

actual = test_data.values
preds = {
    'SARIMA': sarima_forecast.values,
    'Prophet': prophet_forecast.values,
    'LightGBM': lgb_forecast.values,
}

# --- 単純平均 ---
simple_avg = np.mean([preds['SARIMA'], preds['Prophet'], preds['LightGBM']], axis=0)
mape_avg = np.mean(np.abs((actual - simple_avg) / actual)) * 100

# --- 重み付きアンサンブル（グリッドサーチ） ---
best_mape = float('inf')
best_weights = None
for w1, w2 in product(np.arange(0.0, 1.05, 0.1), repeat=2):
    w3 = 1.0 - w1 - w2
    if w3 < 0:
        continue
    weighted = w1 * preds['SARIMA'] + w2 * preds['Prophet'] + w3 * preds['LightGBM']
    mape = np.mean(np.abs((actual - weighted) / actual)) * 100
    if mape < best_mape:
        best_mape = mape
        best_weights = (w1, w2, w3)

print(f"最適重み: SARIMA={best_weights[0]:.1f}, Prophet={best_weights[1]:.1f}, LightGBM={best_weights[2]:.1f}")
print(f"重み付きアンサンブル MAPE: {best_mape:.2f}%")

# --- スタッキング ---
# 訓練期間の交差検証予測でメタ学習器を訓練するのが理想
# ここでは簡易的にテストデータで構築手順を示す
X_stack = np.column_stack([preds['SARIMA'], preds['Prophet'], preds['LightGBM']])
meta_model = Ridge(alpha=1.0)
meta_model.fit(X_stack, actual)
stack_pred = meta_model.predict(X_stack)
mape_stack = np.mean(np.abs((actual - stack_pred) / actual)) * 100

print(f"\nメタ学習器の重み: {dict(zip(['SARIMA', 'Prophet', 'LightGBM'], meta_model.coef_.round(3)))}")

# --- 全モデル比較 ---
all_results = pd.DataFrame([
    {'モデル': 'SARIMA', 'MAPE': np.mean(np.abs((actual - preds['SARIMA']) / actual)) * 100},
    {'モデル': 'Prophet', 'MAPE': np.mean(np.abs((actual - preds['Prophet']) / actual)) * 100},
    {'モデル': 'LightGBM', 'MAPE': np.mean(np.abs((actual - preds['LightGBM']) / actual)) * 100},
    {'モデル': '単純平均', 'MAPE': mape_avg},
    {'モデル': '重み付き', 'MAPE': best_mape},
    {'モデル': 'スタッキング', 'MAPE': mape_stack},
])
print("\n=== 全モデル精度比較 ===")
print(all_results.sort_values('MAPE').to_string(index=False))

# --- 可視化 ---
weighted_pred = (best_weights[0] * preds['SARIMA']
                 + best_weights[1] * preds['Prophet']
                 + best_weights[2] * preds['LightGBM'])

fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(test_data.index, actual, 'ko-', label='実績', linewidth=2)
ax.plot(test_data.index, preds['SARIMA'], 'b--', alpha=0.5, label='SARIMA')
ax.plot(test_data.index, preds['Prophet'], 'r--', alpha=0.5, label='Prophet')
ax.plot(test_data.index, preds['LightGBM'], 'g--', alpha=0.5, label='LightGBM')
ax.plot(test_data.index, weighted_pred, 'm-', linewidth=2, label='重み付きアンサンブル')
ax.legend()
ax.set_title('全モデル予測比較')
ax.set_ylabel('売上')
plt.tight_layout()
plt.show()