住宅価格の格差と影響を与える要因
―空間パネルモデルに基づく研究―

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()

print("=" * 60)
print(f"■ 読み込み完了: {len(df_b)}行 × {df_b.shape[1]}列")
print(f"  年度: {sorted(df_b['年度'].unique())}")
print(f"  都道府県数: {df_b['都道府県'].nunique()}")
print("=" * 60)

# ── 使用列 ─────────────────────────────────────────────────────
TARGET = '標準価格（平均価格）（住宅地）'

FEATURES = [
    '総人口',
    '月間有効求人数（一般）',
    '月間有効求職者数（一般）',
    '消費支出（二人以上の世帯）',
    '65歳以上人口',
    '大学学生数',
    '着工新設住宅戸数',
    '標準価格（平均価格）（商業地）',
]

FEATURE_LABELS = [
    '総人口',
    '有効求人数',
    '有効求職者数',
    '消費支出',
    '65歳以上人口',
    '大学学生数',
    '着工新設住宅戸数',
    '商業地価格',
]

# 数値変換
for c in [TARGET] + FEATURES:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# 欠損除外
df_clean = df_b.dropna(subset=[TARGET] + FEATURES).copy()
print(f"欠損除外後: {len(df_clean)}行")

# 2022年断面データ
df_2022 = df_b[df_b['年度'] == 2022].copy()
df_2022 = df_2022.dropna(subset=[TARGET] + FEATURES).reset_index(drop=True)
print(f"2022年データ: N={len(df_2022)}")

# ── 地域区分（大都市圏 / 地方） ─────────────────────────────────
URBAN = ['東京都', '神奈川県', '埼玉県', '千葉県', '大阪府', '京都府', '愛知県', '兵庫県']

def region_label(pref):
    if pref == '東京都':
        return '東京都'
    elif pref in URBAN:
        return '主要都市圏'
    else:
        return '地方'

df_2022['地域区分'] = df_2022['都道府県'].apply(region_label)

# 地域区分カラーマップ
REGION_COLOR = {'東京都': '#C62828', '主要都市圏': '#1565C0', '地方': '#43A047'}

# ── 時系列用: 代表都道府県の選択 ────────────────────────────────
# 都市圏: 東京都, 大阪府, 愛知県  地方: 北海道, 宮城県, 福岡県, 鹿児島県
PREF_TS = {
    '東京都':   ('#C62828', '-',  'D'),
    '大阪府':   ('#E65100', '--', 's'),
    '愛知県':   ('#FF8F00', ':',  '^'),
    '北海道':   ('#1565C0', '-',  'o'),
    '宮城県':   ('#0277BD', '--', 'v'),
    '福岡県':   ('#1B5E20', ':',  'P'),
    '鹿児島県': ('#6A1B9A', '-.', 'X'),
}

fig1, ax1 = plt.subplots(figsize=(10, 6))

years = sorted(df_clean['年度'].unique())
for pref, (color, ls, marker) in PREF_TS.items():
    ts = df_clean[df_clean['都道府県'] == pref].sort_values('年度')
    ax1.plot(ts['年度'], ts[TARGET] / 1000,
             color=color, linestyle=ls, marker=marker,
             markersize=5, linewidth=1.8, label=pref)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('住宅地 標準価格（千円/㎡）', fontsize=12)
ax1.set_title('住宅地価格の時系列推移（2012〜2023年）\n都市圏と地方の格差',
              fontsize=13, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9, ncol=2, framealpha=0.9)
ax1.set_xticks(years)
ax1.set_xticklabels([str(y) for y in years], rotation=45, ha='right', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.yaxis.set_major_formatter(
    matplotlib.ticker.FuncFormatter(lambda x, _: f'{x:,.0f}')
)

# 都市 / 地方 の区分ラベル
ax1.annotate('都市圏', xy=(2023, df_clean[(df_clean['都道府県']=='東京都') &
              (df_clean['年度']==2023)][TARGET].values[0] / 1000),
             xytext=(2021.5, 360), fontsize=10, color='#C62828',
             arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.2))
ax1.annotate('地方', xy=(2023, df_clean[(df_clean['都道府県']=='鹿児島県') &
              (df_clean['年度']==2023)][TARGET].values[0] / 1000),
             xytext=(2021.0, 60), fontsize=10, color='#6A1B9A',
             arrowprops=dict(arrowstyle='->', color='#6A1B9A', lw=1.2))

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_U5_16_fig1_ts.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2022_U5_16_fig1_ts.png")

df_rank = df_2022[['都道府県', TARGET, '地域区分']].sort_values(TARGET).reset_index(drop=True)

fig2, ax2 = plt.subplots(figsize=(9, 13))

colors2 = [REGION_COLOR[r] for r in df_rank['地域区分']]
bars = ax2.barh(
    range(len(df_rank)),
    df_rank[TARGET] / 1000,
    color=colors2, alpha=0.82, edgecolor='white', height=0.75
)
ax2.set_yticks(range(len(df_rank)))
ax2.set_yticklabels(df_rank['都道府県'], fontsize=9)
ax2.set_xlabel('住宅地 標準価格（千円/㎡）', fontsize=11)
ax2.set_title('47都道府県 住宅地価格ランキング（2022年）\nSSDSE-B 実データ',
              fontsize=12, fontweight='bold')

# 東京都に値ラベル
top_val = df_rank[TARGET].max() / 1000
ax2.text(top_val + 2, len(df_rank) - 1,
         f'{top_val:,.0f}', va='center', fontsize=8, color='#C62828', fontweight='bold')

ax2.axvline(df_2022[TARGET].mean() / 1000, color='gray',
            linestyle='--', linewidth=1.2, label=f'全国平均={df_2022[TARGET].mean()/1000:.0f}千円/㎡')
