"""
2023年度 統計データ分析コンペティション 審査員奨励賞（大学生の部）
論文タイトル: 子ども・子育て支援の充実は合計特殊出生率を高めるか？

目的変数:
  合計特殊出生率 (TFR) ← SSDSE-B「合計特殊出生率」

説明変数 (SSDSE-B 2022年度 都道府県データ, 実データのみ):
  nursery_density  : 保育所等数 / 総人口 × 10000     (保育所数密度)
  capacity_ratio   : 保育所等定員数 / 15歳未満人口    (定員/子供比)
  fill_rate        : 保育所等在所児数 / 保育所等定員数 (充足率; 高いほど需要超)
  waitlist_rate    : 保育所等利用待機児童数 / 保育所等在所児数 × 100 (待機率)
  staff_ratio      : 保育所等保育士数 / 保育所等在所児数 (保育士/在所児比)
  consumption      : 消費支出（二人以上の世帯）  (所得代理)
  land_price       : 標準価格（平均価格）（住宅地）(住宅コスト)
  job_offer_rate   : 月間有効求人数(一般) / 月間有効求職者数(一般) (有効求人倍率)
  aging_rate       : 65歳以上人口 / 総人口 × 100 (高齢化率)

パネルチェック: 2012-2023年 保育所数密度の全国平均推移 vs TFR

分析手法: OLS重回帰, VIF(多重共線性), 標準化係数, 交互作用項検討
データ: SSDSE-B-2026.csv（実データ, np.random 一切使用しない）
"""

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


import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

warnings.filterwarnings('ignore')

# ── パス設定 ─────────────────────────────────────────────
_dir    = os.path.dirname(os.path.abspath(__file__)) if '__file__' in dir() else os.getcwd()
FIG_DIR = os.path.join(_dir, '..', 'html', 'figures')
DATA_B  = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-B-2026.csv')

os.makedirs(FIG_DIR, exist_ok=True)

# ── フォント設定 ─────────────────────────────────────────
plt.rcParams.update({
    'font.family'      : ['Hiragino Sans', 'Hiragino Kaku Gothic ProN',
                           'AppleGothic', 'sans-serif'],
    'axes.unicode_minus': False,
    'figure.dpi'       : 150,
})
DPI = 150

# ── データ読み込み ──────────────────────────────────────
print("データを読み込み中...")
df_raw = pd.read_csv(DATA_B, header=1, encoding='cp932')

# 都道府県レベル (R\d{5}) に絞る
mask = df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)
df_raw = df_raw[mask].copy()

# 全数値列を数値型に変換
for col in df_raw.columns:
    if col not in ['年度', '地域コード', '都道府県']:
        df_raw[col] = pd.to_numeric(df_raw[col], errors='coerce')

# 2022年度データ (47都道府県)
df_2022 = df_raw[df_raw['年度'] == 2022].copy()
print(f"2022年度 都道府県データ: {len(df_2022)}件")

# ── 変数の構築 ─────────────────────────────────────────
df = df_2022[['都道府県']].copy()

# 目的変数: 合計特殊出生率
df['tfr'] = df_2022['合計特殊出生率'].values

# 保育所関連説明変数
df['nursery_density'] = (
    df_2022['保育所等数'].values /
    df_2022['総人口'].values * 10000
)  # 保育所数密度 (人口1万人あたり)

df['capacity_ratio'] = (
    df_2022['保育所等定員数'].values /
    df_2022['15歳未満人口'].replace(0, np.nan).values
)  # 定員/15歳未満人口

df['fill_rate'] = (
    df_2022['保育所等在所児数'].values /
    df_2022['保育所等定員数'].replace(0, np.nan).values
)  # 保育所充足率 (高い = 需要が強い)

df['waitlist_rate'] = (
    df_2022['保育所等利用待機児童数'].values /
    df_2022['保育所等在所児数'].replace(0, np.nan).values * 100
)  # 待機児童率 (供給不足の代理)

df['staff_ratio'] = (
    df_2022['保育所等保育士数'].values /
    df_2022['保育所等在所児数'].replace(0, np.nan).values
)  # 保育士/在所児比 (質指標)

# コントロール変数
df['consumption'] = df_2022['消費支出（二人以上の世帯）'].values  # 所得代理
df['land_price']  = df_2022['標準価格（平均価格）（住宅地）'].values  # 住宅コスト
df['job_offer_rate'] = (
    df_2022['月間有効求人数（一般）'].values /
    df_2022['月間有効求職者数（一般）'].replace(0, np.nan).values
)  # 有効求人倍率
df['aging_rate'] = (
    df_2022['65歳以上人口'].values /
    df_2022['総人口'].replace(0, np.nan).values * 100
)  # 高齢化率

