地域間経済格差の構造分析
ジニ係数と格差決定要因

import os
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

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

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)

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

df_b['地域ブロック'] = df_b['都道府県'].map(region_map)

# ── ジニ係数関数 ───────────────────────────────────────────────────────────────
def gini(arr):
    arr = np.sort(arr[~np.isnan(arr)])
    n = len(arr)
    if n == 0 or np.sum(arr) == 0:
        return np.nan
    idx = np.arange(1, n + 1)
    return (2 * np.sum(idx * arr) - (n + 1) * np.sum(arr)) / (n * np.sum(arr))

# ── 年度ごとの格差指標を計算 ───────────────────────────────────────────────────
COL = '消費支出（二人以上の世帯）'

years = sorted(df_b['年度'].unique())
gini_list, cv_list, q9010_list = [], [], []

for yr in years:
    d = df_b[df_b['年度'] == yr][COL].dropna().values
    g = gini(d)
    cv = d.std() / d.mean() if d.mean() != 0 else np.nan
    q90 = np.percentile(d, 90)
    q10 = np.percentile(d, 10)
    q9010 = q90 / q10 if q10 != 0 else np.nan
    gini_list.append(g)
    cv_list.append(cv)
    q9010_list.append(q9010)

ts = pd.DataFrame({'年度': years, 'ジニ係数': gini_list,
                   '変動係数': cv_list, '90_10比': q9010_list})

# ── 2022年断面データ ───────────────────────────────────────────────────────────
df22 = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
df22['高齢化率'] = df22['65歳以上人口'] / df22['総人口'] * 100
df22['有効求人倍率'] = df22['月間有効求人数（一般）'] / df22['月間有効求職者数（一般）']
df22['高進学率'] = df22['大学学生数'] / df22['総人口'] * 1000
df22['宿泊者数_千人あたり'] = df22['延べ宿泊者数'] / df22['総人口'] * 1000

# ── 統計値の表示 ───────────────────────────────────────────────────────────────
print("=" * 60)
print("格差指標の時系列（2012〜2023年）")
print("=" * 60)
print(ts.to_string(index=False))
print()
print("=" * 60)
print("2022年 主要統計")
print("=" * 60)
gini2022 = gini(df22[COL].dropna().values)
cv2022 = df22[COL].std() / df22[COL].mean()
q90 = np.percentile(df22[COL].dropna(), 90)
q10 = np.percentile(df22[COL].dropna(), 10)
print(f"ジニ係数（消費支出, 2022年）: {gini2022:.4f}")
print(f"変動係数（消費支出, 2022年）: {cv2022:.4f}")
print(f"90/10比（消費支出, 2022年）: {q90/q10:.4f}")
print(f"最大県: {df22.loc[df22[COL].idxmax(),'都道府県']} ({df22[COL].max():,.0f}円)")
print(f"最小県: {df22.loc[df22[COL].idxmin(),'都道府県']} ({df22[COL].min():,.0f}円)")
print()

fig, ax1 = plt.subplots(figsize=(9, 5))
color1, color2 = '#1565C0', '#C62828'

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('ジニ係数', color=color1, fontsize=12)
l1, = ax1.plot(ts['年度'], ts['ジニ係数'], 'o-', color=color1,
               linewidth=2.2, markersize=7, label='ジニ係数')
ax1.tick_params(axis='y', labelcolor=color1)
ax1.set_ylim(0, ts['ジニ係数'].max() * 1.4)
ax1.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x:.3f}'))

ax2 = ax1.twinx()
ax2.set_ylabel('変動係数 (CV)', color=color2, fontsize=12)
l2, = ax2.plot(ts['年度'], ts['変動係数'], 's--', color=color2,
               linewidth=2.2, markersize=7, label='変動係数')
ax2.tick_params(axis='y', labelcolor=color2)
ax2.set_ylim(0, ts['変動係数'].max() * 1.4)
ax2.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x:.3f}'))

