"""
2021年度 統計データ分析コンペティション 審査員奨励賞（大学生の部）
「コロナ禍における観光業への打撃と地域経済への影響」
手法: 差の差分析（DID）、時系列分析、地域比較
データ: SSDSE-B-2026.csv（47都道府県パネルデータ）
"""

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


import os
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

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

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

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)

# ===== 地域マップ =====
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東',
    '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12'
}

# ===== データ加工 =====
# 2018-2023年の47都道府県データ
df = df_b[df_b['年度'].between(2018, 2023)].copy()
df['地方区分'] = df['都道府県'].map(region_map)

# 旅館密度（施設数 per 万人）
df['旅館密度'] = df['旅館営業施設数（ホテルを含む）'] / (df['総人口'] / 10000)

# コロナダミー: 2020,2021年=1, それ以外=0
df['コロナダミー'] = df['年度'].isin([2020, 2021]).astype(int)

# 消費支出（月額・円）
df['消費支出'] = df['消費支出（二人以上の世帯）']

print("=== データ確認 ===")
print(f"観測数: {len(df)} 件（47都道府県×6年）")
print(f"年度: {sorted(df['年度'].unique())}")
print(f"消費支出 記述統計:\n{df['消費支出'].describe().round(0)}")

# ===== 図1: 全国平均消費支出の時系列（2018-2023） =====
ts = df.groupby('年度')['消費支出'].mean()
years = ts.index.tolist()

fig1, ax1 = plt.subplots(figsize=(9, 5))
ax1.axvspan(2019.5, 2021.5, alpha=0.15, color='red', label='COVID-19 影響期（2020-2021年）')
ax1.plot(years, ts.values / 10000, 'o-', color='#1565C0', linewidth=2.5, markersize=7,
         label='全国平均消費支出（二人以上世帯）')
# 変化率注記
for i, yr in enumerate(years):
    val = ts[yr] / 10000
    ax1.annotate(f'{val:.1f}万円', (yr, val), textcoords='offset points',
                 xytext=(0, 10), ha='center', fontsize=9.5, color='#1565C0')

ax1.axvline(2019.5, color='red', linestyle='--', alpha=0.5, linewidth=1)
ax1.axvline(2021.5, color='red', linestyle='--', alpha=0.5, linewidth=1)
ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('消費支出（万円/月）', fontsize=12)
ax1.set_title('図1：全国平均消費支出の時系列変化（2018–2023年）\nコロナ前後の消費動向', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10, loc='lower right')
ax1.set_xticks(years)
ax1.set_xticklabels([f'{y}年' for y in years])
ax1.grid(axis='y', alpha=0.3)
ax1.set_ylim(20, 35)

# コロナ前後の差
pre_avg  = ts[[2018, 2019]].mean() / 10000
post_avg = ts[[2020, 2021]].mean() / 10000
chg_pct  = (post_avg - pre_avg) / pre_avg * 100
print(f"\n=== 図1: 消費支出時系列 ===")
print(f"コロナ前平均: {pre_avg:.2f}万円")
print(f"コロナ期平均: {post_avg:.2f}万円")
print(f"変化率: {chg_pct:.1f}%")

fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2021_U5_3_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("fig1 保存完了")

# ===== 図2: 散布図 旅館密度 vs 2020年消費支出変化率（地域色分け・ラベル） =====
# 2019年と2020年の消費支出変化率を計算（都道府県別）
df_2019 = df[df['年度'] == 2019][['都道府県', '消費支出', '旅館密度', '地方区分']].copy()
df_2019.columns = ['都道府県', '消費支出_2019', '旅館密度_2019', '地方区分']
df_2020 = df[df['年度'] == 2020][['都道府県', '消費支出']].copy()
df_2020.columns = ['都道府県', '消費支出_2020']

scatter_df = pd.merge(df_2019, df_2020, on='都道府県')
scatter_df['変化率2020'] = (scatter_df['消費支出_2020'] - scatter_df['消費支出_2019']) / scatter_df['消費支出_2019'] * 100

print(f"\n=== 図2: 旅館密度 vs 消費支出変化率 ===")
print(f"変化率2020 記述統計:\n{scatter_df['変化率2020'].describe().round(2)}")

# 相関係数
r_val, p_val = stats.pearsonr(scatter_df['旅館密度_2019'], scatter_df['変化率2020'])
print(f"相関係数: r={r_val:.3f}, p={p_val:.4f}")

fig2, ax2 = plt.subplots(figsize=(11, 7))
for region, grp in scatter_df.groupby('地方区分'):
    color = region_colors[region]
    ax2.scatter(grp['旅館密度_2019'], grp['変化率2020'],
                color=color, label=region, s=60, alpha=0.85, zorder=3)
    for _, row in grp.iterrows():
        pref = row['都道府県']
        short = pref.replace('県', '').replace('府', '').replace('都', '').replace('道', '')
        ax2.annotate(short, (row['旅館密度_2019'], row['変化率2020']),
                     textcoords='offset points', xytext=(4, 3),
                     fontsize=7.5, color=color, alpha=0.9)

