"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞（大学生・一般の部）
=================================================================
論文タイトル：日本の電力需要量の影響調査
受賞区分：審査員奨励賞（大学生・一般の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県パネルデータ、2012〜2023年）
  都道府県横断分析：2022年度 N=47都道府県（OLS重回帰）
  時系列分析：2012〜2023年 全国平均の光熱費割合の推移

  代理変数：光熱・水道費（二人以上の世帯）/ 消費支出 × 100
            → 家計支出に占める光熱費の割合を電力需要の代理指標として使用

  Step1. 光熱費割合の時系列トレンド（2012〜2023全国平均）
  Step2. 相関ヒートマップ（光熱費と各説明変数）
  Step3. OLS重回帰（標準化係数）
  Step4. 年平均気温 vs 光熱費割合の散布図

【説明変数】
  ・年平均気温   ：冷暖房需要の代理（低温→暖房需要増）
  ・総人口       ：電力需要の規模要因
  ・消費支出     ：生活水準・所得の代理
  ・高齢化率     ：65歳以上人口比率（在宅時間の増加）
  ・住宅地価格   ：地域の経済・都市度の代理

【手法】
  ・重回帰分析（OLS）
  ・相関分析（ピアソン相関係数）
  ・VIF（分散拡大要因）による多重共線性診断

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ（2012〜2023年度）

【データサイエンス学習ポイント】
  1. 代理変数（proxy variable）の設計と解釈
  2. OLSの標準化係数による説明変数の相対的重要度比較
  3. VIFによる多重共線性の診断
  4. 気温と光熱費の非線形関係（暑くても寒くてもエネルギー消費増）
