"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
==========================================================================
論文タイトル：パネルデータを用いた空き家に関する要因分析
手法：パネルデータ分析（固定効果モデル）、Hausman検定

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県×年度パネルデータ）
          47都道府県 × 12年（2012〜2023年）= 564観測
  目的変数（空き家の代理指標）：
    ・転出率（転出者数 / 総人口）  → 人口流出が空き家増加に寄与
    ・着工新設住宅戸数（人口1万対）→ 新規建設が多い＝空き家リスク上昇の代理
  説明変数：
    ・高齢化率（65歳以上人口 / 総人口）
    ・住宅地地価（円 / 平方メートル）
    ・有効求人倍率（月間有効求人数 / 月間有効求職者数）
    ・消費支出（二人以上世帯、円）

  Step1. 時系列プロット（全国平均）
  Step2. パネルOLS 固定効果モデル（PanelOLS）
  Step3. Hausman 検定（FE vs RE）
  Step4. 散布図（転出率 vs 住宅地地価、都道府県ラベル）
  Step5. 地域別転出率比較（6地域区分）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ（2012〜2023年）

【データサイエンス学習ポイント】
  1. パネルデータの構造（個体×時点）と固定効果モデルの意義
  2. Hausman 検定（FE vs RE）の手順と解釈
  3. クラスタリング標準誤差（都道府県単位）
  4. 地域区分による記述統計の可視化
==========================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2022_U5_1_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import warnings
warnings.filterwarnings('ignore')

from linearmodels.panel import PanelOLS, RandomEffects
import statsmodels.api as sm
from scipy import stats

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
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)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026: 都道府県パネル）
# ================================================================
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()
# 年度を整数に
df_b['年度'] = df_b['年度'].astype(int)
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

print(f"パネルデータ形状: {df_b.shape}  ({df_b['都道府県'].nunique()} 都道府県 × {df_b['年度'].nunique()} 年度)")
print(f"年度範囲: {df_b['年度'].min()} 〜 {df_b['年度'].max()}")

