都道府県別のパネルデータを用いた
合計特殊出生率の決定要因

print("=" * 65)
print("■ データ読み込み・構築（SSDSE-B-2026.csv）")
print("=" * 65)

# SSDSE-B-2026 読み込み
# header=0 でコードベースの列名（A4103 等）を取得し、
# 先頭行（日本語ラベル行）をスキップする
df_b = pd.read_csv(DATA_B, encoding='cp932', header=0)
df_b = df_b.iloc[1:].reset_index(drop=True)  # 日本語ラベル行を除外

# 列名を分析で使いやすい名称にリネーム
df_b = df_b.rename(columns={
    'SSDSE-B-2026': '年度',
    'Code':         '地域コード',
    'Prefecture':   '都道府県',
})

# 47都道府県のみ（地域コード = R + 5桁数字）
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 必要列を数値変換
needed = ['年度', '都道府県', 'A4103', 'A9101', 'J2503', 'A1101',
          'J2506', 'J2505', 'F3103', 'F3102', 'L3221', 'C5401', 'A1303']
for col in needed:
    if col not in ('都道府県',):
        df_b[col] = pd.to_numeric(df_b[col], errors='coerce')

df_b = df_b.dropna(subset=needed)

print(f"読み込み完了: {len(df_b)} 行（47都道府県 × 12年度）")
print(f"年度範囲: {int(df_b['年度'].min())}〜{int(df_b['年度'].max())}")
print(f"都道府県数: {df_b['都道府県'].nunique()}")

# ── 説明変数の計算（実データのみ） ────────────────────────────────
# 婚姻率: 婚姻件数 / 総人口 × 10000
df_b['婚姻率']    = df_b['A9101'] / df_b['A1101'] * 10000
# 保育所密度: 保育所等数 / 総人口 × 10000
df_b['保育所密度'] = df_b['J2503'] / df_b['A1101'] * 10000
# 保育充足率: 在所児数 / 定員数
df_b['保育充足率'] = df_b['J2506'] / df_b['J2505']
# 有効求人倍率: 有効求人数 / 有効求職者数
df_b['有効求人倍率'] = df_b['F3103'] / df_b['F3102']
# 高齢化率: 65歳以上 / 総人口 × 100
df_b['高齢化率']   = df_b['A1303'] / df_b['A1101'] * 100

PREDICTOR_COLS = ['婚姻率', '保育所密度', '保育充足率', '有効求人倍率',
                  'L3221', 'C5401', '高齢化率']

PRED_LABELS = {
    '婚姻率':      '婚姻率（件/万人）',
    '保育所密度':  '保育所密度（施設/万人）',
    '保育充足率':  '保育所充足率（在所/定員）',
    '有効求人倍率': '有効求人倍率',
    'L3221':       '消費支出（円）',
    'C5401':       '住宅地価格（円/㎡）',
    '高齢化率':    '高齢化率（%）',
}

# ── パネルデータ構築 ──────────────────────────────────────────────
df_panel = df_b.set_index(['都道府県', '年度'])
all_needed = ['A4103'] + PREDICTOR_COLS
df_panel = df_panel.dropna(subset=all_needed)

print(f"\nパネルデータ: {len(df_panel)} 観測, "
      f"{df_panel.index.get_level_values(0).nunique()} 都道府県, "
      f"{df_panel.index.get_level_values(1).nunique()} 年度")

print("図1: 全国平均TFRの時系列推移を作成中...")

tfr_year = df_b.groupby('年度')['A4103'].mean()
tfr_q25  = df_b.groupby('年度')['A4103'].quantile(0.25)
tfr_q75  = df_b.groupby('年度')['A4103'].quantile(0.75)

fig1, ax1 = plt.subplots(figsize=(11, 5.5))

ax1.fill_between(tfr_year.index, tfr_q25.values, tfr_q75.values,
                 alpha=0.20, color='#1565C0', label='25〜75パーセンタイル')
ax1.plot(tfr_year.index, tfr_year.values, 'o-',
         color='#1565C0', linewidth=2.2, markersize=7, label='全国平均TFR')
ax1.axhline(2.07, color='#C62828', linestyle=':', linewidth=1.5, alpha=0.7,
            label='人口置換水準（2.07）')
ax1.axhline(1.0, color='#999', linestyle='--', linewidth=1.0, alpha=0.5)

# 主要イベント注記
events = {
    2015: ('少子化社会\n対策大綱', '#2E7D32'),
    2019: ('新型コロナ\n前夜', '#E65100'),
    2020: ('コロナ禍\n始まり', '#C62828'),
    2023: ('過去最低\n1.20→1.09', '#6A1B9A'),
}
for yr, (label, color) in events.items():
    if yr in tfr_year.index:
        val = tfr_year[yr]
        ax1.annotate(label,
                     xy=(yr, val), xytext=(yr, val + 0.06),
                     ha='center', fontsize=9, color=color,
                     arrowprops=dict(arrowstyle='->', color=color, lw=1.2),
                     bbox=dict(boxstyle='round,pad=0.3', fc='white', ec=color, alpha=0.85))

