地方移住の決定要因に関する
パネルデータ分析
―ラグ効果の検証を通じて―

import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────
DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family': 'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi': 150,
    'axes.spines.top': False,
    'axes.spines.right': False,
})

# ── SSDSE-B-2026.csv 読み込み ──────────────────────────────────────
print("=== SSDSE-B 読み込み ===")
df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=1
)

# 47都道府県のみ（地域コード = R + 5桁数字）
df_raw = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_raw['年度'] = pd.to_numeric(df_raw['年度'], errors='coerce')

# 数値型に変換
num_cols = [
    '総人口', '65歳以上人口', '15～64歳人口', '15歳未満人口',
    '転入者数（日本人移動者）', '転入者数（日本人移動者）（男）',
    '転出者数（日本人移動者）',
    '保育所等数', '年平均気温', '婚姻件数',
    '高等学校卒業者数', '高等学校卒業者のうち進学者数',
    '出生数', '死亡数',
    '教育費（二人以上の世帯）',
]
for c in num_cols:
    df_raw[c] = pd.to_numeric(df_raw[c], errors='coerce')

df_raw = df_raw.sort_values(['都道府県', '年度']).reset_index(drop=True)
print(f"SSDSE-B 読み込み完了: {len(df_raw)}行 ({df_raw['都道府県'].nunique()}都道府県, "
      f"年度: {sorted(df_raw['年度'].dropna().astype(int).unique())})")

# ── SSDSE-E-2026.csv 読み込み（面積・県民所得） ─────────────────────
print("=== SSDSE-E 読み込み ===")
df_e_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
    encoding='cp932', header=1
)
df_e = df_e_raw.iloc[1:].copy()
df_e.columns = df_e_raw.iloc[0].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)

# 面積（ha → km²）と1人当たり県民所得
df_e['面積_km2'] = pd.to_numeric(df_e['総面積（北方地域及び竹島を除く）'], errors='coerce') / 100.0
df_e['県民所得_1人'] = pd.to_numeric(df_e['1人当たり県民所得（平成27年基準）'], errors='coerce')
df_cross = df_e[['都道府県', '面積_km2', '県民所得_1人']].copy()
print(f"SSDSE-E 読み込み完了: {len(df_cross)}都道府県")

# ── パネルデータ構築 ──────────────────────────────────────────────
panel = df_raw[['都道府県', '年度',
                '総人口', '65歳以上人口', '15～64歳人口',
                '転入者数（日本人移動者）', '転入者数（日本人移動者）（男）',
                '転出者数（日本人移動者）',
                '保育所等数', '年平均気温', '婚姻件数',
                '高等学校卒業者数', '高等学校卒業者のうち進学者数',
                '教育費（二人以上の世帯）',
                ]].copy()
panel.columns = ['pref', 'year', 'pop', 'pop65', 'pop1564',
                 'inflow', 'inflow_m', 'outflow',
                 'nursery', 'avg_temp', 'marriages',
                 'hs_grads', 'hs_college', 'edu_exp']

# SSDSE-E をマージ（面積のみを使用: 人口密度の計算のため）
panel = panel.merge(df_cross[['都道府県', '面積_km2']].rename(columns={'都道府県': 'pref'}),
                    on='pref', how='left')

# ── 派生変数の計算（実データのみ） ─────────────────────────────────
# 保育所充実度: 保育所等数 / 総人口 × 10000
panel['nursery_rate']   = panel['nursery'] / panel['pop'] * 10000

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

# 婚姻率: 婚姻件数 / 15〜64歳人口 × 1000
panel['marriage_rate']  = panel['marriages'] / panel['pop1564'] * 1000

# 大学進学率: 高校卒業者のうち進学者数 / 高校卒業者数 × 100
panel['college_rate']   = panel['hs_college'] / panel['hs_grads'] * 100