ax1.set_xticks(ts['年度'])
ax1.set_xticklabels([str(y) for y in ts['年度']], rotation=45, fontsize=10)
ax1.grid(axis='y', alpha=0.3)

lines = [l1, l2]
labels = [l.get_label() for l in lines]
ax1.legend(lines, labels, loc='upper right', fontsize=11)

plt.title('都道府県間 消費支出格差の時系列推移（2012〜2023年）', fontsize=13, fontweight='bold', pad=12)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U5_2_fig1.png'), bbox_inches='tight')
plt.close()
print("図1 保存完了")

vals = np.sort(df22[COL].dropna().values)
n = len(vals)
cum_pop = np.linspace(0, 1, n + 1)
cum_income = np.concatenate([[0], np.cumsum(vals) / np.sum(vals)])

fig, ax = plt.subplots(figsize=(7, 7))
ax.plot([0, 1], [0, 1], 'k--', linewidth=1.5, label='完全平等線', alpha=0.6)
ax.fill_between(cum_pop, cum_pop, cum_income, alpha=0.15, color='#1565C0')
ax.plot(cum_pop, cum_income, '-', color='#1565C0', linewidth=2.5,
        label=f'ローレンツ曲線（2022年）')

# ジニ係数注記
ax.text(0.35, 0.72, f'ジニ係数 = {gini2022:.4f}', fontsize=13, color='#1565C0',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#E3F2FD', edgecolor='#1565C0', alpha=0.9))
ax.annotate('', xy=(0.5, 0.4), xytext=(0.35, 0.65),
            arrowprops=dict(arrowstyle='->', color='#1565C0', lw=1.5))
ax.text(0.5, 0.28, 'ジニ面積\n（格差の大きさ）', fontsize=10, color='#1565C0',
        ha='center', alpha=0.8)

ax.set_xlabel('累積人口割合（都道府県数ベース）', fontsize=12)
ax.set_ylabel('累積消費支出割合', fontsize=12)
ax.set_title('ローレンツ曲線：都道府県別消費支出（2022年）', fontsize=13, fontweight='bold')
ax.legend(fontsize=11)
ax.set_xlim(0, 1); ax.set_ylim(0, 1)
ax.grid(alpha=0.3)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U5_2_fig2.png'), bbox_inches='tight')
plt.close()
print("図2 保存完了")

df_rank = df22[['都道府県', COL, '地域ブロック']].dropna().copy()
df_rank = df_rank.sort_values(COL, ascending=True).reset_index(drop=True)
bar_colors = [region_colors.get(r, '#aaa') for r in df_rank['地域ブロック']]

fig, ax = plt.subplots(figsize=(12, 10))
bars = ax.barh(df_rank['都道府県'], df_rank[COL] / 1000, color=bar_colors, edgecolor='white', linewidth=0.5)

# 全国平均線
national_mean = df22[COL].mean() / 1000
ax.axvline(national_mean, color='#333', linewidth=1.5, linestyle='--', label=f'全国平均 {national_mean*1000:,.0f}円')

ax.set_xlabel('消費支出（千円/月）', fontsize=12)
ax.set_title('都道府県別 消費支出ランキング（2022年、地域ブロック別色分け）', fontsize=13, fontweight='bold')
ax.grid(axis='x', alpha=0.3)

# 凡例（地域ブロック）
from matplotlib.patches import Patch
legend_handles = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax.legend(handles=legend_handles, loc='lower right', fontsize=10, title='地域ブロック')
ax.set_xlim(0, df_rank[COL].max() / 1000 * 1.15)

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U5_2_fig3.png'), bbox_inches='tight')
plt.close()
print("図3 保存完了")

reg_vars = {
    '高齢化率(%)': '高齢化率',
    '有効求人倍率': '有効求人倍率',
    '大学生割合(‰)': '高進学率',
    '宿泊者数(千人/千人)': '宿泊者数_千人あたり',
    '合計特殊出生率': '合計特殊出生率',
}

df_reg = df22[['都道府県', COL] + list(reg_vars.values())].dropna().copy()
y_raw = df_reg[COL].values
X_raw = df_reg[list(reg_vars.values())].values

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

X_sm = sm.add_constant(X_std)
model = sm.OLS(y_std, X_sm).fit()
print("=" * 60)
print("OLS回帰結果（標準化変数）")
print("=" * 60)
print(model.summary())

coef = model.params[1:]  # 定数項を除く
pvals = model.pvalues[1:]
ci_arr = model.conf_int()
if hasattr(ci_arr, 'iloc'):
    ci_low = ci_arr.iloc[1:, 0].values
    ci_high = ci_arr.iloc[1:, 1].values
else:
    ci_low = ci_arr[1:, 0]
    ci_high = ci_arr[1:, 1]
var_labels = list(reg_vars.keys())

# 係数の大きい順にソート
order = np.argsort(coef)
coef_s = coef[order]
pvals_s = pvals[order]
ci_low_s = ci_low[order]
ci_high_s = ci_high[order]
labels_s = [var_labels[i] for i in order]

fig, ax = plt.subplots(figsize=(9, 5.5))
colors_bar = ['#C62828' if p < 0.05 else '#BDBDBD' for p in pvals_s]
xerr_low  = coef_s - ci_low_s
xerr_high = ci_high_s - coef_s
bars = ax.barh(labels_s, coef_s, xerr=[xerr_low, xerr_high],
               color=colors_bar, edgecolor='white', linewidth=0.5,
               error_kw=dict(ecolor='#555', capsize=5, linewidth=1.5))
ax.axvline(0, color='black', linewidth=1.2)

# p値注記
for i, (c, p) in enumerate(zip(coef_s, pvals_s)):
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.')
    ax.text(c + (0.01 if c >= 0 else -0.01), i, f' {sig}',
            va='center', ha='left' if c >= 0 else 'right', fontsize=11,
            color='#C62828' if p < 0.05 else '#777')

ax.set_xlabel('標準化偏回帰係数（95% CI）', fontsize=12)
ax.set_title('消費支出格差の決定要因：OLS回帰（標準化係数）\n2022年断面、47都道府県',
             fontsize=12, fontweight='bold')

from matplotlib.patches import Patch
legend_handles2 = [
    Patch(facecolor='#C62828', label='有意（p<0.05）'),
    Patch(facecolor='#BDBDBD', label='非有意'),
]
ax.legend(handles=legend_handles2, loc='lower right', fontsize=10)
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U5_2_fig4.png'), bbox_inches='tight')
plt.close()
print("図4 保存完了")

# ── 最終統計値の出力（HTML用） ─────────────────────────────────────────────────
print()
print("=" * 60)
print("HTML用 統計サマリー")
print("=" * 60)
print(f"2022年ジニ係数: {gini2022:.4f}")
print(f"2022年変動係数: {cv2022:.4f}")
print(f"2022年 90/10比: {q90/q10:.3f}")
print(f"2012年ジニ係数: {ts[ts['年度']==2012]['ジニ係数'].values[0]:.4f}")
print(f"2023年ジニ係数: {ts[ts['年度']==2023]['ジニ係数'].values[0]:.4f}")
print(f"最大消費支出: {df22.loc[df22[COL].idxmax(),'都道府県']} ({df22[COL].max():,.0f}円)")
print(f"最小消費支出: {df22.loc[df22[COL].idxmin(),'都道府県']} ({df22[COL].min():,.0f}円)")
print(f"全国平均消費支出(2022): {df22[COL].mean():,.0f}円")
print(f"OLS R²: {model.rsquared:.3f}")
print(f"OLS 調整済みR²: {model.rsquared_adj:.3f}")
print()
for i, (lab, c, p) in enumerate(zip(var_labels, coef, pvals)):
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.')
    print(f"  {lab}: β={c:.3f}, p={p:.4f} {sig}")

print()
print("全図生成完了")