"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：市町村予算から見る社会増減の時系列要因分析
受賞：審査員奨励賞（大学生・一般の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度）
  対象：全47都道府県 × 12年（2012〜2023）
  目的変数：社会増減率 = (転入者数 - 転出者数) / 総人口 × 1000 (‰)

  分析の流れ
  ・時系列：全国平均 + 東京・大阪・地方代表（秋田・島根）の社会増減率推移
  ・パネルOLS（固定効果モデル）：社会増減の決定要因推定
  ・求人倍率 vs 社会増減率 散布図（2022年）
  ・都道府県別社会増減率ランキング（2022年）

【被説明変数】
  社会増減率（‰）= (転入者数 - 転出者数) / 総人口 × 1000

【説明変数（政策関連）】
  求人倍率         = 月間有効求人数 / 月間有効求職者数
  消費支出（万円）  = 消費支出（二人以上の世帯）/ 10000
  住宅地価格（万円）= 標準価格（平均価格）（住宅地）/ 10000
  大学学生数（千人）= 大学学生数 / 1000
  保育所等数       = 保育所等数

【推定手法】
  linearmodels PanelOLS: entity_effects=True, cov_type='clustered'

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

【データサイエンス学習ポイント】
  1. 時系列パネルデータの整形とインデックス設定
  2. 固定効果モデル（FE）によるパネルOLS推定
  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_U5_7_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from linearmodels import PanelOLS
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: 都道府県別パネルデータ）
# ================================================================
print("=" * 65)
print("■ SSDSE-B-2026 データ読み込み")
print("=" * 65)