=================================================================
"""

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


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
from statsmodels.stats.outliers_influence import variance_inflation_factor
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)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026: 都道府県パネルデータ）
# ================================================================
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()

# 数値変換
num_cols = [
    '光熱・水道費（二人以上の世帯）',
    '消費支出（二人以上の世帯）',
    '年平均気温',
    '最高気温（日最高気温の月平均の最高値）',
    '最低気温（日最低気温の月平均の最低値）',
    '総人口',
    '65歳以上人口',
    '標準価格（平均価格）（住宅地）',
]
for c in num_cols:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# 派生変数
df_b['光熱費割合'] = df_b['光熱・水道費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100
df_b['高齢化率']   = df_b['65歳以上人口'] / df_b['総人口'] * 100

print("=" * 60)
print("■ SSDSE-B 読み込み完了")
print(f"  全レコード数: {len(df_b)}（都道府県×年度）")
print(f"  年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}")
print("=" * 60)

# ================================================================
# ■ 2022年度 横断データ（N=47 都道府県）
# ================================================================
df_2022 = df_b[df_b['年度'] == 2022].copy()

# 説明変数リスト
XVARS = {
    '年平均気温（℃）'    : '年平均気温',
    '総人口（万人）'      : '総人口',
    '消費支出（円）'      : '消費支出（二人以上の世帯）',
    '高齢化率（%）'       : '高齢化率',
    '住宅地価格（円/m²）' : '標準価格（平均価格）（住宅地）',
}
# 目的変数
YVAR = '光熱費割合'

# 分析用データフレーム作成
cols_use = [YVAR, '年平均気温', '総人口', '消費支出（二人以上の世帯）',
            '高齢化率', '標準価格（平均価格）（住宅地）', '都道府県']
df_ols = df_2022[cols_use].dropna().copy().reset_index(drop=True)
N = len(df_ols)

# 変数名マッピング（表示用）
VARNAMES_JP = list(XVARS.keys())
VARNAMES_EN = list(XVARS.values())

print(f"\n■ 2022年度 有効データ: N={N}都道府県")
print(df_ols[[YVAR] + VARNAMES_EN].describe().round(2))

# ── OLS 標準化係数 ─────────────────────────────────────────────
X_raw = df_ols[VARNAMES_EN].values.astype(float)
y_raw = df_ols[YVAR].values.astype(float)

# 標準化
X_std = (X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0, ddof=1)
y_std = (y_raw - y_raw.mean()) / y_raw.std(ddof=1)

model_std = sm.OLS(y_std, sm.add_constant(X_std)).fit()
model_raw = sm.OLS(y_raw, sm.add_constant(X_raw)).fit()

print("\n" + "=" * 60)
print("■ OLS 重回帰（標準化）")
print("=" * 60)
print(f"  N={N}, R²={model_std.rsquared:.4f}, Adj.R²={model_std.rsquared_adj:.4f}")
print(f"\n  {'変数':<22} {'標準化β':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 50)
for name, b, p in zip(VARNAMES_JP,
                      model_std.params[1:],
                      model_std.pvalues[1:]):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<22} {b:>8.4f} {p:>10.4f} {sig:>6}")

# ── VIF 計算 ────────────────────────────────────────────────
print("\n" + "=" * 60)
print("■ VIF（多重共線性診断）")
print("=" * 60)
X_vif = sm.add_constant(X_raw)
vifs = [variance_inflation_factor(X_vif, i + 1) for i in range(X_raw.shape[1])]
print(f"\n  {'変数':<22} {'VIF':>8}")
print("  " + "-" * 32)
for name, v in zip(VARNAMES_JP, vifs):
    flag = ' ← 注意(>5)' if v > 5 else ''
    print(f"  {name:<22} {v:>8.2f}{flag}")

# ── 相関行列 ────────────────────────────────────────────────
print("\n" + "=" * 60)
print("■ 相関行列（ピアソン）")
print("=" * 60)
corr_cols = [YVAR] + VARNAMES_EN
corr_matrix = df_ols[corr_cols].corr()
print(corr_matrix.round(3))

# ================================================================
# ■ 時系列集計（2012〜2023全国平均）
# ================================================================
ts = df_b.groupby('年度').agg(
    光熱費割合_mean=('光熱費割合', 'mean'),
    光熱費割合_std =('光熱費割合', 'std'),
).reset_index()

print("\n■ 光熱費割合 全国平均（年度別）")
print(ts.round(3))

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

# ────────────────────────────────────────────────────────────────
# 図1: 光熱費割合の時系列（全国平均 2012-2023）
# ────────────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(9, 5))

years  = ts['年度'].values
means  = ts['光熱費割合_mean'].values
stds   = ts['光熱費割合_std'].values

ax1.fill_between(years, means - stds, means + stds,
                 alpha=0.18, color='#1565C0', label='±1SD（都道府県間ばらつき）')
ax1.plot(years, means, 'o-', color='#1565C0', linewidth=2.2,
         markersize=7, markerfacecolor='white', markeredgewidth=2, label='全国平均')

# 2020年（コロナ禍）にアノテーション
if 2020 in years:
    idx20 = list(years).index(2020)
    ax1.annotate('2020年\n（コロナ禍）',
                 xy=(2020, means[idx20]),
                 xytext=(2020 - 1.0, means[idx20] + 0.4),
                 arrowprops=dict(arrowstyle='->', color='#C62828'),
                 fontsize=9, color='#C62828')

# 2022年（エネルギー価格高騰）
if 2022 in years:
    idx22 = list(years).index(2022)
    ax1.annotate('2022年\n（エネルギー高騰）',
                 xy=(2022, means[idx22]),
                 xytext=(2021.0, means[idx22] + 0.5),
                 arrowprops=dict(arrowstyle='->', color='#E65100'),
                 fontsize=9, color='#E65100')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('光熱費割合（光熱・水道費 / 消費支出 × 100, %）', fontsize=11)
ax1.set_title('光熱・水道費割合の時系列推移（全国平均）\n2012〜2023年 都道府県データ（SSDSE-B）',
              fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xticks(years)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig1_ts.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2022_U5_5_fig1_ts.png")

# ────────────────────────────────────────────────────────────────
# 図2: 相関ヒートマップ（光熱費と説明変数）
# ────────────────────────────────────────────────────────────────
corr_plot_cols = [YVAR] + VARNAMES_EN
corr_labels = ['光熱費割合'] + VARNAMES_JP
C = df_ols[corr_plot_cols].corr().values

fig2, ax2 = plt.subplots(figsize=(8, 7))
n = len(corr_labels)
im = ax2.imshow(C, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax2, shrink=0.8, label='ピアソン相関係数')

ax2.set_xticks(range(n))
ax2.set_xticklabels(corr_labels, rotation=30, ha='right', fontsize=10)
ax2.set_yticks(range(n))
ax2.set_yticklabels(corr_labels, fontsize=10)

# セル内に相関係数を表示
for i in range(n):
    for j in range(n):
        val = C[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax2.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=9, color=color, fontweight='bold')

ax2.set_title('相関ヒートマップ：光熱費割合と各説明変数\n（2022年度, N=47都道府県）',
              fontsize=12, fontweight='bold')
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_U5_5_fig2_corr.png")

# ────────────────────────────────────────────────────────────────
# 図3: OLS 標準化係数プロット
# ────────────────────────────────────────────────────────────────
betas  = np.asarray(model_std.params[1:])
ses    = np.asarray(model_std.bse[1:])
pvals  = np.asarray(model_std.pvalues[1:])

fig3, ax3 = plt.subplots(figsize=(9, 5))
y_pos = np.arange(len(VARNAMES_JP))
colors3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E' for p in pvals]

ax3.barh(y_pos, betas, color=colors3, alpha=0.78, edgecolor='white', height=0.55)
ax3.errorbar(betas, y_pos, xerr=1.96 * ses, fmt='none',
             color='#333', capsize=4, linewidth=1.3)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(VARNAMES_JP, fontsize=11)
ax3.set_xlabel('標準化偏回帰係数（±1.96SE）', fontsize=11)
ax3.set_title(
    f'OLS重回帰：標準化偏回帰係数\n'
    f'R²={model_std.rsquared:.3f}, Adj.R²={model_std.rsquared_adj:.3f}, N={N}都道府県',
    fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
ax3.legend(handles=[
    Patch(color='#C62828', alpha=0.78, label='p < 0.01'),
    Patch(color='#FF8F00', alpha=0.78, label='p < 0.05'),
    Patch(color='#9E9E9E', alpha=0.78, label='n.s.')
], fontsize=9, loc='lower right')

# 係数値をバーに付記
for i, (b, p) in enumerate(zip(betas, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    ha = 'left' if b >= 0 else 'right'
    ax3.text(b + (0.01 if b >= 0 else -0.01), i,
             f' {b:.3f}{sig}', va='center', ha=ha, fontsize=8.5)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_U5_5_fig3_coef.png")

# ────────────────────────────────────────────────────────────────
# 図4: 年平均気温 vs 光熱費割合 散布図
# ────────────────────────────────────────────────────────────────
x4 = df_ols['年平均気温'].values.astype(float)
y4 = df_ols[YVAR].values.astype(float)
pref = df_ols['都道府県'].values

fig4, ax4 = plt.subplots(figsize=(9, 6))

# 散布図
scatter = ax4.scatter(x4, y4, c=y4, cmap='RdYlBu_r', s=80,
                      edgecolors='gray', linewidths=0.5, alpha=0.85, zorder=3)
plt.colorbar(scatter, ax=ax4, label='光熱費割合（%）', shrink=0.8)

# 回帰直線
z4 = np.polyfit(x4, y4, 1)
xs4 = np.linspace(x4.min() - 0.5, x4.max() + 0.5, 100)
ax4.plot(xs4, np.poly1d(z4)(xs4), 'r--', linewidth=1.8, alpha=0.8, label='線形回帰')
r4, p4 = stats.pearsonr(x4, y4)

# 主要都道府県にラベル
label_prefs = ['北海道', '沖縄県', '東京都', '大阪府', '秋田県', '鹿児島県', '長野県']
for i, pname in enumerate(pref):
    if any(lp in pname for lp in label_prefs):
        ax4.annotate(pname, (x4[i], y4[i]),
                     textcoords='offset points', xytext=(6, 3),
                     fontsize=8.5, color='#333')

ax4.set_xlabel('年平均気温（℃）', fontsize=12)
ax4.set_ylabel('光熱費割合（光熱・水道費 / 消費支出 × 100, %）', fontsize=11)
ax4.set_title(
    f'年平均気温と光熱費割合の関係（2022年度, N={N}都道府県）\n'
    f'r = {r4:.3f}（p = {p4:.4f}）',
    fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2022_U5_5_fig4_scatter.png")

print("\n" + "=" * 60)
print("全図の生成完了（4枚）")
print("  fig1_ts.png      : 光熱費割合の時系列（全国平均 2012-2023）")
print("  fig2_corr.png    : 相関ヒートマップ（光熱費と説明変数）")
print("  fig3_coef.png    : OLS標準化係数プロット")
print("  fig4_scatter.png : 年平均気温 vs 光熱費割合 散布図")
print("=" * 60)