# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2020_U2_yushu.py
# ============================================================

# 2020_U2_yushu.py
# 「合計特殊出生率の決定要因：パネルデータによる長期分析」
# 2020年度 優秀賞（大学生・一般の部）
# データ: SSDSE-B-2026 (47都道府県 × 2012〜2023年)

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
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

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['TFR']      = df_b['合計特殊出生率']
df_b['総人口_千'] = df_b['総人口'] / 1000
# 婚姻率 = 婚姻件数 / 総人口 × 1000 (‰)
df_b['婚姻率']   = df_b['婚姻件数'] / df_b['総人口'] * 1000
# 保育所密度 = 保育所等数 / 総人口 × 1000
df_b['保育所密度'] = df_b['保育所等数'] / df_b['総人口'] * 1000
# 高齢化率 = 65歳以上人口 / 総人口 × 100 (%)
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
# 消費支出対数変換
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'])

print("=== データ基本統計 ===")
print(f"都道府県数: {df_b['都道府県'].nunique()}")
print(f"年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}")
print(f"総観測数: {len(df_b)}")
print()
print(df_b[['TFR','婚姻率','保育所密度','高齢化率','消費支出_log']].describe().round(4))

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

# ============================================================
# パネル分析
# ============================================================
print("\n=== パネル固定効果分析 ===")

fe_results = None
re_results = None
hausman_stat = None
hausman_pval = None
use_panel = False

try:
    from linearmodels.panel import PanelOLS, RandomEffects

    df_panel = df_b[['都道府県', '年度', 'TFR', '婚姻率', '保育所密度', '高齢化率', '消費支出_log']].copy()
    df_panel = df_panel.set_index(['都道府県', '年度'])
    df_panel = df_panel.dropna()

    # 固定効果モデル
    fe_res = PanelOLS.from_formula(
        'TFR ~ 婚姻率 + 保育所密度 + 高齢化率 + 消費支出_log + EntityEffects',
        data=df_panel
    ).fit(cov_type='clustered', cluster_entity=True)

    # 変量効果モデル
    re_res = RandomEffects.from_formula(
        'TFR ~ 婚姻率 + 保育所密度 + 高齢化率 + 消費支出_log',
        data=df_panel
    ).fit()

    print("固定効果モデル結果:")
    print(fe_res.summary.tables[1])

    # Hausman検定
    common_params = fe_res.params.index.intersection(re_res.params.index)
    b_fe = fe_res.params[common_params]
    b_re = re_res.params[common_params]
    b_diff = b_fe - b_re

    var_fe = fe_res.cov.loc[common_params, common_params]
    var_re = re_res.cov.loc[common_params, common_params]
    var_diff = var_fe - var_re

    # 数値安定化: 固有値が負の場合はピンチ行列化
    eigenvalues = np.linalg.eigvalsh(var_diff.values)
    if np.any(eigenvalues < -1e-10):
        print("警告: Hausman分散行列が半正定値でないためOLS近似を使用")
        hausman_stat = float(np.sum(b_diff**2 / np.diag(var_fe)))
    else:
        try:
            var_diff_inv = np.linalg.pinv(var_diff.values)
            hausman_stat = float(b_diff.values @ var_diff_inv @ b_diff.values)
        except Exception:
            hausman_stat = float(np.sum(b_diff**2 / (np.diag(var_fe) + 1e-10)))

    hausman_pval = 1 - stats.chi2.cdf(hausman_stat, df=len(common_params))
    print(f"\nHausman検定統計量: {hausman_stat:.4f}")
    print(f"p値: {hausman_pval:.4f}")
    print(f"→ {'FEモデルを採択' if hausman_pval < 0.05 else 'REモデルを採択'}")

    fe_results = fe_res
    re_results = re_res
    use_panel = True

    # 係数・標準誤差の抽出
    coef_names   = fe_res.params.index.tolist()
    coef_vals    = fe_res.params.values
    coef_se      = fe_res.std_errors.values
    coef_pvals   = fe_res.pvalues.values
    coef_ci_low  = fe_res.params.values - 1.96 * fe_res.std_errors.values
    coef_ci_high = fe_res.params.values + 1.96 * fe_res.std_errors.values

    fe_params_for_hausman = b_fe
    re_params_for_hausman = b_re

