"""
教育用再現コード: 2022年度 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================================
論文タイトル：子育て支援と合計特殊出生率：保育所整備の効果検証

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県データ、2012〜2022年）
  目的変数：合計特殊出生率（TFR）
  分析手法：OLS重回帰、時系列分析、相関分析、散布図

【変数設計】
  - TFR              : 合計特殊出生率（A4103）
  - 保育所密度        : 保育所等数 / 総人口 × 10000（施設数/万人）
  - 保育士密度        : 保育所等保育士数 / 総人口 × 10000（人/万人）
  - 待機児童率        : 保育所等利用待機児童数 / 保育所等定員数 × 100（%）
  - 女性就業代理      : 15〜64歳人口（女）/ 15〜64歳人口 × 100（%）
  - 高齢化率          : 65歳以上人口 / 総人口 × 100（%）

【SSDSE-B使用列】
  A1101: 総人口
  A1302: 15〜64歳人口, A130202: 15〜64歳人口（女）
  A1303: 65歳以上人口
  A4103: 合計特殊出生率
  J2503: 保育所等数
  J2505: 保育所等定員数
  J250502: 保育所等利用待機児童数
  J2506: 保育所等在所児数
  J2526: 保育所等保育士数

【図の構成（4枚）】
  Fig1: TFRと保育所密度の時系列推移（全国平均、2軸グラフ）
  Fig2: 保育所密度 vs TFR 散布図（2022年、都道府県ラベル付き）
  Fig3: OLS回帰係数プロット（TFRの決定要因）
  Fig4: 待機児童率の都道府県別ランキング（2022年）

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

【データサイエンス学習ポイント】
  1. 合計特殊出生率（TFR）の定義と人口置換水準（2.07）
  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_7_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)
# 1行目がコード名（英数字）、2行目（index=0）が日本語ラベル → スキップ
df_all = df_raw.iloc[1:].copy()
df_all.columns = df_raw.columns  # 英数字コード名を保持

# 数値変換対象列
NUM_COLS = [
    'A1101',   # 総人口
    'A1302',   # 15〜64歳人口
    'A130202', # 15〜64歳人口（女）
    'A1303',   # 65歳以上人口
    'A4103',   # 合計特殊出生率
    'J2503',   # 保育所等数
    'J2505',   # 保育所等定員数
    'J250502', # 保育所等利用待機児童数
    'J2506',   # 保育所等在所児数
    'J2526',   # 保育所等保育士数
]
for c in NUM_COLS:
    df_all[c] = pd.to_numeric(df_all[c], errors='coerce')

# 年度列・都道府県列の列名
COL_YEAR = 'SSDSE-B-2026'
COL_PREF = 'Prefecture'

print("=" * 60)
print("■ SSDSE-B-2026 読み込み完了")
print(f"  利用可能年度: {sorted(df_all[COL_YEAR].unique())}")
print(f"  都道府県数（2022年）: {len(df_all[df_all[COL_YEAR]=='2022'])}")
print("=" * 60)

# ================================================================
# ■ 時系列データ構築（2012〜2022年、全国平均）
# ================================================================
# 2023年は保育所等数の集計方法変更があるため2012〜2022を使用
YEARS = list(range(2012, 2023))

ts_records = []
for yr in YEARS:
    ydf = df_all[df_all[COL_YEAR] == str(yr)].copy()
    pop_total  = ydf['A1101'].sum()
    nurseries  = ydf['J2503'].sum()
    nursery_staff = ydf['J2526'].sum()
    waiting    = ydf['J250502'].sum()
    capacity   = ydf['J2505'].sum()

    avg_tfr    = ydf['A4103'].mean()
    density_nursery = nurseries / pop_total * 10000      # 保育所密度
    density_staff   = nursery_staff / pop_total * 10000  # 保育士密度
    waiting_rate    = max(waiting / capacity * 100, 0)   # 待機児童率（≥0）

    ts_records.append({
        'year': yr,
        'TFR': avg_tfr,
        '保育所密度': density_nursery,
        '保育士密度': density_staff,
        '待機児童率': waiting_rate,
    })

df_ts = pd.DataFrame(ts_records)
print("\n■ 時系列データ（全国平均）:")
print(df_ts.to_string(index=False))

# ================================================================
# ■ 2022年 都道府県別クロスセクションデータ
# ================================================================
df22 = df_all[df_all[COL_YEAR] == '2022'].copy().reset_index(drop=True)

pop22   = df22['A1101'].clip(lower=1)
w1564   = df22['A130202']    # 15〜64歳人口（女）
pop1564 = df22['A1302'].clip(lower=1)
pop65p  = df22['A1303']

df22['保育所密度']  = df22['J2503'] / pop22 * 10000
df22['保育士密度']  = df22['J2526'] / pop22 * 10000
df22['待機児童率']  = (df22['J250502'] / df22['J2505'] * 100).clip(lower=0)
df22['女性就業代理'] = w1564 / pop1564 * 100
df22['高齢化率']    = pop65p / pop22 * 100
df22['TFR']        = df22['A4103']

# 欠損値除外
FEATURE_COLS = ['保育所密度', '保育士密度', '待機児童率', '女性就業代理', '高齢化率', 'TFR']
df22_clean = df22.dropna(subset=FEATURE_COLS).copy().reset_index(drop=True)

print(f"\n■ 2022年 都道府県データ（欠損除外後）: N={len(df22_clean)}")
print(df22_clean[['Prefecture'] + FEATURE_COLS].to_string(index=False))

# ================================================================
# ■ OLS重回帰: TFR の決定要因分析（2022年）
# ================================================================
print("\n" + "=" * 60)
print("■ OLS重回帰: TFR の決定要因分析（2022年）")
print("=" * 60)

X_vars = ['保育所密度', '保育士密度', '待機児童率', '女性就業代理', '高齢化率']
y_reg  = df22_clean['TFR'].values
X_reg  = df22_clean[X_vars].values
X_reg_c = sm.add_constant(X_reg)

ols_model = sm.OLS(y_reg, X_reg_c).fit()
print(ols_model.summary())

# 標準化係数（効果量の比較）
X_std = (X_reg - X_reg.mean(axis=0)) / X_reg.std(axis=0)
y_std = (y_reg - y_reg.mean()) / y_reg.std()
ols_std = sm.OLS(y_std, sm.add_constant(X_std)).fit()
beta_std = ols_std.params[1:]  # 定数項除く

print("\n■ 標準化偏回帰係数（β）:")
for name, b in zip(X_vars, beta_std):
    print(f"  {name:<12}: β = {b:+.4f}")

# ================================================================
# ■ 相関分析
# ================================================================
print("\n" + "=" * 60)
print("■ 相関分析（各変数 vs TFR）")
print("=" * 60)
print(f"  {'変数':<12} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 42)
for vname in X_vars:
    col = df22_clean[vname].values
    r, p = stats.pearsonr(col, y_reg)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {vname:<12} {r:>8.4f} {p:>10.4f} {sig:>6}")

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

# ────────────────────────────────────────────────────────────────
# 図1: TFRと保育所密度の時系列推移（2軸グラフ）
# ────────────────────────────────────────────────────────────────
fig1, ax1a = plt.subplots(figsize=(10, 5))
ax1b = ax1a.twinx()

years_arr = df_ts['year'].values
tfr_arr   = df_ts['TFR'].values
den_arr   = df_ts['保育所密度'].values

line1, = ax1a.plot(years_arr, tfr_arr, color='#C62828', marker='o', linewidth=2.2,
                   markersize=6, label='合計特殊出生率（左軸）', zorder=3)
ax1a.axhline(2.07, color='#C62828', linestyle='--', linewidth=1.2, alpha=0.6,
             label='人口置換水準 (2.07)')
ax1a.set_ylabel('合計特殊出生率（TFR）', fontsize=11, color='#C62828')
ax1a.tick_params(axis='y', labelcolor='#C62828')
ax1a.set_ylim(1.1, 2.2)

line2, = ax1b.plot(years_arr, den_arr, color='#1565C0', marker='s', linewidth=2.2,
                   markersize=6, label='保育所密度（右軸）', zorder=3)
ax1b.set_ylabel('保育所密度（施設数/万人）', fontsize=11, color='#1565C0')
ax1b.tick_params(axis='y', labelcolor='#1565C0')
ax1b.set_ylim(1.4, 2.8)

ax1a.set_xlabel('年', fontsize=11)
ax1a.set_xticks(years_arr)
ax1a.set_xticklabels([str(y) for y in years_arr], rotation=45, fontsize=9)
ax1a.grid(True, alpha=0.3, zorder=1)

lines_all = [line1, line2,
             plt.Line2D([0], [0], color='#C62828', linestyle='--', linewidth=1.2, alpha=0.6,
                        label='人口置換水準 (2.07)')]
ax1a.legend(lines_all + [line2],
            ['合計特殊出生率（左軸）', '人口置換水準 (2.07)', '保育所密度（右軸）'],
            fontsize=9, loc='upper right')

ax1a.set_title('TFRと保育所密度の時系列推移（全国47都道府県平均）\n2012〜2022年（SSDSE-B 実データ）',
               fontsize=12, fontweight='bold')
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_H5_7_fig1_timeseries.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2022_H5_7_fig1_timeseries.png")

# ────────────────────────────────────────────────────────────────
# 図2: 保育所密度 vs TFR 散布図（2022年、都道府県ラベル付き）
# ────────────────────────────────────────────────────────────────
fig2, ax2 = plt.subplots(figsize=(11, 8))

x2 = df22_clean['保育所密度'].values
y2 = df22_clean['TFR'].values
prefs = df22_clean['Prefecture'].values

ax2.scatter(x2, y2, color='#1565C0', s=60, alpha=0.7, zorder=3)

# 都道府県ラベル（略称表示）
pref_short = {
    '北海道': '北海', '青森県': '青森', '岩手県': '岩手', '宮城県': '宮城',
    '秋田県': '秋田', '山形県': '山形', '福島県': '福島', '茨城県': '茨城',
    '栃木県': '栃木', '群馬県': '群馬', '埼玉県': '埼玉', '千葉県': '千葉',
    '東京都': '東京', '神奈川県': '神奈', '新潟県': '新潟', '富山県': '富山',
    '石川県': '石川', '福井県': '福井', '山梨県': '山梨', '長野県': '長野',
    '岐阜県': '岐阜', '静岡県': '静岡', '愛知県': '愛知', '三重県': '三重',
    '滋賀県': '滋賀', '京都府': '京都', '大阪府': '大阪', '兵庫県': '兵庫',
    '奈良県': '奈良', '和歌山県': '和歌', '鳥取県': '鳥取', '島根県': '島根',
    '岡山県': '岡山', '広島県': '広島', '山口県': '山口', '徳島県': '徳島',
    '香川県': '香川', '愛媛県': '愛媛', '高知県': '高知', '福岡県': '福岡',
    '佐賀県': '佐賀', '長崎県': '長崎', '熊本県': '熊本', '大分県': '大分',
    '宮崎県': '宮崎', '鹿児島県': '鹿児', '沖縄県': '沖縄',
}
for xi, yi, pn in zip(x2, y2, prefs):
    label = pref_short.get(pn, pn[:2])
    ax2.annotate(label, (xi, yi), textcoords='offset points', xytext=(4, 3),
                 fontsize=7.5, color='#333')

# 回帰直線
z2 = np.polyfit(x2, y2, 1)
xs2 = np.linspace(x2.min(), x2.max(), 100)
ax2.plot(xs2, np.poly1d(z2)(xs2), color='#C62828', linewidth=1.8, alpha=0.7,
         label='回帰直線')

r2, p2 = stats.pearsonr(x2, y2)
ax2.set_xlabel('保育所密度（施設数/万人）', fontsize=12)
ax2.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax2.set_title(f'保育所密度 vs 合計特殊出生率（2022年、47都道府県）\nr = {r2:.3f}  p = {p2:.4f}',
              fontsize=13, fontweight='bold')
ax2.axhline(2.07, color='#888', linestyle=':', linewidth=1.0, label='人口置換水準 (2.07)')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

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

# ────────────────────────────────────────────────────────────────
# 図3: OLS回帰係数プロット（TFRの決定要因）
# ────────────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(9, 5))

coefs3  = np.asarray(ols_model.params[1:])
ses3    = np.asarray(ols_model.bse[1:])
pvals3  = np.asarray(ols_model.pvalues[1:])
y_pos3  = np.arange(len(X_vars))
colors3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E' for p in pvals3]

ax3.barh(y_pos3, coefs3, color=colors3, alpha=0.75, edgecolor='white', height=0.6)
ax3.errorbar(coefs3, y_pos3, xerr=1.96 * ses3, fmt='none', color='#333',
             capsize=5, linewidth=1.5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(X_vars, fontsize=11)
ax3.set_xlabel('OLS回帰係数（±1.96 SE）', fontsize=11)
ax3.set_title(f'TFRの決定要因（OLS重回帰係数プロット）\n'
              f'R² = {ols_model.rsquared:.3f}  N = {len(df22_clean)}都道府県（2022年）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

# 有意水準の凡例
from matplotlib.patches import Patch
legend_patches = [
    Patch(facecolor='#C62828', alpha=0.75, label='p < 0.01'),
    Patch(facecolor='#FF8F00', alpha=0.75, label='p < 0.05'),
    Patch(facecolor='#9E9E9E', alpha=0.75, label='n.s.'),
]
ax3.legend(handles=legend_patches, fontsize=9, loc='lower right')

# 係数と有意水準を注釈
for i, (c, p) in enumerate(zip(coefs3, pvals3)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    offset = 0.003 if c >= 0 else -0.003
    ha = 'left' if c >= 0 else 'right'
    ax3.text(c + offset, i, f' {c:.4f}{sig}', va='center', ha=ha, fontsize=8.5)

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

# ────────────────────────────────────────────────────────────────
# 図4: 待機児童率の都道府県別ランキング（2022年）
# ────────────────────────────────────────────────────────────────
fig4, ax4 = plt.subplots(figsize=(11, 10))

df_rank = df22_clean[['Prefecture', '待機児童率', 'TFR']].copy()
df_rank = df_rank.sort_values('待機児童率', ascending=True).reset_index(drop=True)

y_pos4  = np.arange(len(df_rank))
colors4 = ['#C62828' if v > 0.5 else '#FF8F00' if v > 0.1 else '#1565C0'
           for v in df_rank['待機児童率']]

ax4.barh(y_pos4, df_rank['待機児童率'], color=colors4, alpha=0.75,
         edgecolor='white', height=0.7)
ax4.set_yticks(y_pos4)
ax4.set_yticklabels(df_rank['Prefecture'], fontsize=8.5)
ax4.set_xlabel('待機児童率（待機児童数 / 定員数 × 100, %）', fontsize=11)
ax4.set_title('待機児童率の都道府県別ランキング（2022年）\n（保育所等利用待機児童数 / 保育所等定員数 × 100）',
              fontsize=12, fontweight='bold')
ax4.grid(axis='x', alpha=0.3)
ax4.axvline(0, color='gray', linewidth=0.8)

# TFRを右側に注記（上位・下位）
for i, (_, row) in enumerate(df_rank.iterrows()):
    if row['待機児童率'] > 0.05:
        ax4.text(row['待機児童率'] + 0.01, i, f"TFR {row['TFR']:.2f}",
                 va='center', fontsize=7.5, color='#555')

from matplotlib.patches import Patch
legend4 = [
    Patch(facecolor='#C62828', alpha=0.75, label='0.5%超（高い待機圧力）'),
    Patch(facecolor='#FF8F00', alpha=0.75, label='0.1〜0.5%'),
    Patch(facecolor='#1565C0', alpha=0.75, label='0.1%未満'),
]
ax4.legend(handles=legend4, fontsize=9, loc='lower right')

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

print("\n■ 全図の生成完了（4枚）")
print(f"  fig1_timeseries.png : TFRと保育所密度の時系列推移")
print(f"  fig2_scatter.png    : 保育所密度 vs TFR 散布図（2022年）")
print(f"  fig3_coef.png       : OLS回帰係数プロット")
print(f"  fig4_waiting.png    : 待機児童率の都道府県別ランキング")