df = df.dropna()
print(f"欠損除外後: {len(df)}件 (都道府県)")

# 変数リスト
PRED_COLS = [
    'nursery_density', 'capacity_ratio', 'fill_rate',
    'waitlist_rate', 'staff_ratio',
    'consumption', 'land_price', 'job_offer_rate', 'aging_rate'
]
PRED_LABELS = {
    'nursery_density' : '保育所数密度\n(人口1万対)',
    'capacity_ratio'  : '定員/\n15歳未満人口',
    'fill_rate'       : '充足率\n(在所/定員)',
    'waitlist_rate'   : '待機児童率\n(%)',
    'staff_ratio'     : '保育士/\n在所児比',
    'consumption'     : '消費支出\n(円/世帯)',
    'land_price'      : '住宅地価\n(円/㎡)',
    'job_offer_rate'  : '有効求人\n倍率',
    'aging_rate'      : '高齢化率\n(%)',
}
Y_LABEL = '合計特殊出生率'

X = df[PRED_COLS].astype(float)
y = df['tfr'].astype(float)

print(f"\n目的変数 '{Y_LABEL}' の基本統計:")
print(y.describe().round(4))
print(f"\n説明変数の基本統計:")
print(X.describe().round(4))

# ── VIF 計算 ───────────────────────────────────────────
X_vif = sm.add_constant(X)
vif_vals = []
for i, col in enumerate(X_vif.columns):
    if col == 'const':
        continue
    vif_vals.append({
        'variable': col,
        'label'   : PRED_LABELS[col].replace('\n', ' '),
        'VIF'     : variance_inflation_factor(X_vif.values.astype(float), i)
    })
vif_df = pd.DataFrame(vif_vals)
print("\nVIF:")
print(vif_df[['variable', 'VIF']].round(2).to_string(index=False))

# ── OLS 回帰 ───────────────────────────────────────────
X_ols = sm.add_constant(X)
model = sm.OLS(y, X_ols).fit()
print("\n" + "="*60)
print("OLS 回帰サマリー")
print("="*60)
print(model.summary())
print(f"\nR² = {model.rsquared:.4f}, Adj.R² = {model.rsquared_adj:.4f}")
print(f"F-stat = {model.fvalue:.2f}, p = {model.f_pvalue:.4f}")

# 標準化係数
X_std = (X - X.mean()) / X.std()
y_std = (y - y.mean()) / y.std()
X_std_ols = sm.add_constant(X_std)
model_std = sm.OLS(y_std, X_std_ols).fit()
coef_std = model_std.params.drop('const')
pvals    = model.pvalues.drop('const')

# 交互作用項: 保育所密度 × 高齢化率
df['interact_nursery_aging'] = df['nursery_density'] * df['aging_rate']
X_ia = X.copy()
X_ia['nursery×aging'] = df['interact_nursery_aging'].values
X_ia_std = sm.add_constant(X_ia)
model_ia = sm.OLS(y, X_ia_std).fit()
print("\n--- 交互作用項 (保育所密度 × 高齢化率) 追加モデル ---")
print(f"R² = {model_ia.rsquared:.4f}, Adj.R² = {model_ia.rsquared_adj:.4f}")
ia_coef = model_ia.params.get('nursery×aging', np.nan)
ia_p    = model_ia.pvalues.get('nursery×aging', np.nan)
print(f"交互作用項係数: {ia_coef:.6f}, p = {ia_p:.4f}")

# ── パネルデータ: 2012-2023年 全国平均推移 ────────────────
years_all = sorted(df_raw['年度'].unique())
panel_list = []
for yr in years_all:
    dfy = df_raw[df_raw['年度'] == yr].copy()
    tfr_mean  = dfy['合計特殊出生率'].mean()
    nursery_d = (
        dfy['保育所等数'] / dfy['総人口'].replace(0, np.nan) * 10000
    ).mean()
    panel_list.append({'year': yr, 'tfr': tfr_mean, 'nursery_density': nursery_d})
panel_df = pd.DataFrame(panel_list).sort_values('year')
print("\n--- パネルデータ (全国平均) ---")
print(panel_df.to_string(index=False))

# ────────────────────────────────────────────────────────
# 図1: 散布図行列 (保育所密度・待機率・充足率 vs TFR)
# ────────────────────────────────────────────────────────
fig1, axes = plt.subplots(1, 3, figsize=(14, 5))
fig1.suptitle('子育て支援指標と合計特殊出生率の関係\n（2022年度, 47都道府県）',
               fontsize=13, fontweight='bold', y=1.02)