ax1.set_xlabel('年度', fontsize=11)
ax1.set_ylabel('合計特殊出生率（TFR）', fontsize=11)
ax1.set_title('全国平均 合計特殊出生率の推移（2012〜2023年）\n'
              '出典: SSDSE-B-2026（47都道府県）',
              fontsize=12, fontweight='bold')
ax1.set_ylim(1.1, 1.75)
ax1.legend(fontsize=9, loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xticks(sorted(tfr_year.index))
ax1.set_xticklabels([str(y) for y in sorted(tfr_year.index)],
                    rotation=45, ha='right', fontsize=9)

plt.tight_layout()
out1 = os.path.join(FIG_DIR, '2023_U5_4_fig1_tfr_trend.png')
fig1.savefig(out1, bbox_inches='tight', dpi=150)
plt.close(fig1)
print(f"  -> {os.path.basename(out1)} 保存完了")

print("図2: 固定効果モデルの係数プロットを作成中...")

fig2, ax2 = plt.subplots(figsize=(10, 6))

coefs  = [fe.params[c]      for c in PREDICTOR_COLS]
ses    = [fe.std_errors[c]  for c in PREDICTOR_COLS]
pvals  = [fe.pvalues[c]     for c in PREDICTOR_COLS]
labels = [PRED_LABELS[c]    for c in PREDICTOR_COLS]

bar_colors = [
    '#C62828' if c < 0 and p < 0.05 else
    '#1565C0' if c > 0 and p < 0.05 else
    '#9E9E9E'
    for c, p in zip(coefs, pvals)
]

y_pos = np.arange(len(PREDICTOR_COLS))
ci_95 = [1.96 * s for s in ses]

ax2.barh(y_pos, coefs, xerr=ci_95,
         color=bar_colors, alpha=0.82, edgecolor='white',
         capsize=4, error_kw={'elinewidth': 1.5, 'ecolor': '#444'})
ax2.axvline(0, color='black', linewidth=1.0)
ax2.set_yticks(y_pos)
ax2.set_yticklabels(labels, fontsize=10)
ax2.set_xlabel('固定効果推定量（±1.96 × SE）', fontsize=11)
ax2.set_title('パネル固定効果モデルの回帰係数（95% CI）\n'
              '（被説明変数: 合計特殊出生率, SSDSE-B実データ）',
              fontsize=11, fontweight='bold')
ax2.grid(axis='x', alpha=0.3)
ax2.invert_yaxis()

for i, (c, s, p) in enumerate(zip(coefs, ses, pvals)):
    sig = ('***' if p < 0.001 else '**' if p < 0.01 else
           '*'   if p < 0.05  else '')
    if sig:
        offset = np.sign(c) * (1.96 * s + abs(max(coefs, key=abs)) * 0.03)
        ax2.text(c + offset, i, sig,
                 va='center', ha='left' if c > 0 else 'right',
                 fontsize=11, color='#222', fontweight='bold')

red_p  = mpatches.Patch(color='#C62828', alpha=0.82, label='有意・負の効果（p<0.05）')
blu_p  = mpatches.Patch(color='#1565C0', alpha=0.82, label='有意・正の効果（p<0.05）')
gry_p  = mpatches.Patch(color='#9E9E9E', alpha=0.82, label='非有意（p≥0.05）')
ax2.legend(handles=[blu_p, red_p, gry_p], fontsize=9, loc='lower right')

# Hausman・R2情報をテキストで追記
ax2.text(0.98, 0.02,
         f"Within R²={fe.rsquared:.3f}\nHausman H={H_stat:.1f} (p<0.001)\n→ FE採用",
         transform=ax2.transAxes, ha='right', va='bottom',
         fontsize=9, color='#333',
         bbox=dict(boxstyle='round,pad=0.4', fc='#F3E5F5', ec='#6A1B9A', alpha=0.9))

plt.tight_layout()
out2 = os.path.join(FIG_DIR, '2023_U5_4_fig2_fe_coef.png')
fig2.savefig(out2, bbox_inches='tight', dpi=150)
plt.close(fig2)
print(f"  -> {os.path.basename(out2)} 保存完了")

print("図3: 婚姻率 vs TFR の散布図（2022年横断面）を作成中...")

df_2022 = df_b[df_b['年度'] == 2022].copy()

fig3, ax3 = plt.subplots(figsize=(9, 6.5))

sc3 = ax3.scatter(df_2022['婚姻率'], df_2022['A4103'],
                  c=df_2022['高齢化率'], cmap='RdYlGn_r',
                  s=70, alpha=0.85, edgecolors='white', linewidths=0.5,
                  zorder=3)
cbar3 = plt.colorbar(sc3, ax=ax3)
cbar3.set_label('高齢化率（%）', fontsize=10)

# 回帰直線
x3 = df_2022['婚姻率'].values
y3 = df_2022['A4103'].values
mask3 = ~(np.isnan(x3) | np.isnan(y3))
if mask3.sum() > 2:
    slope3, intercept3, r3, p3, _ = stats.linregress(x3[mask3], y3[mask3])
    xline3 = np.linspace(x3[mask3].min(), x3[mask3].max(), 100)
    ax3.plot(xline3, intercept3 + slope3 * xline3,
             '--', color='#1565C0', linewidth=1.8, label=f'OLS (r={r3:.2f})', zorder=2)

# 代表的な都道府県をラベル表示
highlights3 = ['東京都', '沖縄県', '秋田県', '島根県', '愛知県']
for _, row in df_2022.iterrows():
    if row['都道府県'] in highlights3:
        ax3.annotate(row['都道府県'],
                     xy=(row['婚姻率'], row['A4103']),
                     xytext=(4, 4), textcoords='offset points',
                     fontsize=8.5, color='#333',
                     bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.75))

