"""
教育用再現コード: 2022年 統計データ分析コンペティション 統計数理賞 [高校生の部]
=================================================================
論文タイトル：投票率の地域差：都道府県別社会経済要因の回帰分析
受賞：統計数理賞（高校生の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別データ, 2012〜2023年度）

  ※ 都道府県別投票率の直接データはSSDSE-Bに含まれないため、
    政治参加・市民活動の代理指標として以下を分析：
    - 婚姻率（社会的絆・地域帰属意識の代理）
    - 旅券発行件数（市民活動・社会参加意識の代理）
    - 大学進学率（教育水準→政治関心の代理）

  分析の流れ
  1. 時系列：旅券発行件数（市民活動指標）の推移
  2. 相関分析：社会参加代理指標の相関ヒートマップ
  3. OLS回帰：旅券発行率の決定要因（高齢化率・教育・所得）
  4. 地域別ランキング：旅券発行率（2022年）

【代理変数設計】
  旅券発行率 = 一般旅券発行件数 / 総人口 × 1000（社会活動・市民参加の代理）
  婚姻率 = 婚姻件数 / 総人口 × 1000
  大学進学率 = 高校卒業者のうち進学者数 / 高校卒業者数 × 100
  高齢化率 = 65歳以上人口 / 総人口 × 100
  消費支出_log = log(消費支出（二人以上の世帯）)

【データ出典】
  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_H3_suri.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
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)

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['総人口'] * 1000
df_b['婚姻率'] = df_b['婚姻件数'] / df_b['総人口'] * 1000
df_b['大学進学率'] = df_b['高等学校卒業者のうち進学者数'] / df_b['高等学校卒業者数'].replace(0, np.nan) * 100
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['住宅地価_log'] = np.log(df_b['標準価格（平均価格）（住宅地）'].clip(lower=1))

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

# Fig1: 旅券発行率の時系列推移（地域別平均）
fig, ax = plt.subplots(figsize=(10, 5))
yearly = df_b.groupby(['年度', '地域'])['旅券発行率'].mean().reset_index()
for reg, grp in yearly.groupby('地域'):
    ax.plot(grp['年度'], grp['旅券発行率'], marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('旅券発行率（件/千人）', fontsize=12)
ax.set_title('地域別 旅券発行率の推移（2012〜2023年）\n（社会参加・市民活動の代理指標）', fontsize=13, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H3_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: 相関ヒートマップ
df_2022 = df_b[df_b['年度'] == 2022].copy()
analysis_vars = ['旅券発行率', '婚姻率', '大学進学率', '高齢化率', '消費支出_log', '住宅地価_log']
df_corr = df_2022[analysis_vars].dropna()
corr = df_corr.corr()
fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xticks(range(len(analysis_vars)))
ax.set_yticks(range(len(analysis_vars)))
ax.set_xticklabels(analysis_vars, rotation=45, ha='right', fontsize=9)
ax.set_yticklabels(analysis_vars, fontsize=9)
for i in range(len(analysis_vars)):
    for j in range(len(analysis_vars)):
        val = corr.iloc[i, j]
        ax.text(j, i, f'{val:.2f}', ha='center', va='center',
                fontsize=8, color='white' if abs(val) > 0.5 else 'black')
plt.colorbar(im, ax=ax, label='Pearson相関係数')
ax.set_title('社会参加指標の相関ヒートマップ（2022年, N=47）', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H3_fig2_corr.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: OLS回帰 旅券発行率の決定要因
xvars = ['高齢化率', '大学進学率', '消費支出_log', '婚姻率']
df_reg = df_2022[['旅券発行率'] + xvars].dropna()
X = sm.add_constant(df_reg[xvars])
res = sm.OLS(df_reg['旅券発行率'], X).fit()
print(res.summary())
coefs = res.params.drop('const')
ses = res.bse.drop('const')
pvals = res.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('OLS回帰係数', fontsize=12)
ax.set_title(f'旅券発行率の決定要因（N=47）\nR²={res.rsquared:.3f}（赤=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_H3_fig3_ols.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: 都道府県別旅券発行率ランキング（2022年）
df_sorted = df_2022.sort_values('旅券発行率', ascending=True).dropna(subset=['旅券発行率'])
fig, ax = plt.subplots(figsize=(8, 12))
bar_colors = [region_colors.get(df_sorted[df_sorted['都道府県'] == p]['地域'].values[0], 'gray')
              if len(df_sorted[df_sorted['都道府県'] == p]) > 0 else 'gray'
              for p in df_sorted['都道府県']]
ax.barh(df_sorted['都道府県'], df_sorted['旅券発行率'], color=bar_colors, alpha=0.8)
ax.axvline(df_2022['旅券発行率'].mean(), color='red', linestyle='--', linewidth=1.5,
           label=f'全国平均: {df_2022["旅券発行率"].mean():.2f}')
ax.set_xlabel('旅券発行率（件/千人）', fontsize=12)
ax.set_title('都道府県別 旅券発行率ランキング（2022年）', fontsize=13, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H3_fig4_rank.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")