"""
2022_H5_8_shorei.py
===================
教育用再現コード：気温と光熱・水道費
── 地域気候が家計エネルギー支出に与える影響 ──

2022年度 統計データ分析コンペティション 審査員奨励賞（高校生の部）

分析手法:
  - 相関分析（Pearson）
  - OLS 重回帰（statsmodels）
  - 時系列分析（寒冷地 vs 温暖地の推移比較）
  - 散布図（都道府県ラベル付き）

データ: SSDSE-B-2026.csv（47都道府県・2012〜2023年）
  出典: 統計数理研究所 SSDSE（社会・人口統計体系）

派生変数:
  - 光熱水道費_万円  = 光熱・水道費（二人以上の世帯）/ 10000
  - 光熱水道費割合   = 光熱・水道費 / 消費支出（二人以上の世帯）× 100
  - 高齢化率         = 65歳以上人口 / 総人口 × 100
  - 消費支出_log     = log(消費支出（二人以上の世帯）)

寒冷地（北海道・東北6県）は暖房費が大きく、年平均気温と光熱水道費割合に
負の相関が生じることを可視化・定量化する。

注意: 本スクリプトは実データのみを使用し、合成データ（np.random.seed等）は一切使用しない。
"""

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


import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib import rcParams
import statsmodels.api as sm
from scipy import stats

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────────────
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'

os.makedirs(FIG_DIR, exist_ok=True)

# ── フォント設定 ───────────────────────────────────────────────────────────────
rcParams['font.family'] = ['Hiragino Sans', 'Hiragino Kaku Gothic ProN',
                            'AppleGothic', 'Noto Sans CJK JP', 'sans-serif']
rcParams['axes.unicode_minus'] = False

# ── データ読み込み ─────────────────────────────────────────────────────────────
df_raw = pd.read_csv(DATA_B, header=[0, 1], encoding='cp932')

# カラム名を平坦化（第2レベルを使用）
df_raw.columns = [b if b not in ('', 'Unnamed: 0_level_1') else a
                  for a, b in df_raw.columns]

COL_YEAR  = '年度'
COL_PREF  = '都道府県'
COL_TEMP  = '年平均気温'
COL_TMAX  = '最高気温（日最高気温の月平均の最高値）'
COL_TMIN  = '最低気温（日最低気温の月平均の最低値）'
COL_ENER  = '光熱・水道費（二人以上の世帯）'
COL_CONS  = '消費支出（二人以上の世帯）'
COL_POP65 = '65歳以上人口'
COL_POP   = '総人口'

# 必要列を選択
cols_needed = [COL_YEAR, COL_PREF, COL_TEMP, COL_TMAX, COL_TMIN,
               COL_ENER, COL_CONS, COL_POP65, COL_POP]
df = df_raw[cols_needed].copy()
df = df.dropna(subset=[COL_TEMP, COL_ENER, COL_CONS])
df[COL_YEAR] = df[COL_YEAR].astype(int)

# 派生変数
df['光熱水道費_万円']  = df[COL_ENER] / 10000
df['光熱水道費割合']   = df[COL_ENER] / df[COL_CONS] * 100
df['高齢化率']         = df[COL_POP65] / df[COL_POP] * 100
df['消費支出_log']     = np.log(df[COL_CONS])
df['最高気温']         = df[COL_TMAX]
df['最低気温']         = df[COL_TMIN]

# 2022年のみ抽出（メイン分析用）
df22 = df[df[COL_YEAR] == 2022].copy().reset_index(drop=True)

# 地域区分
KANREI = ['北海道', '青森県', '岩手県', '宮城県', '秋田県', '山形県', '福島県']
ONDAN  = ['沖縄県', '鹿児島県', '宮崎県', '高知県', '愛媛県', '徳島県', '香川県',
          '福岡県', '佐賀県', '長崎県', '熊本県', '大分県']

def get_region(pref):
    if pref in KANREI:
        return '寒冷地'
    elif pref in ONDAN:
        return '温暖地'
    else:
        return 'その他'

df22['地域区分'] = df22[COL_PREF].apply(get_region)
df['地域区分']   = df[COL_PREF].apply(get_region)

