高齢化と地域経済成長の関係：動態的パネル分析

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)

df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['消費支出'] = pd.to_numeric(df_b['消費支出（二人以上の世帯）'], errors='coerce')
df_b['保健医療費_人'] = pd.to_numeric(df_b['保健医療費（二人以上の世帯）'], errors='coerce')
df_b['求人数'] = pd.to_numeric(df_b['月間有効求人数（一般）'], errors='coerce')
df_b['求職者数'] = pd.to_numeric(df_b['月間有効求職者数（一般）'], errors='coerce')
# 求人倍率（有効求人数 / 有効求職者数）
df_b['求人倍率'] = df_b['求人数'] / df_b['求職者数']

# 人口密度の代理変数: 高齢化率の高い地域は過疎地に多い（ここでは総人口/都道府県数で代替）
# SSDSE-Bには面積データがないため、総人口を代理変数として使用
df_b['log_人口'] = np.log(df_b['総人口'])

# 2次項
df_b['高齢化率2'] = df_b['高齢化率'] ** 2

# 消費支出成長率（前年比）
df_b_sorted = df_b.sort_values(['都道府県', '年度']).copy()
df_b_sorted['消費支出成長率'] = df_b_sorted.groupby('都道府県')['消費支出'].pct_change() * 100

# 2012〜2023 絞り込み
df_panel = df_b_sorted[(df_b_sorted['年度'] >= 2012) & (df_b_sorted['年度'] <= 2023)].copy()

print("=== データ概要 ===")
print(f"観測数: {len(df_panel)}")
print(f"都道府県数: {df_panel['都道府県'].nunique()}")
print(f"年度範囲: {df_panel['年度'].min()} - {df_panel['年度'].max()}")
print()

print()
print("=== HTML用統計サマリー ===")
print(f"分析期間: {df_panel['年度'].min()}〜{df_panel['年度'].max()}年")
print(f"サンプル数（延べ）: {len(df_panel.dropna(subset=['消費支出','高齢化率']))} 観測")
print(f"2023年 全国平均高齢化率: {aging_2023['高齢化率'].mean():.1f}%")
print(f"2023年 最高高齢化率: {aging_2023.iloc[0]['都道府県']} {aging_2023.iloc[0]['高齢化率']:.1f}%")
print(f"2023年 最低高齢化率: {aging_2023.iloc[-1]['都道府県']} {aging_2023.iloc[-1]['高齢化率']:.1f}%")

# 2012年との比較
aging_2012 = df_panel[df_panel['年度'] == 2012][['都道府県', '高齢化率']].set_index('都道府県')
aging_2023_idx = aging_2023.set_index('都道府県')
change = (aging_2023_idx['高齢化率'] - aging_2012['高齢化率']).sort_values(ascending=False)
print(f"2012→2023 高齢化率増加幅最大: {change.index[0]} +{change.iloc[0]:.1f}pt")
print(f"2012→2023 高齢化率増加幅最小: {change.index[-1]} +{change.iloc[-1]:.1f}pt")

# 消費支出との相関
corr_data = df_2022[['高齢化率', '消費支出']].dropna()
r2, p2 = stats.pearsonr(corr_data['高齢化率'], corr_data['消費支出'])
print(f"2022年断面 高齢化率×消費支出 相関: r={r2:.3f} p={p2:.4f}")

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

fig1, ax1 = plt.subplots(figsize=(10, 6))

regions_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
df_region_aging = df_panel.groupby(['年度', '地域'])['高齢化率'].mean().reset_index()

for region in regions_order:
    d = df_region_aging[df_region_aging['地域'] == region]
    d = d.sort_values('年度')
    ax1.plot(d['年度'], d['高齢化率'], marker='o', markersize=5,
             color=region_colors[region], linewidth=2, label=region)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('高齢化率（%）', fontsize=12)
