金融資産購入経験の要因分析
―金融教育・損失回避傾向・Digital Capability Indexに注目して―

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import pandas as pd
from scipy import stats
from numpy.linalg import lstsq

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

import os
FIGDIR = os.path.normpath('html/figures') + os.sep
DATA_B = 'data/raw/SSDSE-B-2026.csv'
DATA_E = 'data/raw/SSDSE-E-2026.csv'
os.makedirs(FIGDIR, exist_ok=True)

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
mask_b = (df_b['地域コード'].str.match(r'^R\d{5}$', na=False) &
          (df_b['地域コード'] != 'R00000') &
          (df_b['年度'] == 2022))
df = df_b[mask_b].copy().reset_index(drop=True)

# SSDSE-E
df_e_raw = pd.read_csv(DATA_E, encoding='cp932', header=0)
df_e = df_e_raw.iloc[2:].copy()
df_e.columns = df_e_raw.iloc[1].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)

df['高等学校卒業者数'] = pd.to_numeric(df['高等学校卒業者数'], errors='coerce')
df['高等学校卒業者のうち進学者数'] = pd.to_numeric(df['高等学校卒業者のうち進学者数'], errors='coerce')
df['進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数']

df['総人口'] = pd.to_numeric(df['総人口'], errors='coerce')
df['65歳以上人口'] = pd.to_numeric(df['65歳以上人口'], errors='coerce')
df['保育所等数'] = pd.to_numeric(df['保育所等数'], errors='coerce')
df['年平均気温'] = pd.to_numeric(df['年平均気温'], errors='coerce')

df['高齢化率'] = df['65歳以上人口'] / df['総人口']
df['保育所数千対'] = df['保育所等数'] / df['総人口'] * 10000   # 人口万人対

# SSDSE-E からの変数
def normalize_pref(s):
    return str(s).rstrip('県府都道').strip()

df['pref_short'] = df['都道府県'].apply(normalize_pref)
df_e['pref_short'] = df_e['都道府県'].apply(normalize_pref)

df_e_use = df_e.set_index('pref_short')
for col_e, col_new in [
    ('1人当たり県民所得（平成27年基準）', '1人当たり所得'),
    ('総面積（北方地域及び竹島を除く）', '総面積_ha'),
]:
    df[col_new] = df['pref_short'].map(
        pd.to_numeric(df_e_use[col_e], errors='coerce'))

df['人口密度'] = df['総人口'] / (df['総面積_ha'] / 100)   # 人/km²

# 目的変数: 進学率が全国中央値超 → 1
median_rate = df['進学率'].median()
df['進学率高（binary）'] = (df['進学率'] > median_rate).astype(int)

feature_names = ['1人当たり所得', '高齢化率', '人口密度', '保育所数千対', '年平均気温']
df_model = df[['都道府県', '進学率', '進学率高（binary）'] + feature_names].dropna()
df_model = df_model.reset_index(drop=True)

print(f"サンプル数: {len(df_model)}")
print(f"進学率 (全国): 中央値={median_rate:.3f}, "
      f"binary=1 が {df_model['進学率高（binary）'].sum()} 都道府県")

X_raw = df_model[feature_names].values.astype(float)
y = df_model['進学率高（binary）'].values.astype(float)
n = len(y)

def sigmoid(z):
    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

def logistic_fit(X, y_v, lr=0.1, n_iter=3000):
    n_s, k = X.shape
    b = np.zeros(k)
    for _ in range(n_iter):
        p = sigmoid(X @ b)
        grad = X.T @ (p - y_v) / n_s
        b -= lr * grad
    return b

X_m = X_raw.mean(axis=0)
X_s = X_raw.std(axis=0)
X_s[X_s == 0] = 1
X_std = (X_raw - X_m) / X_s
X_c = np.column_stack([np.ones(n), X_std])

coef = logistic_fit(X_c, y)