ax3.set_xlabel('婚姻率（件/万人）', fontsize=11)
ax3.set_ylabel('合計特殊出生率（TFR）', fontsize=11)
ax3.set_title('婚姻率 vs 合計特殊出生率（2022年 都道府県横断面）\n'
              '色: 高齢化率（緑=低、赤=高）',
              fontsize=11, fontweight='bold')
ax3.legend(fontsize=10, loc='upper left')
ax3.grid(True, alpha=0.3)

plt.tight_layout()
out3 = os.path.join(FIG_DIR, '2023_U5_4_fig3_marriage_tfr.png')
fig3.savefig(out3, bbox_inches='tight', dpi=150)
plt.close(fig3)
print(f"  -> {os.path.basename(out3)} 保存完了")

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 statsmodels.stats.outliers_influence import variance_inflation_factor
from linearmodels import PanelOLS, RandomEffects
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

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

import os
DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

DATA_B = os.path.join(DATA_DIR, 'SSDSE-B-2026.csv')

print("\n" + "=" * 65)
print("■ Step1. パネル固定効果・変量効果モデルの推定")
print("=" * 65)

dep  = df_panel['A4103']
exog = sm.add_constant(df_panel[PREDICTOR_COLS])

# 固定効果モデル（entity effects + clustered SE）
fe = PanelOLS(dep, exog, entity_effects=True).fit(
    cov_type='clustered', cluster_entity=True)

# 変量効果モデル
re = RandomEffects(dep, exog).fit()

print("\n【固定効果モデル（FE）】")
print(f"  Within R²   = {fe.rsquared:.4f}")
print(f"  観測数       = {fe.nobs}")
print(f"  F統計量      = {fe.f_statistic.stat:.3f}  (p={fe.f_statistic.pval:.4f})")
print()
print(f"  {'変数':<16} {'係数':>10} {'SE':>10} {'t値':>8} {'p値':>10}  有意")
print("  " + "-" * 60)
for col in PREDICTOR_COLS + ['const']:
    if col in fe.params.index:
        b  = fe.params[col]
        se = fe.std_errors[col]
        t  = fe.tstats[col]
        p  = fe.pvalues[col]
        sig = ('***' if p < 0.001 else
               '**'  if p < 0.01  else
               '*'   if p < 0.05  else
               '†'   if p < 0.1   else '')
        label = PRED_LABELS.get(col, col)
        print(f"  {label:<16} {b:>10.5f} {se:>10.5f} {t:>8.3f} {p:>10.4f}  {sig}")

print("\n" + "=" * 65)
print("■ Step2. Hausman検定（固定効果 vs 変量効果）")
print("=" * 65)

b_fe = fe.params
b_re = re.params
common = [c for c in b_fe.index if c in b_re.index and c != 'const']

diff    = np.array([b_fe[c] - b_re[c] for c in common])
var_fe  = np.array([fe.cov.loc[c, c]  for c in common])
var_re  = np.array([re.cov.loc[c, c]  for c in common])
var_diff = var_fe - var_re

# 分散差が正の変数のみ使用（標準的なHausman実装）
valid    = var_diff > 0
H_stat   = float(diff[valid] @ np.diag(1.0 / var_diff[valid]) @ diff[valid])
H_df     = int(valid.sum())
H_pval   = 1 - stats.chi2.cdf(H_stat, df=H_df)

