ふるさと納税は地方創生の切り札になりえるか

# ================================================================
# ■ データ読み込み（市区町村パネル）
# ================================================================
df = pd.read_csv(DATA)
df = df.dropna(subset=['直近3年受入額_百万円', '人口増減率', '10万人当たり医師数'])

print("=" * 74)
print("■ 市区町村パネル（ふるさと納税 × 人口動態）")
print("  ふるさと納税: 総務省『ふるさと納税に関する現況調査』団体別一覧（原論文と同じ出典）")
print("  人口動態    : SSDSE-A 2022〜2026 年版（版ごとの収録年次を重ねてパネル化）")
print("=" * 74)
print(f"  市区町村数: {df['地域コード'].nunique()} 団体")
print(f"  年        : {sorted(df['年'].unique())}")
print(f"  観測数    : {len(df)} 行")

# 受入額は分布が大きく歪むので対数を取る（原論文も金額そのものではなく規模で効かせている）
df['log受入額'] = np.log1p(df['直近3年受入額_百万円'])

# ================================================================
# ■ 記述統計（原論文 表3 に対応）
# ================================================================
VARS = ['直近3年受入額_百万円', '自然増減率', '社会増減率', '人口増減率', '10万人当たり医師数']
print("\n【基本統計量】")
desc = df[VARS].describe().T[['count', 'mean', 'std', 'min', '50%', 'max']]
desc.columns = ['度数', '平均', '標準偏差', '最小', '中央', '最大']
print(desc.round(2).to_string())

# ================================================================
# ■ Step1. 固定効果モデルとは何をしているのか（within 変換を手で作る）
# ================================================================
print("\n" + "=" * 74)
print("■ Step1. 固定効果モデル＝『各市区町村の平均からのズレ』で見る")
print("=" * 74)
print("  同じ市区町村の中での動きだけを使うので、")
print("  「もともと魅力的な自治体だから寄付も移住も多い」といった、")
print("  時間で変わらない性質（＝欠落変数）の影響を丸ごと消せる。")


def within(frame, cols, group='地域コード'):
    """各グループ（市区町村）の平均を引く。これが固定効果モデルの中身。"""
    out = frame.copy()
    for c in cols:
        out[c + '_w'] = out[c] - out.groupby(group)[c].transform('mean')
    return out


TARGETS = ['自然増減率', '社会増減率', '人口増減率']
XS = ['log受入額', '10万人当たり医師数']
dw = within(df, TARGETS + XS)

sample = dw[dw['地域コード'] == dw['地域コード'].iloc[0]]
print(f"\n  例：{sample['都道府県'].iloc[0]}{sample['市区町村'].iloc[0]} の人口増減率")
for _, r in sample.iterrows():
    print(f"    {int(r['年'])}年  実測 {r['人口増減率']:+.3f} %"
          f"  平均からのズレ {r['人口増減率_w']:+.3f} %")

# ================================================================
# ■ Step2. 固定効果回帰（原論文 表5 に対応する、人口側の 3 本）
# ================================================================
print("\n" + "=" * 74)
print("■ Step2. 固定効果モデルの推定（被説明変数を 3 通り）")
print("=" * 74)


def fe_ols(frame, y, xs, group='地域コード'):
    """within 変換したデータに最小二乗法をあてる（＝固定効果推定）。

    標準誤差は市区町村ごとにクラスターして、同じ自治体の観測が
    似ていることを考慮する。
    """
    d = frame.dropna(subset=[y + '_w'] + [x + '_w' for x in xs])
    Y = d[y + '_w'].to_numpy()
    X = np.column_stack([d[x + '_w'].to_numpy() for x in xs])
    beta, *_ = np.linalg.lstsq(X, Y, rcond=None)
    resid = Y - X @ beta
    n, k = X.shape
    n_g = d[group].nunique()
    XtX_inv = np.linalg.inv(X.T @ X)
    meat = np.zeros((k, k))
    for _g, idx in d.groupby(group).indices.items():
        Xg, ug = X[idx], resid[idx]
        s = Xg.T @ ug
        meat += np.outer(s, s)
    dof = n - n_g - k
    scale = (n_g / max(n_g - 1, 1)) * ((n - 1) / max(dof, 1))
    cov = XtX_inv @ meat @ XtX_inv * scale
    se = np.sqrt(np.diag(cov))
    t = beta / se
    from scipy import stats as st
    p = 2 * (1 - st.t.cdf(np.abs(t), df=max(dof, 1)))
    ss_res = float(resid @ resid)
    ss_tot = float(((Y - Y.mean()) ** 2).sum())
    return pd.DataFrame({'係数': beta, '標準誤差': se, 't値': t, 'p値': p},
                        index=xs), 1 - ss_res / ss_tot, n, n_g


results = {}
for y in TARGETS:
    tbl, r2, n, ng = fe_ols(dw, y, XS)
    results[y] = tbl
    print(f"\n【被説明変数: {y}】 観測 {n} / 市区町村 {ng} / within R² = {r2:.4f}")
    show = tbl.copy()
    show['判定'] = ['***' if p < .01 else '**' if p < .05 else '*' if p < .1 else 'n.s.'
                    for p in show['p値']]
    print(show.round(4).to_string())

