"""
2018_U3_suri.py
都道府県間経済格差の収束仮説検証：β収束とσ収束の実証分析
統計数理賞（大学生・一般の部）2018年度
データ：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/2018_U3_suri.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)

print("columns:", df_b.columns.tolist())
print("年度range:", df_b['年度'].min(), '-', df_b['年度'].max())
print("都道府県数:", df_b['都道府県'].nunique())

# 消費支出列を特定
cons_col = '消費支出（二人以上の世帯）'
print(f"使用列: {cons_col}")
print(df_b[['都道府県', '年度', cons_col]].head())

# ─── 地域区分 ────────────────────────────────────────────────
# 都道府県列には「北海道」「青森県」等の形式が使われているため接頭語で検索
region_map_raw = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部',
    '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国',
    '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国',
    '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄',
    '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1',
    '関東':         '#e05c5c',
    '中部':         '#f0a500',
    '近畿':         '#5cb85c',
    '中国・四国':   '#9b59b6',
    '九州・沖縄':   '#f39c12'
}

df_b['地域'] = df_b['都道府県'].map(region_map_raw)

# ─── β収束分析（絶対収束 / 絶対β収束） ─────────────────────
start_yr = df_b['年度'].min()   # 2012
end_yr   = df_b['年度'].max()   # 2023
T = end_yr - start_yr           # 11年

start_df = (df_b[df_b['年度'] == start_yr]
            [['都道府県', cons_col]]
            .rename(columns={cons_col: 'start'}))
end_df   = (df_b[df_b['年度'] == end_yr]
            [['都道府県', cons_col]]
            .rename(columns={cons_col: 'end'}))

conv = start_df.merge(end_df, on='都道府県').dropna()
conv['成長率'] = (conv['end'] - conv['start']) / (conv['start'] * T)   # 年率成長率
conv['log初期水準'] = np.log(conv['start'])
conv['地域'] = conv['都道府県'].map(region_map_raw)

X_beta = sm.add_constant(conv['log初期水準'].astype(float))
beta_res = sm.OLS(conv['成長率'].astype(float), X_beta).fit()
beta_coef = beta_res.params['log初期水準']
beta_pval = beta_res.pvalues['log初期水準']
beta_r2   = beta_res.rsquared
beta_n    = int(beta_res.nobs)

print(f"\n=== β収束 ===")
print(f"β係数: {beta_coef:.4f}")
print(f"p値:   {beta_pval:.4f}")
print(f"R²:    {beta_r2:.4f}")
print(f"観測数: {beta_n}")
print(f"収束の判定: {'収束（β<0）' if beta_coef < 0 else '発散（β>0）'}")

# ─── σ収束分析（変動係数の時系列） ───────────────────────────
sigma = (df_b.groupby('年度')[cons_col]
         .apply(lambda x: x.std() / x.mean())
         .reset_index())
sigma.columns = ['年度', 'CV']
sigma = sigma.sort_values('年度')

trend_X = sm.add_constant(sigma['年度'].astype(float))
trend_res = sm.OLS(sigma['CV'].astype(float), trend_X).fit()
trend_coef = trend_res.params['年度']
trend_pval = trend_res.pvalues['年度']

print(f"\n=== σ収束（変動係数） ===")
print(sigma.to_string(index=False))
print(f"トレンド係数: {trend_coef:.6f}")
print(f"p値:          {trend_pval:.4f}")
print(f"収束の判定: {'σ収束（負のトレンド）' if trend_coef < 0 else 'σ発散（正のトレンド）'}")

# ─── 記述統計 ─────────────────────────────────────────────────
print(f"\n=== 消費支出 記述統計 ({start_yr}年・{end_yr}年) ===")
for yr in [start_yr, end_yr]:
    s = df_b[df_b['年度'] == yr][cons_col]
    print(f"{yr}年: 平均={s.mean():.0f}, 標準偏差={s.std():.0f}, "
          f"最小={s.min():.0f}, 最大={s.max():.0f}, CV={s.std()/s.mean():.4f}")

# ═══════════════════════════════════════════════════════════════
# 図1: β収束散布図
# ═══════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(8, 5.5))

for region, grp in conv.groupby('地域'):
    color = region_colors.get(region, '#888888')
    ax.scatter(grp['log初期水準'], grp['成長率'] * 100,
               color=color, s=60, alpha=0.85, label=region, zorder=3)

# 回帰直線
x_line = np.linspace(conv['log初期水準'].min() - 0.05,
                     conv['log初期水準'].max() + 0.05, 200)
y_line = beta_res.params['const'] + beta_coef * x_line
ax.plot(x_line, y_line * 100, color='#1565C0', linewidth=2.0,
        linestyle='--', label='OLS回帰直線', zorder=4)

# 都道府県名ラベル（上位・下位）
extremes = pd.concat([conv.nsmallest(3, '成長率'), conv.nlargest(3, '成長率')])
for _, row in extremes.iterrows():
    label = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax.annotate(label,
                xy=(row['log初期水準'], row['成長率'] * 100),
                xytext=(4, 2), textcoords='offset points', fontsize=8, color='#333333')

# β係数・p値テキスト
sig_str = '***' if beta_pval < 0.001 else ('**' if beta_pval < 0.01 else ('*' if beta_pval < 0.05 else 'n.s.'))
conv_str = '収束' if beta_coef < 0 else '発散'
ax.text(0.05, 0.95,
        f'β = {beta_coef:.4f} {sig_str}\np = {beta_pval:.4f}\nR² = {beta_r2:.3f}\n→ {conv_str}',
        transform=ax.transAxes, fontsize=10, va='top',
        bbox=dict(boxstyle='round,pad=0.5', facecolor='#EFF3FF', alpha=0.9))

ax.set_xlabel(f'log(消費支出{start_yr}年)', fontsize=12)
ax.set_ylabel(f'消費支出の年率成長率 ({start_yr}→{end_yr}年) [%]', fontsize=12)
ax.set_title(f'β収束分析：初期水準と成長率の関係\n（{start_yr}〜{end_yr}年、47都道府県）', fontsize=13, fontweight='bold')
ax.legend(loc='upper right', fontsize=9, framealpha=0.9)
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.8, linestyle=':')

plt.tight_layout()
out_path = os.path.join(FIG_DIR, '2018_U3_fig1.png')
plt.savefig(out_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"\n保存: {out_path}")

# ═══════════════════════════════════════════════════════════════
# 図2: σ収束（変動係数の時系列）
# ═══════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(8, 5))

years_arr = sigma['年度'].values
cv_arr    = sigma['CV'].values

# COVID帯（2019〜2021）グレー塗り
ax.axvspan(2019.5, 2021.5, color='#BDBDBD', alpha=0.3, label='COVID-19 影響期')

# CV折れ線
ax.plot(years_arr, cv_arr, color='#1565C0', linewidth=2.2,
        marker='o', markersize=7, label='変動係数（CV）', zorder=3)

# OLSトレンド線
y_trend = trend_res.params['const'] + trend_coef * years_arr
ax.plot(years_arr, y_trend, color='#C62828', linewidth=1.8,
        linestyle='--', label=f'OLSトレンド (β={trend_coef:.5f})', zorder=4)

# 各年の値をアノテーション
for yr, cv in zip(years_arr, cv_arr):
    ax.annotate(f'{cv:.4f}', xy=(yr, cv), xytext=(0, 8),
                textcoords='offset points', fontsize=7.5,
                ha='center', color='#1565C0')

# 判定テキスト
sig_str2 = '***' if trend_pval < 0.001 else ('**' if trend_pval < 0.01 else ('*' if trend_pval < 0.05 else 'n.s.'))
sigma_str = 'σ収束（格差縮小）' if trend_coef < 0 else 'σ発散（格差拡大）'
ax.text(0.05, 0.10,
        f'トレンド係数 = {trend_coef:.5f} {sig_str2}\np = {trend_pval:.4f}\n→ {sigma_str}',
        transform=ax.transAxes, fontsize=10, va='bottom',
        bbox=dict(boxstyle='round,pad=0.5', facecolor='#FFF9C4', alpha=0.9))

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('変動係数（CV = 標準偏差 / 平均）', fontsize=12)
ax.set_title(f'σ収束分析：都道府県間消費格差の変動係数推移\n（{start_yr}〜{end_yr}年、47都道府県）',
             fontsize=13, fontweight='bold')
ax.legend(loc='upper right', fontsize=9, framealpha=0.9)
ax.grid(True, alpha=0.3)
ax.set_xticks(years_arr)

plt.tight_layout()
out_path = os.path.join(FIG_DIR, '2018_U3_fig2.png')
plt.savefig(out_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"保存: {out_path}")

# ═══════════════════════════════════════════════════════════════
# 図3: 都道府県別消費支出 2012年 vs 最新年の散布図（45度線付き）
# ═══════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(7.5, 6.5))

for region, grp in conv.groupby('地域'):
    color = region_colors.get(region, '#888888')
    ax.scatter(grp['start'] / 1000, grp['end'] / 1000,
               color=color, s=60, alpha=0.85, label=region, zorder=3)

# 45度線
all_vals = pd.concat([conv['start'], conv['end']]) / 1000
vmin, vmax = all_vals.min() * 0.97, all_vals.max() * 1.03
ax.plot([vmin, vmax], [vmin, vmax], color='gray', linewidth=1.5,
        linestyle='--', label='45度線（変化なし）', zorder=2)

# 都道府県名ラベル（上位・下位 消費水準）
top_prefs = conv.nlargest(4, 'end')
bot_prefs = conv.nsmallest(4, 'end')
for _, row in pd.concat([top_prefs, bot_prefs]).iterrows():
    label = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax.annotate(label,
                xy=(row['start'] / 1000, row['end'] / 1000),
                xytext=(4, 2), textcoords='offset points', fontsize=7.5, color='#444')

# 45度線より上 = 改善（消費増加）の注釈
ax.text(0.05, 0.95, '▲ 45度線より上 = 消費増加（成長）',
        transform=ax.transAxes, fontsize=9, va='top', color='#2E7D32',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#E8F5E9', alpha=0.8))
ax.text(0.05, 0.87, '▼ 45度線より下 = 消費減少（後退）',
        transform=ax.transAxes, fontsize=9, va='top', color='#C62828',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#FFEBEE', alpha=0.8))

ax.set_xlabel(f'消費支出 {start_yr}年 [千円/月]', fontsize=12)
ax.set_ylabel(f'消費支出 {end_yr}年 [千円/月]', fontsize=12)
ax.set_title(f'都道府県別消費支出：{start_yr}年 vs {end_yr}年\n（地域色分け、45度線付き）',
             fontsize=13, fontweight='bold')
ax.legend(loc='lower right', fontsize=8.5, framealpha=0.9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
out_path = os.path.join(FIG_DIR, '2018_U3_fig3.png')
plt.savefig(out_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"保存: {out_path}")

# ═══════════════════════════════════════════════════════════════
# 図4: 地域別消費支出の時系列推移（6地域、折れ線）
# ═══════════════════════════════════════════════════════════════
df_b['地域'] = df_b['都道府県'].map(region_map_raw)
region_ts = (df_b.dropna(subset=['地域'])
             .groupby(['年度', '地域'])[cons_col]
             .mean()
             .reset_index())
region_ts.columns = ['年度', '地域', '平均消費支出']

fig, ax = plt.subplots(figsize=(9, 5.5))

region_order = ['関東', '近畿', '中部', '北海道・東北', '中国・四国', '九州・沖縄']
for region in region_order:
    sub = region_ts[region_ts['地域'] == region].sort_values('年度')
    if sub.empty:
        continue
    color = region_colors[region]
    ax.plot(sub['年度'], sub['平均消費支出'] / 1000,
            color=color, linewidth=2.2, marker='o', markersize=5,
            label=region, zorder=3)
    # 最終年ラベル
    last = sub.iloc[-1]
    ax.annotate(region,
                xy=(last['年度'], last['平均消費支出'] / 1000),
                xytext=(4, 0), textcoords='offset points',
                fontsize=8, color=color, va='center', fontweight='bold')

# COVID帯
ax.axvspan(2019.5, 2021.5, color='#BDBDBD', alpha=0.25, label='COVID-19 影響期')

# 全国平均
nat_ts = df_b.groupby('年度')[cons_col].mean().reset_index()
ax.plot(nat_ts['年度'], nat_ts[cons_col] / 1000,
        color='black', linewidth=1.5, linestyle=':', marker='s', markersize=4,
        label='全国平均', zorder=4)

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('消費支出（月平均）[千円]', fontsize=12)
ax.set_title('地域別消費支出の時系列推移\n（6地域別平均、全国平均を含む）',
             fontsize=13, fontweight='bold')
ax.legend(loc='upper left', fontsize=8.5, framealpha=0.9,
          bbox_to_anchor=(0.0, 1.0))
ax.grid(True, alpha=0.3)
all_years = sorted(df_b['年度'].unique())
ax.set_xticks(all_years)
ax.set_xticklabels([str(y) for y in all_years], rotation=45, fontsize=9)

plt.tight_layout()
out_path = os.path.join(FIG_DIR, '2018_U3_fig4.png')
plt.savefig(out_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"保存: {out_path}")

# ─── 図の確認 ─────────────────────────────────────────────────
for i in range(1, 5):
    path = os.path.join(FIG_DIR, f'2018_U3_fig{i}.png')
    size = os.path.getsize(path) / 1024
    print(f"  fig{i}: {size:.1f} KB  {'OK' if size > 10 else 'WARN'}")

print("\nDONE: 2018_U3_suri")