日本の電力需要量の影響調査

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 数値変換
num_cols = [
    '光熱・水道費（二人以上の世帯）',
    '消費支出（二人以上の世帯）',
    '年平均気温',
    '最高気温（日最高気温の月平均の最高値）',
    '最低気温（日最低気温の月平均の最低値）',
    '総人口',
    '65歳以上人口',
    '標準価格（平均価格）（住宅地）',
]
for c in num_cols:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# 派生変数
df_b['光熱費割合'] = df_b['光熱・水道費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100
df_b['高齢化率']   = df_b['65歳以上人口'] / df_b['総人口'] * 100

print("=" * 60)
print("■ SSDSE-B 読み込み完了")
print(f"  全レコード数: {len(df_b)}（都道府県×年度）")
print(f"  年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}")
print("=" * 60)

fig1, ax1 = plt.subplots(figsize=(9, 5))

years  = ts['年度'].values
means  = ts['光熱費割合_mean'].values
stds   = ts['光熱費割合_std'].values

ax1.fill_between(years, means - stds, means + stds,
                 alpha=0.18, color='#1565C0', label='±1SD（都道府県間ばらつき）')
ax1.plot(years, means, 'o-', color='#1565C0', linewidth=2.2,
         markersize=7, markerfacecolor='white', markeredgewidth=2, label='全国平均')

# 2020年（コロナ禍）にアノテーション
if 2020 in years:
    idx20 = list(years).index(2020)
    ax1.annotate('2020年\n（コロナ禍）',
                 xy=(2020, means[idx20]),
                 xytext=(2020 - 1.0, means[idx20] + 0.4),
                 arrowprops=dict(arrowstyle='->', color='#C62828'),
                 fontsize=9, color='#C62828')

# 2022年（エネルギー価格高騰）
if 2022 in years:
    idx22 = list(years).index(2022)
    ax1.annotate('2022年\n（エネルギー高騰）',
                 xy=(2022, means[idx22]),
                 xytext=(2021.0, means[idx22] + 0.5),
                 arrowprops=dict(arrowstyle='->', color='#E65100'),
                 fontsize=9, color='#E65100')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('光熱費割合（光熱・水道費 / 消費支出 × 100, %）', fontsize=11)
ax1.set_title('光熱・水道費割合の時系列推移（全国平均）\n2012〜2023年 都道府県データ（SSDSE-B）',
              fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xticks(years)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig1_ts.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2022_U5_5_fig1_ts.png")

corr_plot_cols = [YVAR] + VARNAMES_EN
corr_labels = ['光熱費割合'] + VARNAMES_JP
C = df_ols[corr_plot_cols].corr().values

fig2, ax2 = plt.subplots(figsize=(8, 7))
n = len(corr_labels)
im = ax2.imshow(C, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax2, shrink=0.8, label='ピアソン相関係数')

ax2.set_xticks(range(n))
ax2.set_xticklabels(corr_labels, rotation=30, ha='right', fontsize=10)
ax2.set_yticks(range(n))
ax2.set_yticklabels(corr_labels, fontsize=10)

# セル内に相関係数を表示
for i in range(n):
    for j in range(n):
        val = C[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax2.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=9, color=color, fontweight='bold')

ax2.set_title('相関ヒートマップ：光熱費割合と各説明変数\n（2022年度, N=47都道府県）',
              fontsize=12, fontweight='bold')
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_U5_5_fig2_corr.png")

betas  = np.asarray(model_std.params[1:])
ses    = np.asarray(model_std.bse[1:])
pvals  = np.asarray(model_std.pvalues[1:])

fig3, ax3 = plt.subplots(figsize=(9, 5))
y_pos = np.arange(len(VARNAMES_JP))
colors3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E' for p in pvals]

ax3.barh(y_pos, betas, color=colors3, alpha=0.78, edgecolor='white', height=0.55)
ax3.errorbar(betas, y_pos, xerr=1.96 * ses, fmt='none',
             color='#333', capsize=4, linewidth=1.3)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(VARNAMES_JP, fontsize=11)
ax3.set_xlabel('標準化偏回帰係数（±1.96SE）', fontsize=11)
ax3.set_title(
    f'OLS重回帰：標準化偏回帰係数\n'
    f'R²={model_std.rsquared:.3f}, Adj.R²={model_std.rsquared_adj:.3f}, N={N}都道府県',
    fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
ax3.legend(handles=[
    Patch(color='#C62828', alpha=0.78, label='p < 0.01'),
    Patch(color='#FF8F00', alpha=0.78, label='p < 0.05'),
    Patch(color='#9E9E9E', alpha=0.78, label='n.s.')
], fontsize=9, loc='lower right')

# 係数値をバーに付記
for i, (b, p) in enumerate(zip(betas, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    ha = 'left' if b >= 0 else 'right'
    ax3.text(b + (0.01 if b >= 0 else -0.01), i,
             f' {b:.3f}{sig}', va='center', ha=ha, fontsize=8.5)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_U5_5_fig3_coef.png")

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats
from statsmodels.stats.outliers_influence import variance_inflation_factor
import warnings
warnings.filterwarnings('ignore')

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

import os
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

df_2022 = df_b[df_b['年度'] == 2022].copy()

# 説明変数リスト
XVARS = {
    '年平均気温（℃）'    : '年平均気温',
    '総人口（万人）'      : '総人口',
    '消費支出（円）'      : '消費支出（二人以上の世帯）',
    '高齢化率（%）'       : '高齢化率',
    '住宅地価格（円/m²）' : '標準価格（平均価格）（住宅地）',
}
# 目的変数
YVAR = '光熱費割合'

# 分析用データフレーム作成
cols_use = [YVAR, '年平均気温', '総人口', '消費支出（二人以上の世帯）',
            '高齢化率', '標準価格（平均価格）（住宅地）', '都道府県']
df_ols = df_2022[cols_use].dropna().copy().reset_index(drop=True)
N = len(df_ols)

# 変数名マッピング（表示用）
VARNAMES_JP = list(XVARS.keys())
VARNAMES_EN = list(XVARS.values())

print(f"\n■ 2022年度 有効データ: N={N}都道府県")
print(df_ols[[YVAR] + VARNAMES_EN].describe().round(2))

# ── OLS 標準化係数 ─────────────────────────────────────────────
X_raw = df_ols[VARNAMES_EN].values.astype(float)
y_raw = df_ols[YVAR].values.astype(float)

# 標準化
X_std = (X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0, ddof=1)
y_std = (y_raw - y_raw.mean()) / y_raw.std(ddof=1)

model_std = sm.OLS(y_std, sm.add_constant(X_std)).fit()
model_raw = sm.OLS(y_raw, sm.add_constant(X_raw)).fit()

print("\n" + "=" * 60)
print("■ OLS 重回帰（標準化）")
print("=" * 60)
print(f"  N={N}, R²={model_std.rsquared:.4f}, Adj.R²={model_std.rsquared_adj:.4f}")
print(f"\n  {'変数':<22} {'標準化β':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 50)
for name, b, p in zip(VARNAMES_JP,
                      model_std.params[1:],
                      model_std.pvalues[1:]):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<22} {b:>8.4f} {p:>10.4f} {sig:>6}")

# ── VIF 計算 ────────────────────────────────────────────────
print("\n" + "=" * 60)
print("■ VIF（多重共線性診断）")
print("=" * 60)
X_vif = sm.add_constant(X_raw)
vifs = [variance_inflation_factor(X_vif, i + 1) for i in range(X_raw.shape[1])]
print(f"\n  {'変数':<22} {'VIF':>8}")
print("  " + "-" * 32)
for name, v in zip(VARNAMES_JP, vifs):
    flag = ' ← 注意(>5)' if v > 5 else ''
    print(f"  {name:<22} {v:>8.2f}{flag}")

# ── 相関行列 ────────────────────────────────────────────────
print("\n" + "=" * 60)
print("■ 相関行列（ピアソン）")
print("=" * 60)
corr_cols = [YVAR] + VARNAMES_EN
corr_matrix = df_ols[corr_cols].corr()
print(corr_matrix.round(3))

ts = df_b.groupby('年度').agg(
    光熱費割合_mean=('光熱費割合', 'mean'),
    光熱費割合_std =('光熱費割合', 'std'),
).reset_index()

print("\n■ 光熱費割合 全国平均（年度別）")
print(ts.round(3))

x4 = df_ols['年平均気温'].values.astype(float)
y4 = df_ols[YVAR].values.astype(float)
pref = df_ols['都道府県'].values

fig4, ax4 = plt.subplots(figsize=(9, 6))

# 散布図
scatter = ax4.scatter(x4, y4, c=y4, cmap='RdYlBu_r', s=80,
                      edgecolors='gray', linewidths=0.5, alpha=0.85, zorder=3)
plt.colorbar(scatter, ax=ax4, label='光熱費割合（%）', shrink=0.8)

# 回帰直線
z4 = np.polyfit(x4, y4, 1)
xs4 = np.linspace(x4.min() - 0.5, x4.max() + 0.5, 100)
ax4.plot(xs4, np.poly1d(z4)(xs4), 'r--', linewidth=1.8, alpha=0.8, label='線形回帰')
r4, p4 = stats.pearsonr(x4, y4)

# 主要都道府県にラベル
label_prefs = ['北海道', '沖縄県', '東京都', '大阪府', '秋田県', '鹿児島県', '長野県']
for i, pname in enumerate(pref):
    if any(lp in pname for lp in label_prefs):
        ax4.annotate(pname, (x4[i], y4[i]),
                     textcoords='offset points', xytext=(6, 3),
                     fontsize=8.5, color='#333')

ax4.set_xlabel('年平均気温（℃）', fontsize=12)
ax4.set_ylabel('光熱費割合（光熱・水道費 / 消費支出 × 100, %）', fontsize=11)
ax4.set_title(
    f'年平均気温と光熱費割合の関係（2022年度, N={N}都道府県）\n'
    f'r = {r4:.3f}（p = {p4:.4f}）',
    fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_U5_5_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2022_U5_5_fig4_scatter.png")

print("\n" + "=" * 60)
print("全図の生成完了（4枚）")
print("  fig1_ts.png      : 光熱費割合の時系列（全国平均 2012-2023）")
print("  fig2_corr.png    : 相関ヒートマップ（光熱費と説明変数）")
print("  fig3_coef.png    : OLS標準化係数プロット")
print("  fig4_scatter.png : 年平均気温 vs 光熱費割合 散布図")
print("=" * 60)

import numpy as np
import pandas as pd
from sklearn.linear_model import LassoCV, Lasso
from sklearn.preprocessing import StandardScaler

df_raw = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)
df = df_raw[(df_raw['年度'] == 2022)
            & df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].reset_index(drop=True)

# 目的変数: 光熱・水道費割合（電力需要の家計側プロキシ）
y = (df['光熱・水道費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100).values

# 説明変数: 気候・人口・経済・住宅など 10 変数
X = pd.DataFrame({
    '年平均気温': pd.to_numeric(df['年平均気温'], errors='coerce'),
    '降水量':     pd.to_numeric(df['降水量（年間）'], errors='coerce'),
    '高齢化率':   df['65歳以上人口'] / df['総人口'] * 100,
    '人口密度代理（住宅地価格）': pd.to_numeric(df['標準価格（平均価格）（住宅地）'], errors='coerce') / 10000,
    '持家延べ面積': pd.to_numeric(df['着工新設持家床面積'], errors='coerce') / pd.to_numeric(df['着工新設持家数'], errors='coerce'),
    '消費支出':   df['消費支出（二人以上の世帯）'] / 100000,
    '合計特殊出生率': pd.to_numeric(df['合計特殊出生率'], errors='coerce'),
    '転入率':     df['転入者数（日本人移動者）'] / df['総人口'] * 100,
    '宿泊施設数': pd.to_numeric(df['旅館営業施設数（ホテルを含む）'], errors='coerce') / df['総人口'] * 10000,
    'ごみ排出量': pd.to_numeric(df['1人1日当たりの排出量'], errors='coerce'),
})

# LASSO は変数のスケールに敏感 → 標準化してから使う（原論文と同じ発想）
scaler = StandardScaler()
Xs = scaler.fit_transform(X)

# λ（正則化の強さ）は 5-fold クロスバリデーションで選ぶ
lasso = LassoCV(cv=5, random_state=0).fit(Xs, y)
print(f'CV で選ばれた λ (alpha) = {lasso.alpha_:.4f}')
print(f'決定係数 R2 = {lasso.score(Xs, y):.3f}')
print()
print('=== LASSO の標準化係数（0 になった変数は「選ばれなかった」） ===')
for name, c in sorted(zip(X.columns, lasso.coef_), key=lambda t: -abs(t[1])):
    mark = '← 選択' if abs(c) > 1e-10 else '（脱落）'
    print(f'  {name:16s} {c:+.4f} {mark}')

# λ を強くすると係数が 0 に潰れていく（LASSO の変数選択機能）
print()
for a in [0.1, 0.3]:
    m = Lasso(alpha=a).fit(Xs, y)
    kept = [(n, c) for n, c in zip(X.columns, m.coef_) if abs(c) > 1e-10]
    print(f'λ = {a}: 残った変数 {len(kept)} 個 →',
          ', '.join(f'{n}({c:+.2f})' for n, c in sorted(kept, key=lambda t: -abs(t[1]))))