"""
2023_U4_katsuyo.py
行動制限下における家計消費の変化に伴う経済波及効果の算出
統計活用奨励賞 [大学生の部] — 2023年度 統計データ分析コンペティション

分析内容:
  1. 消費支出5カテゴリの時系列推移 (2012-2023)、COVID期シェード
  2. 2019年基準 変化率（%）：3期間（行動制限前・行動制限下・緩和後）
  3. OLS回帰：都道府県別 教養娯楽費変化率 → 消費支出変化率
  4. 消費支出削減トップ10都道府県（2020年）

実データ: SSDSE-B-2026.csv（2012-2023年、47都道府県）
合成データ・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_U4_katsuyo.py
# ============================================================


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

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

# ----------------------------------------------------------------
# パス設定
# ----------------------------------------------------------------
_dir = os.path.dirname(os.path.abspath(__file__))
FIG_DIR = os.path.join(_dir, '..', 'html', 'figures')
DATA_B  = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-B-2026.csv')

os.makedirs(FIG_DIR, exist_ok=True)

# ----------------------------------------------------------------
# データ読み込み
# ----------------------------------------------------------------
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)

# 都道府県レベル（R01000〜R47000）のみ抽出
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 消費支出カラム名（SSDSE-B 実際のカラム名）
COL_TOTAL  = '消費支出（二人以上の世帯）'
COL_FOOD   = '食料費（二人以上の世帯）'
COL_ENT    = '教養娯楽費（二人以上の世帯）'
COL_TRANS  = '交通・通信費（二人以上の世帯）'
COL_HOUSE  = '住居費（二人以上の世帯）'
COL_HEALTH = '保健医療費（二人以上の世帯）'

COLS_ALL = [COL_TOTAL, COL_FOOD, COL_ENT, COL_TRANS, COL_HOUSE, COL_HEALTH]
COLS_5   = [COL_FOOD, COL_ENT, COL_TRANS, COL_HOUSE, COL_HEALTH]  # 5カテゴリ

for col in COLS_ALL:
    df_b[col] = pd.to_numeric(df_b[col], errors='coerce')
df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce')

# 全国平均（都道府県平均）を年度別に集計
nat_avg = df_b.groupby('年度')[COLS_ALL].mean()
years   = nat_avg.index.to_numpy()   # 2012-2023

print("データ読み込み完了")
print(f"  年度: {years}")
print(f"  都道府県数: {df_b['地域コード'].nunique()}")
print(nat_avg[[COL_TOTAL, COL_ENT]].to_string())

# ----------------------------------------------------------------
# 色・ラベル設定
# ----------------------------------------------------------------
CAT_LABELS = {
    COL_FOOD:   '食料費',
    COL_ENT:    '教養娯楽費',
    COL_TRANS:  '交通・通信費',
    COL_HOUSE:  '住居費',
    COL_HEALTH: '保健医療費',
}
CAT_COLORS = {
    COL_FOOD:   '#1565C0',
    COL_ENT:    '#C62828',
    COL_TRANS:  '#2E7D32',
    COL_HOUSE:  '#E65100',
    COL_HEALTH: '#6A1B9A',
}

COVID_START = 2020
COVID_END   = 2021

# 2019年基準値（変化率計算用）
base_nat = nat_avg.loc[2019]

# ================================================================
# 図1: 消費5カテゴリ 時系列（2012-2023）、COVID期シェード
# ================================================================
fig, ax = plt.subplots(figsize=(11, 6))

for col in COLS_5:
    vals = nat_avg[col].values / 1000   # 円 → 千円
    ax.plot(years, vals, marker='o', markersize=5,
            color=CAT_COLORS[col], label=CAT_LABELS[col], lw=2)

# 行動制限期間シェード
ax.axvspan(COVID_START - 0.5, COVID_END + 0.5,
           color='#FFCC80', alpha=0.35, zorder=0, label='行動制限期間 (2020-2021)')

# 基準線（2019）
ax.axvline(2019, color='gray', lw=1.2, linestyle='--', alpha=0.7, label='2019年（基準）')

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('月額支出（千円）', fontsize=12)
ax.set_title('図1: 家計消費支出の推移（2012〜2023年、都道府県平均）', fontsize=13)
ax.set_xticks(years)
ax.set_xticklabels([str(y) for y in years], rotation=45, ha='right')
ax.legend(fontsize=10, loc='upper left')
ax.grid(True, axis='y', alpha=0.3)
ax.grid(True, axis='x', alpha=0.15)

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2023_U4_fig1_timeseries.png')
plt.savefig(fig1_path, dpi=150)
plt.close()
print(f"fig1 saved: {fig1_path}")

# ================================================================
# 図2: 2019年基準変化率（%）— 3期間 × 5カテゴリ 集合棒グラフ
# ================================================================
# 3期間: コロナ前平均(2018-2019)、行動制限下(2020-2021)、緩和後(2022-2023)
periods = {
    'コロナ前\n(2018-2019)':  df_b[df_b['年度'].isin([2018, 2019])].groupby('都道府県')[COLS_5].mean(),
    '行動制限下\n(2020-2021)': df_b[df_b['年度'].isin([2020, 2021])].groupby('都道府県')[COLS_5].mean(),
    '緩和後\n(2022-2023)':    df_b[df_b['年度'].isin([2022, 2023])].groupby('都道府県')[COLS_5].mean(),
}

# 2019年 都道府県別基準
base_pref = df_b[df_b['年度'] == 2019].groupby('都道府県')[COLS_5].mean()

# 各期間の全国平均変化率（%）
period_names = list(periods.keys())
change_pct = {}
for pname, pdata in periods.items():
    common_pref = pdata.index.intersection(base_pref.index)
    diff_pct = ((pdata.loc[common_pref] - base_pref.loc[common_pref])
                / base_pref.loc[common_pref] * 100)
    change_pct[pname] = diff_pct.mean()   # 全国平均

n_periods = len(period_names)
n_cats    = len(COLS_5)
x_base    = np.arange(n_periods)
width     = 0.15
offsets   = np.linspace(-(n_cats - 1) / 2, (n_cats - 1) / 2, n_cats) * width

fig, ax = plt.subplots(figsize=(11, 6))

for i, col in enumerate(COLS_5):
    vals = [change_pct[p][col] for p in period_names]
    bars = ax.bar(x_base + offsets[i], vals, width=width * 0.92,
                  color=CAT_COLORS[col], label=CAT_LABELS[col], alpha=0.85)
    for bar, v in zip(bars, vals):
        ypos = bar.get_height() if v >= 0 else bar.get_height() - 0.3
        ax.text(bar.get_x() + bar.get_width() / 2,
                ypos + (0.15 if v >= 0 else -0.8),
                f'{v:+.1f}%', ha='center', va='bottom', fontsize=7.5, rotation=90)

ax.axhline(0, color='black', lw=1.0)
ax.set_xticks(x_base)
ax.set_xticklabels(period_names, fontsize=11)
ax.set_ylabel('2019年比 変化率（%）', fontsize=12)
ax.set_title('図2: 消費カテゴリ別 2019年基準変化率（都道府県平均）', fontsize=13)
ax.legend(fontsize=10, loc='lower right')
ax.grid(True, axis='y', alpha=0.3)

# 行動制限期帯を背景色で区別
ax.axvspan(0.5, 1.5, color='#FFCC80', alpha=0.25, zorder=0)
ax.text(1, ax.get_ylim()[0] + 0.1, '行動制限下', ha='center', fontsize=9, color='#E65100')

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2023_U4_fig2_pct_change.png')
plt.savefig(fig2_path, dpi=150)
plt.close()
print(f"fig2 saved: {fig2_path}")

# ================================================================
# 図3: 散布図 — 都道府県別 教養娯楽費変化 vs 食料費変化（2020 vs 2019）
# ================================================================
df_2019 = df_b[df_b['年度'] == 2019].set_index('都道府県')
df_2020 = df_b[df_b['年度'] == 2020].set_index('都道府県')

common_prefs = df_2019.index.intersection(df_2020.index)

ent_chg  = ((df_2020.loc[common_prefs, COL_ENT]  - df_2019.loc[common_prefs, COL_ENT])
             / df_2019.loc[common_prefs, COL_ENT] * 100)
food_chg = ((df_2020.loc[common_prefs, COL_FOOD] - df_2019.loc[common_prefs, COL_FOOD])
             / df_2019.loc[common_prefs, COL_FOOD] * 100)

# OLS回帰直線
slope, intercept, r_val, p_val, se_slope = stats.linregress(ent_chg.values, food_chg.values)
x_line = np.linspace(ent_chg.min() - 1, ent_chg.max() + 1, 200)
y_line = slope * x_line + intercept

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

ax.scatter(ent_chg.values, food_chg.values, color='#1565C0', alpha=0.75, s=55, zorder=3)
ax.plot(x_line, y_line, color='#C62828', lw=2, linestyle='--',
        label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')

# 都道府県名ラベル（上位・下位5件のみ）
ent_sorted = ent_chg.sort_values()
label_prefs = list(ent_sorted.index[:3]) + list(ent_sorted.index[-3:])
for pref in label_prefs:
    ax.annotate(pref, (ent_chg[pref], food_chg[pref]),
                fontsize=8, xytext=(4, 4), textcoords='offset points',
                color='#333333')

ax.axhline(0, color='gray', lw=0.8, linestyle=':')
ax.axvline(0, color='gray', lw=0.8, linestyle=':')
ax.set_xlabel('教養娯楽費 変化率（2020 vs 2019、%）', fontsize=12)
ax.set_ylabel('食料費 変化率（2020 vs 2019、%）', fontsize=12)
ax.set_title('図3: 教養娯楽費変化 vs 食料費変化（都道府県別、2020年）', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

# 象限ラベル
ax.text(ax.get_xlim()[1] * 0.98, ax.get_ylim()[1] * 0.95,
        '娯楽↑食料↑', ha='right', va='top', fontsize=9, color='#555')
ax.text(ax.get_xlim()[0] * 0.98, ax.get_ylim()[1] * 0.95,
        '娯楽↓食料↑', ha='left', va='top', fontsize=9, color='#555')
ax.text(ax.get_xlim()[1] * 0.98, ax.get_ylim()[0] * 0.95,
        '娯楽↑食料↓', ha='right', va='bottom', fontsize=9, color='#555')
ax.text(ax.get_xlim()[0] * 0.98, ax.get_ylim()[0] * 0.95,
        '娯楽↓食料↓', ha='left', va='bottom', fontsize=9, color='#555')

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2023_U4_fig3_scatter.png')
plt.savefig(fig3_path, dpi=150)
plt.close()
print(f"fig3 saved: {fig3_path}")

# ================================================================
# 図4: 消費支出削減額 トップ10都道府県（2020年）
# ================================================================
total_chg_abs = (df_2020.loc[common_prefs, COL_TOTAL]
                 - df_2019.loc[common_prefs, COL_TOTAL])   # 円（月額、マイナスが削減）
total_chg_pct = total_chg_abs / df_2019.loc[common_prefs, COL_TOTAL] * 100

# 削減が大きい順（変化率の小さい＝最も下がった）トップ10
top10 = total_chg_pct.sort_values().head(10)

fig, ax = plt.subplots(figsize=(10, 6))

colors_bar = ['#C62828' if v < -5 else ('#EF6C00' if v < -2 else '#1565C0')
               for v in top10.values]
bars = ax.barh(top10.index[::-1], top10.values[::-1], color=colors_bar[::-1], alpha=0.85)

for bar, v in zip(bars, top10.values[::-1]):
    ax.text(v - 0.05, bar.get_y() + bar.get_height() / 2,
            f'{v:+.1f}%', va='center', ha='right', fontsize=10, fontweight='600',
            color='white')

ax.axvline(total_chg_pct.mean(), color='navy', lw=1.5, linestyle='--',
           label=f'全国平均 {total_chg_pct.mean():+.1f}%')
ax.set_xlabel('消費支出変化率（2020 vs 2019、%）', fontsize=12)
ax.set_title('図4: 消費支出削減が大きかった都道府県 トップ10（2020年 vs 2019年）', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, axis='x', alpha=0.3)

# 凡例パッチ
p1 = mpatches.Patch(color='#C62828', alpha=0.85, label='−5%超の削減')
p2 = mpatches.Patch(color='#EF6C00', alpha=0.85, label='−2〜−5%の削減')
p3 = mpatches.Patch(color='#1565C0', alpha=0.85, label='−2%未満の削減')
ax.legend(handles=[p1, p2, p3,
                   plt.Line2D([0], [0], color='navy', lw=1.5, linestyle='--',
                              label=f'全国平均 {total_chg_pct.mean():+.1f}%')],
          fontsize=10, loc='lower right')

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2023_U4_fig4_top10.png')
plt.savefig(fig4_path, dpi=150)
plt.close()
print(f"fig4 saved: {fig4_path}")

# ================================================================
# 補足出力: OLS 簡易乗数効果
# ================================================================
print("\n" + "="*60)
print("簡易経済乗数の推計 (全国平均)")
print("="*60)

# 全47都道府県の2020年 消費支出変化率 vs 教養娯楽費変化率
ols_x = ent_chg.values
ols_y = total_chg_pct.loc[common_prefs].values

slope_main, intercept_main, r_main, p_main, se_main = stats.linregress(ols_x, ols_y)

print(f"OLS: 消費支出変化率 = {intercept_main:.3f} + {slope_main:.3f} × 教養娯楽費変化率")
print(f"  R² = {r_main**2:.4f}, p = {p_main:.4f}, SE(slope) = {se_main:.4f}")

# 粗い乗数推計: 全国平均消費支出削減額
nat_2019_total = nat_avg.loc[2019, COL_TOTAL]
nat_2020_total = nat_avg.loc[2020, COL_TOTAL]
delta_monthly  = nat_2020_total - nat_2019_total          # 月額円
delta_annual   = delta_monthly * 12 * 47                  # 47都道府県年額（簡易）

print(f"\n全国平均消費支出: 2019年 {nat_2019_total:,.0f}円/月 → 2020年 {nat_2020_total:,.0f}円/月")
print(f"  月額変化: {delta_monthly:+,.0f}円（{delta_monthly/nat_2019_total*100:+.1f}%）")
print(f"  簡易年間削減額（47都道府県合計参考値）: {delta_annual/1e8:,.0f}億円")
print(f"\n教養娯楽費: 2019年 {nat_avg.loc[2019, COL_ENT]:,.0f}円/月 → 2020年 {nat_avg.loc[2020, COL_ENT]:,.0f}円/月")
ent_delta_pct = (nat_avg.loc[2020, COL_ENT] - nat_avg.loc[2019, COL_ENT]) / nat_avg.loc[2019, COL_ENT] * 100
print(f"  変化率: {ent_delta_pct:+.1f}%（最も大きな削減）")

print("\n全4図の生成が完了しました。")