scatter_vars = [
    ('nursery_density', '保育所数密度\n（人口1万人あたり）', '#1565C0'),
    ('fill_rate',       '保育所充足率\n（在所児/定員）',     '#E65100'),
    ('waitlist_rate',   '待機児童率\n（待機/在所 ×100）',   '#C62828'),
]

for ax, (vcol, vlabel, color) in zip(axes, scatter_vars):
    x_vals = df[vcol].values
    y_vals = y.values

    ax.scatter(x_vals, y_vals, color=color, alpha=0.7, edgecolors='white',
               linewidths=0.5, s=60)

    # 回帰直線
    slope, intercept, r_val, p_val, _ = \
        __import__('scipy').stats.linregress(x_vals, y_vals)
    x_line = np.linspace(x_vals.min(), x_vals.max(), 100)
    ax.plot(x_line, slope * x_line + intercept, color=color, lw=2, linestyle='--')

    # 都道府県ラベル (外れ値のみ)
    threshold = 1.5 * np.std(y_vals)
    y_mean = y_vals.mean()
    for i, (xi, yi, pref) in enumerate(
            zip(x_vals, y_vals, df['都道府県'].values)):
        if abs(yi - y_mean) > threshold:
            ax.annotate(pref, (xi, yi), fontsize=7, ha='left',
                        xytext=(3, 3), textcoords='offset points', color='#333')

    sig_str = '***' if p_val < 0.001 else '**' if p_val < 0.01 \
        else '*' if p_val < 0.05 else 'n.s.'
    ax.set_xlabel(vlabel, fontsize=10)
    ax.set_ylabel(Y_LABEL if ax == axes[0] else '', fontsize=10)
    ax.set_title(f'r = {r_val:.3f} ({sig_str})', fontsize=11)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2023_U5_6_fig1_scatter.png')
fig1.savefig(fig1_path, dpi=DPI, bbox_inches='tight')
plt.close(fig1)
print(f"\n図1保存: {fig1_path}")

# ────────────────────────────────────────────────────────
# 図2: 標準化係数の棒グラフ
# ────────────────────────────────────────────────────────
fig2, ax2 = plt.subplots(figsize=(10, 6))

labels_ordered = [PRED_LABELS[c] for c in PRED_COLS]
coef_vals = [coef_std[c] for c in PRED_COLS]
p_vals_lst = [pvals[c] for c in PRED_COLS]

colors = []
for cv, pv in zip(coef_vals, p_vals_lst):
    if pv < 0.05:
        colors.append('#1565C0' if cv > 0 else '#C62828')
    else:
        colors.append('#90A4AE')

bars = ax2.barh(labels_ordered, coef_vals, color=colors, edgecolor='white',
                linewidth=0.5, height=0.6)

for bar, cv, pv in zip(bars, coef_vals, p_vals_lst):
    sig = '***' if pv < 0.001 else '**' if pv < 0.01 \
        else '*' if pv < 0.05 else 'n.s.'
    xpos = cv + (0.01 if cv >= 0 else -0.01)
    ha   = 'left' if cv >= 0 else 'right'
    ax2.text(xpos, bar.get_y() + bar.get_height() / 2,
             f'p={pv:.3f} {sig}', va='center', ha=ha, fontsize=9, color='#333')

ax2.axvline(0, color='#333', lw=1)
ax2.set_xlabel('標準化回帰係数 (β)', fontsize=11)
ax2.set_title(f'OLS重回帰 標準化係数\n（目的変数: 合計特殊出生率, R²={model.rsquared:.3f}）',
               fontsize=12, fontweight='bold')
ax2.grid(True, axis='x', alpha=0.3)

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#1565C0', label='正の有意効果 (p<0.05)'),
    Patch(facecolor='#C62828', label='負の有意効果 (p<0.05)'),
    Patch(facecolor='#90A4AE', label='非有意 (p≥0.05)'),
]
ax2.legend(handles=legend_elements, loc='lower right', fontsize=9)

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

# ────────────────────────────────────────────────────────
# 図3: VIF棒グラフ
# ────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(9, 5))

vif_labels = [PRED_LABELS[c].replace('\n', ' ') for c in PRED_COLS]
vif_values = [vif_df[vif_df['variable'] == c]['VIF'].values[0] for c in PRED_COLS]
vif_colors = ['#C62828' if v > 10 else '#F9A825' if v > 5 else '#2E7D32'
              for v in vif_values]

bars3 = ax3.barh(vif_labels, vif_values, color=vif_colors, edgecolor='white',
                 linewidth=0.5, height=0.6)

