地域間賃金・雇用格差の構造分析
パネルデータによる長期的要因の解明

import os, numpy as np, pandas as pd
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
from linearmodels.panel import PanelOLS, RandomEffects
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)

print(f"読み込み完了: {df_b.shape[0]}行, {df_b.shape[1]}列")
print(f"年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}")
print(f"都道府県数: {df_b['都道府県'].nunique()}")

df_b['消費支出'] = df_b['消費支出（二人以上の世帯）']

# 求人倍率: 月間有効求人数 / 月間有効求職者数
df_b['求人倍率'] = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）']

# 高齢化率: 65歳以上人口 / 総人口 × 100
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100

# 第三次産業割合の代理: 転入率（移動活発＝都市型サービス産業）
# ただし旅館営業施設数/総人口を用いてサービス業密度を代理
df_b['旅館密度'] = df_b['旅館営業施設数（ホテルを含む）'] / df_b['総人口'] * 1000

# 人口密度の代理: 総人口の対数（都道府県面積が含まれないため）
df_b['log_人口'] = np.log(df_b['総人口'])

# 大学進学率: 高校卒業者のうち進学者数 / 高校卒業者数 × 100
df_b['大学進学率'] = df_b['高等学校卒業者のうち進学者数'] / df_b['高等学校卒業者数'] * 100

# 有効な行のみ使用（欠損値除去）
vars_needed = ['年度', '都道府県', '消費支出', '求人倍率', '高齢化率',
               '旅館密度', 'log_人口', '大学進学率']
df_p = df_b[vars_needed].dropna().copy()

print(f"\nパネルデータ: {df_p.shape[0]}行 ({df_p['都道府県'].nunique()}都道府県 × {df_p['年度'].nunique()}年度)")
print(f"変数の基本統計:")
print(df_p[['消費支出','求人倍率','高齢化率','旅館密度','大学進学率']].describe().round(2))

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

df_p['地域'] = df_p['都道府県'].map(region_map)

print("\n=== 図1: 地域別消費支出推移 ===")
fig1, ax1 = plt.subplots(figsize=(12, 7))

region_groups = df_p.groupby(['地域', '年度'])['消費支出'].mean().reset_index()

for region, color in region_colors.items():
    rg = region_groups[region_groups['地域'] == region].sort_values('年度')
    if len(rg) > 0:
        ax1.plot(rg['年度'], rg['消費支出'] / 1000, color=color,
                 linewidth=2.5, marker='o', markersize=5, label=region)

ax1.set_xlabel('年度', fontsize=13)
ax1.set_ylabel('消費支出（千円/月）', fontsize=13)
ax1.set_title('地域別消費支出の推移（2012〜2023年度）\n47都道府県・地域ブロック別平均', fontsize=14, fontweight='bold')
ax1.legend(loc='upper left', fontsize=11, framealpha=0.9)
ax1.grid(axis='y', alpha=0.4, linestyle='--')
ax1.set_xticks(sorted(df_p['年度'].unique()))
ax1.tick_params(axis='x', rotation=45)

# COVID-19注記
ax1.axvspan(2020, 2021, alpha=0.12, color='gray', label='COVID-19')
ax1.text(2020.15, ax1.get_ylim()[0] + (ax1.get_ylim()[1]-ax1.get_ylim()[0])*0.03,
         'COVID-19', fontsize=10, color='gray')

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

print("\n=== 図2: 散布図（求人倍率 vs 消費支出）===")
# 2019年度のスナップショット
df_2019 = df_p[df_p['年度'] == 2019].copy()

fig2, ax2 = plt.subplots(figsize=(12, 9))

for region, color in region_colors.items():
    sub = df_2019[df_2019['地域'] == region]
    ax2.scatter(sub['求人倍率'], sub['消費支出'] / 1000,
                color=color, s=80, alpha=0.85, label=region, zorder=3)
    for _, row in sub.iterrows():
        ax2.annotate(row['都道府県'].replace('県','').replace('都','').replace('府','').replace('道',''),
                     xy=(row['求人倍率'], row['消費支出'] / 1000),
                     xytext=(3, 3), textcoords='offset points',
                     fontsize=7.5, color='#333333')

