因果推論学習：Phase3

# ローカル環境の場合
# pip install statsmodels

# 必要なライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')

# 日本語フォント設定（matplotlib）
plt.rcParams['font.sans-serif'] = ['Hiragino Sans', 'Yu Gothic', 'Meiryo', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("✅ statsmodels環境構築完了！")

# サンプルデータ: 週次問い合わせ数（52週 × 2年 = 104週）
np.random.seed(42)
weeks = pd.date_range(start='2023-01-01', periods=104, freq='W')

# トレンド成分（徐々に増加）
trend = np.linspace(100, 150, 104)

# 季節性成分（春に増加、冬に減少）
seasonal = 20 * np.sin(2 * np.pi * np.arange(104) / 52)

# ランダム成分
noise = np.random.normal(0, 5, 104)

# 問い合わせ数 = トレンド + 季節性 + ノイズ
inquiries = trend + seasonal + noise

# DataFrameに変換
df = pd.DataFrame({
'date': weeks,
'inquiries': inquiries
})
df.set_index('date', inplace=True)

# 可視化
fig, axes = plt.subplots(4, 1, figsize=(12, 10))

# 1. 元データ
axes[0].plot(df.index, df['inquiries'], color='#9c27b0', linewidth=2)
axes[0].set_title('週次問い合わせ数（元データ）', fontsize=14, fontweight='bold')
axes[0].set_ylabel('問い合わせ数')
axes[0].grid(alpha=0.3)

# 2. トレンド成分
axes[1].plot(weeks, trend, color='#2196f3', linewidth=2)
axes[1].set_title('トレンド成分', fontsize=12)
axes[1].set_ylabel('トレンド')
axes[1].grid(alpha=0.3)

# 3. 季節性成分
axes[2].plot(weeks, seasonal, color='#ff9800', linewidth=2)
axes[2].set_title('季節性成分', fontsize=12)
axes[2].set_ylabel('季節性')
axes[2].grid(alpha=0.3)

# 4. ランダム成分
axes[3].plot(weeks, noise, color='#4caf50', linewidth=2)
axes[3].set_title('ランダム成分（ノイズ）', fontsize=12)
axes[3].set_ylabel('ノイズ')
axes[3].set_xlabel('日付')
axes[3].grid(alpha=0.3)

plt.tight_layout()
plt.show()

print(f"データ期間: {df.index[0]} 〜 {df.index[-1]}")
print(f"サンプル数: {len(df)}")
print(f"平均問い合わせ数: {df['inquiries'].mean():.1f}")
print(f"標準偏差: {df['inquiries'].std():.1f}")

def adf_test(series, name=''):
"""
ADF検定を実行し、結果を分かりやすく出力する関数

帰無仮説（H0）: 単位根が存在する（= 非定常）
対立仮説（H1）: 単位根が存在しない（= 定常）

p値 < 0.05 なら H0を棄却 → 定常 p値 >= 0.05 なら H0を棄却できない → 非定常
"""
result = adfuller(series, autolag='AIC')

print(f'=== ADF検定結果: {name} ===')
print(f'ADF統計量: {result[0]:.4f}')
print(f'p値: {result[1]:.4f}')
print(f'ラグ次数: {result[2]}')
print(f'サンプルサイズ: {result[3]}')
print('臨界値:')
for key, value in result[4].items():
print(f' {key}: {value:.3f}')

# 判定
if result[1] < 0.05: print(f'✅ 結論: 定常である（p={result[1]:.4f} < 0.05）') else: print(f'⚠️ 結論: 非定常である（p={result[1]:.4f} >= 0.05）')
print()

return result

# 元データのADF検定
adf_test(df['inquiries'], '元データ')

# 1階差分を取る（非定常 → 定常化）
df['inquiries_diff'] = df['inquiries'].diff().dropna()

# 差分データのADF検定
adf_test(df['inquiries_diff'].dropna(), '1階差分データ')

# 可視化
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 元データ
axes[0].plot(df.index, df['inquiries'], color='#9c27b0', linewidth=2)
axes[0].set_title('元データ（非定常の可能性）', fontsize=14, fontweight='bold')
axes[0].set_ylabel('問い合わせ数')
axes[0].grid(alpha=0.3)

# 1階差分
axes[1].plot(df.index[1:], df['inquiries_diff'].dropna(), color='#2196f3', linewidth=2)
axes[1].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[1].set_title('1階差分データ（定常化済み）', fontsize=14, fontweight='bold')
axes[1].set_ylabel('変化量')
axes[1].set_xlabel('日付')
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

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

# 元データのACF
plot_acf(df['inquiries'], lags=30, ax=axes[0, 0], color='#9c27b0')
axes[0, 0].set_title('元データのACF', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('ラグ')
axes[0, 0].set_ylabel('自己相関')

# 元データのPACF
plot_pacf(df['inquiries'], lags=30, ax=axes[0, 1], color='#9c27b0', method='ywm')
axes[0, 1].set_title('元データのPACF', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('ラグ')
axes[0, 1].set_ylabel('偏自己相関')

# 差分データのACF
plot_acf(df['inquiries_diff'].dropna(), lags=30, ax=axes[1, 0], color='#2196f3')
axes[1, 0].set_title('差分データのACF', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('ラグ')
axes[1, 0].set_ylabel('自己相関')

# 差分データのPACF
plot_pacf(df['inquiries_diff'].dropna(), lags=30, ax=axes[1, 1], color='#2196f3', method='ywm')
axes[1, 1].set_title('差分データのPACF', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('ラグ')
axes[1, 1].set_ylabel('偏自己相関')

plt.tight_layout()
plt.show()

# ACF/PACFの数値確認
acf_values = acf(df['inquiries_diff'].dropna(), nlags=10)
pacf_values = pacf(df['inquiries_diff'].dropna(), nlags=10, method='ywm')

print("=== ACF値（ラグ0-10）===")
for i, val in enumerate(acf_values):
print(f"ラグ {i}: {val:.3f}")

print("\n=== PACF値（ラグ0-10）===")
for i, val in enumerate(pacf_values):
print(f"ラグ {i}: {val:.3f}")

# データの準備（訓練データとテストデータに分割）
train_size = int(len(df) * 0.8) # 80%を訓練、20%をテスト
train = df['inquiries'][:train_size]
test = df['inquiries'][train_size:]

print(f"訓練データ: {len(train)}週")
print(f"テストデータ: {len(test)}週")

# ARIMAパラメータのグリッドサーチ（簡易版）
# p, d, q の候補
p_range = range(0, 3)
d_range = range(0, 2)
q_range = range(0, 3)

best_aic = np.inf
best_params = None
best_model = None

print("\n=== ARIMAパラメータの探索 ===")
for p, d, q in itertools.product(p_range, d_range, q_range):
try:
# モデルの推定
model = ARIMA(train, order=(p, d, q))
fitted_model = model.fit()

# AICを比較
if fitted_model.aic < best_aic: best_aic = fitted_model.aic best_params = (p, d, q) best_model = fitted_model print(f"ARIMA({p},{d},{q}) - AIC: {fitted_model.aic:.2f}") except: continue print(f"\n✅ 最良モデル: ARIMA{best_params}") print(f"AIC: {best_aic:.2f}") # モデルサマリー print("\n" + "="*60) print(best_model.summary()) print("="*60) # 残差診断 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 1. 残差プロット residuals = best_model.resid axes[0, 0].plot(residuals, color='#9c27b0') axes[0, 0].axhline(y=0, color='red', linestyle='--') axes[0, 0].set_title('残差プロット', fontsize=12, fontweight='bold') axes[0, 0].set_ylabel('残差') axes[0, 0].grid(alpha=0.3) # 2. 残差のヒストグラム axes[0, 1].hist(residuals, bins=20, color='#9c27b0', alpha=0.7, edgecolor='black') axes[0, 1].set_title('残差の分布', fontsize=12, fontweight='bold') axes[0, 1].set_xlabel('残差') axes[0, 1].set_ylabel('頻度') axes[0, 1].grid(alpha=0.3) # 3. 残差のACF from statsmodels.graphics.tsaplots import plot_acf plot_acf(residuals, lags=20, ax=axes[1, 0], color='#9c27b0') axes[1, 0].set_title('残差のACF', fontsize=12, fontweight='bold') # 4. Q-Qプロット（正規性の確認） from scipy import stats stats.probplot(residuals, dist="norm", plot=axes[1, 1]) axes[1, 1].set_title('Q-Qプロット（正規性）', fontsize=12, fontweight='bold') axes[1, 1].grid(alpha=0.3) plt.tight_layout() plt.show() # Ljung-Box検定（残差の独立性） from statsmodels.stats.diagnostic import acorr_ljungbox lb_test = acorr_ljungbox(residuals, lags=10, return_df=True) print("\n=== Ljung-Box検定（残差の独立性）===") print(lb_test) print("\np値 > 0.05 なら残差は独立（良いモデル）")

# 予測の実行
forecast_steps = len(test)
forecast = best_model.forecast(steps=forecast_steps)

# 予測区間の計算
forecast_obj = best_model.get_forecast(steps=forecast_steps)
forecast_df = forecast_obj.summary_frame(alpha=0.05) # 95%予測区間

# 可視化
plt.figure(figsize=(14, 6))

# 訓練データ
plt.plot(train.index, train.values, label='訓練データ', color='#9c27b0', linewidth=2)

# テストデータ（実測値）
plt.plot(test.index, test.values, label='テストデータ（実測）',
color='#2196f3', linewidth=2, linestyle='--')

# 予測値
plt.plot(test.index, forecast.values, label='予測値',
color='#ff9800', linewidth=2)

# 95%予測区間
plt.fill_between(test.index,
forecast_df['mean_ci_lower'],
forecast_df['mean_ci_upper'],
color='#ff9800', alpha=0.2, label='95%予測区間')

plt.title(f'ARIMA{best_params}による問い合わせ数予測', fontsize=16, fontweight='bold')
plt.xlabel('日付')
plt.ylabel('問い合わせ数')
plt.legend(loc='best')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# 予測精度の評価
from sklearn.metrics import mean_absolute_error, mean_squared_error

mae = mean_absolute_error(test, forecast)
rmse = np.sqrt(mean_squared_error(test, forecast))
mape = np.mean(np.abs((test - forecast) / test)) * 100

print("\n=== 予測精度 ===")
print(f"MAE（平均絶対誤差）: {mae:.2f}")
print(f"RMSE（平均二乗誤差平方根）: {rmse:.2f}")
print(f"MAPE（平均絶対パーセント誤差）: {mape:.2f}%")

# トレンド
trend_seasonal = np.linspace(100, 180, 156)

# 年次季節性（52週周期で春に急増）
seasonal_yearly = 30 * np.sin(2 * np.pi * np.arange(156) / 52 - np.pi/2)

# ノイズ
noise_seasonal = np.random.normal(0, 8, 156)

# 問い合わせ数
inquiries_seasonal = trend_seasonal + seasonal_yearly + noise_seasonal

df_seasonal = pd.DataFrame({
'date': weeks_seasonal,
'inquiries': inquiries_seasonal
})
df_seasonal.set_index('date', inplace=True)

# 訓練/テストデータ分割
train_size_s = int(len(df_seasonal) * 0.8)
train_s = df_seasonal['inquiries'][:train_size_s]
test_s = df_seasonal['inquiries'][train_size_s:]

# SARIMA実装
# SARIMA(p, d, q)(P, D, Q, s)
# p, d, q: 通常のARIMAパラメータ
# P, D, Q: 季節性パラメータ
# s: 季節周期（週次データで年次季節性なら52）

from statsmodels.tsa.statespace.sarimax import SARIMAX

# モデル1: 季節性なし（比較用）
model_no_season = ARIMA(train_s, order=(1, 1, 1))
fitted_no_season = model_no_season.fit()

# モデル2: 季節性あり（SARIMA）
# (1,1,1)(1,1,1,52) = ARIMA(1,1,1) + 季節性ARIMA(1,1,1,52)
model_sarima = SARIMAX(train_s, order=(1, 1, 1),
seasonal_order=(1, 1, 1, 52))
fitted_sarima = model_sarima.fit(disp=False)

print("=== モデル比較 ===")
print(f"ARIMA(1,1,1) AIC: {fitted_no_season.aic:.2f}")
print(f"SARIMA(1,1,1)(1,1,1,52) AIC: {fitted_sarima.aic:.2f}")
print(f"\n✅ AICが小さいモデル: {'SARIMA' if fitted_sarima.aic < fitted_no_season.aic else 'ARIMA'}") # 予測 forecast_no_season = fitted_no_season.forecast(steps=len(test_s)) forecast_sarima = fitted_sarima.forecast(steps=len(test_s)) # 可視化 plt.figure(figsize=(14, 6)) plt.plot(train_s.index, train_s.values, label='訓練データ', color='#9c27b0', linewidth=2) plt.plot(test_s.index, test_s.values, label='テストデータ（実測）', color='#2196f3', linewidth=2, linestyle='--') plt.plot(test_s.index, forecast_no_season.values, label='ARIMA予測（季節性なし）', color='#ff9800', linewidth=2, alpha=0.7) plt.plot(test_s.index, forecast_sarima, label='SARIMA予測（季節性あり）', color='#4caf50', linewidth=2) # 春の時期を強調表示 for year in [2021, 2022, 2023, 2024]: spring_start = pd.Timestamp(f'{year}-03-01') spring_end = pd.Timestamp(f'{year}-05-01') plt.axvspan(spring_start, spring_end, alpha=0.1, color='orange', label='春シーズン' if year == 2021 else '') plt.title('SARIMA: 季節性を考慮した予測', fontsize=16, fontweight='bold') plt.xlabel('日付') plt.ylabel('問い合わせ数') plt.legend(loc='best') plt.grid(alpha=0.3) plt.tight_layout() plt.show() # 予測精度比較 mae_no_season = mean_absolute_error(test_s, forecast_no_season) mae_sarima = mean_absolute_error(test_s, forecast_sarima) print("\n=== 予測精度比較 ===") print(f"ARIMA MAE: {mae_no_season:.2f}") print(f"SARIMA MAE: {mae_sarima:.2f}") print(f"改善率: {(mae_no_season - mae_sarima) / mae_no_season * 100:.1f}%")

# リスティング広告費（万円/日）
listing_ad = 50 + 10 * np.sin(2 * np.pi * np.arange(365) / 365) + np.random.normal(0, 5, 365)

# SNS広告費（万円/日）- リスティングより少し遅れて変動
sns_ad = 30 + 8 * np.sin(2 * np.pi * np.arange(365) / 365 - 0.5) + np.random.normal(0, 3, 365)

# CV数 - 広告費の影響を3日遅れで受ける
cv_base = 100 + 20 * np.sin(2 * np.pi * np.arange(365) / 365)
cv = np.zeros(365)
cv[:3] = cv_base[:3] + np.random.normal(0, 5, 3)

for i in range(3, 365):
# 3日前の広告費の影響
ad_effect = 0.5 * listing_ad[i-3] + 0.3 * sns_ad[i-3]
cv[i] = cv_base[i] + ad_effect + np.random.normal(0, 5)

# DataFrameに変換
df_var = pd.DataFrame({
'listing_ad': listing_ad,
'sns_ad': sns_ad,
'cv': cv
}, index=days)

# 可視化
fig, axes = plt.subplots(3, 1, figsize=(14, 10))

axes[0].plot(df_var.index, df_var['listing_ad'], color='#9c27b0', linewidth=2)
axes[0].set_title('リスティング広告費', fontsize=12, fontweight='bold')
axes[0].set_ylabel('広告費（万円/日）')
axes[0].grid(alpha=0.3)

axes[1].plot(df_var.index, df_var['sns_ad'], color='#2196f3', linewidth=2)
axes[1].set_title('SNS広告費', fontsize=12, fontweight='bold')
axes[1].set_ylabel('広告費（万円/日）')
axes[1].grid(alpha=0.3)

axes[2].plot(df_var.index, df_var['cv'], color='#4caf50', linewidth=2)
axes[2].set_title('CV数', fontsize=12, fontweight='bold')
axes[2].set_ylabel('CV数/日')
axes[2].set_xlabel('日付')
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.show()

# 定常性の確認（各変数でADF検定）
print("=== 定常性検定（ADF） ===")
for col in df_var.columns:
result = adfuller(df_var[col])
print(f"\n{col}:")
print(f" ADF統計量: {result[0]:.4f}")
print(f" p値: {result[1]:.4f}")
print(f" 判定: {'定常' if result[1] < 0.05 else '非定常'}") # 相関分析 print("\n=== 相関行列 ===") print(df_var.corr()) # VARモデルの推定 from statsmodels.tsa.api import VAR # 訓練/テストデータ分割 train_size_var = int(len(df_var) * 0.8) train_var = df_var[:train_size_var] test_var = df_var[train_size_var:] # VARモデルの作成 model_var = VAR(train_var) # 最適ラグ次数の選択（AIC、BIC基準） lag_order = model_var.select_order(maxlags=10) print("\n=== 最適ラグ次数の選択 ===") print(lag_order.summary()) # 推奨ラグ次数（AIC基準） optimal_lag = lag_order.aic print(f"\n✅ AIC基準の最適ラグ: {optimal_lag}") # VARモデルの推定 fitted_var = model_var.fit(optimal_lag) print("\n" + "="*60) print(fitted_var.summary()) print("="*60)

# 各ペアでグレンジャー因果性を検定
# 検定の方向: X → Y （Xが Yを グレンジャー因果するか）

print("=== グレンジャー因果性検定 ===")
print("\n1. リスティング広告費 → CV数")
result_listing_cv = grangercausalitytests(df_var[['cv', 'listing_ad']],
maxlag=7, verbose=True)

print("\n" + "="*60)
print("\n2. SNS広告費 → CV数")
result_sns_cv = grangercausalitytests(df_var[['cv', 'sns_ad']],
maxlag=7, verbose=True)

print("\n" + "="*60)
print("\n3. CV数 → リスティング広告費（逆方向）")
result_cv_listing = grangercausalitytests(df_var[['listing_ad', 'cv']],
maxlag=7, verbose=True)

# 結果の整理
def summarize_granger(test_result, max_lag=7):
"""グレンジャー因果性検定結果を整理"""
summary = []
for lag in range(1, max_lag + 1):
# F検定のp値を取得
p_value = test_result[lag][0]['ssr_ftest'][1]
summary.append({
'lag': lag,
'p_value': p_value,
'significant': '✅' if p_value < 0.05 else '❌' }) return pd.DataFrame(summary) print("\n=== グレンジャー因果性検定サマリー ===") print("\nリスティング広告費 → CV数:") print(summarize_granger(result_listing_cv)) print("\nSNS広告費 → CV数:") print(summarize_granger(result_sns_cv)) print("\nCV数 → リスティング広告費:") print(summarize_granger(result_cv_listing))

# 各変数ペアでグレンジャー因果性を検定（ラグ3で評価）
variables = ['listing_ad', 'sns_ad', 'cv']
granger_matrix = pd.DataFrame(np.zeros((3, 3)),
index=variables,
columns=variables)

# listing_ad → cv
result = grangercausalitytests(df_var[['cv', 'listing_ad']], maxlag=3, verbose=False)
granger_matrix.loc['listing_ad', 'cv'] = result[3][0]['ssr_ftest'][1] # p値

# sns_ad → cv
result = grangercausalitytests(df_var[['cv', 'sns_ad']], maxlag=3, verbose=False)
granger_matrix.loc['sns_ad', 'cv'] = result[3][0]['ssr_ftest'][1]

# listing_ad → sns_ad
result = grangercausalitytests(df_var[['sns_ad', 'listing_ad']], maxlag=3, verbose=False)
granger_matrix.loc['listing_ad', 'sns_ad'] = result[3][0]['ssr_ftest'][1]

# sns_ad → listing_ad
result = grangercausalitytests(df_var[['listing_ad', 'sns_ad']], maxlag=3, verbose=False)
granger_matrix.loc['sns_ad', 'listing_ad'] = result[3][0]['ssr_ftest'][1]

# cv → listing_ad
result = grangercausalitytests(df_var[['listing_ad', 'cv']], maxlag=3, verbose=False)
granger_matrix.loc['cv', 'listing_ad'] = result[3][0]['ssr_ftest'][1]

# cv → sns_ad
result = grangercausalitytests(df_var[['sns_ad', 'cv']], maxlag=3, verbose=False)
granger_matrix.loc['cv', 'sns_ad'] = result[3][0]['ssr_ftest'][1]

print("=== グレンジャー因果性行列（p値, ラグ3）===")
print(granger_matrix)
print("\np < 0.05 で有意な因果関係") # ネットワーク図の作成 G = nx.DiGraph() G.add_nodes_from(variables) # 有意な因果関係のみエッジを追加（p < 0.05） for i, var1 in enumerate(variables): for j, var2 in enumerate(variables): if i != j and granger_matrix.loc[var1, var2] < 0.05: G.add_edge(var1, var2, weight=1 - granger_matrix.loc[var1, var2]) # 可視化 plt.figure(figsize=(10, 8)) pos = nx.spring_layout(G, seed=42) # ノードの色分け node_colors = {'listing_ad': '#9c27b0', 'sns_ad': '#2196f3', 'cv': '#4caf50'} colors = [node_colors[node] for node in G.nodes()] # 描画 nx.draw_networkx_nodes(G, pos, node_color=colors, node_size=3000, alpha=0.9) nx.draw_networkx_labels(G, pos, font_size=12, font_weight='bold', font_color='white') nx.draw_networkx_edges(G, pos, edge_color='gray', arrows=True, arrowsize=20, width=2, arrowstyle='->')

plt.title('グレンジャー因果性ネットワーク（ラグ3日）', fontsize=16, fontweight='bold')
plt.axis('off')
plt.tight_layout()
plt.show()

# 有意な因果関係のサマリー
print("\n=== 有意なグレンジャー因果関係（p < 0.05）===") for i, var1 in enumerate(variables): for j, var2 in enumerate(variables): if i != j and granger_matrix.loc[var1, var2] < 0.05: print(f"✅ {var1} → {var2} (p={granger_matrix.loc[var1, var2]:.4f})")

# IRFのプロット
fig = irf.plot(orth=False, impulse='listing_ad', response='cv',
figsize=(10, 6))
plt.suptitle('リスティング広告費のショックに対するCV数の応答',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

fig = irf.plot(orth=False, impulse='sns_ad', response='cv',
figsize=(10, 6))
plt.suptitle('SNS広告費のショックに対するCV数の応答',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# 累積IRF（総効果）
irf_cum = irf.cum_effects
print("=== 累積インパルス応答（10日後まで）===")
print("\nリスティング広告費 1万円増加の累積効果:")
print(f"CV数の累積増加: {irf_cum[9, 2, 0]:.2f} 件")

print("\nSNS広告費 1万円増加の累積効果:")
print(f"CV数の累積増加: {irf_cum[9, 2, 1]:.2f} 件")

# IRFを手動でプロット（カスタマイズ版）
plt.figure(figsize=(14, 6))

# リスティング広告費 → CV数
plt.subplot(1, 2, 1)
irf_listing_cv = irf.irfs[:, 2, 0] # CV(index=2), listing_ad(index=0)
plt.plot(range(len(irf_listing_cv)), irf_listing_cv,
marker='o', color='#9c27b0', linewidth=2, markersize=8)
plt.axhline(y=0, color='red', linestyle='--', alpha=0.5)
plt.title('リスティング広告費 → CV数', fontsize=12, fontweight='bold')
plt.xlabel('日数（ラグ）')
plt.ylabel('CV数の変化')
plt.grid(alpha=0.3)

# SNS広告費 → CV数
plt.subplot(1, 2, 2)
irf_sns_cv = irf.irfs[:, 2, 1] # CV(index=2), sns_ad(index=1)
plt.plot(range(len(irf_sns_cv)), irf_sns_cv,
marker='o', color='#2196f3', linewidth=2, markersize=8)
plt.axhline(y=0, color='red', linestyle='--', alpha=0.5)
plt.title('SNS広告費 → CV数', fontsize=12, fontweight='bold')
plt.xlabel('日数（ラグ）')
plt.ylabel('CV数の変化')
plt.grid(alpha=0.3)

plt.tight_layout()
plt.show()

# 実務的解釈
print("\n=== 実務的解釈 ===")
print(f"リスティング広告費を1万円増やすと:")
print(f" - 当日のCV増加: {irf_listing_cv[0]:.2f} 件")
print(f" - 3日後のCV増加: {irf_listing_cv[3]:.2f} 件")
print(f" - 10日間の累積CV増加: {irf_cum[9, 2, 0]:.2f} 件")
print(f"\nSNS広告費を1万円増やすと:")
print(f" - 当日のCV増加: {irf_sns_cv[0]:.2f} 件")
print(f" - 3日後のCV増加: {irf_sns_cv[3]:.2f} 件")
print(f" - 10日間の累積CV増加: {irf_cum[9, 2, 1]:.2f} 件")

# 予測結果をDataFrameに変換
forecast_df = pd.DataFrame(forecast_var,
index=test_var.index,
columns=train_var.columns)

# 可視化
fig, axes = plt.subplots(3, 1, figsize=(14, 12))

for i, col in enumerate(df_var.columns):
# 訓練データ
axes[i].plot(train_var.index, train_var[col],
label='訓練データ', color='#9c27b0', linewidth=2)

# テストデータ（実測）
axes[i].plot(test_var.index, test_var[col],
label='テストデータ（実測）', color='#2196f3',
linewidth=2, linestyle='--')

# 予測値
axes[i].plot(forecast_df.index, forecast_df[col],
label='予測値', color='#ff9800', linewidth=2)

axes[i].set_title(f'{col} の予測', fontsize=12, fontweight='bold')
axes[i].set_ylabel(col)
axes[i].legend(loc='best')
axes[i].grid(alpha=0.3)

axes[2].set_xlabel('日付')
plt.tight_layout()
plt.show()

# 予測精度の評価
print("=== VAR予測精度 ===")
for col in df_var.columns:
mae = mean_absolute_error(test_var[col], forecast_df[col])
rmse = np.sqrt(mean_squared_error(test_var[col], forecast_df[col]))
mape = np.mean(np.abs((test_var[col] - forecast_df[col]) / test_var[col])) * 100

print(f"\n{col}:")
print(f" MAE: {mae:.2f}")
print(f" RMSE: {rmse:.2f}")
print(f" MAPE: {mape:.2f}%")

# シナリオ分析: 広告費を20%増やしたら
print("\n=== シナリオ分析 ===")
print("リスティング広告費を20%増やした場合のCV予測:")

# 新しいシナリオデータ
scenario_data = train_var.copy()
scenario_data['listing_ad'] = scenario_data['listing_ad'] * 1.2

# シナリオ予測
forecast_scenario = fitted_var.forecast(scenario_data.values[-optimal_lag:],
steps=30)
forecast_scenario_df = pd.DataFrame(forecast_scenario,
columns=train_var.columns)

print(f"通常時の30日間平均CV: {forecast_var[:30, 2].mean():.1f} 件/日")
print(f"広告費20%増の30日間平均CV: {forecast_scenario_df['cv'].mean():.1f} 件/日")
print(f"CV増加率: {(forecast_scenario_df['cv'].mean() / forecast_var[:30, 2].mean() - 1) * 100:.1f}%")

# 季節性のあるデータを使用（Week 1で作成したdf_seasonal）
data_ucm = df_seasonal['inquiries']

# モデル1: トレンド + 不規則成分
model_trend = UnobservedComponents(data_ucm, level='local linear trend')
fitted_trend = model_trend.fit(disp=False)

# モデル2: トレンド + 季節性（52週周期）
model_seasonal = UnobservedComponents(
data_ucm,
level='local linear trend',
seasonal=52 # 52週周期の季節性
)
fitted_seasonal = model_seasonal.fit(disp=False)

# モデル比較
print("=== モデル比較（AIC） ===")
print(f"トレンドのみ: {fitted_trend.aic:.2f}")
print(f"トレンド + 季節性: {fitted_seasonal.aic:.2f}")
print(f"\n✅ 最良モデル: {'季節性あり' if fitted_seasonal.aic < fitted_trend.aic else 'トレンドのみ'}") # 成分の抽出 trend = fitted_seasonal.level.smoothed seasonal = fitted_seasonal.seasonal.smoothed # 可視化 fig, axes = plt.subplots(4, 1, figsize=(14, 12)) # 元データ axes[0].plot(data_ucm.index, data_ucm.values, color='#9c27b0', linewidth=2) axes[0].set_title('元データ', fontsize=12, fontweight='bold') axes[0].set_ylabel('問い合わせ数') axes[0].grid(alpha=0.3) # トレンド成分 axes[1].plot(data_ucm.index, trend, color='#2196f3', linewidth=2) axes[1].set_title('トレンド成分（抽出）', fontsize=12, fontweight='bold') axes[1].set_ylabel('トレンド') axes[1].grid(alpha=0.3) # 季節性成分 axes[2].plot(data_ucm.index, seasonal, color='#ff9800', linewidth=2) axes[2].set_title('季節性成分（抽出）', fontsize=12, fontweight='bold') axes[2].set_ylabel('季節性') axes[2].grid(alpha=0.3) # 残差 residuals = data_ucm - trend - seasonal axes[3].plot(data_ucm.index, residuals, color='#4caf50', linewidth=2) axes[3].axhline(y=0, color='red', linestyle='--', alpha=0.5) axes[3].set_title('残差（不規則成分）', fontsize=12, fontweight='bold') axes[3].set_ylabel('残差') axes[3].set_xlabel('日付') axes[3].grid(alpha=0.3) plt.tight_layout() plt.show() # サマリー print("\n" + "="*60) print(fitted_seasonal.summary()) print("="*60) # 予測 forecast_ucm = fitted_seasonal.forecast(steps=52) # 1年先まで予測 forecast_ci = fitted_seasonal.get_forecast(steps=52).summary_frame(alpha=0.05) # 予測の可視化 plt.figure(figsize=(14, 6)) plt.plot(data_ucm.index, data_ucm.values, label='観測値', color='#9c27b0', linewidth=2) plt.plot(forecast_ucm.index, forecast_ucm.values, label='予測値', color='#ff9800', linewidth=2) plt.fill_between(forecast_ucm.index, forecast_ci['mean_ci_lower'], forecast_ci['mean_ci_upper'], color='#ff9800', alpha=0.2, label='95%予測区間') plt.title('UCMによる1年先予測', fontsize=16, fontweight='bold') plt.xlabel('日付') plt.ylabel('問い合わせ数') plt.legend(loc='best') plt.grid(alpha=0.3) plt.tight_layout() plt.show()

# 構造変化を含むデータの生成
np.random.seed(42)
n = 200
time = np.arange(n)

# 100時点目で構造変化（傾きが変わる）
y = np.zeros(n)
y[:100] = 50 + 0.5 * time[:100] + np.random.normal(0, 3, 100)
y[100:] = 50 + 1.5 * time[100:] + np.random.normal(0, 3, 100)

df_break = pd.DataFrame({
'time': time,
'y': y
})

# 可視化
plt.figure(figsize=(12, 6))
plt.plot(df_break['time'], df_break['y'], color='#9c27b0', linewidth=2)
plt.axvline(x=100, color='red', linestyle='--', linewidth=2,
label='構造変化点（真値）')
plt.title('構造変化を含む時系列データ', fontsize=14, fontweight='bold')
plt.xlabel('時点')
plt.ylabel('値')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

# Chow検定（特定の時点での構造変化を検定）
# 仮説: 100時点目で構造変化があるか

# 期間1: 0-99
# 期間2: 100-199
period1 = df_break[:100]
period2 = df_break[100:]

# 全期間でのOLS
model_full = OLS(df_break['y'], sm.add_constant(df_break['time'])).fit()

# 期間1でのOLS
model_p1 = OLS(period1['y'], sm.add_constant(period1['time'])).fit()

# 期間2でのOLS
model_p2 = OLS(period2['y'], sm.add_constant(period2['time'])).fit()

# Chow検定統計量
rss_full = model_full.ssr
rss_p1 = model_p1.ssr
rss_p2 = model_p2.ssr
rss_restricted = rss_p1 + rss_p2

k = 2 # パラメータ数（切片+傾き）
n1, n2 = len(period1), len(period2)

chow_stat = ((rss_full - rss_restricted) / k) / (rss_restricted / (n1 + n2 - 2*k))

# F分布でのp値計算
from scipy.stats import f as f_dist
p_value = 1 - f_dist.cdf(chow_stat, k, n1 + n2 - 2*k)

print("=== Chow検定結果 ===")
print(f"Chow統計量: {chow_stat:.4f}")
print(f"p値: {p_value:.6f}")
print(f"判定: {'構造変化あり ✅' if p_value < 0.05 else '構造変化なし'}") # 期間別の傾きを比較 print(f"\n=== 期間別パラメータ ===") print(f"期間1（0-99）の傾き: {model_p1.params[1]:.4f}") print(f"期間2（100-199）の傾き: {model_p2.params[1]:.4f}") # CUSUM検定（構造変化点を探索） import statsmodels.api as sm # トレンド回帰 X = sm.add_constant(df_break['time']) model_cusum = OLS(df_break['y'], X).fit() resids = model_cusum.resid # CUSUMプロット from statsmodels.stats.diagnostic import recursive_ols # 再帰的残差のプロット res_ols = recursive_ols(model_cusum) fig, ax = plt.subplots(figsize=(12, 6)) ax.plot(df_break['time'], res_ols.cusum, color='#9c27b0', linewidth=2) ax.axhline(y=0, color='black', linestyle='-', alpha=0.3) # 5%有意水準の信頼区間 n_points = len(df_break) critical_value = 0.948 # 5%水準のCUSUM臨界値 upper = critical_value * np.sqrt(n_points) lower = -critical_value * np.sqrt(n_points) ax.plot(df_break['time'], [upper]*n_points, 'r--', label='5%信頼区間') ax.plot(df_break['time'], [lower]*n_points, 'r--') ax.axvline(x=100, color='orange', linestyle='--', linewidth=2, label='真の構造変化点') ax.set_title('CUSUM検定（構造変化の検出）', fontsize=14, fontweight='bold') ax.set_xlabel('時点') ax.set_ylabel('CUSUM統計量') ax.legend() ax.grid(alpha=0.3) plt.tight_layout() plt.show() print("\n=== CUSUM検定の解釈 ===") print("CUSUM統計量が信頼区間を超えた時点 → 構造変化の可能性")

import lingam
from lingam.utils import make_dot

# シンプルな因果構造のデータ生成
# 真の因果関係: X1 → X2 → X3
np.random.seed(42)
n = 1000

# X1は外生変数
X1 = np.random.laplace(0, 1, n)

# X2 = 2*X1 + ノイズ
X2 = 2 * X1 + np.random.laplace(0, 1, n)

# X3 = -1.5*X2 + ノイズ
X3 = -1.5 * X2 + np.random.laplace(0, 1, n)

# データフレーム化
data_lingam = pd.DataFrame({
'X1': X1,
'X2': X2,
'X3': X3
})

# DirectLiNGAMで因果関係を推定
model_lingam = lingam.DirectLiNGAM()
model_lingam.fit(data_lingam)

# 因果順序の推定結果
print("=== 推定された因果順序 ===")
print(f"因果順序: {model_lingam.causal_order_}")
print("（0=X1, 1=X2, 2=X3）")

# 隣接行列（因果係数）
print("\n=== 隣接行列（因果係数）===")
print(pd.DataFrame(model_lingam.adjacency_matrix_,
columns=data_lingam.columns,
index=data_lingam.columns))

print("\n解釈:")
print(" X1 → X2 の係数: {:.2f} （真値: 2.0）".format(
model_lingam.adjacency_matrix_[1, 0]))
print(" X2 → X3 の係数: {:.2f} （真値: -1.5）".format(
model_lingam.adjacency_matrix_[2, 1]))

# 因果グラフの可視化（graphviz必要）
try:
dot = make_dot(model_lingam.adjacency_matrix_, labels=data_lingam.columns)
dot.render('lingam_graph', format='png', cleanup=True)
print("\n✅ 因果グラフを 'lingam_graph.png' に保存しました")
except:
print("\n⚠️ graphvizがインストールされていないため、グラフ可視化はスキップ")
print(" グラフ可視化が必要な場合: pip install graphviz")

print("\n=== LiNGAMの注意点 ===")
print("1. 非ガウス性が前提（データがほぼ正規分布だと失敗）")
print("2. 線形関係が前提（非線形な因果関係は検出できない）")
print("3. 非巡回グラフが前提（フィードバックループがあると失敗）")
print("4. Phase 4-5で詳しく学習します")

# ステップ1: データ準備
print("=== ステップ1: データ準備 ===")
# （実際は自社データを読み込む）
# df = pd.read_csv('marketing_data.csv', parse_dates=['date'])

# ステップ2: 探索的データ分析
print("\n=== ステップ2: EDA ===")
# - 時系列プロット
# - 基本統計量
# - 相関分析
# - 定常性チェック

# ステップ3: 単変量時系列分析（ARIMA/SARIMA）
print("\n=== ステップ3: CV数のSARIMA予測 ===")
# - ADF検定
# - ACF/PACF分析
# - パラメータ選択
# - モデル推定と予測

# ステップ4: 多変量時系列分析（VAR）
print("\n=== ステップ4: VAR分析 ===")
# - 最適ラグ選択
# - VAR推定
# - グレンジャー因果性検定
# - IRF分析

# ステップ5: 成分分解（UCM）
print("\n=== ステップ5: トレンド・季節性分解 ===")
# - UCMでトレンド抽出
# - 季節性パターンの確認
# - 構造変化の検出

# ステップ6: ビジネス示唆
print("\n=== ステップ6: ビジネス示唆 ===")
print("1. 広告効果の遅れ: リスティングは3日後、SNSは5日後にピーク")
print("2. 効果測定期間: 広告効果は10日間持続するため、最低2週間で評価")
print("3. 予算配分: リスティングの累積効果が高いため、予算シェア60%推奨")
print("4. 季節性: 3-4月は通常の1.5倍の広告費投下が効果的")
print("5. 長期トレンド: 市場全体が年率5%成長、追い風状態")