for bar, v in zip(bars3, vif_values):
    ax3.text(v + 0.15, bar.get_y() + bar.get_height() / 2,
             f'{v:.2f}', va='center', ha='left', fontsize=9, color='#333')

ax3.axvline(5,  color='#F9A825', lw=1.5, linestyle='--', alpha=0.8, label='VIF=5 (警告)')
ax3.axvline(10, color='#C62828', lw=1.5, linestyle='--', alpha=0.8, label='VIF=10 (問題)')
ax3.set_xlabel('VIF（分散膨張係数）', fontsize=11)
ax3.set_title('多重共線性の診断：VIF\n（VIF<5: 問題なし, 5≤VIF<10: 要注意, VIF≥10: 問題）',
               fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, axis='x', alpha=0.3)
ax3.set_xlim(0, max(vif_values) * 1.25 + 2)

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

# ────────────────────────────────────────────────────────
# 図4: 時系列 (2012-2023) 保育所密度 vs TFR (双軸)
# ────────────────────────────────────────────────────────
fig4, ax4a = plt.subplots(figsize=(11, 5))
ax4b = ax4a.twinx()

years_plot   = panel_df['year'].values
tfr_plot     = panel_df['tfr'].values
density_plot = panel_df['nursery_density'].values

line1, = ax4a.plot(years_plot, density_plot, color='#1565C0', lw=2.5,
                   marker='o', markersize=6, label='保育所数密度\n（人口1万対, 左軸）')
line2, = ax4b.plot(years_plot, tfr_plot, color='#E65100', lw=2.5,
                   marker='s', markersize=6, label='合計特殊出生率\n（右軸）',
                   linestyle='--')

ax4a.set_xlabel('年度', fontsize=11)
ax4a.set_ylabel('保育所数密度（人口1万人あたり）', fontsize=11, color='#1565C0')
ax4b.set_ylabel('合計特殊出生率', fontsize=11, color='#E65100')
ax4a.tick_params(axis='y', labelcolor='#1565C0')
ax4b.tick_params(axis='y', labelcolor='#E65100')

ax4a.set_title('全国平均：保育所数密度と合計特殊出生率の推移（2012–2023年）\n'
                '（47都道府県平均）',
                fontsize=12, fontweight='bold')

# 注釈: 2015年ピーク
ax4b.annotate('TFR 局所最高\n(2015: 1.53)',
               xy=(2015, panel_df[panel_df['year'] == 2015]['tfr'].values[0]),
               xytext=(2016.5, 1.50),
               fontsize=9, color='#E65100',
               arrowprops=dict(arrowstyle='->', color='#E65100', lw=1.2))

# 凡例を統合
lines = [line1, line2]
labels = [l.get_label() for l in lines]
ax4a.legend(lines, labels, loc='upper left', fontsize=9)
ax4a.grid(True, alpha=0.3)
ax4a.set_xticks(years_plot)
ax4a.set_xticklabels([str(y) for y in years_plot], rotation=45, ha='right', fontsize=9)

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

# ── 結果サマリー ────────────────────────────────────────
print("\n" + "="*60)
print("分析結果サマリー")
print("="*60)
print(f"サンプルサイズ : {len(df)} 都道府県 (2022年度)")
print(f"R²             : {model.rsquared:.4f}")
print(f"Adj.R²         : {model.rsquared_adj:.4f}")
print(f"F統計量        : {model.fvalue:.2f}  (p={model.f_pvalue:.4f})")
print()
print("標準化係数と有意性:")
for col in PRED_COLS:
    lbl = PRED_LABELS[col].replace('\n', ' ')
    beta = coef_std[col]
    pv   = pvals[col]
    sig  = '***' if pv < 0.001 else '**' if pv < 0.01 \
           else '*' if pv < 0.05 else 'n.s.'
    print(f"  {lbl:25s}: β={beta:+.4f}  p={pv:.4f} {sig}")
print()
print("VIF (多重共線性診断):")
for col in PRED_COLS:
    vif_v = vif_df[vif_df['variable'] == col]['VIF'].values[0]
    flag  = '[問題]' if vif_v > 10 else '[要注意]' if vif_v > 5 else '[OK]'
    print(f"  {PRED_LABELS[col].replace(chr(10), ' '):25s}: VIF={vif_v:.2f} {flag}")
print()
print(f"交互作用項 (保育所密度 × 高齢化率): β={ia_coef:.6f}, p={ia_p:.4f}")
print()
print("生成図:")
for path in [fig1_path, fig2_path, fig3_path, fig4_path]:
    print(f"  {path}")