ax2.legend(fontsize=9, loc='lower right')
ax2.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
legend_handles = [
    Patch(color=REGION_COLOR['東京都'],     alpha=0.82, label='東京都'),
    Patch(color=REGION_COLOR['主要都市圏'], alpha=0.82, label='主要都市圏'),
    Patch(color=REGION_COLOR['地方'],       alpha=0.82, label='地方'),
]
ax2.legend(handles=legend_handles + [
    plt.Line2D([0], [0], color='gray', linestyle='--',
               label=f'全国平均={df_2022[TARGET].mean()/1000:.0f}千円/㎡')
], fontsize=9, loc='lower right')

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

n_feat = len(FEATURE_LABELS)
x_pos = np.arange(n_feat)
width = 0.25

fig3, ax3 = plt.subplots(figsize=(12, 6))

bars_ols   = ax3.bar(x_pos - width, ols_coefs   / 1000, width, label='OLS',
                      color='#1565C0', alpha=0.80, edgecolor='white')
bars_ridge = ax3.bar(x_pos,          ridge_coefs / 1000, width, label='Ridge (α=1.0)',
                      color='#2E7D32', alpha=0.80, edgecolor='white')
bars_lasso = ax3.bar(x_pos + width,  lasso_coefs / 1000, width, label='Lasso (α=0.01)',
                      color='#E65100', alpha=0.80, edgecolor='white')

ax3.axhline(0, color='gray', linewidth=0.8, linestyle='--')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(FEATURE_LABELS, rotation=30, ha='right', fontsize=10)
ax3.set_ylabel('標準化係数（×10³ 円/㎡ 単位）', fontsize=11)
ax3.set_title('OLS・Ridge・Lasso の係数比較（2022年断面, 標準化後）\n'
              '正則化による係数の収縮効果',
              fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(axis='y', alpha=0.3)

# Lasso でゼロになった変数にマーク
for i, (lc, name) in enumerate(zip(lasso_coefs, FEATURE_LABELS)):
    if abs(lc) < 1e-6:
        ax3.text(i + width, 0.5, 'Lasso=0', ha='center', va='bottom',
                 fontsize=7.5, color='#E65100', rotation=0)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_U5_16_fig3_model.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_U5_16_fig3_model.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 sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from scipy import stats
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)

y_2022 = df_2022[TARGET].values
X_2022 = df_2022[FEATURES].values

# OLS（標準化あり → 標準化係数）
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_2022)

ols_model = sm.OLS(y_2022, sm.add_constant(X_scaled)).fit()
ols_coefs = np.array(ols_model.params[1:])

# Ridge (alpha=1.0, random_state=42)
ridge = Ridge(alpha=1.0, random_state=42)
ridge.fit(X_scaled, y_2022)
ridge_coefs = ridge.coef_

# Lasso (alpha=0.01, random_state=42)
lasso = Lasso(alpha=0.01, random_state=42, max_iter=10000)
lasso.fit(X_scaled, y_2022)
lasso_coefs = lasso.coef_

print("\n" + "=" * 60)
print("■ OLS / Ridge / Lasso 係数比較（標準化後）")
print("=" * 60)
print(f"{'変数':<20} {'OLS':>12} {'Ridge':>12} {'Lasso':>12}")
print("-" * 60)
for name, o, r, la in zip(FEATURE_LABELS, ols_coefs, ridge_coefs, lasso_coefs):
    print(f"{name:<20} {o:>12.1f} {r:>12.1f} {la:>12.1f}")
print(f"\nOLS R² = {ols_model.rsquared:.4f}")
print(ols_model.summary2())

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

for region, color in REGION_COLOR.items():
    sub = df_2022[df_2022['地域区分'] == region]
    ax4.scatter(sub['総人口'] / 1e6,
                sub[TARGET] / 1000,
                c=color, s=60, alpha=0.85, label=region, zorder=3)
    # 都道府県名ラベル（東京・大阪・秋田など主要県）
    for _, row in sub.iterrows():
        if row[TARGET] > 100000 or row['都道府県'] in ['秋田県', '青森県', '鳥取県',
                                                         '北海道', '大阪府']:
            ax4.annotate(row['都道府県'],
                         xy=(row['総人口']/1e6, row[TARGET]/1000),
                         xytext=(3, 3), textcoords='offset points',
                         fontsize=7.5, color=color)

# 回帰直線（全体）
x4 = df_2022['総人口'].values / 1e6
y4 = df_2022[TARGET].values / 1000
z4 = np.polyfit(x4, y4, 1)
xs4 = np.linspace(x4.min(), x4.max(), 200)
ax4.plot(xs4, np.poly1d(z4)(xs4), 'k--', linewidth=1.2, alpha=0.6, label='回帰直線（全体）')

r4, p4 = stats.pearsonr(x4, y4)
ax4.set_xlabel('総人口（百万人）', fontsize=12)
ax4.set_ylabel('住宅地 標準価格（千円/㎡）', fontsize=12)
ax4.set_title(f'総人口と住宅地価格の関係（2022年）\n'
              f'r = {r4:.3f}  (p = {p4:.4f})  N=47都道府県',
              fontsize=12, fontweight='bold')
ax4.legend(fontsize=10, loc='upper left')
ax4.grid(True, alpha=0.3)

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

print("\n" + "=" * 60)
print("■ 全図の生成完了（4枚）")
print("  2022_U5_16_fig1_ts.png     : 住宅地価格 時系列")
print("  2022_U5_16_fig2_rank.png   : 47都道府県 価格ランキング")
print("  2022_U5_16_fig3_model.png  : OLS vs Ridge vs Lasso 係数比較")
print("  2022_U5_16_fig4_scatter.png: 総人口 vs 住宅地価格 散布図")
print("=" * 60)