演習：Store Salesデータを徹底分析しよう

import pandas as pd
import matplotlib.pyplot as plt

train = pd.read_csv('train.csv', parse_dates=['date'])
stores = pd.read_csv('stores.csv')
oil = pd.read_csv('oil.csv', parse_dates=['date'])
holidays = pd.read_csv('holidays_events.csv', parse_dates=['date'])
transactions = pd.read_csv('transactions.csv', parse_dates=['date'])

# 基本情報
for name, df in [('train', train), ('stores', stores), ('oil', oil),
                  ('holidays', holidays), ('transactions', transactions)]:
    print(f"\n{name}: {df.shape[0]:,}行 x {df.shape[1]}列")
    print(f"  欠損: {df.isnull().sum().to_dict()}")

# 日次売上合計の可視化
daily_total = train.groupby('date')['sales'].sum()
fig, ax = plt.subplots(figsize=(16, 5))
daily_total.plot(ax=ax)
ax.set_title('全店舗合計日次売上推移')
ax.set_ylabel('Sales')
plt.show()

# Top 5カテゴリ
top5 = train.groupby('family')['sales'].sum().sort_values(ascending=False).head(5)
print("\n売上Top 5カテゴリ:")
print(top5)
# GROCERY I, BEVERAGES, PRODUCE, CLEANING, DAIRY

# 店舗タイプ別
merged = train.merge(stores, on='store_nbr')
type_sales = merged.groupby('type')['sales'].mean().sort_values(ascending=False)
print("\n店舗タイプ別平均売上:")
print(type_sales)

from statsmodels.tsa.seasonal import STL
import numpy as np

top3_families = ['GROCERY I', 'BEVERAGES', 'PRODUCE']

fig, axes = plt.subplots(3, 4, figsize=(20, 12))
for i, family in enumerate(top3_families):
    series = train[train['family'] == family].groupby('date')['sales'].sum()
    series = series.asfreq('D').fillna(method='ffill')

    stl = STL(series, period=7, robust=True)
    result = stl.fit()

    axes[i, 0].plot(series.index, series.values, linewidth=0.5)
    axes[i, 0].set_title(f'{family}: 原系列')
    axes[i, 1].plot(result.trend.index, result.trend.values)
    axes[i, 1].set_title(f'{family}: トレンド')
    axes[i, 2].plot(result.seasonal.index, result.seasonal.values, linewidth=0.5)
    axes[i, 2].set_title(f'{family}: 季節性')
    axes[i, 3].plot(result.resid.index, result.resid.values, linewidth=0.5)
    axes[i, 3].set_title(f'{family}: 残差')

plt.tight_layout()
plt.show()

# 曜日別売上
weekday = train.groupby(train['date'].dt.dayofweek)['sales'].mean()
weekday.index = ['月', '火', '水', '木', '金', '土', '日']
weekday.plot(kind='bar', title='曜日別平均売上')
plt.show()

# 月別売上
monthly = train.groupby(train['date'].dt.month)['sales'].mean()
monthly.plot(kind='bar', title='月別平均売上')
plt.show()

# 外れ値の特定（残差が3σ以上）
for family in top3_families:
    series = train[train['family'] == family].groupby('date')['sales'].sum()
    series = series.asfreq('D').fillna(method='ffill')
    stl = STL(series, period=7, robust=True)
    result = stl.fit()
    threshold = result.resid.std() * 3
    outliers = result.resid[abs(result.resid) > threshold]
    print(f"\n{family}: 外れ値 {len(outliers)}件")
    print(outliers.sort_values(ascending=False).head(5))

# 石油価格の補間
oil['dcoilwtico'] = oil['dcoilwtico'].interpolate(method='linear')
fig, ax1 = plt.subplots(figsize=(16, 5))
ax1.plot(oil['date'], oil['dcoilwtico'], color='orange', label='石油価格')
ax2 = ax1.twinx()
ax2.plot(daily_total.index, daily_total.values, color='blue', alpha=0.3, label='売上')
ax1.set_ylabel('Oil Price (USD)')
ax2.set_ylabel('Sales')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.title('石油価格と売上の推移')
plt.show()

# ラグ付き相関
merged = daily_total.reset_index()
merged.columns = ['date', 'sales']
merged = merged.merge(oil, on='date', how='left')
merged['dcoilwtico'] = merged['dcoilwtico'].interpolate()

for lag in [0, 7, 14, 30]:
    corr = merged['sales'].corr(merged['dcoilwtico'].shift(lag))
    print(f"ラグ{lag:2d}日: 相関 = {corr:.3f}")

# 祝日影響分析
holiday_types = holidays['type'].unique()
for htype in holiday_types:
    hol_dates = holidays[holidays['type'] == htype]['date'].unique()
    on_holiday = daily_total[daily_total.index.isin(hol_dates)].mean()
    off_holiday = daily_total[~daily_total.index.isin(hol_dates)].mean()
    effect = (on_holiday / off_holiday - 1) * 100
    print(f"{htype:12s}: 売上影響 = {effect:+.1f}%")

# 地震前後の分析
earthquake = pd.Timestamp('2016-04-16')
window = 60
eq_period = daily_total[
    (daily_total.index >= earthquake - pd.Timedelta(days=window)) &
    (daily_total.index <= earthquake + pd.Timedelta(days=window))
]
fig, ax = plt.subplots(figsize=(14, 5))
eq_period.plot(ax=ax)
ax.axvline(earthquake, color='red', linestyle='--', label='地震発生')
ax.legend()
ax.set_title('地震前後の売上推移')
plt.show()

# プロモーション効果
for family in top3_families:
    fd = train[train['family'] == family]
    promo = fd[fd['onpromotion'] > 0]['sales'].mean()
    no_promo = fd[fd['onpromotion'] == 0]['sales'].mean()
    lift = (promo / no_promo - 1) * 100 if no_promo > 0 else 0
    print(f"{family:15s}: プロモリフト = {lift:+.1f}%")

主要パターン:
1. 強い週次季節性（日曜が最大、月〜木は低め）
2. 年次季節性（12月ピーク、1月低下）
3. 2015年以降の成長トレンド
4. 祝日前日の買い込み効果
5. 石油価格下落期の売上低迷

組み込むべき特徴量:
- ラグ特徴量: lag_1, lag_7, lag_14, lag_28
- ローリング統計: rolling_mean_7, rolling_mean_28, rolling_std_7
- カレンダー: dayofweek, month, day, is_weekend
- 祝日: is_holiday, days_to_holiday, holiday_type
- 石油: oil_price, oil_ma7, oil_change_7d
- プロモーション: onpromotion, promo_ratio
- 店舗: store_type, cluster

予測困難なケース:
- 地震のような突発イベント
- 新商品や終売品（過去データなし）
- プロモーション情報が欠損する2014年以前

方針:
1. ベースライン: 移動平均（7日/28日）
2. SARIMA: 週次季節性を組み込み
3. Prophet: 祝日効果を自動モデル化
4. LightGBM: 全特徴量を投入
5. アンサンブル: 上記モデルの加重平均