"""
教育用再現コード: 2023年 統計データ分析コンペティション 審査員奨励賞（大学生）
=================================================================================
論文タイトル：市区町村ごとの失業率の要因分析

【分析概要】
  データ：SSDSE-B-2026（47都道府県, 2022年）
  目的変数：求職圧力指数（Jobseeker Pressure Index）
             = 月間有効求職者数 / (月間有効求職者数 + 月間有効求人数) × 100
             → 値が高いほど求職者が求人を上回る「就職難」状態を示す
             （SSDSE-Bは市区町村レベルの失業率を直接提供しないため、
               都道府県別の労働市場データからプロキシを構築する）

  Step1. データ読み込み・変数構築（SSDSE-B 2022年）
  Step2. 記述統計・分布確認（ヒストグラム + KDE）
  Step3. 相関係数ヒートマップ
  Step4. OLS重回帰分析（標準化係数・VIF）
  Step5. 散布図（消費支出 vs 求職圧力）

【変数】
  目的変数:
    求職圧力指数 = 月間有効求職者数 / (月間有効求職者数 + 月間有効求人数) × 100

  説明変数:
    高齢化率    = 65歳以上人口 / 総人口 × 100
    男性率      = 15〜64歳人口（男）/ 15〜64歳人口 × 100
    消費支出    = 消費支出（二人以上の世帯）[円/月]
    住宅地価格  = 標準価格（平均価格）（住宅地）[円/m²]
    TFR         = 合計特殊出生率
    大学進学率  = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
    年平均気温  = 年平均気温 [℃]
    宿泊者per capita = 延べ宿泊者数 / 総人口

【手法】
  重回帰分析（OLS）、VIF（分散拡大係数）、標準化回帰係数

【図の出力】
  html/figures/2023_U5_5_fig1_dist.png   ... 求職圧力指数の分布（ヒストグラム+KDE）
  html/figures/2023_U5_5_fig2_corr.png   ... 相関係数ヒートマップ
  html/figures/2023_U5_5_fig3_coef.png   ... OLS標準化係数棒グラフ（エラーバー付き）
  html/figures/2023_U5_5_fig4_scatter.png ... 消費支出 vs 求職圧力指数 散布図

【データ出典】
  独立行政法人統計センター「SSDSE（教育用標準データセット）」
  https://www.nstac.go.jp/use/literacy/ssdse/
  ※ 合成データ・乱数（np.random等）は一切使用しない
=================================================================================
"""

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


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

from scipy import stats as scipy_stats
from scipy.stats import gaussian_kde
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

# ── パス設定 ─────────────────────────────────────────────────────────────────
FIG_DIR  = 'html/figures'
DATA_DIR = 'data/raw'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family':        'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi':         150,
    'axes.spines.top':    False,
    'axes.spines.right':  False,
})

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step1. データ読み込み・変数構築
# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 65)
print("■ データ読み込み（SSDSE-B-2026 実データのみ）")
print("=" * 65)

df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=1
)

# 47都道府県のみ（地域コード R + 5桁数字、かつ全国合計 R00000 を除く）
df_raw = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_raw['年度'] = pd.to_numeric(df_raw['年度'], errors='coerce')

# 2022年データを使用（最新の完全データ年）
YEAR = 2022
df = df_raw[df_raw['年度'] == YEAR].copy().reset_index(drop=True)
print(f"SSDSE-B {YEAR}年: {len(df)}都道府県")

# 数値変換
NUM_COLS = [
    '月間有効求職者数（一般）', '月間有効求人数（一般）',
    '総人口', '15～64歳人口', '15～64歳人口（男）', '65歳以上人口',
    '消費支出（二人以上の世帯）',
    '標準価格（平均価格）（住宅地）',
    '合計特殊出生率',
    '高等学校卒業者数', '高等学校卒業者のうち進学者数',
    '年平均気温',
    '延べ宿泊者数',
]
for c in NUM_COLS:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# ── 変数の構築 ────────────────────────────────────────────────────────────────
# 目的変数：求職圧力指数（Jobseeker Pressure Index）
df['求職圧力指数'] = (
    df['月間有効求職者数（一般）'] /
    (df['月間有効求職者数（一般）'] + df['月間有効求人数（一般）'])
    * 100
)

