空き家を減らすために

# ============================================================
# STEP1: SSDSE-E（都道府県データ）を読み込み、47都道府県を抽出
#   高齢者率 = 65歳以上人口 / 総人口   × 100
#   空き家率 = 空き家数     / 総住宅数 × 100  （＝全空き家の率）
# ============================================================
print("STEP1: SSDSE-E を読み込み中 ...")
raw_e = pd.read_csv(DATA_E, encoding='cp932', header=0)
codes = list(raw_e.columns)
names = raw_e.iloc[1].tolist()            # 2行目 = 日本語の変数名
code_of = dict(zip(names, codes))         # 変数名 -> 変数コード

e = raw_e.iloc[2:].copy()                 # 3行目以降が実データ
e = e[e['SSDSE-E-2026'].astype(str).str.match(r'^R\d{5}$')]
e = e[e['SSDSE-E-2026'] != 'R00000']      # 全国計を除外

def col_e(jp):
    return pd.to_numeric(e[code_of[jp]], errors='coerce')

df_e = pd.DataFrame({
    '都道府県':  e['Prefecture'].values,
    '高齢者率':  (col_e('65歳以上人口') / col_e('総人口')   * 100).values,
    '空き家率':  (col_e('空き家数')     / col_e('総住宅数') * 100).values,
})
print(f"  抽出: {len(df_e)} 都道府県")
print(df_e.head(5).round(2).to_string(index=False))

# ============================================================
# STEP2: SSDSE-B から住宅地の地価（標準価格）を取得
# ============================================================
print("STEP2: SSDSE-B（地価）を読み込み中 ...")
b = pd.read_csv(DATA_B, encoding='cp932', header=1)
b = b[b['地域コード'].astype(str).str.match(r'^R\d{5}')].copy()
b['年度'] = pd.to_numeric(b['年度'], errors='coerce')
latest = int(b['年度'].max())
b = b[b['年度'] == latest]
b['地価'] = pd.to_numeric(b['標準価格（平均価格）（住宅地）'], errors='coerce')
df_b = b[['都道府県', '地価']].reset_index(drop=True)
print(f"  対象年度: {latest} 年度 / {len(df_b)} 都道府県")
print(f"  地価レンジ: {df_b['地価'].min():,.0f} 〜 {df_b['地価'].max():,.0f} 円/m^2")

# ============================================================
# STEP3: SSDSE-A（市区町村）を都道府県ごとに集計し、
#        「高齢者のみ世帯」割合を算出
#   高齢者のみ世帯 = 65歳以上の単独世帯 + 夫65歳以上・妻60歳以上の夫婦のみ世帯
# ============================================================
print("STEP3: SSDSE-A から高齢者のみ世帯割合を集計中 ...")
raw_a  = pd.read_csv(DATA_A, encoding='cp932', header=0)
a_name = raw_a.iloc[1].tolist()
a_code = dict(zip(a_name, list(raw_a.columns)))
a = raw_a.iloc[2:].copy()

def col_a(jp):
    return pd.to_numeric(a[a_code[jp]], errors='coerce')

a['一般世帯数'] = col_a('一般世帯数')
a['高齢単独']   = col_a('65歳以上世帯員の単独世帯数')
a['高齢夫婦']   = col_a('夫65歳以上、妻60歳以上の夫婦のみの世帯数')
g = a.groupby('Prefecture').agg(
        一般世帯数=('一般世帯数', 'sum'),
        高齢単独=('高齢単独', 'sum'),
        高齢夫婦=('高齢夫婦', 'sum')).reset_index()
g['高齢者のみ世帯割合'] = (g['高齢単独'] + g['高齢夫婦']) / g['一般世帯数'] * 100
df_a = g[['Prefecture', '高齢者のみ世帯割合']].rename(columns={'Prefecture': '都道府県'})
print(f"  {len(df_a)} 都道府県 / 高齢者のみ世帯割合 レンジ: "
      f"{df_a['高齢者のみ世帯割合'].min():.1f}% 〜 {df_a['高齢者のみ世帯割合'].max():.1f}%")

# ============================================================
# STEP4: 3つのデータを都道府県名で結合（47都道府県のパネル）
# ============================================================
df = df_e.merge(df_b, on='都道府県').merge(df_a, on='都道府県')
print(f"STEP4: 結合完了 -> {len(df)} 都道府県 × {df.shape[1]} 変数")
print(df.head().round(2).to_string(index=False))

# ============================================================
# STEP5: 相関分析（全空き家率ベース）
# ============================================================
def report(x, y, label):
    r, p = stats.pearsonr(df[x], df[y])
    print(f"  {label:26s}: r = {r:+.3f}  (p = {p:.4f})")