ax1.set_title('地域別 高齢化率の時系列推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax1.legend(loc='upper left', fontsize=10)
ax1.set_xticks(range(2012, 2024))
ax1.grid(True, alpha=0.3)
ax1.set_ylim(15, 40)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_U2_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("図1 保存完了")

df_2022 = df_panel[df_panel['年度'] == 2022].copy()
df_2022 = df_2022.dropna(subset=['高齢化率', '消費支出'])

fig2, ax2 = plt.subplots(figsize=(11, 8))

for region in regions_order:
    d = df_2022[df_2022['地域'] == region]
    ax2.scatter(d['高齢化率'], d['消費支出'] / 1000, color=region_colors[region],
                s=60, alpha=0.8, label=region, zorder=3)

# 都道府県ラベル
for _, row in df_2022.iterrows():
    ax2.annotate(row['都道府県'].replace('県', '').replace('都', '').replace('府', '').replace('道', ''),
                 (row['高齢化率'], row['消費支出'] / 1000),
                 fontsize=7, ha='center', va='bottom', color='#333333',
                 xytext=(0, 4), textcoords='offset points')

# 2次曲線フィット
x_fit = df_2022['高齢化率'].values
y_fit = df_2022['消費支出'].values / 1000
coeffs = np.polyfit(x_fit, y_fit, 2)
x_range = np.linspace(x_fit.min(), x_fit.max(), 200)
y_range = np.polyval(coeffs, x_range)
ax2.plot(x_range, y_range, 'k--', linewidth=2, label=f'2次フィット\n({coeffs[0]:.1f}x²+{coeffs[1]:.1f}x+{coeffs[2]:.0f})')

# 相関係数
r, p = stats.pearsonr(x_fit, y_fit)
ax2.text(0.05, 0.95, f'r = {r:.3f} (p = {p:.3f})', transform=ax2.transAxes,
         fontsize=11, va='top', bbox=dict(boxstyle='round,pad=0.4', facecolor='white', alpha=0.8))

ax2.set_xlabel('高齢化率（%）', fontsize=12)
ax2.set_ylabel('消費支出（千円）', fontsize=12)
ax2.set_title('高齢化率と消費支出の関係（2022年断面、47都道府県）', fontsize=14, fontweight='bold')
ax2.legend(loc='upper right', fontsize=9)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_U2_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("図2 保存完了")

fig3, ax3 = plt.subplots(figsize=(9, 5))

if fe_params is not None:
    # const以外の係数のみ
    plot_keys = [k for k in fe_params.index if k != 'const']
    plot_labels = [param_names_en.get(k, k) for k in plot_keys]
    plot_vals = [fe_params[k] for k in plot_keys]
    # 標準化（係数の絶対値）で色分け
    colors_bar = ['#e05c5c' if v < 0 else '#4e9af1' for v in plot_vals]

    # 信頼区間
    # conf_int の列名を動的に取得（linearmodels: 'lower'/'upper', statsmodels: 0/1）
    ci_cols = list(fe_conf_int.columns)
    lo_col, hi_col = ci_cols[0], ci_cols[1]
    lo = [fe_conf_int.loc[k, lo_col] for k in plot_keys]
    hi = [fe_conf_int.loc[k, hi_col] for k in plot_keys]
    xerr_lo = [abs(v - l) for v, l in zip(plot_vals, lo)]
    xerr_hi = [abs(h - v) for v, h in zip(plot_vals, hi)]

    bars = ax3.barh(plot_labels, plot_vals, color=colors_bar, alpha=0.8, height=0.6)
    ax3.errorbar(plot_vals, plot_labels, xerr=[xerr_lo, xerr_hi],
                 fmt='none', color='black', capsize=5, linewidth=1.5)
    ax3.axvline(0, color='black', linewidth=0.8, linestyle='-')
    ax3.set_xlabel('推定係数', fontsize=12)
    ax3.set_title('固定効果モデル（FE）の推定係数と95%信頼区間', fontsize=13, fontweight='bold')

ax3.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_U2_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("図3 保存完了")

df_fe = df_panel.dropna(subset=['消費支出', '高齢化率', '求人倍率', '保健医療費_人', 'log_人口', '高齢化率2']).copy()

fe_results = None
fe_params = None
fe_conf_int = None

try:
    from linearmodels.panel import PanelOLS
    df_fe = df_fe.set_index(['都道府県', '年度'])
    y = df_fe['消費支出']
    X = sm.add_constant(df_fe[['高齢化率', '高齢化率2', '求人倍率', '保健医療費_人', 'log_人口']])
    model = PanelOLS(y, X, entity_effects=True)
    fe_results = model.fit(cov_type='clustered', cluster_entity=True)
    print("=== PanelOLS 固定効果モデル結果 ===")
    print(fe_results.summary)
    fe_params = fe_results.params
    fe_conf_int = fe_results.conf_int()
    print()
    print("係数:")
    print(fe_params)
except Exception as e:
    print(f"PanelOLS 失敗 ({e})、OLS にフォールバック")
    df_fe_ols = df_panel.dropna(subset=['消費支出', '高齢化率', '求人倍率', '保健医療費_人', 'log_人口', '高齢化率2']).copy()
    X_ols = sm.add_constant(df_fe_ols[['高齢化率', '高齢化率2', '求人倍率', '保健医療費_人', 'log_人口']])
    y_ols = df_fe_ols['消費支出']
    ols_model = sm.OLS(y_ols, X_ols).fit()
    print("=== OLS フォールバック結果 ===")
    print(ols_model.summary())
    fe_params = ols_model.params
    fe_conf_int = ols_model.conf_int()

# 係数の情報を取得
if fe_params is not None:
    param_names_en = {'const': '定数項', '高齢化率': '高齢化率', '高齢化率2': '高齢化率²',
                      '求人倍率': '求人倍率', '保健医療費_人': '保健医療費/人', 'log_人口': 'log(人口)'}
    param_display = {param_names_en.get(k, k): v for k, v in fe_params.items() if k != 'const'}

aging_2023 = df_panel[df_panel['年度'] == 2023][['都道府県', '高齢化率']].sort_values('高齢化率', ascending=False)
top5_prefs = aging_2023.head(5)['都道府県'].tolist()
bot5_prefs = aging_2023.tail(5)['都道府県'].tolist()
print(f"高齢化率上位5: {top5_prefs}")
print(f"高齢化率下位5: {bot5_prefs}")

df_top = df_panel[df_panel['都道府県'].isin(top5_prefs)].groupby('年度')['消費支出'].mean() / 1000
df_bot = df_panel[df_panel['都道府県'].isin(bot5_prefs)].groupby('年度')['消費支出'].mean() / 1000

fig4, ax4 = plt.subplots(figsize=(10, 6))
ax4.plot(df_top.index, df_top.values, marker='o', markersize=6,
         color='#e05c5c', linewidth=2.5,
         label=f'高齢化率上位5県\n（{", ".join([p.replace("県","").replace("都","").replace("府","").replace("道","") for p in top5_prefs])}）')
ax4.plot(df_bot.index, df_bot.values, marker='s', markersize=6,
         color='#4e9af1', linewidth=2.5,
         label=f'高齢化率下位5県\n（{", ".join([p.replace("県","").replace("都","").replace("府","").replace("道","") for p in bot5_prefs])}）')

ax4.fill_between(df_top.index, df_top.values, df_bot.values, alpha=0.08, color='gray')
ax4.set_xlabel('年度', fontsize=12)
ax4.set_ylabel('平均消費支出（千円）', fontsize=12)
ax4.set_title('高齢化率 上位群 vs 下位群 の消費支出比較（2012〜2023年）', fontsize=13, fontweight='bold')
ax4.legend(loc='upper left', fontsize=9)
ax4.set_xticks(range(2012, 2024))
ax4.grid(True, alpha=0.3)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2019_U2_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("図4 保存完了")