# 回帰直線
x_all = df_2019['求人倍率'].values
y_all = df_2019['消費支出'].values / 1000
slope, intercept, r_val, p_val, se = stats.linregress(x_all, y_all)
x_line = np.linspace(x_all.min(), x_all.max(), 100)
ax2.plot(x_line, slope * x_line + intercept, 'k--', linewidth=1.8,
         label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})')

ax2.set_xlabel('求人倍率（月間有効求人数/求職者数）', fontsize=13)
ax2.set_ylabel('消費支出（千円/月）', fontsize=13)
ax2.set_title('求人倍率と消費支出の関係（2019年度・47都道府県）', fontsize=14, fontweight='bold')
ax2.legend(loc='upper left', fontsize=10, framealpha=0.9)
ax2.grid(alpha=0.3, linestyle='--')

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2019_U1_fig2.png')
fig2.savefig(fig2_path, dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f"保存: {fig2_path}")
print(f"相関係数: r={r_val:.4f}, p={p_val:.4f}")

print("\n=== 図3: FE係数プロット ===")
fig3, ax3 = plt.subplots(figsize=(9, 6))

if use_panel and fe_result is not None:
    params = fe_result.params.drop('Intercept', errors='ignore')
    std_errors = fe_result.std_errors.drop('Intercept', errors='ignore')
    pvalues = fe_result.pvalues.drop('Intercept', errors='ignore')

    # 95% CI
    ci_lower = params - 1.96 * std_errors
    ci_upper = params + 1.96 * std_errors

    var_labels = {
        '求人倍率': '求人倍率',
        '高齢化率': '高齢化率 (%)',
        '旅館密度': 'サービス業密度\n(旅館/千人)',
        '大学進学率': '大学進学率 (%)'
    }

    ys = np.arange(len(params))
    colors_bar = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvalues]

    bars = ax3.barh(ys, params.values, xerr=1.96*std_errors.values,
                    color=colors_bar, alpha=0.85, capsize=5, height=0.5,
                    error_kw={'linewidth': 2, 'ecolor': '#333'})

    ax3.set_yticks(ys)
    ax3.set_yticklabels([var_labels.get(v, v) for v in params.index], fontsize=12)
    ax3.axvline(0, color='black', linewidth=1.0, linestyle='-')
    ax3.set_xlabel('推定係数（95% CI）', fontsize=12)
    ax3.set_title('固定効果モデル（FE）：推定係数と95%信頼区間\n目的変数：消費支出（円/月）',
                  fontsize=13, fontweight='bold')

    # 凡例
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor='#e05c5c', alpha=0.85, label='有意（p < 0.05）'),
        Patch(facecolor='#aaaaaa', alpha=0.85, label='非有意（p >= 0.05）')
    ]
    ax3.legend(handles=legend_elements, loc='lower right', fontsize=11)
    ax3.grid(axis='x', alpha=0.3, linestyle='--')

else:
    # OLS fallback
    params = ols_result.params.drop('const', errors='ignore')
    std_errors = ols_result.bse.drop('const', errors='ignore')
    pvalues = ols_result.pvalues.drop('const', errors='ignore')

    ci_lower = params - 1.96 * std_errors
    ci_upper = params + 1.96 * std_errors

    ys = np.arange(len(params))
    colors_bar = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvalues]

    ax3.barh(ys, params.values, xerr=1.96*std_errors.values,
             color=colors_bar, alpha=0.85, capsize=5, height=0.5,
             error_kw={'linewidth': 2, 'ecolor': '#333'})

    ax3.set_yticks(ys)
    ax3.set_yticklabels(params.index, fontsize=12)
    ax3.axvline(0, color='black', linewidth=1.0)
    ax3.set_xlabel('推定係数（95% CI）', fontsize=12)
    ax3.set_title('OLSモデル：推定係数と95%信頼区間\n目的変数：消費支出（円/月）',
                  fontsize=13, fontweight='bold')
    ax3.grid(axis='x', alpha=0.3, linestyle='--')

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