COLOR_MAP = {'寒冷地': '#1565C0', '温暖地': '#E65100', 'その他': '#78909C'}

print(f"2022年データ: {len(df22)}都道府県")
print(f"年平均気温 範囲: {df22[COL_TEMP].min():.1f}〜{df22[COL_TEMP].max():.1f}℃")
print(f"光熱水道費割合 範囲: {df22['光熱水道費割合'].min():.1f}〜{df22['光熱水道費割合'].max():.1f}%")

# 相関係数（参考出力）
r, p = stats.pearsonr(df22[COL_TEMP], df22['光熱水道費割合'])
print(f"\n年平均気温 × 光熱水道費割合: r={r:.3f}, p={p:.4f}")

# ══════════════════════════════════════════════════════════════════════════════
# 図1：年平均気温の都道府県別ランキング（2022年）
# ══════════════════════════════════════════════════════════════════════════════
df22_sorted = df22.sort_values(COL_TEMP, ascending=True).reset_index(drop=True)
colors1 = [COLOR_MAP[g] for g in df22_sorted['地域区分']]

fig1, ax1 = plt.subplots(figsize=(10, 11), dpi=150)
bars = ax1.barh(range(len(df22_sorted)), df22_sorted[COL_TEMP], color=colors1,
                edgecolor='white', linewidth=0.5, height=0.8)
ax1.set_yticks(range(len(df22_sorted)))
ax1.set_yticklabels(df22_sorted[COL_PREF], fontsize=8.5)
ax1.set_xlabel('年平均気温（℃）', fontsize=11)
ax1.set_title('図1　都道府県別 年平均気温ランキング（2022年）\n寒冷地（青）は低温、温暖地（橙）は高温', fontsize=12, fontweight='bold', pad=12)
ax1.axvline(df22[COL_TEMP].mean(), color='#333', linestyle='--', linewidth=1.2, alpha=0.7, label=f'全国平均 {df22[COL_TEMP].mean():.1f}℃')
ax1.set_xlim(0, df22[COL_TEMP].max() + 2)

# 気温の値ラベル
for i, (val, pref) in enumerate(zip(df22_sorted[COL_TEMP], df22_sorted[COL_PREF])):
    ax1.text(val + 0.15, i, f'{val:.1f}', va='center', fontsize=7.5, color='#333')

patch_k = mpatches.Patch(color='#1565C0', label='寒冷地（北海道・東北）')
patch_o = mpatches.Patch(color='#E65100', label='温暖地（九州・四国南部）')
patch_s = mpatches.Patch(color='#78909C', label='その他')
ax1.legend(handles=[patch_k, patch_o, patch_s], loc='lower right', fontsize=9)
ax1.grid(axis='x', alpha=0.3)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

fig1.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2022_H5_8_fig1_temp_ranking.png')
fig1.savefig(fig1_path, dpi=150, bbox_inches='tight')
plt.close(fig1)
print(f"図1 保存: {fig1_path}")

# ══════════════════════════════════════════════════════════════════════════════
# 図2：年平均気温 vs 光熱水道費割合 散布図（2022年）
# ══════════════════════════════════════════════════════════════════════════════
fig2, ax2 = plt.subplots(figsize=(9, 7), dpi=150)

for region, color in COLOR_MAP.items():
    mask = df22['地域区分'] == region
    ax2.scatter(df22.loc[mask, COL_TEMP], df22.loc[mask, '光熱水道費割合'],
                color=color, s=70, zorder=3, label=region, alpha=0.85, edgecolors='white', linewidth=0.8)

# 都道府県ラベル（全都道府県）
for _, row in df22.iterrows():
    pref_short = row[COL_PREF].replace('県', '').replace('都', '').replace('道', '').replace('府', '')
    ax2.annotate(pref_short,
                 (row[COL_TEMP], row['光熱水道費割合']),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, color='#333', zorder=4)

