"""
教育用再現コード: 2022年度 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=============================================================================
論文タイトル：女性就労と少子化：女性就業率・育児環境と合計特殊出生率の関係
受　　賞：審査員奨励賞（高校生の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県パネルデータ, 2012〜2023年）
  目的変数：合計特殊出生率（TFR, A4103）

  使用変数:
    TFR              = 合計特殊出生率 (A4103)
    女性就業率_代理  = 15〜64歳人口（女）/ 15〜64歳人口 × 100 (A130202 / A1302)
    保育所密度       = 保育所等数 / 総人口 × 10000 (J2503 / A1101)
    高齢化率         = 65歳以上人口 / 総人口 × 100 (A1303 / A1101)
    婚姻率           = 婚姻件数 / 総人口 × 1000 (A9101 / A1101)
    消費支出_log     = log(消費支出（二人以上の世帯）) (log(L3221))

  分析内容:
    Fig1: TFRの時系列推移（地域別平均, 2012〜2023年）
    Fig2: 女性就業率代理 vs TFR 散布図（2022年, 都道府県ラベル付き）
    Fig3: OLS回帰係数プロット（TFRの決定要因）
    Fig4: 相関ヒートマップ（TFR・女性就業率・保育所密度・婚姻率・高齢化率）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ（2012〜2023年度版）

【データサイエンス学習ポイント】
  1. 女性就業率の逆U字仮説（Nordic paradox — 就業率↑でTFR↑の実証）
  2. 相関行列の読み方（ヒートマップの色と値の解釈）
  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/2022_H5_10_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
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
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026: 都道府県パネルデータ）
# ================================================================
df_raw = pd.read_csv(DATA_B, encoding='cp932', header=0)
# 行0はコード行（英字）、行1はラベル行（日本語）、行2以降がデータ
data = df_raw.iloc[1:].copy()
data.columns = df_raw.columns

# 列名マッピング（コード → 意味）
# SSDSE-B-2026  = 年度
# Code          = 地域コード
# Prefecture    = 都道府県
# A4103         = 合計特殊出生率
# A1302         = 15〜64歳人口（総数）
# A130202       = 15〜64歳人口（女）
# A1101         = 総人口
# A1303         = 65歳以上人口
# A9101         = 婚姻件数
# J2503         = 保育所等数
# L3221         = 消費支出（二人以上の世帯）

num_cols = ['A4103', 'A1302', 'A130202', 'A1101', 'A1303', 'A9101', 'J2503', 'L3221']
for c in num_cols:
    data[c] = pd.to_numeric(data[c], errors='coerce')
data['SSDSE-B-2026'] = pd.to_numeric(data['SSDSE-B-2026'], errors='coerce')

# 年度・都道府県・国全体コードを整理
data = data[data['Code'] != 'Z00000'].copy()   # 全国計除外
data = data.dropna(subset=num_cols + ['SSDSE-B-2026']).copy()
data = data[(data['SSDSE-B-2026'] >= 2012) & (data['SSDSE-B-2026'] <= 2023)].copy()

print("=" * 60)
print(f"■ 読み込み完了: {len(data)}件（{data['SSDSE-B-2026'].nunique()}年 × {data['Prefecture'].nunique()}都道府県）")
print("=" * 60)

# ── 特徴量エンジニアリング ────────────────────────────────────────
data['TFR']             = data['A4103'].astype(float)
data['女性就業率_代理'] = data['A130202'] / data['A1302'] * 100   # 15-64歳女性/15-64歳総人口
data['保育所密度']      = data['J2503'] / data['A1101'] * 10000   # 保育所等数/総人口×10000
data['高齢化率']        = data['A1303'] / data['A1101'] * 100     # 65歳以上/総人口×100
data['婚姻率']          = data['A9101'] / data['A1101'] * 1000    # 婚姻件数/総人口×1000
data['消費支出_log']    = np.log(data['L3221'].clip(lower=1))     # log(消費支出)

data['年度'] = data['SSDSE-B-2026'].astype(int)
data['都道府県'] = data['Prefecture']

vars_stat = ['TFR', '女性就業率_代理', '保育所密度', '高齢化率', '婚姻率', '消費支出_log']

# ================================================================
# ■ 地域区分の定義（北海道・東北 / 関東甲信越 / 東海・北陸 / 近畿 / 中国・四国 / 九州・沖縄）
# ================================================================
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東甲信越', '栃木県': '関東甲信越', '群馬県': '関東甲信越',
    '埼玉県': '関東甲信越', '千葉県': '関東甲信越', '東京都': '関東甲信越',
    '神奈川県': '関東甲信越', '新潟県': '関東甲信越', '山梨県': '関東甲信越',
    '長野県': '関東甲信越',
    '富山県': '東海・北陸', '石川県': '東海・北陸', '福井県': '東海・北陸',
    '静岡県': '東海・北陸', '愛知県': '東海・北陸', '三重県': '東海・北陸',
    '岐阜県': '東海・北陸',
    '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
data['地域'] = data['都道府県'].map(region_map).fillna('その他')

region_order  = ['北海道・東北', '関東甲信越', '東海・北陸', '近畿', '中国・四国', '九州・沖縄']
region_colors = ['#1565C0', '#C62828', '#2E7D32', '#E65100', '#6A1B9A', '#00695C']

# 2022年断面データ（地域列を含む）
df22 = data[data['年度'] == 2022].copy()

print("\n【変数の記述統計（2022年断面）】")
print(df22[vars_stat].describe().round(3))

# ================================================================
# ■ OLS回帰（2022年断面）
# ================================================================
print("\n" + "=" * 60)
print("■ OLS回帰（2022年断面, N=47都道府県）")
print("=" * 60)

df22_clean = df22.dropna(subset=vars_stat).copy()
X_vars = ['女性就業率_代理', '保育所密度', '高齢化率', '婚姻率', '消費支出_log']
X_ols = sm.add_constant(df22_clean[X_vars].astype(float))
y_ols = df22_clean['TFR'].astype(float)

ols_model = sm.OLS(y_ols, X_ols).fit()
print(ols_model.summary2())

# 相関行列（2022年）
print("\n【相関行列（2022年, N=47）】")
corr_vars = ['TFR', '女性就業率_代理', '保育所密度', '婚姻率', '高齢化率']
corr_df = df22_clean[corr_vars].astype(float).corr()
print(corr_df.round(3))

# ================================================================
# ■ 図の生成（4枚）
# ================================================================

# ────────────────────────────────────────────────────────────────
# 図1: TFRの時系列推移（地域別平均, 2012〜2023年）
# ────────────────────────────────────────────────────────────────
print("\n図1を生成中 ...")

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

tfr_region = (
    data[data['地域'].isin(region_order)]
    .groupby(['年度', '地域'])['TFR']
    .mean()
    .reset_index()
)

for reg, col in zip(region_order, region_colors):
    subset = tfr_region[tfr_region['地域'] == reg].sort_values('年度')
    ax1.plot(subset['年度'], subset['TFR'], marker='o', markersize=4,
             linewidth=1.8, color=col, label=reg, alpha=0.9)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('合計特殊出生率の時系列推移（地域別平均）\nSSDSE-B 2012〜2023年 実データ',
              fontsize=13, fontweight='bold')
ax1.legend(fontsize=9, loc='upper right', ncol=2)
ax1.set_xlim(2011.5, 2023.5)
ax1.set_xticks(range(2012, 2024))
ax1.xaxis.set_tick_params(rotation=45)
ax1.grid(True, alpha=0.3)
ax1.axvline(2020, color='gray', linestyle=':', linewidth=1.2, alpha=0.6)
ax1.text(2020.1, ax1.get_ylim()[0] + 0.03, 'COVID-19', fontsize=8, color='gray')

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2022_H5_10_fig1_timeseries.png')
fig1.savefig(fig1_path, bbox_inches='tight', dpi=150)
plt.close(fig1)
print(f"図1保存: {fig1_path}")

# ────────────────────────────────────────────────────────────────
# 図2: 女性就業率代理 vs TFR 散布図（2022年, 都道府県ラベル付き）
# ────────────────────────────────────────────────────────────────
print("図2を生成中 ...")

fig2, ax2 = plt.subplots(figsize=(11, 7))

scatter_colors = [region_colors[region_order.index(reg)] if reg in region_order else '#888888'
                  for reg in df22_clean['地域']]

ax2.scatter(df22_clean['女性就業率_代理'], df22_clean['TFR'],
            c=scatter_colors, s=60, alpha=0.8, edgecolors='white', linewidths=0.5, zorder=3)

# 都道府県ラベル（略称）
label_map = {
    '北海道': '北海', '青森県': '青森', '岩手県': '岩手', '宮城県': '宮城', '秋田県': '秋田',
    '山形県': '山形', '福島県': '福島', '茨城県': '茨城', '栃木県': '栃木', '群馬県': '群馬',
    '埼玉県': '埼玉', '千葉県': '千葉', '東京都': '東京', '神奈川県': '神奈', '新潟県': '新潟',
    '富山県': '富山', '石川県': '石川', '福井県': '福井', '山梨県': '山梨', '長野県': '長野',
    '岐阜県': '岐阜', '静岡県': '静岡', '愛知県': '愛知', '三重県': '三重', '滋賀県': '滋賀',
    '京都府': '京都', '大阪府': '大阪', '兵庫県': '兵庫', '奈良県': '奈良', '和歌山県': '和歌',
    '鳥取県': '鳥取', '島根県': '島根', '岡山県': '岡山', '広島県': '広島', '山口県': '山口',
    '徳島県': '徳島', '香川県': '香川', '愛媛県': '愛媛', '高知県': '高知', '福岡県': '福岡',
    '佐賀県': '佐賀', '長崎県': '長崎', '熊本県': '熊本', '大分県': '大分', '宮崎県': '宮崎',
    '鹿児島県': '鹿児', '沖縄県': '沖縄',
}
for _, row in df22_clean.iterrows():
    short = label_map.get(row['都道府県'], row['都道府県'][:2])
    ax2.annotate(short,
                 xy=(row['女性就業率_代理'], row['TFR']),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, alpha=0.85)

# 回帰直線
x_emp = df22_clean['女性就業率_代理'].astype(float).values
y_tfr = df22_clean['TFR'].astype(float).values
r2, p2 = stats.pearsonr(x_emp, y_tfr)
z2 = np.polyfit(x_emp, y_tfr, 1)
xs2 = np.linspace(x_emp.min(), x_emp.max(), 100)
ax2.plot(xs2, np.poly1d(z2)(xs2), color='#C62828', linewidth=1.5,
         linestyle='--', alpha=0.8, label=f'回帰直線 r={r2:.3f} (p={p2:.3f})')

# 凡例（地域別）
from matplotlib.lines import Line2D
handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor=col,
                  markersize=8, label=reg)
           for reg, col in zip(region_order, region_colors)]
handles.append(Line2D([0], [0], color='#C62828', linestyle='--', linewidth=1.5, label='回帰直線'))
ax2.legend(handles=handles, fontsize=8, loc='upper left')

ax2.set_xlabel('女性就業率代理（15〜64歳女性/15〜64歳総人口 ×100, %）', fontsize=11)
ax2.set_ylabel('合計特殊出生率（TFR）', fontsize=11)
ax2.set_title('女性就業率代理 vs 合計特殊出生率（2022年, 47都道府県）\nSSDSE-B 実データ',
              fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2022_H5_10_fig2_scatter.png')
fig2.savefig(fig2_path, bbox_inches='tight', dpi=150)
plt.close(fig2)
print(f"図2保存: {fig2_path}")

# ────────────────────────────────────────────────────────────────
# 図3: OLS回帰係数プロット（TFRの決定要因）
# ────────────────────────────────────────────────────────────────
print("図3を生成中 ...")

coef_names = ['女性就業率_代理', '保育所密度', '高齢化率', '婚姻率', '消費支出_log']
coef_labels = [
    '女性就業率代理\n(15-64歳女性比率)',
    '保育所密度\n(保育所数/万人)',
    '高齢化率\n(65歳以上/総人口)',
    '婚姻率\n(婚姻件数/千人)',
    '消費支出(log)\n(二人以上世帯)',
]

coefs3 = np.asarray(ols_model.params[1:])   # 定数項を除く
ses3   = np.asarray(ols_model.bse[1:])
pvals3 = np.asarray(ols_model.pvalues[1:])

bar_colors3 = []
for p in pvals3:
    if p < 0.01:
        bar_colors3.append('#C62828')    # 赤: p<0.01
    elif p < 0.05:
        bar_colors3.append('#FF8F00')    # オレンジ: p<0.05
    elif p < 0.1:
        bar_colors3.append('#1565C0')    # 青: p<0.10
    else:
        bar_colors3.append('#9E9E9E')    # 灰: n.s.

fig3, ax3 = plt.subplots(figsize=(9, 5))
y_pos3 = np.arange(len(coef_names))
bars3 = ax3.barh(y_pos3, coefs3, color=bar_colors3, alpha=0.80,
                 edgecolor='white', height=0.6)
ax3.errorbar(coefs3, y_pos3, xerr=1.96 * ses3,
             fmt='none', color='#333333', capsize=4, linewidth=1.2, zorder=5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(coef_labels, fontsize=10)
ax3.set_xlabel('標準化なし回帰係数（±1.96 SE）', fontsize=11)
ax3.set_title(f'OLS回帰係数プロット（目的変数: TFR）\nR²={ols_model.rsquared:.3f}, N=47都道府県（2022年）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#C62828', alpha=0.80, label='p < 0.01'),
    Patch(facecolor='#FF8F00', alpha=0.80, label='p < 0.05'),
    Patch(facecolor='#1565C0', alpha=0.80, label='p < 0.10'),
    Patch(facecolor='#9E9E9E', alpha=0.80, label='n.s.'),
]
ax3.legend(handles=legend_elements, fontsize=9, loc='lower right')

# p値注記
for i, (coef, se, pval) in enumerate(zip(coefs3, ses3, pvals3)):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else '†' if pval < 0.10 else ''
    offset = 0.003 if coef >= 0 else -0.003
    ax3.text(coef + 1.96 * se + offset, i, f' {sig}',
             va='center', ha='left', fontsize=10, color='black')

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2022_H5_10_fig3_coef.png')
fig3.savefig(fig3_path, bbox_inches='tight', dpi=150)
plt.close(fig3)
print(f"図3保存: {fig3_path}")

# ────────────────────────────────────────────────────────────────
# 図4: 相関ヒートマップ（TFR・女性就業率・保育所密度・婚姻率・高齢化率）
# ────────────────────────────────────────────────────────────────
print("図4を生成中 ...")

heat_vars  = ['TFR', '女性就業率_代理', '保育所密度', '婚姻率', '高齢化率']
heat_labels = [
    'TFR\n(合計特殊出生率)',
    '女性就業率代理\n(15-64歳)',
    '保育所密度\n(保育所/万人)',
    '婚姻率\n(件/千人)',
    '高齢化率\n(65歳以上%)',
]

corr_mat = df22_clean[heat_vars].astype(float).corr()

fig4, ax4 = plt.subplots(figsize=(8, 6.5))
im = ax4.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')

cbar = plt.colorbar(im, ax=ax4, fraction=0.046, pad=0.04)
cbar.set_label('相関係数', fontsize=10)

n_vars4 = len(heat_vars)
ax4.set_xticks(range(n_vars4))
ax4.set_yticks(range(n_vars4))
ax4.set_xticklabels(heat_labels, fontsize=8.5, rotation=0, ha='center')
ax4.set_yticklabels(heat_labels, fontsize=8.5)

# 相関係数と有意星を表示
for i in range(n_vars4):
    for j in range(n_vars4):
        val = corr_mat.iloc[i, j]
        if i == j:
            ax4.text(j, i, '1.00', ha='center', va='center',
                     fontsize=10, fontweight='bold', color='black')
        else:
            # p値計算
            x_c = df22_clean[heat_vars[j]].astype(float).values
            y_c = df22_clean[heat_vars[i]].astype(float).values
            mask_c = np.isfinite(x_c) & np.isfinite(y_c)
            _, pval_c = stats.pearsonr(x_c[mask_c], y_c[mask_c])
            sig_c = '***' if pval_c < 0.001 else '**' if pval_c < 0.01 else '*' if pval_c < 0.05 else ''
            text_color = 'white' if abs(val) > 0.6 else 'black'
            ax4.text(j, i, f'{val:.2f}{sig_c}',
                     ha='center', va='center', fontsize=9,
                     color=text_color, fontweight='bold' if abs(val) > 0.3 else 'normal')

ax4.set_title('相関ヒートマップ（2022年, N=47都道府県）\n* p<0.05, ** p<0.01, *** p<0.001',
              fontsize=12, fontweight='bold', pad=14)

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2022_H5_10_fig4_heatmap.png')
fig4.savefig(fig4_path, bbox_inches='tight', dpi=150)
plt.close(fig4)
print(f"図4保存: {fig4_path}")

# ================================================================
# ■ 分析結果サマリー
# ================================================================
print("\n" + "=" * 60)
print("■ 分析結果サマリー")
print("=" * 60)
print(f"  使用データ: SSDSE-B-2026（実データ）")
print(f"  対象期間: 2012〜2023年, 47都道府県パネル")
print()
print("  OLS回帰結果（2022年断面, N=47）:")
print(f"  {'変数':<20} {'係数':>8} {'SE':>8} {'p値':>10} {'有意'}")
print("  " + "-" * 55)
for name, coef, se, pval in zip(X_vars, coefs3, ses3, pvals3):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else '†' if pval < 0.10 else 'n.s.'
    print(f"  {name:<20} {coef:>8.4f} {se:>8.4f} {pval:>10.4f}  {sig}")
print(f"  R² = {ols_model.rsquared:.4f}, N=47")
print()
print("  相関分析（2022年, TFRとの相関）:")
for var in heat_vars[1:]:
    x_v = df22_clean[var].astype(float).values
    y_v = df22_clean['TFR'].astype(float).values
    mask_v = np.isfinite(x_v) & np.isfinite(y_v)
    r_v, p_v = stats.pearsonr(x_v[mask_v], y_v[mask_v])
    sig_v = '***' if p_v < 0.001 else '**' if p_v < 0.01 else '*' if p_v < 0.05 else 'n.s.'
    print(f"  TFR ↔ {var:<18}: r={r_v:+.3f} ({sig_v})")

print("\n全図の生成完了（4枚）")
print(f"  fig1: 2022_H5_10_fig1_timeseries.png  (TFR時系列推移)")
print(f"  fig2: 2022_H5_10_fig2_scatter.png     (女性就業率 vs TFR散布図)")
print(f"  fig3: 2022_H5_10_fig3_coef.png        (OLS回帰係数プロット)")
print(f"  fig4: 2022_H5_10_fig4_heatmap.png     (相関ヒートマップ)")