# 回帰直線
x_fit = np.linspace(scatter_df['旅館密度_2019'].min(), scatter_df['旅館密度_2019'].max(), 100)
slope, intercept = np.polyfit(scatter_df['旅館密度_2019'], scatter_df['変化率2020'], 1)
ax2.plot(x_fit, slope * x_fit + intercept, '--', color='#555', linewidth=1.5,
         label=f'回帰直線 (r={r_val:.2f}, p={p_val:.3f})')
ax2.axhline(0, color='gray', linewidth=0.8, linestyle='-', alpha=0.5)

ax2.set_xlabel('旅館密度（施設数 / 万人）', fontsize=12)
ax2.set_ylabel('消費支出変化率（2019→2020年, %）', fontsize=12)
ax2.set_title('図2：旅館密度と2020年消費支出変化率の関係\n観光依存度が高い県ほど消費が落ち込んだか？', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='upper right', ncol=2)
ax2.grid(alpha=0.3)

fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2021_U5_3_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("fig2 保存完了")

# ===== 図3: DID回帰係数プロット =====
# 観光依存度: 2019年旅館密度で分類（上位47%=高依存、それ以外=低依存）
prefs_2019 = df[df['年度'] == 2019][['都道府県', '旅館密度']].copy()
median_density = prefs_2019['旅館密度'].median()
high_tourism = set(prefs_2019[prefs_2019['旅館密度'] >= median_density]['都道府県'])

df['treat'] = df['都道府県'].isin(high_tourism).astype(int)
df['post']  = df['コロナダミー']
df['treat_post'] = df['treat'] * df['post']

# DID回帰: 消費支出 ~ treat + post + treat*post + 固定効果(簡易)
# 都道府県固定効果のためのダミー変数
X_vars = ['treat', 'post', 'treat_post']
X = sm.add_constant(df[X_vars])
y_did = df['消費支出']
model_did = sm.OLS(y_did, X).fit()

print(f"\n=== 図3: DID回帰結果 ===")
print(model_did.summary())

# 係数と95%CI
coef_names = ['定数項（Intercept）', '観光高依存（treat）', 'コロナ期（post）', 'DID交差項（treat×post）']
coefs = model_did.params.values
ci_low = model_did.conf_int()[0].values
ci_high = model_did.conf_int()[1].values
pvals_did = model_did.pvalues.values

# DID交差項のみハイライト（index=3）
did_coef  = coefs[3]
did_ci_lo = ci_low[3]
did_ci_hi = ci_high[3]
did_pval  = pvals_did[3]
print(f"DID交差項係数: {did_coef:.1f}円, 95%CI=[{did_ci_lo:.1f}, {did_ci_hi:.1f}], p={did_pval:.4f}")

fig3, ax3 = plt.subplots(figsize=(9, 5))
colors_coef = ['#888888', '#4e9af1', '#e05c5c', '#e65c00']
y_pos = np.arange(len(coef_names))

for i, (c, lo, hi, pv, nm) in enumerate(zip(coefs, ci_low, ci_high, pvals_did, coef_names)):
    color = colors_coef[i]
    lw = 3 if i == 3 else 1.5
    ax3.plot([lo, hi], [y_pos[i], y_pos[i]], '-', color=color, linewidth=lw, zorder=2)
    mk = '*' if pv < 0.05 else 'o'
    ms = 10 if i == 3 else 8
    ax3.plot(c, y_pos[i], mk, color=color, markersize=ms, zorder=3)
    sig_str = f"p={pv:.3f}" + (' *' if pv < 0.05 else ' n.s.')
    ax3.annotate(f'{c:+,.0f}円 ({sig_str})', (c, y_pos[i]),
                 textcoords='offset points', xytext=(8, 0), fontsize=9.5, color=color, va='center')

ax3.axvline(0, color='black', linewidth=1, linestyle='--', alpha=0.6)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(coef_names, fontsize=11)
ax3.set_xlabel('回帰係数（円/月）', fontsize=12)
ax3.set_title('図3：差の差（DID）分析の回帰係数と95%信頼区間\n観光高依存×コロナ期の交差項で打撃の差を推定', fontsize=13, fontweight='bold')
ax3.grid(axis='x', alpha=0.3)

# DID交差項を囲む矩形
ax3.axhspan(2.6, 3.4, alpha=0.08, color='#e65c00')

fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2021_U5_3_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("fig3 保存完了")

