"""
2021年度 統計データ分析コンペティション 審査員奨励賞（高校生の部）
「都道府県別家計所得格差と社会保障支出の効果分析」
教育用再現スクリプト

使用データ: 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/2021_H5_5_shorei.py
# ============================================================


import os
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

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

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)

# 2022年断面データ
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)

# ── 変数の定義と派生変数計算 ──────────────────────────────────────────
# 所得代理変数: 消費支出（二人以上の世帯）
df['消費支出'] = df['消費支出（二人以上の世帯）']

# 大学進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数
df['大学進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数']

# 高齢化率 = 65歳以上人口 / 総人口
df['高齢化率'] = df['65歳以上人口'] / df['総人口']

# 有効求人倍率（労働市場の活況 = 第三次産業・経済の代理）
df['有効求人倍率'] = df['月間有効求人数（一般）'] / df['月間有効求職者数（一般）']

# 転入超過率（人口移動の活発さ）
df['転入超過率'] = (df['転入者数（日本人移動者）'] - df['転出者数（日本人移動者）']) / df['総人口']

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

df['地域'] = df['都道府県'].map(region_map)
df['color'] = df['地域'].map(region_colors)

# ── 格差指標の計算 ─────────────────────────────────────────────────
income = df['消費支出'].values
cv = income.std() / income.mean()  # 変動係数
p90p10 = np.percentile(income, 90) / np.percentile(income, 10)  # 90/10分位数比

# ジニ係数
def gini(x):
    x = np.sort(x)
    n = len(x)
    cumx = np.cumsum(x)
    return (n + 1 - 2 * np.sum(cumx) / cumx[-1]) / n

gini_coef = gini(income)

print("=== 格差指標（2022年断面） ===")
print(f"平均消費支出: {income.mean():,.0f} 円/月")
print(f"変動係数 (CV): {cv:.4f}")
print(f"90/10分位数比: {p90p10:.4f}")
print(f"ジニ係数: {gini_coef:.4f}")
print(f"最大 - {df.loc[df['消費支出'].idxmax(), '都道府県']}: {income.max():,.0f} 円")
print(f"最小 - {df.loc[df['消費支出'].idxmin(), '都道府県']}: {income.min():,.0f} 円")

# ── 重回帰分析 ─────────────────────────────────────────────────────
X_cols = ['大学進学率', '高齢化率', '有効求人倍率', '転入超過率']
X_raw = df[X_cols].copy()
X_std = (X_raw - X_raw.mean()) / X_raw.std()
Y_raw = df['消費支出']
Y_std = (Y_raw - Y_raw.mean()) / Y_raw.std()

X_sm = sm.add_constant(X_std)
model = sm.OLS(Y_std, X_sm).fit()

print("\n=== 重回帰分析結果（標準化変数） ===")
print(model.summary())

coefs   = model.params[1:]
pvalues = model.pvalues[1:]
conf_int = model.conf_int().iloc[1:]

print("\n標準化偏回帰係数:")
for col, b, p in zip(X_cols, coefs, pvalues):
    sig = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "n.s."
    print(f"  {col}: β={b:.4f}, p={p:.4f} {sig}")
print(f"\nR²={model.rsquared:.3f}, 調整済みR²={model.rsquared_adj:.3f}")

# 相関分析
print("\n=== 相関分析 ===")
for col in X_cols:
    r, p = stats.pearsonr(df[col].dropna(), df.loc[df[col].notna(), '消費支出'])
    print(f"  {col} vs 消費支出: r={r:.4f}, p={p:.4f}")

# ── 図1: 都道府県別消費支出ランキング（棒グラフ） ──────────────────
df_sorted = df.sort_values('消費支出', ascending=True).reset_index(drop=True)

fig, ax = plt.subplots(figsize=(14, 10))
bars = ax.barh(
    df_sorted['都道府県'],
    df_sorted['消費支出'] / 1000,
    color=[region_colors[r] for r in df_sorted['地域']],
    edgecolor='white', linewidth=0.3
)
ax.axvline(income.mean() / 1000, color='#333', linestyle='--', linewidth=1.5, label=f'全国平均 {income.mean()/1000:.1f}千円')
ax.set_xlabel('消費支出（千円/月, 二人以上の世帯）', fontsize=12)
ax.set_title('都道府県別 消費支出ランキング（2022年）\n─ 家計所得格差の可視化 ─', fontsize=14, fontweight='bold')
ax.legend(fontsize=11, loc='lower right')

# 凡例パッチ
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax.legend(handles=legend_elements + [
    plt.Line2D([0], [0], color='#333', linestyle='--', linewidth=1.5, label=f'全国平均 {income.mean()/1000:.1f}千円')
], fontsize=10, loc='lower right', ncol=2)

ax.tick_params(axis='y', labelsize=9)
ax.tick_params(axis='x', labelsize=10)
ax.set_xlim(200, 350)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_5_fig1.png'), bbox_inches='tight')
plt.close(fig)
print("\n図1 保存完了: 2021_H5_5_fig1.png")

# ── 図2: ローレンツ曲線 ──────────────────────────────────────────
income_sorted = np.sort(income)
n = len(income_sorted)
lorenz_x = np.concatenate([[0], np.arange(1, n + 1) / n])
lorenz_y = np.concatenate([[0], np.cumsum(income_sorted) / income_sorted.sum()])

fig, ax = plt.subplots(figsize=(8, 7))
ax.plot([0, 1], [0, 1], 'k--', linewidth=1.5, label='完全均等分布線')
ax.plot(lorenz_x, lorenz_y, color='#e05c5c', linewidth=2.5, label=f'ローレンツ曲線（ジニ係数={gini_coef:.4f}）')
ax.fill_between(lorenz_x, lorenz_x, lorenz_y, alpha=0.15, color='#e05c5c', label='格差面積')

# 注目点を追加
idx10 = int(0.1 * n)
idx90 = int(0.9 * n)
ax.annotate(
    f'下位10%が得る所得シェア\n= {lorenz_y[idx10]:.1%}',
    xy=(lorenz_x[idx10], lorenz_y[idx10]),
    xytext=(0.2, 0.05),
    arrowprops=dict(arrowstyle='->', color='#333', lw=1.2),
    fontsize=10, color='#333'
)

ax.set_xlabel('都道府県の累積比率（所得の低い順）', fontsize=12)
ax.set_ylabel('消費支出の累積シェア', fontsize=12)
ax.set_title('ローレンツ曲線：都道府県間消費支出格差（2022年）', fontsize=13, fontweight='bold')
ax.legend(fontsize=11, loc='upper left')
ax.set_xlim(0, 1); ax.set_ylim(0, 1)
ax.grid(True, alpha=0.3)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_5_fig2.png'), bbox_inches='tight')
plt.close(fig)
print("図2 保存完了: 2021_H5_5_fig2.png")

# ── 図3: 有効求人倍率 vs 消費支出（散布図・地域色分け） ──────────────
fig, ax = plt.subplots(figsize=(11, 8))

for region, color in region_colors.items():
    mask = df['地域'] == region
    ax.scatter(
        df.loc[mask, '有効求人倍率'],
        df.loc[mask, '消費支出'] / 1000,
        c=color, s=80, alpha=0.85, edgecolors='white', linewidth=0.5, label=region
    )

# 都道府県ラベル
for _, row in df.iterrows():
    ax.annotate(
        row['都道府県'].replace('都', '').replace('道', '').replace('府', '').replace('県', ''),
        (row['有効求人倍率'], row['消費支出'] / 1000),
        fontsize=7, ha='center', va='bottom',
        xytext=(0, 4), textcoords='offset points', color='#333'
    )

# 回帰直線
x_reg = df['有効求人倍率'].values
y_reg = df['消費支出'].values / 1000
slope, intercept, r_val, p_val, _ = stats.linregress(x_reg, y_reg)
x_line = np.linspace(x_reg.min(), x_reg.max(), 100)
ax.plot(x_line, slope * x_line + intercept, color='#333', linewidth=1.8,
        linestyle='--', label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})')

ax.set_xlabel('有効求人倍率（労働市場の活況度）', fontsize=12)
ax.set_ylabel('消費支出（千円/月）', fontsize=12)
ax.set_title('有効求人倍率 vs 消費支出（2022年，47都道府県）\n─ 労働市場の強さと家計所得の関係 ─',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=10, loc='upper left', ncol=2)
ax.grid(True, alpha=0.25)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_5_fig3.png'), bbox_inches='tight')
plt.close(fig)
print("図3 保存完了: 2021_H5_5_fig3.png")

# ── 図4: 標準化偏回帰係数プロット ────────────────────────────────
var_labels = ['大学進学率', '高齢化率', '有効求人倍率', '転入超過率']
betas  = coefs.values
ci_lo  = (betas - conf_int.iloc[:, 0].values)
ci_hi  = (conf_int.iloc[:, 1].values - betas)

colors_bar = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvalues]

fig, ax = plt.subplots(figsize=(9, 5))
y_pos = range(len(var_labels))
bars2 = ax.barh(y_pos, betas, xerr=[ci_lo, ci_hi], color=colors_bar,
                align='center', edgecolor='white', linewidth=0.5,
                error_kw=dict(ecolor='#555', elinewidth=1.5, capsize=5))
ax.axvline(0, color='#333', linewidth=1.0)
ax.set_yticks(y_pos)
ax.set_yticklabels(var_labels, fontsize=12)
ax.set_xlabel('標準化偏回帰係数 (β) ± 95%CI', fontsize=12)
ax.set_title('消費支出（所得）格差の規定要因\n─ 重回帰分析・標準化偏回帰係数（2022年, N=47）─',
             fontsize=13, fontweight='bold')
ax.set_xlim(-0.4, 0.9)
ax.grid(True, axis='x', alpha=0.3)

# β値をバー上に表示
for i, (b, p) in enumerate(zip(betas, pvalues)):
    sig_str = '*' if p < 0.05 else ''
    ax.text(b + (0.04 if b >= 0 else -0.04), i,
            f'β={b:.3f}{sig_str}', va='center',
            ha='left' if b >= 0 else 'right', fontsize=10, fontweight='bold')

# 凡例
from matplotlib.patches import Patch as MPatch
legend_handles = [
    MPatch(facecolor='#e05c5c', label='有意 (p<0.05)'),
    MPatch(facecolor='#aaaaaa', label='非有意')
]
ax.legend(handles=legend_handles, fontsize=10, loc='lower right')
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_5_fig4.png'), bbox_inches='tight')
plt.close(fig)
print("図4 保存完了: 2021_H5_5_fig4.png")

print("\n=== 全図の生成完了 ===")
print(f"  図1: 都道府県別消費支出ランキング")
print(f"  図2: ローレンツ曲線（ジニ係数={gini_coef:.4f}）")
print(f"  図3: 有効求人倍率 vs 消費支出（散布図）")
print(f"  図4: 標準化偏回帰係数プロット（R²={model.rsquared:.3f}）")
print(f"\n格差サマリー:")
print(f"  変動係数: {cv:.4f}")
print(f"  90/10分位数比: {p90p10:.4f}")
print(f"  ジニ係数: {gini_coef:.4f}")
print(f"  最大: {df.loc[df['消費支出'].idxmax(), '都道府県']} ({income.max():,.0f}円)")
print(f"  最小: {df.loc[df['消費支出'].idxmin(), '都道府県']} ({income.min():,.0f}円)")