# 回帰直線
x_reg = df22[COL_TEMP].values
y_reg = df22['光熱水道費割合'].values
slope, intercept, r_val, p_val, se = stats.linregress(x_reg, y_reg)
x_line = np.linspace(x_reg.min(), x_reg.max(), 100)
ax2.plot(x_line, slope * x_line + intercept, color='#C62828', linewidth=1.8,
         linestyle='-', zorder=2, label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')

ax2.set_xlabel('年平均気温（℃）', fontsize=12)
ax2.set_ylabel('光熱水道費割合（%）', fontsize=12)
ax2.set_title('図2　年平均気温 vs 光熱水道費割合（2022年、47都道府県）\n'
              '気温が低い地域ほど家計に占める光熱費の割合が高い', fontsize=11, fontweight='bold', pad=12)

ax2.text(0.97, 0.97, f'r = {r_val:.3f}\np = {p_val:.4f}',
         transform=ax2.transAxes, ha='right', va='top',
         fontsize=10, bbox=dict(boxstyle='round,pad=0.4', facecolor='#FFF9C4', alpha=0.9))

ax2.legend(fontsize=9, loc='upper right')
ax2.grid(alpha=0.3)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

fig2.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2022_H5_8_fig2_scatter.png')
fig2.savefig(fig2_path, dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f"図2 保存: {fig2_path}")

# ══════════════════════════════════════════════════════════════════════════════
# 図3：OLS 回帰係数プロット（光熱水道費の決定要因）
# ══════════════════════════════════════════════════════════════════════════════
# 説明変数: 年平均気温、最高気温、最低気温、高齢化率、消費支出_log
var_labels = {
    COL_TEMP:    '年平均気温（℃）',
    '最高気温':  '最高気温（℃）',
    '最低気温':  '最低気温（℃）',
    '高齢化率':  '高齢化率（%）',
    '消費支出_log': '消費支出（対数）',
}
X_cols = list(var_labels.keys())
X_data = df22[X_cols].copy()
y_data = df22['光熱水道費割合'].copy()

# 欠損除去
valid = X_data.notna().all(axis=1) & y_data.notna()
X_data = X_data[valid]
y_data = y_data[valid]

# 標準化（β係数のプロット用）
X_std = (X_data - X_data.mean()) / X_data.std()
X_const = sm.add_constant(X_std)
model = sm.OLS(y_data, X_const).fit()
print("\nOLS 回帰結果:")
print(model.summary())

coefs  = model.params[1:]   # 定数項除く
ci_low = model.conf_int().iloc[1:, 0]
ci_hi  = model.conf_int().iloc[1:, 1]
pvals  = model.pvalues[1:]
labels = [var_labels[c] for c in X_cols]

# 係数の大きい順にソート
order  = np.argsort(np.abs(coefs.values))[::-1]
coefs_s  = coefs.values[order]
ci_low_s = ci_low.values[order]
ci_hi_s  = ci_hi.values[order]
pvals_s  = pvals.values[order]
labels_s = [labels[i] for i in order]

fig3, ax3 = plt.subplots(figsize=(8, 5), dpi=150)
colors3 = ['#1565C0' if c < 0 else '#E65100' for c in coefs_s]
sig_marker = ['*' if p < 0.05 else '' for p in pvals_s]

y_pos = np.arange(len(labels_s))
ax3.barh(y_pos, coefs_s, color=colors3, alpha=0.75, height=0.55, zorder=3)
ax3.errorbar(coefs_s, y_pos,
             xerr=[coefs_s - ci_low_s, ci_hi_s - coefs_s],
             fmt='none', color='#333', capsize=4, linewidth=1.2, zorder=4)
ax3.axvline(0, color='#333', linewidth=1, linestyle='-')
ax3.set_yticks(y_pos)
ax3.set_yticklabels([f'{lbl} {mk}' for lbl, mk in zip(labels_s, sig_marker)], fontsize=10)
ax3.set_xlabel('標準化回帰係数（β）', fontsize=11)
ax3.set_title(f'図3　OLS 回帰係数プロット（光熱水道費割合の決定要因）\n'
              f'目的変数: 光熱水道費割合（%）　R² = {model.rsquared:.3f}　n=47', fontsize=10, fontweight='bold', pad=10)

# * 凡例
ax3.text(0.98, 0.02, '* p<0.05　青=負の効果　橙=正の効果\n95%信頼区間を横棒で表示',
         transform=ax3.transAxes, ha='right', va='bottom', fontsize=9,
         bbox=dict(boxstyle='round,pad=0.3', facecolor='#F8F9FA', alpha=0.9))

ax3.grid(axis='x', alpha=0.3)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)

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