# 教育費（二人以上の世帯）: 時間変動する所得代理変数（SSDSE-B実データ）
# 注: 1人当たり県民所得（SSDSE-E）は時点固定のため固定効果推定では識別不能
panel['edu_expense']    = panel['edu_exp']

# 人口密度: 総人口 / 面積(km²)
panel['pop_density']    = panel['pop'] / panel['面積_km2']

# 目的変数: 転入者数（人）※原論文に合わせた線形
panel = panel.sort_values(['pref', 'year']).reset_index(drop=True)

# ── ラグ変数の作成 ──────────────────────────────────────────────────
LAG_VARS = ['nursery_rate', 'aging_rate', 'avg_temp', 'marriage_rate',
            'college_rate', 'edu_expense', 'pop_density']

for var in LAG_VARS:
    panel[f'{var}_lag1'] = panel.groupby('pref')[var].shift(1)
    panel[f'{var}_lag2'] = panel.groupby('pref')[var].shift(2)

# ラグ2期まで使えるデータ（2016年〜）
panel_est = panel[panel['year'] >= 2016].dropna(
    subset=['inflow_m'] + LAG_VARS +
           [f'{v}_lag1' for v in LAG_VARS] +
           [f'{v}_lag2' for v in LAG_VARS]
).copy()
print(f"\n推定用サンプル: {len(panel_est)}観測 "
      f"({panel_est['pref'].nunique()}都道府県, "
      f"年度: {sorted(panel_est['year'].unique())})")

# ── 固定効果パネル回帰（within変換） ───────────────────────────────
def within_transform(df, dep_var, expl_vars, id_col='pref'):
    """within推定量（固定効果モデル）: 個体内平均を除去"""
    df2 = df.copy()
    all_vars = [dep_var] + expl_vars
    group_means = df2.groupby(id_col)[all_vars].transform('mean')
    for v in all_vars:
        df2[v + '_dm'] = df2[v] - group_means[v]
    y = df2[dep_var + '_dm'].values
    X = sm.add_constant(df2[[v + '_dm' for v in expl_vars]].values)
    res = sm.OLS(y, X).fit(cov_type='HC1')
    return res

EXPL_CURRENT = LAG_VARS
EXPL_LAG1    = [v + '_lag1' for v in LAG_VARS]
EXPL_LAG2    = [v + '_lag2' for v in LAG_VARS]

res_t0 = within_transform(panel_est, 'inflow_m', EXPL_CURRENT)
res_t1 = within_transform(panel_est, 'inflow_m', EXPL_LAG1)
res_t2 = within_transform(panel_est, 'inflow_m', EXPL_LAG2)

print("\n=== 固定効果モデル（当期, 男性転入者数）===")
tbl = res_t0.summary2().tables[1]
cols_needed = [c for c in ['Coef.', 'Std.Err.', 't', 'P>|t|'] if c in tbl.columns]
print(tbl[cols_needed].to_string())
print(f"\nR² (within) = {res_t0.rsquared:.4f}")

# ── Figure 1: 転入者数の時系列推移 ───────────────────────────────
print("\nFigure 1: 転入者数の時系列推移...")
YEARS_ALL = sorted(panel['year'].dropna().astype(int).unique())

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 全国合計転入者数
nat = panel.groupby('year')['inflow'].sum() / 1e4  # 万人

ax = axes[0]
ax.plot(nat.index, nat.values, 'o-', color='#1565C0', lw=2.5)
ax.fill_between(nat.index, nat.values, nat.values.min(),
                alpha=0.15, color='#1565C0')
ax.set_xlabel('年度')
ax.set_ylabel('全国転入者数（万人）')
ax.set_title('全国転入者数（日本人）の推移', fontsize=12, fontweight='bold')
ax.set_xticks(YEARS_ALL)
ax.set_xticklabels([str(y) for y in YEARS_ALL], rotation=45)