# ================================================================
# ■ 変数生成
# ================================================================
df_b['転出率']      = df_b['転出者数（日本人移動者）'] / df_b['総人口'] * 1000   # 人口千対
df_b['転入率']      = df_b['転入者数（日本人移動者）'] / df_b['総人口'] * 1000
df_b['着工率']      = df_b['着工新設住宅戸数'] / df_b['総人口'] * 10000          # 人口1万対
df_b['高齢化率']    = df_b['65歳以上人口'] / df_b['総人口'] * 100                # %
df_b['地価_log']    = np.log(df_b['標準価格（平均価格）（住宅地）'].clip(lower=1))
df_b['求人倍率']    = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'].replace(0, np.nan)
df_b['消費支出_log']= np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['純移動率']    = (df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']) / df_b['総人口'] * 1000

# ================================================================
# ■ 地域区分（6地域）
# ================================================================
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部',
    '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
df_b['地域'] = df_b['都道府県'].map(region_map)
region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']

# ================================================================
# ■ 図1: 転出率と着工率の時系列（全国平均）
# ================================================================
print("\n[図1] 時系列プロット作成中...")
ts = df_b.groupby('年度').agg(
    転出率_mean=('転出率', 'mean'),
    着工率_mean=('着工率', 'mean'),
    転入率_mean=('転入率', 'mean'),
).reset_index()

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
fig.suptitle('転出率・着工新設住宅の全国平均推移（2012〜2023年）', fontsize=14, fontweight='bold', y=1.01)

# 左：転出率
ax = axes[0]
ax.plot(ts['年度'], ts['転出率_mean'], 'o-', color='#C62828', linewidth=2.2,
        markersize=6, label='転出率（千人対）')
ax.plot(ts['年度'], ts['転入率_mean'], 's--', color='#1565C0', linewidth=2,
        markersize=5, label='転入率（千人対）', alpha=0.8)
ax.fill_between(ts['年度'], ts['転出率_mean'], ts['転入率_mean'],
                where=ts['転出率_mean'] > ts['転入率_mean'],
                alpha=0.15, color='#C62828', label='転出超過')
ax.fill_between(ts['年度'], ts['転出率_mean'], ts['転入率_mean'],
                where=ts['転出率_mean'] <= ts['転入率_mean'],
                alpha=0.15, color='#1565C0')
ax.set_xlabel('年度', fontsize=11)
ax.set_ylabel('率（人口千対）', fontsize=11)
ax.set_title('転出率・転入率', fontsize=12)
ax.legend(fontsize=10)
ax.grid(axis='y', alpha=0.4)
ax.xaxis.set_major_locator(mticker.MultipleLocator(2))

# 右：着工率
ax2 = axes[1]
ax2.plot(ts['年度'], ts['着工率_mean'], 'D-', color='#2E7D32', linewidth=2.2,
         markersize=6, label='着工新設住宅戸数（人口1万対）')
ax2.set_xlabel('年度', fontsize=11)
ax2.set_ylabel('戸数（人口1万対）', fontsize=11)
ax2.set_title('着工新設住宅戸数', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(axis='y', alpha=0.4)
ax2.xaxis.set_major_locator(mticker.MultipleLocator(2))

# 補足テキスト
axes[0].annotate('コロナ禍\n2020年', xy=(2020, ts.loc[ts['年度']==2020,'転出率_mean'].values[0]),
                 xytext=(2021, ts.loc[ts['年度']==2020,'転出率_mean'].values[0]+0.5),
                 fontsize=9, color='gray',
                 arrowprops=dict(arrowstyle='->', color='gray', lw=1))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig1_ts.png")

# ================================================================
# ■ パネルデータ準備（PanelOLS 用）
# ================================================================
print("\n[パネル分析] 固定効果モデル推定中...")

EXOG_VARS = ['高齢化率', '地価_log', '求人倍率', '消費支出_log']
ENDOG_VAR  = '転出率'

df_panel = df_b[['都道府県', '年度', ENDOG_VAR] + EXOG_VARS].dropna().copy()
print(f"  分析対象: {len(df_panel)} 観測 ({df_panel['都道府県'].nunique()} 都道府県)")

# MultiIndex 設定
df_panel = df_panel.set_index(['都道府県', '年度'])

# ── 固定効果モデル（FE）
fe_model = PanelOLS(
    df_panel[ENDOG_VAR],
    sm.add_constant(df_panel[EXOG_VARS]),
    entity_effects=True,
    time_effects=False,
)
fe_res = fe_model.fit(cov_type='clustered', cluster_entity=True)
print("\n[固定効果モデル推定結果]")
print(fe_res.summary.tables[1])

# ── ランダム効果モデル（RE）— Hausman 検定用
re_model = RandomEffects(
    df_panel[ENDOG_VAR],
    sm.add_constant(df_panel[EXOG_VARS]),
)
re_res = re_model.fit(cov_type='robust')

# ================================================================
# ■ Hausman 検定（FE vs RE）
# ================================================================
print("\n[Hausman検定] FE vs RE 係数差の検定...")
fe_params = fe_res.params.drop('const', errors='ignore')
re_params = re_res.params.drop('const', errors='ignore')
common_vars = fe_params.index.intersection(re_params.index)

b_diff = fe_params[common_vars].values - re_params[common_vars].values

# 共分散行列差
fe_cov = fe_res.cov.loc[common_vars, common_vars].values
re_cov = re_res.cov.loc[common_vars, common_vars].values
cov_diff = fe_cov - re_cov

try:
    cov_inv = np.linalg.pinv(cov_diff)
    H_stat = float(b_diff @ cov_inv @ b_diff)
    H_df   = len(common_vars)
    H_pval = 1 - stats.chi2.cdf(H_stat, H_df)
    hausman_ok = True
except Exception:
    H_stat, H_df, H_pval = np.nan, len(common_vars), np.nan
    hausman_ok = False

print(f"  Hausman 統計量 H = {H_stat:.4f}")
print(f"  自由度 = {H_df}")
print(f"  p値 = {H_pval:.4f}")
if hausman_ok:
    if H_pval < 0.05:
        print("  → p<0.05: 固定効果モデルが適切（個体固有の偏りを制御）")
    else:
        print("  → p≥0.05: ランダム効果モデルも許容")

# ================================================================
# ■ 図2: 固定効果係数プロット
# ================================================================
print("\n[図2] 固定効果係数プロット作成中...")

coef_df = pd.DataFrame({
    'coef': fe_res.params.drop('const', errors='ignore'),
    'se':   fe_res.std_errors.drop('const', errors='ignore'),
    'pval': fe_res.pvalues.drop('const', errors='ignore'),
})
coef_df['ci_lo'] = coef_df['coef'] - 1.96 * coef_df['se']
coef_df['ci_hi'] = coef_df['coef'] + 1.96 * coef_df['se']
coef_df = coef_df.reindex(EXOG_VARS)

var_labels = {
    '高齢化率':    '高齢化率（%）',
    '地価_log':    '住宅地地価\n（log, 円/㎡）',
    '求人倍率':    '有効求人倍率',
    '消費支出_log':'消費支出\n（log, 円）',
}
labels = [var_labels.get(v, v) for v in coef_df.index]
colors = ['#C62828' if p < 0.05 else '#90A4AE' for p in coef_df['pval']]
sig_markers = ['★' if p < 0.01 else '✓' if p < 0.05 else '' for p in coef_df['pval']]

fig, ax = plt.subplots(figsize=(8, 5))
y_pos = np.arange(len(coef_df))
ax.barh(y_pos, coef_df['coef'], color=colors, height=0.55, alpha=0.85)
ax.errorbar(coef_df['coef'], y_pos,
            xerr=1.96 * coef_df['se'],
            fmt='none', color='#333', elinewidth=1.5, capsize=5)
ax.axvline(0, color='black', linewidth=0.9, linestyle='-')

for i, (coef, sig, pval) in enumerate(zip(coef_df['coef'], sig_markers, coef_df['pval'])):
    xoff = coef_df['ci_hi'].iloc[i] + abs(coef_df['ci_hi'].max()) * 0.03
    pstr = f'p={pval:.3f}' if pval >= 0.001 else 'p<0.001'
    ax.text(xoff, i, f'{sig} {pstr}', va='center', fontsize=9,
            color='#C62828' if pval < 0.05 else '#666')

ax.set_yticks(y_pos)
ax.set_yticklabels(labels, fontsize=11)
ax.set_xlabel('推定係数（従属変数：転出率 人口千対）', fontsize=11)
ax.set_title('固定効果モデル（FE）の推定係数\n（エラーバー：95%信頼区間、標準誤差はクラスタリング）',
             fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.4, linestyle='--')

# 凡例
from matplotlib.patches import Patch
legend_els = [
    Patch(facecolor='#C62828', label='p < 0.05（統計的有意）'),
    Patch(facecolor='#90A4AE', label='p ≥ 0.05（非有意）'),
]
ax.legend(handles=legend_els, fontsize=10, loc='lower right')

# Hausman 結果をテキスト注釈
h_text = (f'Hausman検定: H={H_stat:.2f}, p={H_pval:.3f}\n'
          f'→ 固定効果モデル{"適切" if H_pval < 0.05 else "（RE許容）"}')
ax.text(0.02, 0.98, h_text, transform=ax.transAxes,
        fontsize=9, va='top', color='#333',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#EFF3FF', alpha=0.8))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig2_fe_coef.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig2_fe_coef.png")

# ================================================================
# ■ 図3: 転出率 vs 住宅地地価 散布図（都道府県ラベル）
# ================================================================
print("\n[図3] 散布図（転出率 vs 住宅地地価）作成中...")

# 2012〜2023 年の平均で集計
df_avg = df_b.groupby('都道府県').agg(
    転出率_avg=('転出率', 'mean'),
    地価_avg=('標準価格（平均価格）（住宅地）', 'mean'),
    地域=('地域', 'first'),
).reset_index()

region_colors = {
    '北海道・東北': '#1565C0',
    '関東':        '#C62828',
    '中部':        '#2E7D32',
    '近畿':        '#E65100',
    '中国・四国':  '#6A1B9A',
    '九州・沖縄':  '#00838F',
}

fig, ax = plt.subplots(figsize=(10, 7))
for region, grp in df_avg.groupby('地域'):
    ax.scatter(grp['地価_avg'], grp['転出率_avg'],
               color=region_colors.get(region, 'gray'),
               s=70, alpha=0.85, label=region, zorder=3)

# 都道府県ラベル（一部省略せず全て）
for _, row in df_avg.iterrows():
    pref_short = row['都道府県'].replace('県','').replace('府','').replace('都','').replace('道','')
    ax.annotate(pref_short, (row['地価_avg'], row['転出率_avg']),
                textcoords='offset points', xytext=(4, 2),
                fontsize=7, alpha=0.85, color='#333')

# OLS 回帰線
x_all = df_avg['地価_avg'].values
y_all = df_avg['転出率_avg'].values
slope, intercept, r_val, p_val, _ = stats.linregress(x_all, y_all)
xline = np.linspace(x_all.min(), x_all.max(), 100)
ax.plot(xline, slope * xline + intercept, '--', color='gray', linewidth=1.5,
        label=f'回帰線 (r={r_val:.2f}, p={p_val:.3f})')

ax.set_xlabel('住宅地平均地価（円/㎡）\n2012〜2023年平均', fontsize=12)
ax.set_ylabel('転出率（人口千対）\n2012〜2023年平均', fontsize=12)
ax.set_title('転出率 vs 住宅地地価（都道府県別）', fontsize=13, fontweight='bold')
ax.legend(fontsize=10, ncol=2, loc='upper right')
ax.grid(alpha=0.3, linestyle='--')
ax.xaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'{int(x):,}'))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig3_scatter.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig3_scatter.png")

# ================================================================
# ■ 図4: 地域別転出率比較（箱ひげ図）
# ================================================================
print("\n[図4] 地域別転出率比較作成中...")

# 2022年（最新の完全年）で比較（全期間の平均も補助で）
df_region_data = []
for region in region_order:
    vals = df_b[df_b['地域'] == region]['転出率'].dropna().values
    df_region_data.append(vals)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
fig.suptitle('地域別 転出率の分布（2012〜2023年）', fontsize=14, fontweight='bold')

# 左：箱ひげ図
ax = axes[0]
bp = ax.boxplot(df_region_data, patch_artist=True, notch=False,
                medianprops=dict(color='white', linewidth=2.5),
                whiskerprops=dict(linewidth=1.5),
                flierprops=dict(marker='o', markersize=4, alpha=0.5))
for patch, region in zip(bp['boxes'], region_order):
    patch.set_facecolor(region_colors[region])
    patch.set_alpha(0.8)
ax.set_xticks(range(1, len(region_order)+1))
ax.set_xticklabels([r.replace('・', '\n') for r in region_order], fontsize=9)
ax.set_ylabel('転出率（人口千対）', fontsize=11)
ax.set_title('地域別転出率（箱ひげ図）', fontsize=12)
ax.grid(axis='y', alpha=0.4)

# 右：年別推移（地域別折れ線）
ax2 = axes[1]
for region in region_order:
    ts_r = df_b[df_b['地域'] == region].groupby('年度')['転出率'].mean()
    ax2.plot(ts_r.index, ts_r.values, 'o-', color=region_colors[region],
             linewidth=1.8, markersize=4, label=region, alpha=0.9)

ax2.set_xlabel('年度', fontsize=11)
ax2.set_ylabel('転出率（人口千対）', fontsize=11)
ax2.set_title('地域別 転出率の年次推移', fontsize=12)
ax2.legend(fontsize=9, loc='upper right', ncol=2)
ax2.grid(alpha=0.4)
ax2.xaxis.set_major_locator(mticker.MultipleLocator(2))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig4_region.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig4_region.png")

# ================================================================
# ■ 推定結果サマリー出力
# ================================================================
print("\n" + "="*60)
print("[固定効果モデル 主要指標]")
print(f"  R² (within)  : {fe_res.rsquared:.4f}")
print(f"  R² (between) : {fe_res.rsquared_between:.4f}")
print(f"  観測数       : {fe_res.nobs}")
print(f"  エンティティ : {fe_res.entity_info['total']}")
print()
print("[係数（転出率への影響）]")
for var in EXOG_VARS:
    if var in fe_res.params:
        c = fe_res.params[var]
        p = fe_res.pvalues[var]
        sig = '***' if p < 0.001 else '** ' if p < 0.01 else '*  ' if p < 0.05 else '   '
        print(f"  {var:18s}: {c:+.4f}  {sig}  (p={p:.4f})")
print()
print("[Hausman 検定]")
print(f"  H = {H_stat:.4f},  df = {H_df},  p = {H_pval:.4f}")
if H_pval < 0.05:
    print("  → 固定効果モデルを採用（ランダム効果と有意差あり）")
else:
    print("  → p≥0.05（ランダム効果モデルも統計的に許容）")
print()
print("[生成図一覧]")
for f in ['2022_U5_1_fig1_ts.png', '2022_U5_1_fig2_fe_coef.png',
          '2022_U5_1_fig3_scatter.png', '2022_U5_1_fig4_region.png']:
    path = os.path.join(FIG_DIR, f)
    size = os.path.getsize(path) // 1024 if os.path.exists(path) else -1
    print(f"  {f} ({size} KB)")
print("="*60)
print("\n完了: 全4図を html/figures/ に保存しました。")