# 標準誤差（Fisher情報行列逆）
p_hat = sigmoid(X_c @ coef)
W = p_hat * (1 - p_hat)
W = np.maximum(W, 1e-6)
H = X_c.T @ (W[:, None] * X_c)
try:
    cov_mat = np.linalg.pinv(H)
except Exception:
    cov_mat = np.eye(len(coef)) * 0.01
se_coef = np.sqrt(np.abs(np.diag(cov_mat)))

z_stats = coef / (se_coef + 1e-12)
p_vals = 2 * (1 - stats.norm.cdf(np.abs(z_stats)))
or_vals = np.exp(coef)
ci_low = np.exp(coef - 1.96 * se_coef)
ci_high = np.exp(coef + 1.96 * se_coef)

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

ax = axes[0]
counts = df_model['進学率高（binary）'].value_counts().sort_index()
labels_pie = ['進学率 低（0）', '進学率 高（1）']
colors_pie = ['#90A4AE', '#2E7D32']
wedges, texts, autotexts = ax.pie(
    counts.values, labels=labels_pie, colors=colors_pie,
    autopct='%1.1f%%', startangle=90,
    textprops={'fontsize': 11})
ax.set_title("大学進学率（中央値超=高）の割合", fontsize=12)

ax = axes[1]
groups_sorted = df_model.sort_values('進学率')
bars = ax.bar(range(len(groups_sorted)), groups_sorted['進学率'] * 100,
              color=['#C62828' if v == 1 else '#1565C0'
                     for v in groups_sorted['進学率高（binary）']],
              alpha=0.8)
ax.axhline(median_rate * 100, color='black', lw=2, linestyle='--',
           label=f'中央値 {median_rate*100:.1f}%')
ax.set_xlabel("都道府県（進学率昇順）")
ax.set_ylabel("大学進学率（%）")
ax.set_title("都道府県別大学進学率と中央値閾値", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, axis='y', alpha=0.3)

red_p = mpatches.Patch(color='#C62828', alpha=0.8, label='高（中央値超）')
blue_p = mpatches.Patch(color='#1565C0', alpha=0.8, label='低（中央値以下）')
axes[1].legend(handles=[red_p, blue_p, plt.Line2D([0], [0], color='black',
               lw=2, linestyle='--', label=f'中央値={median_rate*100:.1f}%')],
               fontsize=9)

fig.suptitle("図1: 目的変数（大学進学率 高/低）の分布", fontsize=13)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig1_dist.png", dpi=150)
plt.close()
print("fig1 saved")

def calc_vif(X):
    n_s, k_s = X.shape
    vifs = []
    for j in range(k_s):
        Xj = X[:, j]
        Xi = np.delete(X, j, axis=1)
        Xi_c = np.column_stack([np.ones(n_s), Xi])
        b_j, _, _, _ = lstsq(Xi_c, Xj, rcond=None)
        pred_j = Xi_c @ b_j
        r2_j = 1 - np.var(Xj - pred_j) / max(np.var(Xj), 1e-10)
        vifs.append(1 / max(1 - r2_j, 1e-10))
    return vifs

vif_vals = calc_vif(X_raw)

fig, ax = plt.subplots(figsize=(8, 5))
colors_vif = ['#C62828' if v > 10 else ('#F9A825' if v > 5 else '#2E7D32') for v in vif_vals]
bars = ax.barh(feature_names, vif_vals, color=colors_vif, alpha=0.85)
ax.axvline(5, color='#F9A825', lw=1.5, linestyle='--', label='VIF=5（警戒水準）')
ax.axvline(10, color='#C62828', lw=1.5, linestyle='--', label='VIF=10（問題水準）')

for bar, v in zip(bars, vif_vals):
    ax.text(v + 0.05, bar.get_y() + bar.get_height() / 2,
            f'{v:.2f}', va='center', fontsize=10)