print(f"\n  Hausman検定統計量 H = {H_stat:.3f}")
print(f"  自由度              = {H_df}")
print(f"  p値                 = {H_pval:.4f}")
if H_pval < 0.05:
    print("  → 帰無仮説（変量効果が一致推定量）を棄却")
    print("  → 固定効果モデル（FE）を採用")
else:
    print("  → 帰無仮説を棄却できない → 変量効果モデル（RE）を採用")

print("\n" + "=" * 65)
print("■ Step3. VIF（分散膨張係数）による多重共線性診断")
print("=" * 65)

X_vif = df_b[PREDICTOR_COLS].dropna().values
vif_values = {}
print(f"\n  {'変数':<16} {'VIF':>8}  判定")
print("  " + "-" * 40)
for i, col in enumerate(PREDICTOR_COLS):
    vif = variance_inflation_factor(X_vif, i)
    vif_values[col] = vif
    flag = ('⚠ 高' if vif > 10 else '○ 良好')
    label = PRED_LABELS.get(col, col)
    print(f"  {label:<16} {vif:>8.2f}  {flag}")

print("\n  ※ VIF > 10 は多重共線性の懸念（パネルFEでは都道府県固定効果が吸収するため参考値）")

print("\n" + "=" * 65)
print("■ 図の生成（4枚）")
print("=" * 65)

print("図4: 保育所密度 vs TFR の散布図（複数年）を作成中...")

TARGET_YEARS = [2016, 2019, 2022]
YEAR_COLORS  = {2016: '#1565C0', 2019: '#2E7D32', 2022: '#C62828'}
YEAR_MARKERS = {2016: 'o', 2019: 's', 2022: '^'}

fig4, ax4 = plt.subplots(figsize=(9, 6.5))

for yr in TARGET_YEARS:
    sub = df_b[df_b['年度'] == yr].copy()
    x4 = sub['保育所密度'].values
    y4 = sub['A4103'].values
    mask4 = ~(np.isnan(x4) | np.isnan(y4))
    c4 = YEAR_COLORS[yr]
    m4 = YEAR_MARKERS[yr]

    ax4.scatter(x4[mask4], y4[mask4],
                color=c4, marker=m4, s=60, alpha=0.72,
                edgecolors='white', linewidths=0.4,
                label=f'{yr}年', zorder=3)

    if mask4.sum() > 2:
        slope4, intercept4, r4, p4, _ = stats.linregress(x4[mask4], y4[mask4])
        xline4 = np.linspace(x4[mask4].min(), x4[mask4].max(), 80)
        ax4.plot(xline4, intercept4 + slope4 * xline4,
                 '--', color=c4, linewidth=1.4,
                 label=f'{yr}年 OLS (r={r4:.2f})', alpha=0.85, zorder=2)

ax4.set_xlabel('保育所密度（施設/万人）', fontsize=11)
ax4.set_ylabel('合計特殊出生率（TFR）', fontsize=11)
ax4.set_title('保育所密度 vs 合計特殊出生率（2016・2019・2022年）\n'
              '47都道府県別 横断面散布図',
              fontsize=11, fontweight='bold')
ax4.legend(fontsize=9, loc='upper right', ncol=2)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
out4 = os.path.join(FIG_DIR, '2023_U5_4_fig4_childcare_tfr.png')
fig4.savefig(out4, bbox_inches='tight', dpi=150)
plt.close(fig4)
print(f"  -> {os.path.basename(out4)} 保存完了")

print("\n" + "=" * 65)
print("■ 分析結果サマリ")
print("=" * 65)
print(f"  パネルデータ: 47都道府県 × 12年度 = {fe.nobs} 観測")
print(f"  固定効果モデル Within R² = {fe.rsquared:.4f}")
print(f"  Hausman検定: H={H_stat:.3f}, df={H_df}, p={H_pval:.4f} → FE採用")
print()
print("  主要な決定要因（FE推定）:")
for col in PREDICTOR_COLS:
    b = fe.params[col]
    p = fe.pvalues[col]
    sig = ('***' if p < 0.001 else '**' if p < 0.01 else
           '*'   if p < 0.05  else '†'  if p < 0.1  else ' ns')
    direction = '正' if b > 0 else '負'
    print(f"    {PRED_LABELS[col]:<20}: {b:>10.5f} ({direction}, {sig})")

print()
print("  VIF診断:")
for col, vif in vif_values.items():
    flag = '⚠ 高' if vif > 10 else '○ 良好'
    print(f"    {PRED_LABELS[col]:<20}: VIF={vif:.1f}  {flag}")

print()
print("生成図:")
print(f"  {os.path.basename(out1)}")
print(f"  {os.path.basename(out2)}")
print(f"  {os.path.basename(out3)}")
print(f"  {os.path.basename(out4)}")
print()
print("※ 使用データ: SSDSE-B-2026.csv のみ（合成データなし）")