# 地方主要県の推移
ax2 = axes[1]
highlight_prefs = ['鳥取県', '島根県', '長野県', '宮崎県']
colors_h = ['#E53935', '#8E24AA', '#1E88E5', '#43A047']
for pref, col in zip(highlight_prefs, colors_h):
    sub = panel[panel['pref'] == pref].sort_values('year')
    ax2.plot(sub['year'], sub['inflow'] / 1000, 'o-',
             color=col, lw=2, label=pref)
ax2.set_xlabel('年度')
ax2.set_ylabel('転入者数（千人）')
ax2.set_title('地方主要県の転入者数推移', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10, ncol=2)
ax2.set_xticks(YEARS_ALL)
ax2.set_xticklabels([str(y) for y in YEARS_ALL], rotation=45)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H3_fig1_trend.png'), bbox_inches='tight')
plt.close()

print("Figure 1 saved.")

# ── Figure 2: 固定効果モデル係数（当期） ────────────────────────────
print("Figure 2: 固定効果モデル係数（当期）...")
fig, ax = plt.subplots(figsize=(10, 5))

var_labels = ['保育所充実度\n(保育所数/人口万)', '高齢化率\n(%)',
              '年平均気温\n(℃)', '婚姻率\n(婚姻数/15-64歳人口千)',
              '大学進学率\n(%)', '教育費\n(二人以上世帯, 円)',
              '人口密度\n(人/km²)']
coefs = res_t0.params[1:]
cis   = np.asarray(res_t0.conf_int())[1:]
pvals = res_t0.pvalues[1:]

colors_bar = ['#E53935' if c < 0 else '#1E88E5' for c in coefs]
sig_marks  = ['***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
              for p in pvals]

y_pos = np.arange(len(var_labels))
ax.barh(y_pos, coefs, color=colors_bar, alpha=0.8, height=0.6)
ax.errorbar(coefs, y_pos, xerr=[coefs - cis[:, 0], cis[:, 1] - coefs],
            fmt='none', color='black', capsize=5, lw=1.5)
ax.axvline(0, color='black', lw=1)
ax.set_yticks(y_pos)
ax.set_yticklabels([f'{l} {m}' for l, m in zip(var_labels, sig_marks)], fontsize=10)
ax.set_xlabel('固定効果推定量（男性転入者数, 人）')
ax.set_title('固定効果モデル推定結果（当期）\n男性転入者数への影響',
             fontsize=12, fontweight='bold')

neg_patch = mpatches.Patch(color='#E53935', alpha=0.8, label='負の効果')
pos_patch = mpatches.Patch(color='#1E88E5', alpha=0.8, label='正の効果')
ax.legend(handles=[neg_patch, pos_patch], fontsize=10)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H3_fig2_fe.png'), bbox_inches='tight')
plt.close()

print("Figure 2 saved.")

# ── Figure 3: ラグ効果の係数推移 ────────────────────────────────────
print("Figure 3: ラグ効果の係数推移...")
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

lag_labels = ['当期 (t)', '1期ラグ (t-1)', '2期ラグ (t-2)']

def plot_lag_coef(ax, idx, title_var, title_text, results_list):
    """ラグ別係数をバーチャートで描画"""
    coef_list = [r.params[idx] for r in results_list]
    ci_list   = [np.asarray(r.conf_int())[idx] for r in results_list]
    pval_list = [r.pvalues[idx] for r in results_list]
    bar_colors = ['#E53935' if c < 0 else '#1E88E5' for c in coef_list]
    ax.bar(range(3), coef_list, color=bar_colors, alpha=0.85, width=0.5)
    for i, (ci, pv) in enumerate(zip(ci_list, pval_list)):
        ax.errorbar(i, coef_list[i],
                    yerr=[[coef_list[i] - ci[0]], [ci[1] - coef_list[i]]],
                    fmt='none', color='black', capsize=6, lw=2)
        if pv < 0.05:
            ax.text(i, max(coef_list[i], 0) + abs(ci[1] - coef_list[i]) + abs(coef_list[i]) * 0.05,
                    '***' if pv < 0.001 else ('**' if pv < 0.01 else '*'),
                    ha='center', fontsize=12, color='red')
    ax.axhline(0, color='black', lw=1.2)
    ax.set_xticks(range(3))
    ax.set_xticklabels(lag_labels, fontsize=11)
    ax.set_ylabel('推定係数（男性転入者数, 人）')
    ax.set_title(title_text, fontsize=12, fontweight='bold')
    return coef_list