# 説明変数
df['高齢化率']   = df['65歳以上人口'] / df['総人口'] * 100
df['男性率']     = df['15～64歳人口（男）'] / df['15～64歳人口'] * 100
df['消費支出']   = df['消費支出（二人以上の世帯）']          # 円/月
df['住宅地価格'] = df['標準価格（平均価格）（住宅地）']        # 円/m²
df['TFR']        = df['合計特殊出生率']
df['大学進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100
df['年平均気温'] = df['年平均気温']
df['宿泊者per_capita'] = df['延べ宿泊者数'] / df['総人口']

# 分析変数名
Y_NAME = '求職圧力指数'
X_NAMES = ['高齢化率', '男性率', '消費支出', '住宅地価格', 'TFR', '大学進学率', '年平均気温', '宿泊者per_capita']

# 表示用変数名（短縮）
X_LABELS = {
    '高齢化率':        '高齢化率\n（%）',
    '男性率':          '男性率\n（生産年齢，%）',
    '消費支出':        '消費支出\n（円/月）',
    '住宅地価格':      '住宅地価格\n（円/m²）',
    'TFR':             '合計特殊\n出生率',
    '大学進学率':      '大学進学率\n（%）',
    '年平均気温':      '年平均気温\n（℃）',
    '宿泊者per_capita':'宿泊者\nper capita',
}

# 欠損除外
VARS_ALL = [Y_NAME] + X_NAMES + ['都道府県']
df_ana = df[VARS_ALL].dropna().reset_index(drop=True)
N = len(df_ana)
prefs = df_ana['都道府県'].tolist()

y = df_ana[Y_NAME].values
X = df_ana[X_NAMES].values

print(f"\n分析対象: {N}都道府県")
print(f"\n{Y_NAME} 記述統計:")
print(f"  平均   = {y.mean():.2f}%")
print(f"  標準偏差 = {y.std():.2f}%")
print(f"  最小   = {y.min():.2f}%  ({prefs[y.argmin()]})")
print(f"  最大   = {y.max():.2f}%  ({prefs[y.argmax()]})")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step2. 相関分析
# ═══════════════════════════════════════════════════════════════════════════════
df_corr = df_ana[X_NAMES + [Y_NAME]].copy()
corr_matrix = df_corr.corr()

print("\n【説明変数と求職圧力指数の相関】")
for xn in X_NAMES:
    r, p = scipy_stats.pearsonr(df_ana[xn], y)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'ns'
    print(f"  {xn:<16} r = {r:+.3f}  p = {p:.4f}  {sig}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step3. VIF 計算
# ═══════════════════════════════════════════════════════════════════════════════
X_reg = sm.add_constant(X)
vif_vals = [variance_inflation_factor(X_reg, i + 1) for i in range(len(X_NAMES))]
vif_df = pd.DataFrame({'変数': X_NAMES, 'VIF': vif_vals})

print("\n【VIF（分散拡大係数）】")
for i, (xn, v) in enumerate(zip(X_NAMES, vif_vals)):
    flag = '  ★多重共線性の疑い' if v > 5 else ''
    print(f"  {xn:<16} VIF = {v:.2f}{flag}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step4. OLS 重回帰（非標準化・標準化）
# ═══════════════════════════════════════════════════════════════════════════════
# 非標準化 OLS
ols_result = sm.OLS(y, X_reg).fit()
print("\n【OLS 重回帰結果（非標準化）】")
print(f"  R² = {ols_result.rsquared:.3f}  AdjR² = {ols_result.rsquared_adj:.3f}")
print(f"  F統計量 = {ols_result.fvalue:.2f}  p = {ols_result.f_pvalue:.4f}")
print(ols_result.summary().tables[1])

# 標準化 OLS（Zスコア変換）
y_z = (y - y.mean()) / y.std()
X_z = (X - X.mean(axis=0)) / X.std(axis=0)
X_z_const = sm.add_constant(X_z)
ols_std = sm.OLS(y_z, X_z_const).fit()

std_coefs = ols_std.params[1:]        # 定数項を除く
std_pvals = ols_std.pvalues[1:]
std_bse   = ols_std.bse[1:]           # 標準誤差（95%CI用）

print("\n【標準化回帰係数】")
for xn, sc, pv in zip(X_NAMES, std_coefs, std_pvals):
    sig = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'ns'
    print(f"  {xn:<16} β = {sc:+.3f}  p = {pv:.4f}  {sig}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ 図の生成（4枚）
# ═══════════════════════════════════════════════════════════════════════════════

# ────────────────────────────────────────────────────────────────────────────
# 図1：求職圧力指数の分布（ヒストグラム + KDE）
# ────────────────────────────────────────────────────────────────────────────
print("\n図1: 分布図（ヒストグラム + KDE）を作成中...")

fig1, ax1 = plt.subplots(figsize=(9, 5))

# ヒストグラム
n_bins = 10
counts, bin_edges, _ = ax1.hist(
    y, bins=n_bins, color='#1565C0', alpha=0.65, edgecolor='white',
    label='度数', density=True
)

# KDE（カーネル密度推定）
kde = gaussian_kde(y, bw_method='scott')
x_kde = np.linspace(y.min() - 1, y.max() + 1, 300)
ax1.plot(x_kde, kde(x_kde), color='#E65100', linewidth=2.5, label='KDE（カーネル密度）')

# 平均・中央値の縦線
ax1.axvline(y.mean(),   color='#1565C0', linestyle='--', linewidth=1.8,
            label=f'平均 = {y.mean():.2f}%')
ax1.axvline(np.median(y), color='#2E7D32', linestyle=':', linewidth=1.8,
            label=f'中央値 = {np.median(y):.2f}%')

ax1.set_xlabel('求職圧力指数（%）', fontsize=12)
ax1.set_ylabel('密度', fontsize=12)
ax1.set_title(
    f'求職圧力指数の分布 — 47都道府県（{YEAR}年）\n'
    f'（月間有効求職者数 / (求職者数 + 求人数) × 100）\nデータ：SSDSE-B-2026 実データ',
    fontsize=12, fontweight='bold'
)
ax1.legend(fontsize=10)
ax1.grid(axis='y', alpha=0.3)

# 上位・下位都道府県ラベル
highlight = ['北海道', '沖縄県', '東京都', '愛知県', '秋田県']
for pref in highlight:
    if pref in prefs:
        idx = prefs.index(pref)
        kde_val = float(kde(np.array([y[idx]]))[0])
        ax1.annotate(
            pref.replace('道','').replace('府','').replace('都','').replace('県',''),
            xy=(y[idx], 0.005),
            xytext=(y[idx], kde_val * 0.55),
            fontsize=8, color='#C62828',
            arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.0),
            ha='center'
        )

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2023_U5_5_fig1_dist.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  -> 2023_U5_5_fig1_dist.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図2：相関係数ヒートマップ
# ────────────────────────────────────────────────────────────────────────────
print("図2: 相関係数ヒートマップを作成中...")

fig2, ax2 = plt.subplots(figsize=(10, 9))
corr_vals = corr_matrix.values
n_all = len(corr_matrix.columns)
col_labels_short = [X_LABELS.get(c, c) if c != Y_NAME else '求職圧力\n指数' for c in corr_matrix.columns]

im = ax2.imshow(corr_vals, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
cbar = plt.colorbar(im, ax=ax2, fraction=0.046, pad=0.04)
cbar.set_label('Pearson相関係数', fontsize=10)

ax2.set_xticks(range(n_all))
ax2.set_yticks(range(n_all))
ax2.set_xticklabels(col_labels_short, fontsize=8, rotation=35, ha='right')
ax2.set_yticklabels(col_labels_short, fontsize=8)
ax2.set_title(
    f'Pearson相関係数行列（n={N}都道府県）\n* p<0.05  ** p<0.01  *** p<0.001\nデータ：SSDSE-B-2026',
    fontsize=12, fontweight='bold'
)

for i in range(n_all):
    for j in range(n_all):
        val = corr_vals[i, j]
        sig_str = ''
        if i != j:
            try:
                _, pv = scipy_stats.pearsonr(
                    df_corr.iloc[:, i].dropna(),
                    df_corr.iloc[:, j].dropna()
                )
                sig_str = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else ''
            except Exception:
                pass
        text_color = 'white' if abs(val) > 0.6 else 'black'
        ax2.text(j, i, f'{val:.2f}{sig_str}', ha='center', va='center',
                 fontsize=7.5, fontweight='bold', color=text_color)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_U5_5_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  -> 2023_U5_5_fig2_corr.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図3：OLS 標準化係数棒グラフ（エラーバー付き）
# ────────────────────────────────────────────────────────────────────────────
print("図3: 標準化係数棒グラフを作成中...")

# 95% CI（±1.96×SE）
ci95 = 1.96 * std_bse

# 有意性フラグで色分け
bar_colors = []
for pv in std_pvals:
    if pv < 0.01:
        bar_colors.append('#C62828')   # 深紅：高有意
    elif pv < 0.05:
        bar_colors.append('#E65100')   # オレンジ：有意
    else:
        bar_colors.append('#90CAF9')   # 薄青：非有意

short_labels = [X_LABELS.get(xn, xn) for xn in X_NAMES]

# 標準化係数の降順ソート
sort_idx = np.argsort(std_coefs)
sorted_labels  = [short_labels[i] for i in sort_idx]
sorted_coefs   = std_coefs[sort_idx]
sorted_ci      = ci95[sort_idx]
sorted_pvals   = std_pvals[sort_idx]
sorted_colors  = [bar_colors[i] for i in sort_idx]

fig3, ax3 = plt.subplots(figsize=(10, 6))

bars3 = ax3.barh(
    np.arange(len(X_NAMES)), sorted_coefs,
    color=sorted_colors, alpha=0.85, edgecolor='white',
    xerr=sorted_ci, error_kw=dict(ecolor='#333', capsize=4, linewidth=1.2)
)
ax3.axvline(0, color='black', linewidth=1.0)
ax3.set_yticks(np.arange(len(X_NAMES)))
ax3.set_yticklabels(sorted_labels, fontsize=10)
ax3.set_xlabel('標準化回帰係数（β）', fontsize=12)
ax3.set_title(
    f'OLS 重回帰：標準化係数（β）と95%信頼区間\n'
    f'目的変数：求職圧力指数　R²={ols_result.rsquared:.3f}  AdjR²={ols_result.rsquared_adj:.3f}  n={N}',
    fontsize=12, fontweight='bold'
)
ax3.grid(axis='x', alpha=0.3)
ax3.invert_yaxis()

# 有意性ラベル
for i, (bar, pv) in enumerate(zip(bars3, sorted_pvals)):
    w = bar.get_width()
    sig_str = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'ns'
    offset = sorted_ci[i] + 0.02
    ax3.text(w + offset if w >= 0 else w - offset - 0.05,
             bar.get_y() + bar.get_height() / 2,
             sig_str, va='center', fontsize=10,
             color='#C62828' if sig_str != 'ns' else '#999')

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#C62828', alpha=0.85, label='p < 0.01（高度有意）'),
    Patch(facecolor='#E65100', alpha=0.85, label='p < 0.05（有意）'),
    Patch(facecolor='#90CAF9', alpha=0.85, label='p ≥ 0.05（非有意）'),
]
ax3.legend(handles=legend_elements, fontsize=9, loc='lower right')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2023_U5_5_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  -> 2023_U5_5_fig3_coef.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図4：消費支出 vs 求職圧力指数 散布図（回帰直線 + 都道府県ラベル）
# ────────────────────────────────────────────────────────────────────────────
print("図4: 散布図（消費支出 vs 求職圧力指数）を作成中...")

x_scatter = df_ana['消費支出'].values
y_scatter = y

r_s, p_s = scipy_stats.pearsonr(x_scatter, y_scatter)
coef_fit  = np.polyfit(x_scatter, y_scatter, 1)
x_fit     = np.linspace(x_scatter.min(), x_scatter.max(), 200)
y_fit     = np.polyval(coef_fit, x_fit)

fig4, ax4 = plt.subplots(figsize=(10, 7))

scatter = ax4.scatter(
    x_scatter / 1000, y_scatter,    # 千円単位に変換
    c='#1565C0', alpha=0.75, s=70,
    edgecolors='white', linewidth=0.6, zorder=3
)
ax4.plot(x_fit / 1000, y_fit, color='#E65100', linewidth=2.2,
         linestyle='--', label=f'回帰直線 (r={r_s:.3f})', zorder=2)

# 都道府県ラベル（代表的な都道府県）
LABEL_PREFS = ['北海道', '東京都', '大阪府', '愛知県', '沖縄県',
               '秋田県', '青森県', '鹿児島県', '福岡県', '神奈川県']
for i, pref in enumerate(prefs):
    if pref in LABEL_PREFS:
        short = pref.replace('都','').replace('道','').replace('府','').replace('県','')
        ax4.annotate(
            short, (x_scatter[i] / 1000, y_scatter[i]),
            textcoords='offset points', xytext=(5, 3),
            fontsize=8.5, color='#333',
            bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor='none')
        )

ax4.set_xlabel('消費支出（二人以上の世帯）[千円/月]', fontsize=12)
ax4.set_ylabel('求職圧力指数（%）', fontsize=12)
ax4.set_title(
    f'消費支出と求職圧力指数の関係（都道府県別, {YEAR}年）\n'
    f'r = {r_s:.3f}  {"p < 0.05" if p_s < 0.05 else f"p = {p_s:.3f}"}  '
    f'回帰係数 = {coef_fit[0]*1000:.4f} (%/千円)',
    fontsize=12, fontweight='bold'
)
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.25)

# テキストボックス：解釈
ax4.text(
    0.03, 0.95,
    f'r = {r_s:.3f}\n{"p < 0.05" if p_s < 0.05 else "p = {:.3f}".format(p_s)}\n'
    f'{"消費支出が高いほど\n求職圧力が低い傾向" if r_s < 0 else "消費支出が高いほど\n求職圧力が高い傾向"}',
    transform=ax4.transAxes, fontsize=10, va='top',
    bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.85)
)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_U5_5_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  -> 2023_U5_5_fig4_scatter.png 保存完了")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ 最終サマリー
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 65)
print("✓ 全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {os.path.abspath(FIG_DIR)}")
print("  2023_U5_5_fig1_dist.png    - 求職圧力指数の分布（ヒストグラム+KDE）")
print("  2023_U5_5_fig2_corr.png    - 相関係数ヒートマップ")
print("  2023_U5_5_fig3_coef.png    - OLS標準化係数棒グラフ（エラーバー付き）")
print("  2023_U5_5_fig4_scatter.png - 消費支出 vs 求職圧力指数散布図")
print()
print("【主要知見】")
print(f"  分析対象: {N}都道府県（SSDSE-B {YEAR}年, 実データ）")
print(f"  求職圧力指数: 平均 {y.mean():.2f}%  SD={y.std():.2f}%")
print(f"  OLS  R² = {ols_result.rsquared:.3f}  AdjR² = {ols_result.rsquared_adj:.3f}")
print(f"  VIF 最大値: {max(vif_vals):.2f}  最小値: {min(vif_vals):.2f}")
sig_vars = [X_NAMES[i] for i, p in enumerate(std_pvals) if p < 0.05]
print(f"  有意な説明変数（p<0.05）: {sig_vars if sig_vars else '—（なし）'}")
print(f"  r（消費支出 vs 求職圧力）: {r_s:.3f}  p={p_s:.4f}")
print()
print("  ※ 合成データ・乱数は一切使用していません")
print("  ※ データ出典: 統計数理研究所 SSDSE-B-2026")