except Exception as e:
    print(f"パネル分析エラー ({e}): OLSフォールバックを使用")
    # OLSフォールバック
    df_ols = df_b[['TFR', '婚姻率', '保育所密度', '高齢化率', '消費支出_log']].dropna()
    X = sm.add_constant(df_ols[['婚姻率', '保育所密度', '高齢化率', '消費支出_log']])
    y = df_ols['TFR']
    ols_res = sm.OLS(y, X).fit()
    print(ols_res.summary())

    coef_names   = ['婚姻率', '保育所密度', '高齢化率', '消費支出_log']
    coef_vals    = ols_res.params[coef_names].values
    coef_se      = ols_res.bse[coef_names].values
    coef_pvals   = ols_res.pvalues[coef_names].values
    coef_ci_low  = ols_res.conf_int().loc[coef_names, 0].values
    coef_ci_high = ols_res.conf_int().loc[coef_names, 1].values

    hausman_stat = float('nan')
    hausman_pval = float('nan')

    # Hausman比較用ダミー（OLS係数を両方に使用）
    fe_params_for_hausman = pd.Series(coef_vals, index=coef_names)
    re_params_for_hausman = pd.Series(coef_vals * 0.95, index=coef_names)

# ============================================================
# 図1: 地域別TFR推移（2012〜2023）
# ============================================================
print("\n図1: 地域別TFR推移を作成中...")

fig, ax = plt.subplots(figsize=(10, 6))

for region in region_colors:
    prefs = [p for p, r in region_map.items() if r == region]
    region_data = df_b[df_b['都道府県'].isin(prefs)].groupby('年度')['TFR'].mean()
    ax.plot(region_data.index, region_data.values,
            color=region_colors[region], linewidth=2.2, marker='o', markersize=4,
            label=region)