# 左: 保育所充実度のラグ別係数（idx=1: 第1説明変数）
c1 = plot_lag_coef(axes[0], 1, '保育所充実度',
                   '保育所充実度のラグ別係数\n（保育所数/人口万）',
                   [res_t0, res_t1, res_t2])

# 右: 高齢化率のラグ別係数（idx=2）
c2 = plot_lag_coef(axes[1], 2, '高齢化率',
                   '高齢化率のラグ別係数\n（65歳以上人口/総人口 × 100）',
                   [res_t0, res_t1, res_t2])

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H3_fig3_lag.png'), bbox_inches='tight')
plt.close()

print("Figure 3 saved.")

# ── Figure 4: VIF確認 ─────────────────────────────────────────────
print("Figure 4: VIF確認...")
fig, ax = plt.subplots(figsize=(10, 5))

expl_data = panel_est[EXPL_CURRENT].dropna()
X_vif = expl_data.values
vif_values = []
for i in range(X_vif.shape[1]):
    vif_values.append(variance_inflation_factor(X_vif, i))

var_labels_vif = ['保育所充実度', '高齢化率', '年平均気温',
                  '婚姻率', '大学進学率', '教育費(二人以上世帯)', '人口密度']
colors_vif = ['#E53935' if v >= 10 else ('#FB8C00' if v >= 5 else '#43A047')
              for v in vif_values]

bars = ax.barh(range(len(var_labels_vif)), vif_values, color=colors_vif, alpha=0.85)
ax.axvline(10, color='#E53935', ls='--', lw=1.5, label='VIF=10 (問題あり)')
ax.axvline(5,  color='#FB8C00', ls='--', lw=1.5, label='VIF=5 (注意)')
ax.set_yticks(range(len(var_labels_vif)))
ax.set_yticklabels(var_labels_vif, fontsize=11)
ax.set_xlabel('VIF値')
ax.set_title('多重共線性の確認（VIF）', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)

for bar, val in zip(bars, vif_values):
    ax.text(bar.get_width() + 0.1, bar.get_y() + bar.get_height() / 2,
            f'{val:.2f}', va='center', fontsize=10)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2025_H3_fig4_vif.png'), bbox_inches='tight')
plt.close()

print("Figure 4 saved.")

print("\n=== ラグ効果まとめ（保育所充実度係数） ===")
for i, (res, lag) in enumerate(zip([res_t0, res_t1, res_t2], ['当期', 't-1期', 't-2期'])):
    coef = res.params[1]
    pval = res.pvalues[1]
    sig = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.'))
    print(f"  {lag}: beta={coef:.2f}, p={pval:.3f} {sig}")

print("\n=== ラグ効果まとめ（高齢化率係数） ===")
for i, (res, lag) in enumerate(zip([res_t0, res_t1, res_t2], ['当期', 't-1期', 't-2期'])):
    coef = res.params[2]
    pval = res.pvalues[2]
    sig = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.'))
    print(f"  {lag}: beta={coef:.2f}, p={pval:.3f} {sig}")

print("\n分析完了。html/figures/ に図を保存しました。")
print(f"  2025_H3_fig1_trend.png - 転入者数の時系列推移")
print(f"  2025_H3_fig2_fe.png    - 固定効果モデル推定結果")
print(f"  2025_H3_fig3_lag.png   - ラグ効果の係数推移")
print(f"  2025_H3_fig4_vif.png   - VIF確認")