print("\n=== パネル分析 ===")
fe_result = None
re_result = None
ols_result = None
use_panel = False

try:
    df_panel = df_p.set_index(['都道府県', '年度'])
    fe_res = PanelOLS.from_formula(
        '消費支出 ~ 求人倍率 + 高齢化率 + 旅館密度 + 大学進学率 + EntityEffects',
        data=df_panel
    ).fit(cov_type='clustered', cluster_entity=True)

    re_res = RandomEffects.from_formula(
        '消費支出 ~ 求人倍率 + 高齢化率 + 旅館密度 + 大学進学率',
        data=df_panel
    ).fit()

    fe_result = fe_res
    re_result = re_res
    use_panel = True

    print("固定効果モデル (FE):")
    print(fe_res.summary.tables[1])

    print("\nランダム効果モデル (RE):")
    print(re_res.summary.tables[1])

    # Hausman検定 (手動実装)
    # 共通変数の係数差を利用
    common_vars = ['求人倍率', '高齢化率', '旅館密度', '大学進学率']
    b_fe = fe_res.params[common_vars].values
    b_re = re_res.params[common_vars].values
    diff = b_fe - b_re

    # 共分散行列の差
    cov_fe = fe_res.cov.loc[common_vars, common_vars].values
    cov_re = re_res.cov.loc[common_vars, common_vars].values
    cov_diff = cov_fe - cov_re

    try:
        hausman_stat = float(diff @ np.linalg.inv(cov_diff) @ diff)
        hausman_df = len(common_vars)
        hausman_pval = 1 - stats.chi2.cdf(hausman_stat, df=hausman_df)
        print(f"\nHausman検定: χ²={hausman_stat:.4f}, df={hausman_df}, p={hausman_pval:.4f}")
        print("→ FE採択" if hausman_pval < 0.05 else "→ RE採択")
    except np.linalg.LinAlgError:
        # 擬似逆行列で代替
        hausman_stat = float(diff @ np.linalg.pinv(cov_diff) @ diff)
        hausman_df = len(common_vars)
        hausman_pval = 1 - stats.chi2.cdf(hausman_stat, df=hausman_df)
        print(f"\nHausman検定 (擬似逆行列): χ²={hausman_stat:.4f}, df={hausman_df}, p={hausman_pval:.4f}")

except Exception as e:
    print(f"パネル分析エラー: {e}")
    print("OLSフォールバックを使用")

    X_vars = ['求人倍率', '高齢化率', '旅館密度', '大学進学率']
    X_ols = sm.add_constant(df_p[X_vars])
    ols_result = sm.OLS(df_p['消費支出'], X_ols).fit()
    print(ols_result.summary())
    use_panel = False

    hausman_stat = None
    hausman_pval = None

print("\n=== 図4: FE vs RE 係数比較 ===")
fig4, ax4 = plt.subplots(figsize=(10, 6))

if use_panel and fe_result is not None and re_result is not None:
    common_vars = ['求人倍率', '高齢化率', '旅館密度', '大学進学率']
    var_labels = {
        '求人倍率': '求人倍率',
        '高齢化率': '高齢化率',
        '旅館密度': 'サービス業密度',
        '大学進学率': '大学進学率'
    }

    fe_coefs = fe_result.params[common_vars].values
    re_coefs = re_result.params[common_vars].values
    fe_se    = fe_result.std_errors[common_vars].values
    re_se    = re_result.std_errors[common_vars].values

    x_pos = np.arange(len(common_vars))
    width = 0.35

    bars_fe = ax4.bar(x_pos - width/2, fe_coefs, width,
                      yerr=1.96*fe_se, capsize=6,
                      label='固定効果 (FE)', color='#1565C0', alpha=0.85,
                      error_kw={'linewidth': 2})
    bars_re = ax4.bar(x_pos + width/2, re_coefs, width,
                      yerr=1.96*re_se, capsize=6,
                      label='変量効果 (RE)', color='#E65100', alpha=0.85,
                      error_kw={'linewidth': 2})

    ax4.set_xticks(x_pos)
    ax4.set_xticklabels([var_labels.get(v, v) for v in common_vars], fontsize=12)
    ax4.axhline(0, color='black', linewidth=1.0)
    ax4.set_ylabel('推定係数', fontsize=12)

    try:
        hausman_stat_v = hausman_stat
        hausman_pval_v = hausman_pval
        title_str = (f'Hausman検定によるFE・RE係数比較\n'
                     f'Hausman χ²={hausman_stat_v:.3f}, p={hausman_pval_v:.4f} '
                     f'→ {"FE採択" if hausman_pval_v < 0.05 else "RE採択"}')
    except NameError:
        title_str = 'Hausman検定によるFE・RE係数比較'

    ax4.set_title(title_str, fontsize=13, fontweight='bold')
    ax4.legend(fontsize=12, loc='upper right')
    ax4.grid(axis='y', alpha=0.3, linestyle='--')