# ===== 図4: 箱ひげ図 観光高依存県 vs 低依存県のコロナ前後消費変化 =====
# 上位10県（旅館密度）vs 下位10県
top10_prefs    = prefs_2019.nlargest(10, '旅館密度')['都道府県'].tolist()
bottom10_prefs = prefs_2019.nsmallest(10, '旅館密度')['都道府県'].tolist()

print(f"\n=== 図4: 観光高依存TOP10 ===")
for p in top10_prefs:
    d = prefs_2019[prefs_2019['都道府県'] == p]['旅館密度'].values[0]
    print(f"  {p}: {d:.2f}")
print(f"観光低依存BOTTOM10: {bottom10_prefs}")

# コロナ前後の消費変化（2019→2020,2021の変化率）
df_pre   = df[df['年度'] == 2019][['都道府県', '消費支出']].rename(columns={'消費支出': '消費支出_pre'})
df_post2 = df[df['年度'].isin([2020, 2021])].groupby('都道府県')['消費支出'].mean().reset_index()
df_post2.columns = ['都道府県', '消費支出_post']
chg_df = pd.merge(df_pre, df_post2, on='都道府県')
chg_df['消費変化率'] = (chg_df['消費支出_post'] - chg_df['消費支出_pre']) / chg_df['消費支出_pre'] * 100

top10_chg    = chg_df[chg_df['都道府県'].isin(top10_prefs)]['消費変化率']
bottom10_chg = chg_df[chg_df['都道府県'].isin(bottom10_prefs)]['消費変化率']

print(f"\n高依存TOP10   消費変化率: mean={top10_chg.mean():.2f}%, median={top10_chg.median():.2f}%")
print(f"低依存BOTTOM10 消費変化率: mean={bottom10_chg.mean():.2f}%, median={bottom10_chg.median():.2f}%")

# t検定
t_stat, t_pval = stats.ttest_ind(top10_chg, bottom10_chg)
print(f"t検定: t={t_stat:.3f}, p={t_pval:.4f}")

fig4, ax4 = plt.subplots(figsize=(8, 6))
bp = ax4.boxplot(
    [bottom10_chg.values, top10_chg.values],
    labels=['観光低依存\n(旅館密度 下位10県)', '観光高依存\n(旅館密度 上位10県)'],
    patch_artist=True,
    medianprops=dict(color='black', linewidth=2)
)
bp['boxes'][0].set_facecolor('#AED6F1')
bp['boxes'][1].set_facecolor('#F1948A')
bp['boxes'][0].set_alpha(0.8)
bp['boxes'][1].set_alpha(0.8)

# 個別データ点
ax4.scatter([1]*len(bottom10_chg), bottom10_chg.values, color='#1565C0', s=50, zorder=5, alpha=0.7)
ax4.scatter([2]*len(top10_chg),    top10_chg.values,    color='#C0392B', s=50, zorder=5, alpha=0.7)

ax4.axhline(0, color='gray', linewidth=1, linestyle='--', alpha=0.6)
ax4.set_ylabel('消費支出変化率（2019→2020-21年平均, %）', fontsize=11)
ax4.set_title(f'図4：観光依存度別のコロナ禍における消費支出変化\nt検定: t={t_stat:.2f}, p={t_pval:.3f}', fontsize=13, fontweight='bold')
ax4.grid(axis='y', alpha=0.3)

# 平均値マーク
ax4.plot(1, bottom10_chg.mean(), 'D', color='#1565C0', markersize=10, zorder=6, label='平均値')
ax4.plot(2, top10_chg.mean(),    'D', color='#C0392B', markersize=10, zorder=6)
ax4.legend(fontsize=10)

fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2021_U5_3_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("fig4 保存完了")

# ===== 統計値サマリー出力 =====
print("\n" + "="*60)
print("【HTML作成用 統計値サマリー】")
print("="*60)
print(f"[時系列] コロナ前(2018-19)全国平均消費支出: {pre_avg:.2f}万円/月")
print(f"[時系列] コロナ期(2020-21)全国平均消費支出: {post_avg:.2f}万円/月")
print(f"[時系列] コロナ前後変化率: {chg_pct:+.1f}%")
print(f"[散布図] 旅館密度×消費変化率 相関係数: r={r_val:.3f}, p={p_val:.4f}")
print(f"[DID]   DID交差項係数: {did_coef:+,.0f}円/月")
print(f"[DID]   95%CI: [{did_ci_lo:,.0f}, {did_ci_hi:,.0f}]")
print(f"[DID]   p値: {did_pval:.4f}")
print(f"[箱ひげ] 高依存TOP10 平均変化率: {top10_chg.mean():+.2f}%")
print(f"[箱ひげ] 低依存BTM10 平均変化率: {bottom10_chg.mean():+.2f}%")
print(f"[箱ひげ] t検定: t={t_stat:.3f}, p={t_pval:.4f}")
print(f"[観光TOP10] {', '.join(top10_prefs)}")