print("STEP5: 相関分析")
report('高齢者率',           '空き家率', '高齢者率 × 空き家率')
report('地価',               '空き家率', '地価 × 空き家率')
report('高齢者のみ世帯割合', '空き家率', '高齢者のみ世帯 × 空き家率')
report('高齢者率',           '地価',     '高齢者率 × 地価')

# 原論文が外れ値扱いした7大都市圏を除外した地価×空き家率
outliers = ['東京都', '神奈川県', '大阪府', '京都府', '愛知県', '埼玉県', '兵庫県']
sub = df[~df['都道府県'].isin(outliers)]
r_o, p_o = stats.pearsonr(sub['地価'], sub['空き家率'])
print(f"  地価 × 空き家率（7大都市圏を除外, n={len(sub)}）: r = {r_o:+.3f} (p = {p_o:.4f})")

# ============================================================
# 図1: 高齢者率 × 空き家率（散布図＋回帰直線）
# ============================================================
print("図1 を生成中 ...")
fig, ax = plt.subplots(figsize=(9, 6.5))
x = df['高齢者率'].values
y = df['空き家率'].values
ax.scatter(x, y, s=60, c='#6aa84f', alpha=0.85,
           edgecolors='white', linewidth=0.6, zorder=3)
sl, ic, r, p, se = stats.linregress(x, y)
xs = np.linspace(x.min(), x.max(), 100)
ax.plot(xs, sl * xs + ic, color='#e69138', lw=2.2, zorder=2,
        label=f'回帰直線 y={sl:.3f}x{ic:+.3f}  (r={r:.3f})')
for _, row in df.iterrows():
    ax.annotate(short(row['都道府県']), (row['高齢者率'], row['空き家率']),
                xytext=(3, 3), textcoords='offset points',
                fontsize=7, color='#555')
ax.set_xlabel('高齢者率（65歳以上人口の割合, %）', fontsize=12)
ax.set_ylabel('空き家率（空き家数/総住宅数, %）', fontsize=12)
ax.set_title('図1  高齢者率 × 空き家率（47都道府県, SSDSE-E）',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=10, loc='upper left')
ax.grid(alpha=0.3, linestyle='--')
ax.set_axisbelow(True)
fig.text(0.5, -0.03,
    '※SSDSE-E（全空き家率）による再現。原論文 図1 は住宅・土地統計調査の「その他の住宅」割合を用いており、\n'
    '  傾き（原論文 y=0.025x-0.3504）や個別数値は原論文の報告値とは一致しない（符号＝向きの比較用）。',
    ha='center', fontsize=8, color='#777')
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig1.png'), bbox_inches='tight')
plt.close(fig)
print(f"  -> fig1 保存 (r={r:.3f}, p={p:.4f})")

# ============================================================
# 図2: 地価 × 空き家率（大都市圏を外れ値として区別／回帰直線2本）
# ============================================================
print("図2 を生成中 ...")
fig, ax = plt.subplots(figsize=(9, 6.5))
is_out = df['都道府県'].isin(outliers)
ax.scatter(df.loc[~is_out, '地価'], df.loc[~is_out, '空き家率'],
           s=60, c='#3d85c6', alpha=0.85, edgecolors='white',
           linewidth=0.6, label='通常の県', zorder=3)
ax.scatter(df.loc[is_out, '地価'], df.loc[is_out, '空き家率'],
           s=95, c='#cc0000', alpha=0.9, marker='D', edgecolors='white',
           linewidth=0.6, label='大都市圏（外れ値）', zorder=4)
for _, row in df[is_out].iterrows():
    ax.annotate(short(row['都道府県']), (row['地価'], row['空き家率']),
                xytext=(4, 2), textcoords='offset points',
                fontsize=8, color='#990000')
sl1, ic1, r1, p1, _ = stats.linregress(df['地価'], df['空き家率'])
sl2, ic2, r2, p2, _ = stats.linregress(sub['地価'], sub['空き家率'])
xs1 = np.linspace(df['地価'].min(), df['地価'].max(), 100)
xs2 = np.linspace(sub['地価'].min(), sub['地価'].max(), 100)
ax.plot(xs1, sl1 * xs1 + ic1, color='#999', lw=1.8, ls='--',
        label=f'全47件 (r={r1:.3f})', zorder=2)
ax.plot(xs2, sl2 * xs2 + ic2, color='#e69138', lw=2.2,
        label=f'外れ値除外 (r={r2:.3f})', zorder=2)