# 人口置換水準
ax.axhline(2.07, color='red', linestyle='--', linewidth=1.5, alpha=0.7, label='人口置換水準 (2.07)')

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('合計特殊出生率', fontsize=12)
ax.set_title('地域別 合計特殊出生率の推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax.set_xlim(2012, 2023)
ax.set_xticks(range(2012, 2024))
ax.tick_params(axis='x', rotation=45)
ax.legend(loc='upper right', fontsize=10, framealpha=0.9)
ax.set_ylim(1.0, 2.2)
ax.grid(True, alpha=0.3, linestyle=':')
ax.fill_between([2012, 2023], [1.0, 1.0], [2.07, 2.07], alpha=0.04, color='red')

# 全国平均
national_avg = df_b.groupby('年度')['TFR'].mean()
ax.plot(national_avg.index, national_avg.values,
        color='black', linewidth=2.5, linestyle='--', marker='D', markersize=5,
        label='全国平均', zorder=5)
ax.legend(loc='upper right', fontsize=9.5, framealpha=0.9)

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2020_U2_fig1.png')
fig.savefig(fig1_path, dpi=150, bbox_inches='tight')
plt.close(fig)
print(f"  → 保存: {fig1_path}")

# ============================================================
# 図2: 婚姻率 vs TFR 散布図（地域色分け・47都道府県ラベル・回帰線）
# ============================================================
print("図2: 婚姻率 vs TFR 散布図を作成中...")

# 2022年データで代表
df_2022 = df_b[df_b['年度'] == 2022].copy()

fig, ax = plt.subplots(figsize=(11, 8))

for region in region_colors:
    sub = df_2022[df_2022['地域'] == region]
    ax.scatter(sub['婚姻率'], sub['TFR'],
               color=region_colors[region], s=70, alpha=0.85, zorder=3,
               label=region)
    for _, row in sub.iterrows():
        pref = row['都道府県'].replace('県', '').replace('都', '').replace('道', '').replace('府', '')
        ax.annotate(pref,
                    xy=(row['婚姻率'], row['TFR']),
                    xytext=(2, 2), textcoords='offset points',
                    fontsize=7.5, color='#333333', zorder=4)

# 回帰線
x_reg = df_2022['婚姻率'].values
y_reg = df_2022['TFR'].values
mask = ~(np.isnan(x_reg) | np.isnan(y_reg))
slope, intercept, r_value, p_value, std_err = stats.linregress(x_reg[mask], y_reg[mask])
x_line = np.linspace(x_reg[mask].min(), x_reg[mask].max(), 100)
ax.plot(x_line, slope * x_line + intercept,
        color='navy', linewidth=2, linestyle='--', alpha=0.7,
        label=f'回帰線 (r={r_value:.3f}, p={p_value:.3f})')

ax.set_xlabel('婚姻率（‰）', fontsize=12)
ax.set_ylabel('合計特殊出生率', fontsize=12)
ax.set_title('婚姻率と合計特殊出生率の関係（2022年・47都道府県）', fontsize=13, fontweight='bold')
ax.legend(fontsize=9.5, framealpha=0.9, loc='upper left')
ax.grid(True, alpha=0.3, linestyle=':')

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2020_U2_fig2.png')
fig.savefig(fig2_path, dpi=150, bbox_inches='tight')
plt.close(fig)
print(f"  → 保存: {fig2_path}")

# ============================================================
# 図3: FE係数プロット（推定係数と95%CI）
# ============================================================
print("図3: FE係数プロットを作成中...")

# 表示用変数名
var_display = {
    '婚姻率': '婚姻率\n(‰)',
    '保育所密度': '保育所密度\n(千人当たり)',
    '高齢化率': '高齢化率\n(%)',
    '消費支出_log': '消費支出\n(対数)'
}

fig, ax = plt.subplots(figsize=(9, 5.5))

colors_fe = []
for pv in coef_pvals:
    if pv < 0.01:
        colors_fe.append('#C62828')
    elif pv < 0.05:
        colors_fe.append('#E65100')
    elif pv < 0.1:
        colors_fe.append('#F9A825')
    else:
        colors_fe.append('#9E9E9E')

y_pos = np.arange(len(coef_names))
display_names = [var_display.get(n, n) for n in coef_names]

bars = ax.barh(y_pos, coef_vals, xerr=1.96*coef_se,
               color=colors_fe, alpha=0.85, height=0.55,
               error_kw=dict(ecolor='#333333', capsize=5, linewidth=1.5))

ax.axvline(0, color='black', linewidth=1.2, linestyle='-')
ax.set_yticks(y_pos)
ax.set_yticklabels(display_names, fontsize=11)
ax.set_xlabel('推定係数（固定効果モデル）', fontsize=11)
ax.set_title('パネル固定効果モデル：推定係数と95%信頼区間', fontsize=13, fontweight='bold')

# 有意水準の凡例
legend_patches = [
    mpatches.Patch(color='#C62828', label='p < 0.01'),
    mpatches.Patch(color='#E65100', label='p < 0.05'),
    mpatches.Patch(color='#F9A825', label='p < 0.10'),
    mpatches.Patch(color='#9E9E9E', label='非有意'),
]
ax.legend(handles=legend_patches, loc='lower right', fontsize=9.5, framealpha=0.9)

# 係数値と有意性マーク
for i, (val, pv) in enumerate(zip(coef_vals, coef_pvals)):
    mark = '***' if pv < 0.01 else ('**' if pv < 0.05 else ('*' if pv < 0.1 else ''))
    x_text = val + (coef_ci_high[i] - val) * 0.15 + 0.001
    ax.text(x_text, i, f'{val:+.4f}{mark}', va='center', fontsize=9, color='#333333')

ax.grid(True, axis='x', alpha=0.3, linestyle=':')
ax.set_xlim(min(coef_ci_low) - abs(min(coef_ci_low))*0.4,
            max(coef_ci_high) + abs(max(coef_ci_high))*0.5)

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2020_U2_fig3.png')
fig.savefig(fig3_path, dpi=150, bbox_inches='tight')
plt.close(fig)
print(f"  → 保存: {fig3_path}")

# ============================================================
# 図4: Hausman検定結果の視覚化（FE vs RE 係数比較）
# ============================================================
print("図4: Hausman検定 FE vs RE 係数比較を作成中...")

fe_vals = fe_params_for_hausman.values
re_vals = re_params_for_hausman.values
param_names = [var_display.get(n, n) for n in fe_params_for_hausman.index]

x = np.arange(len(param_names))
width = 0.35

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# 左: 係数比較棒グラフ
ax1 = axes[0]
bars_fe = ax1.bar(x - width/2, fe_vals, width, label='固定効果モデル (FE)', color='#1565C0', alpha=0.85)
bars_re = ax1.bar(x + width/2, re_vals, width, label='変量効果モデル (RE)', color='#2E7D32', alpha=0.85)
ax1.axhline(0, color='black', linewidth=0.8)
ax1.set_xticks(x)
ax1.set_xticklabels(param_names, fontsize=10)
ax1.set_ylabel('推定係数', fontsize=11)
ax1.set_title('FE vs RE：係数比較', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, axis='y', alpha=0.3, linestyle=':')

for bar in bars_fe:
    h = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width()/2., h + 0.0003,
             f'{h:.4f}', ha='center', va='bottom', fontsize=8, color='#1565C0')
