合計特殊出生率の決定要因：
パネルデータによる長期分析

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=== 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}")

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}")

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}")

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}")

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)

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}")