df_b_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b_raw[df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

print(f"都道府県数: {df_b['都道府県'].nunique()}")
print(f"年度範囲  : {df_b['年度'].min()}〜{df_b['年度'].max()}")
print(f"総行数    : {len(df_b)} (47都道府県 × 12年)")

# ================================================================
# ■ 変数の構築
# ================================================================
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

# 目的変数：社会増減率（‰）
df_b['社会増減']   = df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']
df_b['社会増減率'] = df_b['社会増減'] / df_b['総人口'] * 1000

# 説明変数
df_b['求人倍率']         = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'].clip(lower=1)
df_b['消費支出_万円']     = df_b['消費支出（二人以上の世帯）'] / 10000
df_b['住宅地価格_万円']   = df_b['標準価格（平均価格）（住宅地）'] / 10000
df_b['大学学生数_千人']   = df_b['大学学生数'] / 1000
df_b['保育所等数']        = df_b['保育所等数']

print("\n■ 社会増減率 基本統計（‰）")
print(df_b['社会増減率'].describe().round(3))

# ================================================================
# ■ 都道府県区分（地域グループ）
# ================================================================
metro_prefs  = ['東京都', '大阪府']
rural_prefs  = ['秋田県', '島根県']
highlight_prefs = metro_prefs + rural_prefs

region_colors = {
    '東京都': '#C62828',
    '大阪府': '#1565C0',
    '秋田県': '#2E7D32',
    '島根県': '#6A1B9A',
}

# ================================================================
# ■ パネルOLS（固定効果モデル）
# ================================================================
print("\n" + "=" * 65)
print("■ パネルOLS 固定効果モデル推定")
print("=" * 65)

# パネル用インデックス設定
df_panel = df_b.copy()
df_panel = df_panel.set_index(['都道府県', '年度'])

XVARS = ['求人倍率', '消費支出_万円', '住宅地価格_万円', '大学学生数_千人', '保育所等数']
YVAR  = '社会増減率'

panel_data = df_panel[[YVAR] + XVARS].dropna()
print(f"推定サンプル: {len(panel_data)} obs")

mod = PanelOLS(
    dependent=panel_data[YVAR],
    exog=panel_data[XVARS],
    entity_effects=True,
    time_effects=False,
)
res = mod.fit(cov_type='clustered', cluster_entity=True)
print(res.summary)

# 係数・標準誤差・p値を整理
fe_params   = res.params
fe_stderr   = res.std_errors
fe_pvalues  = res.pvalues
fe_rsq      = res.rsquared

print("\n■ 推定結果まとめ")
print(f"  Within R²: {fe_rsq:.4f}")
for var in XVARS:
    sig = '***' if fe_pvalues[var] < 0.001 else '**' if fe_pvalues[var] < 0.01 else '*' if fe_pvalues[var] < 0.05 else 'n.s.'
    print(f"  {var:<18} β={fe_params[var]:+.4f}  SE={fe_stderr[var]:.4f}  p={fe_pvalues[var]:.4f}  {sig}")

# ================================================================
# ■ 図1：社会増減率 時系列（全国 + 都市・地方代表）
# ================================================================
print("\n" + "=" * 65)
print("■ 図1：社会増減率 時系列")
print("=" * 65)

# 全国平均（47都道府県の単純平均）
nat_avg = df_b.groupby('年度')['社会増減率'].mean()

fig1, ax1 = plt.subplots(figsize=(10, 5.5))
ax1.axhline(0, color='gray', linewidth=0.8, linestyle='--', alpha=0.6)

# 全国平均
ax1.plot(nat_avg.index, nat_avg.values, color='#455A64', linewidth=2.2,
         linestyle='--', label='全国平均', zorder=3)

# ハイライト都道府県
for pref in highlight_prefs:
    tmp = df_b[df_b['都道府県'] == pref].sort_values('年度')
    ax1.plot(tmp['年度'], tmp['社会増減率'],
             color=region_colors[pref], linewidth=2.0,
             marker='o', markersize=4, label=pref, zorder=4)

# 背景グレー（全47都道府県）
for pref in df_b['都道府県'].unique():
    if pref not in highlight_prefs:
        tmp = df_b[df_b['都道府県'] == pref].sort_values('年度')
        ax1.plot(tmp['年度'], tmp['社会増減率'],
                 color='#BDBDBD', linewidth=0.6, alpha=0.45, zorder=1)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('社会増減率（‰）', fontsize=12)
ax1.set_title('社会増減率（転入-転出）の時系列推移\n全47都道府県・2012〜2023年度（SSDSE-B実データ）',
              fontsize=13, fontweight='bold')
ax1.legend(loc='upper right', fontsize=9, framealpha=0.9)
ax1.grid(True, alpha=0.3)
ax1.set_xticks(sorted(df_b['年度'].unique()))
ax1.tick_params(axis='x', rotation=45)

# 注釈：COVID-19
ax1.axvspan(2019.5, 2021.5, color='#FFECB3', alpha=0.4, zorder=0)
ax1.text(2020.5, ax1.get_ylim()[0] * 0.85 if ax1.get_ylim()[0] < 0 else -0.5,
         'COVID-19', ha='center', fontsize=8, color='#795548', style='italic')

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

# ================================================================
# ■ 図2：FE係数プロット（社会増減の決定要因）
# ================================================================
print("\n■ 図2：FE係数プロット")

var_labels = {
    '求人倍率'       : '求人倍率',
    '消費支出_万円'   : '消費支出（万円）',
    '住宅地価格_万円' : '住宅地価格（万円）',
    '大学学生数_千人' : '大学学生数（千人）',
    '保育所等数'      : '保育所等数',
}

coefs  = [fe_params[v] for v in XVARS]
ses    = [fe_stderr[v] for v in XVARS]
pvals  = [fe_pvalues[v] for v in XVARS]
labels = [var_labels[v] for v in XVARS]

bar_colors = []
for p in pvals:
    if p < 0.01:
        bar_colors.append('#C62828')
    elif p < 0.05:
        bar_colors.append('#FF8F00')
    else:
        bar_colors.append('#9E9E9E')

fig2, ax2 = plt.subplots(figsize=(9, 5))
ypos = np.arange(len(XVARS))

bars = ax2.barh(ypos, coefs, color=bar_colors, alpha=0.8,
                edgecolor='white', height=0.55)
ax2.errorbar(coefs, ypos, xerr=1.96 * np.array(ses),
             fmt='none', color='#333', capsize=5, linewidth=1.5)
ax2.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax2.set_yticks(ypos)
ax2.set_yticklabels(labels, fontsize=11)
ax2.set_xlabel('固定効果モデル 推定係数（±1.96 SE）', fontsize=11)
ax2.set_title(f'社会増減率の決定要因（パネルOLS 固定効果モデル）\n'
              f'Within R²={fe_rsq:.3f}、N=47都道府県×12年（クラスタリング標準誤差）',
              fontsize=12, fontweight='bold')
ax2.invert_yaxis()
ax2.grid(axis='x', alpha=0.3)

# p値注釈
for i, (c, p) in enumerate(zip(coefs, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    offset = 1.96 * ses[i] + abs(max(coefs) - min(coefs)) * 0.02
    ha = 'left' if c >= 0 else 'right'
    x_pos = c + (1.96 * ses[i] + abs(max(coefs) - min(coefs)) * 0.01) if c >= 0 else c - (1.96 * ses[i] + abs(max(coefs) - min(coefs)) * 0.01)
    ax2.text(x_pos, i, sig, va='center', ha=ha, fontsize=9, color=bar_colors[i], fontweight='bold')

legend_patches = [
    mpatches.Patch(color='#C62828', alpha=0.8, label='p < 0.01'),
    mpatches.Patch(color='#FF8F00', alpha=0.8, label='p < 0.05'),
    mpatches.Patch(color='#9E9E9E', alpha=0.8, label='n.s. (p ≥ 0.05)'),
]
ax2.legend(handles=legend_patches, fontsize=9, loc='lower right')

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

# ================================================================
# ■ 図3：求人倍率 vs 社会増減率 散布図（2022年）
# ================================================================
print("\n■ 図3：求人倍率 vs 社会増減率 散布図（2022年）")

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

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

# 全都道府県プロット（グレー）
ax3.scatter(df_2022['求人倍率'], df_2022['社会増減率'],
            color='#78909C', s=50, alpha=0.6, zorder=2, label='その他都道府県')

# ハイライト都道府県
for pref in highlight_prefs:
    row = df_2022[df_2022['都道府県'] == pref]
    if len(row) > 0:
        ax3.scatter(row['求人倍率'], row['社会増減率'],
                    color=region_colors[pref], s=130, zorder=5,
                    edgecolors='white', linewidths=1.2, label=pref)
        ax3.annotate(pref,
                     xy=(row['求人倍率'].values[0], row['社会増減率'].values[0]),
                     xytext=(6, 4), textcoords='offset points',
                     fontsize=9, color=region_colors[pref], fontweight='bold')

# 回帰直線
x_vals = df_2022['求人倍率'].values
y_vals = df_2022['社会増減率'].values
z = np.polyfit(x_vals, y_vals, 1)
xs = np.linspace(x_vals.min(), x_vals.max(), 100)
ax3.plot(xs, np.poly1d(z)(xs), color='#C62828', linewidth=1.8,
         linestyle='--', alpha=0.8, label=f'回帰直線', zorder=3)

from scipy import stats as sp_stats
r_val, p_val = sp_stats.pearsonr(x_vals, y_vals)

ax3.axhline(0, color='gray', linewidth=0.8, linestyle=':', alpha=0.7)
ax3.set_xlabel('求人倍率（月間有効求人数 / 月間有効求職者数）', fontsize=12)
ax3.set_ylabel('社会増減率（‰）', fontsize=12)
ax3.set_title(f'求人倍率と社会増減率の関係（2022年度・47都道府県）\n'
              f'ピアソン相関: r = {r_val:.3f}（p = {p_val:.4f}）',
              fontsize=12, fontweight='bold')
ax3.legend(fontsize=9, loc='upper left')
ax3.grid(True, alpha=0.3)

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

# ================================================================
# ■ 図4：社会増減率 都道府県ランキング（2022年）
# ================================================================
print("\n■ 図4：社会増減率 都道府県ランキング（2022年）")

df_rank = df_2022[['都道府県', '社会増減率']].sort_values('社会増減率', ascending=True).reset_index(drop=True)

# 色付け：正（流入超）= 暖色、負（流出超）= 寒色
rank_colors = ['#C62828' if v >= 0 else '#1565C0' for v in df_rank['社会増減率']]

fig4, ax4 = plt.subplots(figsize=(9, 13))
ypos4 = np.arange(len(df_rank))

bars4 = ax4.barh(ypos4, df_rank['社会増減率'], color=rank_colors,
                 alpha=0.8, edgecolor='white', height=0.7)
ax4.axvline(0, color='#333', linewidth=0.9, linestyle='-')
ax4.set_yticks(ypos4)
ax4.set_yticklabels(df_rank['都道府県'], fontsize=9)
ax4.set_xlabel('社会増減率（‰）', fontsize=12)
ax4.set_title('社会増減率 都道府県ランキング（2022年度）\n正値＝転入超過（流入）、負値＝転出超過（流出）',
              fontsize=12, fontweight='bold')
ax4.grid(axis='x', alpha=0.3)

# 値ラベル
for i, (v, pref) in enumerate(zip(df_rank['社会増減率'], df_rank['都道府県'])):
    ha = 'left' if v >= 0 else 'right'
    offset = 0.05 if v >= 0 else -0.05
    ax4.text(v + offset, i, f'{v:+.2f}', va='center', ha=ha, fontsize=7.5,
             color='#333')

# ハイライト：全国平均線
nat_avg_2022 = df_2022['社会増減率'].mean()
ax4.axvline(nat_avg_2022, color='#FF8F00', linewidth=1.5, linestyle='--',
            alpha=0.9, label=f'全国平均 {nat_avg_2022:+.2f}‰')
ax4.legend(fontsize=10, loc='lower right')

legend_patches2 = [
    mpatches.Patch(color='#C62828', alpha=0.8, label='転入超過（流入）'),
    mpatches.Patch(color='#1565C0', alpha=0.8, label='転出超過（流出）'),
]
ax4.legend(handles=legend_patches2 + [
    mpatches.Patch(color='#FF8F00', alpha=0.9, label=f'全国平均 {nat_avg_2022:+.2f}‰')
], fontsize=9, loc='lower right')

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

# ================================================================
# ■ 完了サマリ
# ================================================================
print("\n" + "=" * 65)
print("■ 全図の生成完了（4枚）")
print("=" * 65)
print(f"  fig1_ts.png      : 社会増減率 時系列（全国 + 都市・地方代表）")
print(f"  fig2_fe_coef.png : FE係数プロット（社会増減の決定要因）")
print(f"  fig3_scatter.png : 求人倍率 vs 社会増減率 散布図（2022年）")
print(f"  fig4_rank.png    : 社会増減率 都道府県ランキング（2022年）")
print(f"\n  パネルOLS結果（固定効果モデル、クラスタリング標準誤差）：")
for var in XVARS:
    sig = '***' if fe_pvalues[var] < 0.001 else '**' if fe_pvalues[var] < 0.01 else '*' if fe_pvalues[var] < 0.05 else 'n.s.'
    print(f"    {var_labels[var]:<18} β={fe_params[var]:+.4f}  p={fe_pvalues[var]:.4f}  {sig}")
print(f"\n  Within R² = {fe_rsq:.4f}")