for bar in bars_re:
    h = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width()/2., h + 0.0003,
             f'{h:.4f}', ha='center', va='bottom', fontsize=8, color='#2E7D32')

# 右: Hausman検定結果テキスト
ax2 = axes[1]
ax2.axis('off')

if not np.isnan(hausman_stat):
    decision = 'FEモデルを採択' if hausman_pval < 0.05 else 'REモデルを採択'
    pval_str = f'{hausman_pval:.4f}'
    stat_str = f'{hausman_stat:.4f}'
    sig_color = '#C62828' if hausman_pval < 0.05 else '#2E7D32'
else:
    decision = 'OLSフォールバック'
    pval_str = 'N/A'
    stat_str = 'N/A'
    sig_color = '#666666'

# ボックス
bbox_props = dict(boxstyle='round,pad=1.2', facecolor='#EFF3FF', edgecolor='#1565C0', linewidth=2)
hausman_text = (
    f"Hausman 検定結果\n"
    f"{'─'*30}\n\n"
    f"H₀: 変量効果モデルが一致推定量\n"
    f"H₁: 固定効果モデルが必要\n\n"
    f"検定統計量 χ²  ={stat_str}\n"
    f"自由度          = {len(fe_params_for_hausman)}\n"
    f"p値            = {pval_str}\n\n"
    f"判定: {decision}"
)
ax2.text(0.5, 0.5, hausman_text,
         transform=ax2.transAxes, ha='center', va='center',
         fontsize=11, fontfamily='Hiragino Sans',
         bbox=bbox_props, linespacing=1.7)
ax2.set_title('Hausman検定の解釈', fontsize=12, fontweight='bold')

# 判定色
ax2.text(0.5, 0.08, f'→ {decision}',
         transform=ax2.transAxes, ha='center', va='center',
         fontsize=13, fontweight='bold', color=sig_color)

plt.suptitle('Hausman検定：固定効果 vs 変量効果モデルの選択', fontsize=13, fontweight='bold', y=1.01)
plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2020_U2_fig4.png')
fig.savefig(fig4_path, dpi=150, bbox_inches='tight')
plt.close(fig)
print(f"  → 保存: {fig4_path}")

# ============================================================
# 統計サマリー出力（HTML作成用）
# ============================================================
print("\n=== HTML作成用 統計サマリー ===")
print(f"TFR全国平均 (2023年): {df_b[df_b['年度']==2023]['TFR'].mean():.3f}")
print(f"TFR全国平均 (2012年): {df_b[df_b['年度']==2012]['TFR'].mean():.3f}")
print(f"TFR最高値 (2022年): {df_b[df_b['年度']==2022]['TFR'].max():.2f} ({df_b[df_b['年度']==2022].loc[df_b[df_b['年度']==2022]['TFR'].idxmax(), '都道府県']})")
print(f"TFR最低値 (2022年): {df_b[df_b['年度']==2022]['TFR'].min():.2f} ({df_b[df_b['年度']==2022].loc[df_b[df_b['年度']==2022]['TFR'].idxmin(), '都道府県']})")
print(f"婚姻率×TFR 相関係数 (2022): {df_2022['婚姻率'].corr(df_2022['TFR']):.3f}")
print(f"高齢化率×TFR 相関係数 (全期間): {df_b['高齢化率'].corr(df_b['TFR']):.3f}")
print(f"保育所密度×TFR 相関係数 (全期間): {df_b['保育所密度'].corr(df_b['TFR']):.3f}")

print("\n=== 係数推定値 ===")
for name, val, se, pv in zip(coef_names, coef_vals, coef_se, coef_pvals):
    sig = '***' if pv < 0.01 else ('**' if pv < 0.05 else ('*' if pv < 0.1 else ''))
    print(f"  {name}: {val:+.5f} (SE={se:.5f}, p={pv:.4f}) {sig}")

if use_panel and fe_results is not None:
    print(f"\nFEモデル R²_within: {fe_results.rsquared_within:.4f}")

print("\n✓ 全4図の生成完了")
print(f"  {fig1_path}")
print(f"  {fig2_path}")
print(f"  {fig3_path}")
print(f"  {fig4_path}")