"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：周辺環境からみた新型コロナウイルス感染症の要因分析
賞：審査員奨励賞（大学生・一般の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ 2012〜2023）
  分析方針：延べ宿泊者数（人口1人あたり）をモビリティ／接触密度の代理指標として，
           COVID-19 拡大との関連要因を都道府県別クロスセクションで探索する。
           2019年（コロナ前）断面での重回帰分析と
           2019→2020 変化率との比較を中心に展開する。

  Step1. 時系列グラフ（fig1）：延べ宿泊者数（人口比）2012〜2023
  Step2. 相関ヒートマップ（fig2）：宿泊密度と環境変数の相関
  Step3. OLS標準化係数（fig3）：2019年断面の重回帰分析
  Step4. 変化率横棒グラフ（fig4）：2019→2020 都道府県別宿泊者変化率

【説明変数】
  人口密度       : 総人口 / 面積（代理：総人口）
  高齢化率       : 65歳以上人口 / 総人口
  観光客数代理   : 旅館営業施設客室数（ホテル含む）
  交通・通信費   : 交通・通信費（二人以上の世帯）
  消費支出       : 消費支出（二人以上の世帯）
  求人数         : 月間有効求人数（一般）
  地価（住宅地） : 標準価格（平均価格）（住宅地）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計数理研究所 SSDSE より公表の実データ