# ================================================================
# ■ Step3. 原論文の結論と照らし合わせる
# ================================================================
print("\n" + "=" * 74)
print("■ Step3. 原論文の結論と、この再現の結果")
print("=" * 74)
b_soc = results['社会増減率'].loc['log受入額']
b_nat = results['自然増減率'].loc['log受入額']
b_pop = results['人口増減率'].loc['log受入額']
print("  原論文（2012/2014/2016）:")
print("    ・社会増減率  … 関係が確認されなかった")
print("    ・自然増減率  … 負の効果")
print("    ・結論        … ふるさと納税は人口増加をもたらさない")
print("\n  この再現（2020〜2023）:")
for name, b in [('社会増減率', b_soc), ('自然増減率', b_nat), ('人口増減率', b_pop)]:
    verdict = ('有意な関係は見られない' if b['p値'] >= .05
               else ('正の効果' if b['係数'] > 0 else '負の効果'))
    print(f"    ・{name}  … 係数 {b['係数']:+.4f} (p={b['p値']:.3f}) → {verdict}")
print("\n  → 『寄付が集まっても人口増加には結びつかない』という原論文の主旨が、")
print("     時点を変えても同じ向きで確認できるかどうかが読みどころ。")

# ================================================================
# ■ 図1: ふるさと納税受入額の推移（原論文 図1 に対応）
# ================================================================
fig, ax = plt.subplots(figsize=(9, 5))
tot = (df.groupby('年')['ふるさと納税受入額_千円'].sum() / 1e6)
ax.bar(tot.index.astype(int), tot.values, color='#1565C0')
for x, v in zip(tot.index.astype(int), tot.values):
    ax.text(x, v, f'{v:,.0f}', ha='center', va='bottom', fontsize=11)
ax.set_xlabel('年度')
ax.set_ylabel('ふるさと納税受入額の合計（十億円）')
ax.set_title('図1: ふるさと納税受入額の推移（分析対象の市区町村合計）\n'
             '出典: 総務省「ふるさと納税に関する現況調査」', fontsize=13)
ax.grid(axis='y', alpha=.3)
fig.savefig(os.path.join(FIG_DIR, '2020_U2_furusato_fig1_trend.png'),
            dpi=DPI, bbox_inches='tight')
plt.close(fig)
print("\n  → html/figures/2020_U2_furusato_fig1_trend.png 保存完了")

# ================================================================
# ■ 図2: 固定効果の考え方（実測値と、平均からのズレ）
# ================================================================
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
top = (df.groupby('地域コード')['直近3年受入額_百万円'].mean()
       .sort_values(ascending=False).head(6).index)
sub = df[df['地域コード'].isin(top)]
subw = within(sub, ['人口増減率'])
for code, g in subw.groupby('地域コード'):
    label = f"{g['都道府県'].iloc[0]}{g['市区町村'].iloc[0]}"
    axes[0].plot(g['年'], g['人口増減率'], marker='o', label=label)
    axes[1].plot(g['年'], g['人口増減率_w'], marker='o', label=label)
axes[0].set_title('実測の人口増減率\n（自治体ごとの水準の違いが大きい）', fontsize=12)
axes[1].set_title('自治体の平均を引いたあと\n（固定効果モデルが見ているのはこちら）', fontsize=12)
for ax in axes:
    ax.set_xlabel('年')
    ax.set_ylabel('人口増減率（％）')
    ax.grid(alpha=.3)
    ax.legend(fontsize=9)
axes[1].axhline(0, color='#C62828', lw=1.2, ls='--')
fig.suptitle('図2: 固定効果モデルは「その自治体の中での動き」だけを使う', fontsize=14)
fig.savefig(os.path.join(FIG_DIR, '2020_U2_furusato_fig2_within.png'),
            dpi=DPI, bbox_inches='tight')
plt.close(fig)
print("  → html/figures/2020_U2_furusato_fig2_within.png 保存完了")

# ================================================================
# ■ 図3: 推定結果（係数と 95% 信頼区間）
# ================================================================
fig, ax = plt.subplots(figsize=(9, 5))
ys = list(results.keys())
coef = [results[y].loc['log受入額', '係数'] for y in ys]
se = [results[y].loc['log受入額', '標準誤差'] for y in ys]
pos = np.arange(len(ys))
ax.errorbar(coef, pos, xerr=[1.96 * s for s in se], fmt='o', color='#1565C0',
            capsize=6, markersize=9, lw=2)
ax.axvline(0, color='#C62828', ls='--', lw=1.4)
ax.set_yticks(pos)
ax.set_yticklabels(ys, fontsize=12)
ax.set_xlabel('ふるさと納税受入額（対数）の係数と 95% 信頼区間')
ax.set_title('図3: ふるさと納税は人口を増やしたか（固定効果モデル）\n'
             '横棒が 0 をまたいでいれば「効果があるとは言えない」', fontsize=13)
ax.grid(axis='x', alpha=.3)
fig.savefig(os.path.join(FIG_DIR, '2020_U2_furusato_fig3_coef.png'),
            dpi=DPI, bbox_inches='tight')
plt.close(fig)
print("  → html/figures/2020_U2_furusato_fig3_coef.png 保存完了")

# ================================================================
# ■ 再現できていない部分（隠さずに書く）
# ================================================================
print("\n" + "=" * 74)
print("■ この再現でできていないこと")
print("=" * 74)
print("  1. 仮説1（経済効果）の被説明変数『納税義務者1人当たり課税対象所得』は")
print("     SSDSE に収録が無く、この教材では検証できていない。")
print("     取得先: 総務省「市町村税課税状況等の調」")
print("  2. 時点が原論文（2012/2014/2016）ではなく 2020〜2023 年。")
print("     SSDSE-A の版を重ねて作れるのがこの範囲のため。")
print("  3. 大学進学率・保育施設利用者比率・高齢者施設定員比率は、")
print("     市区町村パネルとして揃わないため統制変数から外している。")
print("=" * 74)