# ══════════════════════════════════════════════════════════════════════════════
# 図4：光熱水道費の時系列推移（寒冷地 vs 温暖地の比較）
# ══════════════════════════════════════════════════════════════════════════════
years_avail = sorted(df[COL_YEAR].unique())

# 地域区分ごとの年次平均を計算
ts = df.groupby([COL_YEAR, '地域区分']).agg(
    光熱水道費_万円=('光熱水道費_万円', 'mean'),
    光熱水道費割合=('光熱水道費割合', 'mean')
).reset_index()

ts_k = ts[ts['地域区分'] == '寒冷地'].sort_values(COL_YEAR)
ts_o = ts[ts['地域区分'] == '温暖地'].sort_values(COL_YEAR)
ts_s = ts[ts['地域区分'] == 'その他'].sort_values(COL_YEAR)

fig4, (ax4a, ax4b) = plt.subplots(1, 2, figsize=(12, 5), dpi=150)

# 左：光熱水道費_万円
ax4a.plot(ts_k[COL_YEAR], ts_k['光熱水道費_万円'], 'o-', color='#1565C0', linewidth=2,
          markersize=5, label='寒冷地（北海道・東北）')
ax4a.plot(ts_o[COL_YEAR], ts_o['光熱水道費_万円'], 's-', color='#E65100', linewidth=2,
          markersize=5, label='温暖地（九州・四国南部）')
ax4a.plot(ts_s[COL_YEAR], ts_s['光熱水道費_万円'], '^--', color='#78909C', linewidth=1.5,
          markersize=4, label='その他')
ax4a.set_xlabel('年度', fontsize=10)
ax4a.set_ylabel('光熱水道費（万円）', fontsize=10)
ax4a.set_title('光熱水道費の推移（万円）', fontsize=11, fontweight='bold')
ax4a.legend(fontsize=8)
ax4a.grid(alpha=0.3)
ax4a.spines['top'].set_visible(False)
ax4a.spines['right'].set_visible(False)
ax4a.set_xticks(years_avail)
ax4a.set_xticklabels([str(y) for y in years_avail], rotation=45, fontsize=8)

# 右：光熱水道費割合
ax4b.plot(ts_k[COL_YEAR], ts_k['光熱水道費割合'], 'o-', color='#1565C0', linewidth=2,
          markersize=5, label='寒冷地（北海道・東北）')
ax4b.plot(ts_o[COL_YEAR], ts_o['光熱水道費割合'], 's-', color='#E65100', linewidth=2,
          markersize=5, label='温暖地（九州・四国南部）')
ax4b.plot(ts_s[COL_YEAR], ts_s['光熱水道費割合'], '^--', color='#78909C', linewidth=1.5,
          markersize=4, label='その他')
ax4b.set_xlabel('年度', fontsize=10)
ax4b.set_ylabel('光熱水道費割合（%）', fontsize=10)
ax4b.set_title('消費支出に占める光熱水道費割合の推移', fontsize=11, fontweight='bold')
ax4b.legend(fontsize=8)
ax4b.grid(alpha=0.3)
ax4b.spines['top'].set_visible(False)
ax4b.spines['right'].set_visible(False)
ax4b.set_xticks(years_avail)
ax4b.set_xticklabels([str(y) for y in years_avail], rotation=45, fontsize=8)

fig4.suptitle('図4　光熱水道費の時系列推移：寒冷地 vs 温暖地（2012〜2023年）',
              fontsize=12, fontweight='bold', y=1.02)
fig4.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2022_H5_8_fig4_timeseries.png')
fig4.savefig(fig4_path, dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f"図4 保存: {fig4_path}")

print("\n全図の保存完了。")
print(f"出力先: {os.path.abspath(FIG_DIR)}")