ax.set_xlabel("VIF（分散膨張係数）")
ax.set_title("図2: 説明変数のVIF（多重共線性チェック）", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig2_vif.png", dpi=150)
plt.close()
print("fig2 saved")

fig, ax = plt.subplots(figsize=(9, 6))

or_display = or_vals[1:]
ci_low_d = ci_low[1:]
ci_high_d = ci_high[1:]
p_display = p_vals[1:]

y_pos = np.arange(len(feature_names))
colors_or = ['#C62828' if or_v > 1 else '#1565C0' for or_v in or_display]

ax.axvline(1.0, color='black', lw=1.5, linestyle='-')

x_max = min(max(ci_high_d) * 1.1, 20)
for i, (or_v, cl, ch, p_v) in enumerate(zip(or_display, ci_low_d, ci_high_d, p_display)):
    cl_plot = max(cl, 0.05)
    ch_plot = min(ch, x_max - 0.5)
    ax.plot([cl_plot, ch_plot], [i, i], color=colors_or[i], lw=2.5, alpha=0.8)
    ax.scatter([or_v], [i], color=colors_or[i], s=80, zorder=5)
    sig = "***" if p_v < 0.001 else ("**" if p_v < 0.01 else ("*" if p_v < 0.05 else "n.s."))
    ax.text(min(ch_plot + 0.1, x_max - 0.2), i,
            f'OR={or_v:.2f} {sig}', va='center', fontsize=10)

ax.set_yticks(y_pos)
ax.set_yticklabels(feature_names, fontsize=11)
ax.set_xlabel("オッズ比（95% CI）")
ax.set_title("図3: ロジスティック回帰 オッズ比 Forest Plot\n（目的変数: 大学進学率 高/低）", fontsize=13)
ax.set_xlim(0.0, x_max)

red_patch = mpatches.Patch(color='#C62828', alpha=0.8, label='OR>1（進学率高を増加）')
blue_patch = mpatches.Patch(color='#1565C0', alpha=0.8, label='OR<1（進学率高を減少）')
ax.legend(handles=[red_patch, blue_patch], fontsize=10)
ax.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig3_odds.png", dpi=150)
plt.close()
print("fig3 saved")

fig, axes = plt.subplots(1, 3, figsize=(13, 5))

# 各変数の平均値（他変数固定用）
x_mean_std = np.zeros(len(feature_names))   # standardized mean = 0

def predict_prob(var_idx, var_range_raw, coef, X_m, X_s):
    """1変数を変化させたときの予測確率"""
    var_range_std = (var_range_raw - X_m[var_idx]) / X_s[var_idx]
    probs = []
    for v_std in var_range_std:
        x_pts = x_mean_std.copy()
        x_pts[var_idx] = v_std
        z = coef[0] + coef[1:] @ x_pts
        probs.append(sigmoid(z) * 100)
    return np.array(probs)

var_configs = [
    (0, '1人当たり所得', '1人当たり県民所得（万円）', None),
    (1, '高齢化率', '高齢化率', None),
    (3, '保育所数千対', '保育所数（人口万人あたり）', None),
]

for ax, (vi, vname, xlabel, _) in zip(axes, var_configs):
    col_data = X_raw[:, vi]
    var_range = np.linspace(col_data.min(), col_data.max(), 200)
    prob_line = predict_prob(vi, var_range, coef, X_m, X_s)

    ax.plot(var_range, prob_line, color='#1565C0', lw=2.5)
    ax.scatter(col_data, sigmoid(X_c @ coef) * 100, color='gray', alpha=0.5, s=20)
    ax.set_xlabel(xlabel)
    ax.set_ylabel("大学進学率高の予測確率（%）")
    ax.set_title(f"{vname}\nと進学率予測確率")
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 100)

fig.suptitle("図4: 主要変数別 大学進学率高の予測確率", fontsize=13)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig4_prob.png", dpi=150)
plt.close()
print("fig4 saved")
print("All figures saved.")