少子化問題に向き合う

raw = pd.read_csv(DATA_B, encoding='cp932', skiprows=[1])

def make_frame(year):
    """指定年の47都道府県について、少子化分析用の指標を計算して返す。"""
    d = raw[raw['SSDSE-B-2026'] == year].copy()
    d['births_p1000']   = d['A4101'] / d['A1101'] * 1000          # 人口千人あたり出生数
    d['marr_p1000']     = d['A9101'] / d['A1101'] * 1000          # 人口千人あたり婚姻件数
    d['netmarr_p1000']  = (d['A9101'] - d['A9201']) / d['A1101'] * 1000  # 千人あたり純結婚数
    d['tfr']            = d['A4103']                              # 合計特殊出生率(TFR)
    d['edu']            = d['L322108']                           # 教育費(二人以上の世帯)
    d['childcare_p1000'] = (d['E1101'] + d['J2503']) / d['A1101'] * 1000  # 千人あたり(幼稚園+保育所等)
    return d

d = make_frame(YEAR)
print(f'SSDSE-B-2026 読み込み: N = {len(d)} 都道府県（断面年: {YEAR}年 ＝ 原論文と同じ）')
print(f'人口千人あたり出生数: 平均 {d["births_p1000"].mean():.2f} 人'
      f'（{d["births_p1000"].min():.2f}〜{d["births_p1000"].max():.2f}）')
print(f'合計特殊出生率(TFR): 平均 {d["tfr"].mean():.3f}'
      f'（{d["tfr"].min():.2f}〜{d["tfr"].max():.2f}）')
print('原論文の報告値（参考）: 人口置換水準 TFR=2.07、2019年の全国TFR=1.34')

# ══════════════════════════════════════════════════════════════
# 図1【実再現】千人あたり婚姻件数 × TFR / 千人あたり出生数（原論文 図[4][5]）
#   原論文の主張: TFRは婚姻件数とあまり相関しない(報告値 r≒0.26)が、

r_tfr,  p_tfr  = stats.pearsonr(d['marr_p1000'], d['tfr'])
r_bir,  p_bir  = stats.pearsonr(d['marr_p1000'], d['births_p1000'])
print(f'千人あたり婚姻件数 × TFR          : r = {r_tfr:+.4f} (p = {p_tfr:.3f})'
      f'   ← 原論文の報告値 r≒0.26')
print(f'千人あたり婚姻件数 × 千人あたり出生数: r = {r_bir:+.4f} (p = {p_bir:.1e})'
      f'   ← 原論文の報告値: 強い正の相関')
print('→ TFRとの相関は弱く、出生数の方が婚姻件数と強く相関する、という原論文の判断を再現。')

years_main = list(range(2012, 2019))          # 原論文 図[10] と同じ 2012〜2018
years_ext  = list(range(2019, 2024))          # 発展: 収録データで直近まで延長
def corr_year(year):
    dy = make_frame(year)
    return dy['netmarr_p1000'].corr(dy['births_p1000'])
r_main = [corr_year(y) for y in years_main]
r_ext  = [corr_year(y) for y in years_ext]
print(f"{'年度':<6}{'相関係数(実再現)':>16}")
for y, r in zip(years_main, r_main):
    print(f'{y:<6}{r:>16.4f}')
print('→ 原論文の報告どおり 0.5前後の正の相関で推移し、全体として弱まる傾向を再現。')
print('（発展）2019〜2023: ' + ', '.join(f'{y}:{r:.3f}' for y, r in zip(years_ext, r_ext)))

r_cc,  p_cc  = stats.pearsonr(d['childcare_p1000'], d['births_p1000'])
d_ex = d[~d['Prefecture'].isin(EXCLUDE)]
r_ccx, _     = stats.pearsonr(d_ex['childcare_p1000'], d_ex['births_p1000'])

REPORTED_EDU = {
    '私立中学校数 × 教育費（平均未満）':  0.079,   # 図[6] 相関なし
    'TFR × 教育費（平均未満）':          -0.098,   # 図[7] 相関なし
    '私立中学校数 × 教育費（平均以上）':  0.45,    # 図[8] 正の相関
    'TFR × 教育費（平均以上）':          -0.53,    # 図[9] 負の相関
}
for k, v in REPORTED_EDU.items():
    print(f'  {k:<28} r = {v:+.3f}')
print('  ※ 私立中学校数の平均は全国で約16.6校。平均以上の県では教育費が高く TFR は下がる。')

print()

r_eduT, p_eduT = stats.pearsonr(d['edu'], d['tfr'])
r_eduB, p_eduB = stats.pearsonr(d['edu'], d['births_p1000'])
print(f'教育費 × TFR            : r = {r_eduT:+.4f} (p = {p_eduT:.3f})')
print(f'教育費 × 千人あたり出生数: r = {r_eduB:+.4f} (p = {p_eduB:.3f})')
print('→ 全県でも教育費とTFRは負の方向。原論文の「教育費が高いほどTFRが下がる」と整合。')

REPORTED_CARE = {'育児時間（男女合計）': 0.528, '育児時間（女性）': 0.38, '育児時間（男性）': 0.02}
for k, v in REPORTED_CARE.items():
    print(f'  {k:<16} r = {v:+.3f}')
print('  ※ 男女合計は弱い正の相関。内訳は女性が弱い正・男性はほとんど相関なし。')
print('  → 男女間の育児時間の格差是正が課題、という原論文の主張の根拠。')