"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：地域格差と人口移動：固定効果パネルモデルによる転出入要因の分析
受賞：審査員奨励賞（大学生・一般の部）

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

  分析の流れ
  1. 時系列：純移動率の推移（東京圏・地方）
  2. 地域間格差の可視化（都道府県別純移動率 2022年）
  3. パネルOLS固定効果モデル：純移動率の決定要因
  4. 散布図：消費支出 vs 純移動率（2022年）

【被説明変数】
  純移動率 = (転入者数 - 転出者数) / 総人口 × 1000 (‰)

【説明変数】
  求人倍率 = 月間有効求人数 / 月間有効求職者数
  消費支出_log = log(消費支出)
  住宅地価_log = log(標準価格（平均価格）（住宅地）)
  高齢化率 = 65歳以上人口 / 総人口 × 100
  保育所密度 = 保育所等数 / 総人口 × 10000

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

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from linearmodels.panel import PanelOLS, RandomEffects
import statsmodels.api as sm
from scipy import stats

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)

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

# 変数生成
df_b['純移動率'] = (df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']) / df_b['総人口'] * 1000
df_b['求人倍率'] = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'].replace(0, np.nan)
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['住宅地価_log'] = np.log(df_b['標準価格（平均価格）（住宅地）'].clip(lower=1))
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['保育所密度'] = df_b['保育所等数'] / df_b['総人口'] * 10000

region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
df_b['地域'] = df_b['都道府県'].map(region_map)
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12',
}

# Fig1: 純移動率の時系列（代表都道府県）
target_prefs = ['東京都', '大阪府', '愛知県', '秋田県', '島根県', '沖縄県']
fig, ax = plt.subplots(figsize=(10, 5))
for pref in target_prefs:
    grp = df_b[df_b['都道府県'] == pref]
    ax.plot(grp['年度'], grp['純移動率'], marker='o', markersize=4, label=pref)
ax.axhline(0, color='black', linewidth=0.8, linestyle='--')
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('純移動率（‰）', fontsize=12)
ax.set_title('代表都道府県の純移動率推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_11_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: 都道府県別純移動率ランキング（2022年）
df_2022 = df_b[df_b['年度'] == 2022].copy().dropna(subset=['純移動率'])
df_sorted = df_2022.sort_values('純移動率', ascending=True)
fig, ax = plt.subplots(figsize=(8, 12))
bar_colors = ['#e05c5c' if v > 0 else '#4e9af1' for v in df_sorted['純移動率']]
ax.barh(df_sorted['都道府県'], df_sorted['純移動率'], color=bar_colors, alpha=0.8)
ax.axvline(0, color='black', linewidth=1)
ax.set_xlabel('純移動率（‰）', fontsize=12)
ax.set_title('都道府県別 純移動率ランキング（2022年）\n（赤=流入超過, 青=流出超過）', fontsize=13, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_11_fig2_rank.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: パネルOLS固定効果
panel_vars = ['求人倍率', '消費支出_log', '住宅地価_log', '高齢化率', '保育所密度']
df_panel = df_b[['年度', '都道府県', '純移動率'] + panel_vars].dropna()
df_panel = df_panel.set_index(['都道府県', '年度'])
y = df_panel['純移動率']
X = sm.add_constant(df_panel[panel_vars])
try:
    mod = PanelOLS(y, X, entity_effects=True, drop_absorbed=True)
    res = mod.fit(cov_type='clustered', cluster_entity=True)
    coefs = res.params.drop('const', errors='ignore')
    ses = res.std_errors.drop('const', errors='ignore')
    pvals = res.pvalues.drop('const', errors='ignore')
except Exception as e:
    print(f"FE error: {e}")
    df_reg = df_b[df_b['年度'] == 2022][['純移動率'] + panel_vars].dropna()
    res_ols = sm.OLS(df_reg['純移動率'], sm.add_constant(df_reg[panel_vars])).fit()
    coefs = res_ols.params.drop('const')
    ses = res_ols.bse.drop('const')
    pvals = res_ols.pvalues.drop('const')

fig, ax = plt.subplots(figsize=(8, 5))
colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
        error_kw={'elinewidth': 1.5, 'capsize': 4})
ax.set_yticks(range(len(coefs)))
ax.set_yticklabels(coefs.index, fontsize=10)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('回帰係数', fontsize=12)
ax.set_title('純移動率の決定要因 — FEパネルOLS係数\n（赤=p<0.05, 灰=非有意）', fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_11_fig3_fe_coef.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: 散布図 消費支出 vs 純移動率（2022年）
df_sc = df_b[df_b['年度'] == 2022][['都道府県', '純移動率', '消費支出_log', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_sc.groupby('地域'):
    ax.scatter(grp['消費支出_log'], grp['純移動率'],
               color=region_colors.get(reg, 'gray'), label=reg, s=60, alpha=0.8)
for _, row in df_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['消費支出_log'], row['純移動率']),
                fontsize=6, alpha=0.6, xytext=(2, 2), textcoords='offset points')
x_vals = df_sc['消費支出_log']
y_vals = df_sc['純移動率']
slope, intercept, r, p, _ = stats.linregress(x_vals, y_vals)
xr = np.linspace(x_vals.min(), x_vals.max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線 (r={r:.2f}, p={p:.3f})')
ax.axhline(0, color='gray', linewidth=0.8, linestyle=':')
ax.set_xlabel('消費支出（log）', fontsize=12)
ax.set_ylabel('純移動率（‰）', fontsize=12)
ax.set_title('消費支出 vs 純移動率（2022年, N=47）', fontsize=13, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_11_fig4_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")