else:
    # OLS only: show coefficient bar chart with CI
    if ols_result is not None:
        params = ols_result.params.drop('const', errors='ignore')
        std_errors = ols_result.bse.drop('const', errors='ignore')
        pvalues = ols_result.pvalues.drop('const', errors='ignore')

        ys = np.arange(len(params))
        ax4.bar(ys, params.values, color='#1565C0', alpha=0.8)
        ax4.set_xticks(ys)
        ax4.set_xticklabels(params.index, fontsize=11)
        ax4.axhline(0, color='black', linewidth=1.0)
        ax4.set_title('OLS推定係数（パネル分析フォールバック）', fontsize=13, fontweight='bold')
        ax4.grid(axis='y', alpha=0.3, linestyle='--')

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2019_U1_fig4.png')
fig4.savefig(fig4_path, dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f"保存: {fig4_path}")

print("\n" + "="*60)
print("統計値サマリー（HTML用）")
print("="*60)

print(f"\n■ データ規模")
print(f"  観測数: {df_p.shape[0]} ({df_p['都道府県'].nunique()}都道府県 × {df_p['年度'].nunique()}年)")
print(f"  期間: {df_p['年度'].min()}〜{df_p['年度'].max()}年度")

print(f"\n■ 消費支出の記述統計")
cs = df_p['消費支出']
print(f"  平均: {cs.mean():,.0f}円/月")
print(f"  標準偏差: {cs.std():,.0f}円/月")
print(f"  最大: {cs.max():,.0f}円/月 ({df_p.loc[cs.idxmax(),'都道府県']})")
print(f"  最小: {cs.min():,.0f}円/月 ({df_p.loc[cs.idxmin(),'都道府県']})")

print(f"\n■ 求人倍率の記述統計")
kj = df_p['求人倍率']
print(f"  平均: {kj.mean():.3f}")
print(f"  最大: {kj.max():.3f}")
print(f"  最小: {kj.min():.3f}")

print(f"\n■ 相関分析（2019年度スナップショット）")
print(f"  求人倍率 vs 消費支出: r={r_val:.4f}, p={p_val:.4f}")

if use_panel and fe_result is not None:
    print(f"\n■ 固定効果モデル（FE）結果")
    for var in fe_result.params.index:
        if var == 'Intercept':
            continue
        coef = fe_result.params[var]
        se_v = fe_result.std_errors[var]
        pv   = fe_result.pvalues[var]
        sig  = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'n.s.'
        print(f"  {var}: β={coef:+.2f}, SE={se_v:.2f}, p={pv:.4f} {sig}")

    try:
        print(f"\n■ Hausman検定")
        print(f"  χ²={hausman_stat:.4f}, df={len(common_vars)}, p={hausman_pval:.4f}")
        print(f"  判断: {'固定効果モデルを採択' if hausman_pval < 0.05 else '変量効果モデルを採択'}")
    except NameError:
        pass

print("\n=== 完了 ===")
print(f"図1: {fig1_path}")
print(f"図2: {fig2_path}")
print(f"図3: {fig3_path}")
print(f"図4: {fig4_path}")