【データサイエンス学習ポイント】
  1. パネルデータの年次断面抽出と時系列可視化
  2. 相関ヒートマップによる多変数の俯瞰
  3. 標準化係数（Beta）による変数重要度の比較
  4. COVID前後の変化率分析
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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_6_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
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_DIR = 'data/raw'
DATA_B = os.path.join(DATA_DIR, '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()

# 数値列を変換
num_cols = [
    '総人口', '65歳以上人口', '延べ宿泊者数', '外国人延べ宿泊者数',
    '旅館営業施設客室数（ホテルを含む）', '旅館営業施設数（ホテルを含む）',
    '年平均気温', '交通・通信費（二人以上の世帯）', '消費支出（二人以上の世帯）',
    '月間有効求人数（一般）', '標準価格（平均価格）（住宅地）',
]
for c in num_cols:
    if c in df_b.columns:
        df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

print("=" * 60)
print(f"■ SSDSE-B 都道府県データ: {df_b.shape[0]}行 × {df_b.shape[1]}列")
print(f"  年度範囲: {df_b['年度'].min()} 〜 {df_b['年度'].max()}")
print(f"  都道府県数: {df_b['都道府県'].nunique()}")
print("=" * 60)

# ================================================================
# ■ 宿泊密度（一人あたり延べ宿泊者数）を計算
# ================================================================
df_b['宿泊密度'] = df_b['延べ宿泊者数'] / df_b['総人口'].clip(1)
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'].clip(1) * 100  # %
df_b['外国人宿泊率'] = df_b['外国人延べ宿泊者数'] / df_b['延べ宿泊者数'].clip(1) * 100  # %

# ================================================================
# ■ 時系列データ（全国集計）
# ================================================================
# 全国平均（都道府県平均）の年次推移
ts = df_b.groupby('年度')['宿泊密度'].mean().reset_index()
ts.columns = ['年度', '宿泊密度_平均']

print("\n■ 宿泊密度の年次推移（全国平均）:")
print(ts.to_string(index=False))

# ================================================================
# ■ 2019年断面データ（COVID前）
# ================================================================
df_2019 = df_b[df_b['年度'] == 2019].copy()
df_2020 = df_b[df_b['年度'] == 2020].copy()

print(f"\n■ 2019年断面: N={len(df_2019)}")
print(f"  宿泊密度 平均={df_2019['宿泊密度'].mean():.2f}, 中央値={df_2019['宿泊密度'].median():.2f}")

# ================================================================
# ■ 分析変数（2019年断面）
# ================================================================
# 説明変数
ENV_VARS = {
    '高齢化率(%)':      '高齢化率',
    '宿泊施設客室数':   '旅館営業施設客室数（ホテルを含む）',
    '交通通信費(万円)': '交通・通信費（二人以上の世帯）',
    '消費支出(万円)':   '消費支出（二人以上の世帯）',
    '有効求人数(千)':   '月間有効求人数（一般）',
    '地価住宅地':       '標準価格（平均価格）（住宅地）',
}

df_2019 = df_2019.copy()
df_2019['高齢化率'] = df_2019['高齢化率']
df_2019['宿泊施設客室数']   = pd.to_numeric(df_2019['旅館営業施設客室数（ホテルを含む）'], errors='coerce')
df_2019['交通通信費(万円)'] = pd.to_numeric(df_2019['交通・通信費（二人以上の世帯）'], errors='coerce') / 10000
df_2019['消費支出(万円)']   = pd.to_numeric(df_2019['消費支出（二人以上の世帯）'], errors='coerce') / 10000
df_2019['有効求人数(千)']   = pd.to_numeric(df_2019['月間有効求人数（一般）'], errors='coerce') / 1000
df_2019['地価住宅地']       = pd.to_numeric(df_2019['標準価格（平均価格）（住宅地）'], errors='coerce')

XVARS = ['高齢化率', '宿泊施設客室数', '交通通信費(万円)', '消費支出(万円)', '有効求人数(千)', '地価住宅地']
XVAR_LABELS = {
    '高齢化率':         '高齢化率（%）',
    '宿泊施設客室数':   '宿泊施設客室数',
    '交通通信費(万円)': '交通・通信費（万円）',
    '消費支出(万円)':   '消費支出（万円）',
    '有効求人数(千)':   '月間有効求人数（千件）',
    '地価住宅地':       '住宅地地価（千円/m²）',
}

df_2019_clean = df_2019[['都道府県', '宿泊密度'] + XVARS].dropna()
N = len(df_2019_clean)
print(f"\n■ 2019年 分析用データ: N={N}（欠損除外後）")

# ================================================================
# ■ 相関分析
# ================================================================
corr_cols = ['宿泊密度'] + XVARS
corr_mat = df_2019_clean[corr_cols].corr()
print("\n■ 相関行列（2019年断面）:")
print(corr_mat.round(3))

# ================================================================
# ■ OLS重回帰（2019年断面）
# ================================================================
y_reg = df_2019_clean['宿泊密度'].values
X_reg = df_2019_clean[XVARS].values

# 標準化（Betaコefのため）
y_std = (y_reg - y_reg.mean()) / y_reg.std()
X_std = (X_reg - X_reg.mean(axis=0)) / X_reg.std(axis=0)

model_std = sm.OLS(y_std, sm.add_constant(X_std)).fit()
model_raw = sm.OLS(y_reg, sm.add_constant(X_reg)).fit()

print("\n■ OLS重回帰（標準化係数、2019年断面）:")
print(f"  R² = {model_std.rsquared:.4f}")
print(f"  {'変数':<20} {'Beta':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 50)
for i, var in enumerate(XVARS):
    beta = model_std.params[i + 1]
    p    = model_std.pvalues[i + 1]
    sig  = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {var:<20} {beta:>8.4f} {p:>10.4f} {sig:>6}")

# ================================================================
# ■ 2019→2020 変化率
# ================================================================
df_merge = pd.merge(
    df_2019[['都道府県', '宿泊密度']].rename(columns={'宿泊密度': '宿泊密度_2019'}),
    df_2020[['都道府県', '宿泊密度']].rename(columns={'宿泊密度': '宿泊密度_2020'}),
    on='都道府県'
)
df_merge['変化率'] = (df_merge['宿泊密度_2020'] - df_merge['宿泊密度_2019']) / df_merge['宿泊密度_2019'].clip(0.001) * 100
df_merge = df_merge.sort_values('変化率').reset_index(drop=True)

print("\n■ 2019→2020 宿泊密度変化率 上位・下位5件:")
print(df_merge[['都道府県', '変化率']].tail(5).to_string())
print("...")
print(df_merge[['都道府県', '変化率']].head(5).to_string())


# ================================================================
# ■ 図の生成（4枚）
# ================================================================

# ────────────────────────────────────────────────────────────────
# 図1: 延べ宿泊者数（一人あたり）の時系列 2012〜2023 — COVID期シェード
# ────────────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(10, 5))

years  = ts['年度'].values
values = ts['宿泊密度_平均'].values

# COVID期（2020〜2021）シェード
ax1.axvspan(2019.5, 2021.5, alpha=0.18, color='#C62828', label='COVID-19 流行期（2020〜2021）')

ax1.plot(years, values, 'o-', color='#1565C0', linewidth=2.2, markersize=6,
         label='延べ宿泊者数（人口一人あたり）全国平均', zorder=3)

# 2019→2020 矢印注釈
y19 = ts[ts['年度'] == 2019]['宿泊密度_平均'].values
y20 = ts[ts['年度'] == 2020]['宿泊密度_平均'].values
if len(y19) and len(y20):
    ax1.annotate('', xy=(2020, y20[0]), xytext=(2019, y19[0]),
                 arrowprops=dict(arrowstyle='->', color='#C62828', lw=2.0))
    ax1.text(2019.5, (y19[0] + y20[0]) / 2 + 0.05, '▼COVID衝撃', color='#C62828',
             fontsize=9, ha='center')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('延べ宿泊者数（人口一人あたり：泊/人）', fontsize=12)
ax1.set_title('延べ宿泊者数（人口あたり）の推移 2012〜2023\n（都道府県平均・SSDSE-B実データ）',
              fontsize=13, fontweight='bold')
ax1.set_xticks(years)
ax1.set_xticklabels([str(y) for y in years], rotation=45, ha='right')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_U5_6_fig1_ts.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2022_U5_6_fig1_ts.png")

# ────────────────────────────────────────────────────────────────
# 図2: 相関ヒートマップ（宿泊密度と環境変数）
# ────────────────────────────────────────────────────────────────
fig2, ax2 = plt.subplots(figsize=(8, 7))

# 列ラベルを短くする
label_map = {
    '宿泊密度':         '宿泊密度\n(一人あたり)',
    '高齢化率':         '高齢化率\n(%)',
    '宿泊施設客室数':   '宿泊施設\n客室数',
    '交通通信費(万円)': '交通・\n通信費',
    '消費支出(万円)':   '消費支出',
    '有効求人数(千)':   '有効求人数',
    '地価住宅地':       '住宅地地価',
}
cm_data = df_2019_clean[corr_cols].rename(columns=label_map)
cm = cm_data.corr()

im = ax2.imshow(cm.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax2, fraction=0.046, pad=0.04, label='相関係数')

labels = list(cm.columns)
ax2.set_xticks(range(len(labels)))
ax2.set_yticks(range(len(labels)))
ax2.set_xticklabels(labels, fontsize=9)
ax2.set_yticklabels(labels, fontsize=9)

for i in range(len(labels)):
    for j in range(len(labels)):
        v = cm.values[i, j]
        tc = 'white' if abs(v) > 0.6 else 'black'
        ax2.text(j, i, f'{v:.2f}', ha='center', va='center', fontsize=9, color=tc)

ax2.set_title('相関ヒートマップ：宿泊密度と環境変数\n（2019年断面, 都道府県 N=47）',
              fontsize=12, fontweight='bold')
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_U5_6_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_U5_6_fig2_corr.png")

# ────────────────────────────────────────────────────────────────
# 図3: OLS標準化係数プロット（2019年断面）
# ────────────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(8, 5))

betas = model_std.params[1:]
ses   = model_std.bse[1:]
pvals = model_std.pvalues[1:]

xlabels = [XVAR_LABELS.get(v, v) for v in XVARS]
colors3 = ['#C62828' if p < 0.05 else '#90A4AE' for p in pvals]
y_pos   = np.arange(len(XVARS))

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, zorder=5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(xlabels, fontsize=10)
ax3.set_xlabel('標準化回帰係数 Beta（±1.96SE）', fontsize=11)
ax3.set_title(f'OLS標準化係数（2019年断面, N={N}都道府県）\n'
              f'R²={model_std.rsquared:.3f}  目的変数：延べ宿泊者数（人口あたり）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

sig_patch   = mpatches.Patch(color='#C62828', alpha=0.78, label='p < 0.05')
insig_patch = mpatches.Patch(color='#90A4AE', alpha=0.78, label='n.s.')
ax3.legend(handles=[sig_patch, insig_patch], 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 ''
    ax3.text(b + (0.01 if b >= 0 else -0.01), i,
             f' {b:.3f}{sig}',
             va='center', ha='left' if b >= 0 else 'right', fontsize=8.5)

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

# ────────────────────────────────────────────────────────────────
# 図4: 都道府県別宿泊者変化率 2019→2020（横棒グラフ）
# ────────────────────────────────────────────────────────────────
fig4, ax4 = plt.subplots(figsize=(9, 13))

colors4 = ['#C62828' if v < 0 else '#1565C0' for v in df_merge['変化率']]
y_pos4  = np.arange(len(df_merge))

ax4.barh(y_pos4, df_merge['変化率'], color=colors4, alpha=0.75,
         edgecolor='white', height=0.75)
ax4.axvline(0, color='black', linewidth=0.8)
ax4.set_yticks(y_pos4)
ax4.set_yticklabels(df_merge['都道府県'], fontsize=9)
ax4.set_xlabel('変化率（%）', fontsize=11)
ax4.set_title('延べ宿泊者数（人口あたり）の変化率：2019年→2020年\n（都道府県別, SSDSE-B実データ）',
              fontsize=12, fontweight='bold')
ax4.grid(axis='x', alpha=0.3)

dec_patch = mpatches.Patch(color='#C62828', alpha=0.75, label='減少（赤）')
inc_patch = mpatches.Patch(color='#1565C0', alpha=0.75, label='増加（青）')
ax4.legend(handles=[dec_patch, inc_patch], fontsize=9, loc='lower right')

# 国全体平均線
nat_mean = df_merge['変化率'].mean()
ax4.axvline(nat_mean, color='#FF8F00', linestyle='--', linewidth=1.5,
            label=f'全国平均 {nat_mean:.1f}%')
ax4.legend(handles=[dec_patch, inc_patch,
                    mpatches.Patch(color='#FF8F00', label=f'全国平均 {nat_mean:.1f}%')],
           fontsize=9, loc='lower right')

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

print("\n全図の生成完了（4枚）")
print("  2022_U5_6_fig1_ts.png     : 宿泊密度の時系列（COVID期シェード）")
print("  2022_U5_6_fig2_corr.png   : 相関ヒートマップ")
print("  2022_U5_6_fig3_coef.png   : OLS標準化係数プロット（2019年断面）")
print("  2022_U5_6_fig4_change.png : 都道府県別変化率 2019→2020（横棒）")