地域別賃金水準の決定要因
消費支出・雇用指標を用いたパネル分析

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

from linearmodels.panel import PanelOLS
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

import os
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)
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

# 変数生成
df_b['消費支出_万円'] = df_b['消費支出（二人以上の世帯）'] / 10000
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['求人倍率'] = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'].replace(0, np.nan)
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['大学学生_log'] = np.log(df_b['大学学生数'].clip(lower=1))
df_b['純移動率'] = (df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']) / df_b['総人口'] * 1000

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

fig, ax = plt.subplots(figsize=(10, 5))
yearly = df_b.groupby(['年度', '地域'])['消費支出_万円'].mean().reset_index()
for reg, grp in yearly.groupby('地域'):
    ax.plot(grp['年度'], grp['消費支出_万円'], marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('消費支出（万円/月）', fontsize=12)
ax.set_title('地域別 消費支出（賃金代理）の推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_14_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

df_box = df_b[['都道府県', '地域', '求人倍率']].dropna()
region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
fig, ax = plt.subplots(figsize=(10, 5))
bxdata = [df_box[df_box['地域'] == reg]['求人倍率'].values for reg in region_order]
bp = ax.boxplot(bxdata, patch_artist=True, labels=region_order)
for patch, reg in zip(bp['boxes'], region_order):
    patch.set_facecolor(region_colors.get(reg, 'gray'))
    patch.set_alpha(0.7)
ax.set_xlabel('地域', fontsize=12)
ax.set_ylabel('求人倍率', fontsize=12)
ax.set_title('地域別 求人倍率の分布（2012〜2023年）', fontsize=13, fontweight='bold')
ax.grid(axis='y', alpha=0.3)
ax.axhline(1.0, color='red', linestyle='--', linewidth=1.5, label='倍率=1.0')
ax.legend(fontsize=10)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_14_fig2_boxplot.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

panel_vars = ['求人倍率', '高齢化率', '大学学生_log', '純移動率']
df_panel = df_b[['年度', '都道府県', '消費支出_log'] + panel_vars].dropna()
df_panel = df_panel.set_index(['都道府県', '年度'])
y = df_panel['消費支出_log']
X = sm.add_constant(df_panel[panel_vars])
try:
    mod = PanelOLS(y, X, entity_effects=True, drop_absorbed=True)
    res = mod.fit(cov_type='clustered', cluster_entity=True)
    coefs = res.params.drop('const', errors='ignore')
    ses = res.std_errors.drop('const', errors='ignore')
    pvals = res.pvalues.drop('const', errors='ignore')
    print(res.summary)
except Exception as e:
    print(f"FE error: {e}")
    df_ols = df_b[df_b['年度'] == 2022][['消費支出_log'] + panel_vars].dropna()
    res_ols = sm.OLS(df_ols['消費支出_log'], sm.add_constant(df_ols[panel_vars])).fit()
    coefs = res_ols.params.drop('const')
    ses = res_ols.bse.drop('const')
    pvals = res_ols.pvalues.drop('const')

fig, ax = plt.subplots(figsize=(8, 5))
colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
        error_kw={'elinewidth': 1.5, 'capsize': 4})
ax.set_yticks(range(len(coefs)))
ax.set_yticklabels(coefs.index, fontsize=10)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('回帰係数', fontsize=12)
ax.set_title('消費支出（賃金代理）の決定要因 — FEパネルOLS\n（赤=p<0.05, 灰=非有意）', fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_14_fig3_fe_coef.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

df_2022 = df_b[df_b['年度'] == 2022].copy()
df_sc = df_2022[['都道府県', '求人倍率', '消費支出_万円', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_sc.groupby('地域'):
    ax.scatter(grp['求人倍率'], grp['消費支出_万円'],
               color=region_colors.get(reg, 'gray'), label=reg, s=60, alpha=0.8)
for _, row in df_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['求人倍率'], row['消費支出_万円']),
                fontsize=6, alpha=0.6, xytext=(2, 2), textcoords='offset points')
slope, intercept, r, p, _ = stats.linregress(df_sc['求人倍率'], df_sc['消費支出_万円'])
xr = np.linspace(df_sc['求人倍率'].min(), df_sc['求人倍率'].max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線 (r={r:.2f}, p={p:.3f})')
ax.set_xlabel('求人倍率', fontsize=12)
ax.set_ylabel('消費支出（万円/月）', fontsize=12)
ax.set_title('求人倍率 vs 消費支出（2022年, N=47）', fontsize=13, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_14_fig4_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")