"""
教育用再現コード: 2021年 統計データ分析コンペティション 総務大臣賞 [大学生・一般の部]
=================================================================
論文タイトル：COVID-19が地域経済に与えた影響：消費支出・雇用指標のパネル分析
受賞：総務大臣賞（大学生・一般の部）— 最優秀賞

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度）
  対象：全47都道府県 × 12年（2012〜2023）

  COVID-19ダミー変数: 2020年・2021年 = 1（COVID期）, それ以外 = 0

  分析の流れ
  1. 時系列：COVID前後（2019→2020）の消費支出急落を可視化
  2. 差の差（DiD）比較：観光依存度が高い地域 vs 低い地域の消費変動
  3. パネルOLS固定効果：COVID期ダミーが消費支出・求人倍率に与えた効果
  4. 地域別COVID影響スコア：都道府県別 2019→2020年の消費変化率ランキング

【被説明変数】
  消費支出（二人以上の世帯）（log変換）

【主要変数】
  COVID期ダミー = (年度 == 2020 or 2021).astype(int)
  求人倍率 = 月間有効求人数（一般）/ 月間有効求職者数（一般）
  宿泊密度 = 延べ宿泊者数 / 総人口（観光依存度の代理）
  高齢化率 = 65歳以上人口 / 総人口 × 100

【推定手法】
  linearmodels PanelOLS: entity_effects=True, time_effects=True, cov_type='clustered'

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）

【データサイエンス学習ポイント】
  1. 自然実験としてのCOVID-19（ダミー変数を用いたパネル分析）
  2. 差の差推定（DiD）の考え方と並行トレンド仮定
  3. 個体効果と時間効果の両方を制御した双方向固定効果モデル
  4. 観光依存度指標の設計と政策的示唆
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2021_U1_daijin.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from linearmodels.panel import PanelOLS
import statsmodels.api as sm
from scipy import stats

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

import os
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

# 変数生成
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['消費支出_万円'] = df_b['消費支出（二人以上の世帯）'] / 10000
df_b['求人倍率'] = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'].replace(0, np.nan)
df_b['宿泊密度'] = df_b['延べ宿泊者数'] / df_b['総人口'].replace(0, np.nan)
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['COVID期'] = df_b['年度'].isin([2020, 2021]).astype(int)

region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
df_b['地域'] = df_b['都道府県'].map(region_map)
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12',
}

# Fig1: 消費支出の時系列（COVID前後の急落）
fig, ax = plt.subplots(figsize=(10, 5))
yearly = df_b.groupby(['年度', '地域'])['消費支出_万円'].mean().reset_index()
for reg, grp in yearly.groupby('地域'):
    ax.plot(grp['年度'], grp['消費支出_万円'], marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.axvspan(2019.5, 2021.5, alpha=0.1, color='red', label='COVID期（2020-2021年）')
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('消費支出（万円/月）', fontsize=12)
ax.set_title('地域別 消費支出の推移（2012〜2023年）\nCOVID-19の影響（2020-2021年の急落）', fontsize=13, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_U1_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: 求人倍率の時系列（COVID期の急落）
fig, ax = plt.subplots(figsize=(10, 5))
yearly_k = df_b.groupby(['年度', '地域'])['求人倍率'].mean().reset_index()
for reg, grp in yearly_k.groupby('地域'):
    ax.plot(grp['年度'], grp['求人倍率'], marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.axvspan(2019.5, 2021.5, alpha=0.1, color='red', label='COVID期')
ax.axhline(1.0, color='black', linestyle='--', linewidth=0.8, label='倍率=1.0')
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('有効求人倍率', fontsize=12)
ax.set_title('地域別 有効求人倍率の推移（2012〜2023年）\nCOVID-19の雇用への影響', fontsize=13, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_U1_fig2_koyou_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: パネルOLS固定効果（双方向）— COVID期ダミーの効果
panel_vars = ['COVID期', '求人倍率', '宿泊密度', '高齢化率']
df_panel = df_b[['年度', '都道府県', '消費支出_log'] + panel_vars].dropna()
df_panel = df_panel.set_index(['都道府県', '年度'])
y = df_panel['消費支出_log']
X = sm.add_constant(df_panel[panel_vars])
try:
    mod = PanelOLS(y, X, entity_effects=True, time_effects=True, drop_absorbed=True)
    res = mod.fit(cov_type='clustered', cluster_entity=True)
    coefs = res.params.drop('const', errors='ignore')
    ses = res.std_errors.drop('const', errors='ignore')
    pvals = res.pvalues.drop('const', errors='ignore')
    print(res.summary)
except Exception as e:
    print(f"2-way FE error: {e}")
    mod = PanelOLS(y, X, entity_effects=True, drop_absorbed=True)
    res = mod.fit(cov_type='clustered', cluster_entity=True)
    coefs = res.params.drop('const', errors='ignore')
    ses = res.std_errors.drop('const', errors='ignore')
    pvals = res.pvalues.drop('const', errors='ignore')

fig, ax = plt.subplots(figsize=(8, 5))
colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
        error_kw={'elinewidth': 1.5, 'capsize': 4})
ax.set_yticks(range(len(coefs)))
ax.set_yticklabels(coefs.index, fontsize=10)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('双方向固定効果推定係数', fontsize=12)
ax.set_title('消費支出への影響 — 双方向FEパネルOLS\n（赤=p<0.05, 灰=非有意）', fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_U1_fig3_fe_coef.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: 都道府県別 COVID影響スコア（2019→2020 消費変化率ランキング）
df_2019 = df_b[df_b['年度'] == 2019][['都道府県', '消費支出_万円', '地域']].copy()
df_2020 = df_b[df_b['年度'] == 2020][['都道府県', '消費支出_万円']].copy()
df_covid = df_2019.merge(df_2020, on='都道府県', suffixes=('_2019', '_2020'))
df_covid['変化率'] = (df_covid['消費支出_万円_2020'] - df_covid['消費支出_万円_2019']) / df_covid['消費支出_万円_2019'] * 100
df_covid = df_covid.sort_values('変化率', ascending=True)

fig, ax = plt.subplots(figsize=(8, 12))
bar_colors = [region_colors.get(r, 'gray') for r in df_covid['地域']]
ax.barh(df_covid['都道府県'], df_covid['変化率'], color=bar_colors, alpha=0.8)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('消費支出変化率（2019→2020年, %）', fontsize=12)
ax.set_title('COVID-19による消費支出変化率ランキング\n（2019→2020年、都道府県別）', fontsize=13, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_U1_fig4_covid_rank.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")