ax.set_xlabel('住宅地の地価（標準価格, 円/m^2）', fontsize=12)
ax.set_ylabel('空き家率（%）', fontsize=12)
ax.set_title('図2  地価 × 空き家率（47都道府県, SSDSE-E+B）',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=9, loc='upper right')
ax.grid(alpha=0.3, linestyle='--')
ax.set_axisbelow(True)
fig.text(0.5, -0.03,
    '※SSDSE（全空き家率）による再現。原論文 図2・3 は「その他の住宅」割合を用いており数値は一致しない。',
    ha='center', fontsize=8, color='#777')
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig2.png'), bbox_inches='tight')
plt.close(fig)
print(f"  -> fig2 保存 (全47 r={r1:.3f}, 除外後 r={r2:.3f})")

# ============================================================
# 図3: 高齢者のみ世帯割合 × 空き家率（散布図＋回帰直線）
# ============================================================
print("図3 を生成中 ...")
fig, ax = plt.subplots(figsize=(9, 6.5))
x = df['高齢者のみ世帯割合'].values
y = df['空き家率'].values
ax.scatter(x, y, s=60, c='#8e7cc3', alpha=0.85,
           edgecolors='white', linewidth=0.6, zorder=3)
sl, ic, r, p, se = stats.linregress(x, y)
xs = np.linspace(x.min(), x.max(), 100)
ax.plot(xs, sl * xs + ic, color='#e69138', lw=2.2, zorder=2,
        label=f'回帰直線 y={sl:.3f}x{ic:+.3f}  (r={r:.3f})')
for _, row in df.iterrows():
    ax.annotate(short(row['都道府県']),
                (row['高齢者のみ世帯割合'], row['空き家率']),
                xytext=(3, 3), textcoords='offset points',
                fontsize=7, color='#555')
ax.set_xlabel('高齢者のみ世帯の割合（%, SSDSE-A集計）', fontsize=12)
ax.set_ylabel('空き家率（%）', fontsize=12)
ax.set_title('図3  高齢者のみ世帯割合 × 空き家率（47都道府県）',
             fontsize=13, fontweight='bold')
ax.legend(fontsize=10, loc='upper left')
ax.grid(alpha=0.3, linestyle='--')
ax.set_axisbelow(True)
fig.text(0.5, -0.03,
    '※SSDSE（全空き家率／高齢者のみ世帯はSSDSE-A集計）による再現。原論文 図7 とは変数定義が異なり数値は一致しない。',
    ha='center', fontsize=8, color='#777')
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig3.png'), bbox_inches='tight')
plt.close(fig)
print(f"  -> fig3 保存 (r={r:.3f}, p={p:.4f})")

# ============================================================
# 図4: 空き家化の理由（国土交通省アンケートの報告値を可視化）
#   ※SSDSE未収録。原論文 図8 の報告値を再掲（再計算ではない）
# ============================================================
print("図4（報告値の可視化）を生成中 ...")
reasons = ['居住者の死亡のため', '居住者の転出のため',
           '家屋の建て替え等のため', 'その他', '入院・施設入居のため']
vals   = [68.1, 23.4, 4.3, 4.3, 0.0]
colors = ['#cc0000', '#e69138', '#6aa84f', '#999999', '#3d85c6']
fig, ax = plt.subplots(figsize=(8.5, 6))
wedges, _, _ = ax.pie(
    vals, colors=colors, startangle=90, counterclock=False,
    autopct=lambda v: f'{v:.1f}%' if v > 0.5 else '',
    pctdistance=0.75, wedgeprops=dict(edgecolor='white', linewidth=1.2))
ax.legend(wedges, [f'{n}（{v:.1f}%）' for n, v in zip(reasons, vals)],
          loc='center left', bbox_to_anchor=(0.98, 0.5), fontsize=10)
ax.set_title('図4  空き家化の理由（国土交通省アンケート）\n'
             '※原論文の報告値の可視化（再計算ではない）',
             fontsize=12, fontweight='bold')
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig4.png'), bbox_inches='tight')
plt.close(fig)
print("  -> fig4 保存（報告値の可視化）")

# ============================================================
# 参考: 空き家処理費用の試算（原論文の報告値, 再掲）
#   ※SSDSE未収録の単価・区分を用いた原論文独自の試算
# ============================================================
print("空き家処理費用（原論文の報告値, 再掲）:")
for name, oku in [('解体', 3270), ('売却(改修)', 134), ('賃貸(管理)', 123)]:
    print(f"  {name:9s}: 約 {oku:>4d} 億円")