"""
2019_U5_3_shorei.py
都道府県別観光消費の地域経済波及効果：宿泊・消費支出データによる重回帰分析
審査員奨励賞（大学生部門）2019年度
使用データ：SSDSE-B-2026.csv（実データのみ）
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2019_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
from statsmodels.stats.outliers_influence import variance_inflation_factor
from matplotlib.patches import Patch

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)

print("=== SSDSE-B カラム一覧 ===")
print(df_b.columns.tolist())
print()

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

# ── 2019年度クロスセクションデータ ─────────────────────────────────────────
df19 = df_b[df_b['年度'] == 2019].copy()
df19 = df19.dropna(subset=[
    '消費支出（二人以上の世帯）', '旅館営業施設数（ホテルを含む）',
    '延べ宿泊者数', '総人口', '65歳以上人口', '外国人延べ宿泊者数'
]).reset_index(drop=True)

# 代理変数の作成
df19['旅館密度']  = (df19['旅館営業施設数（ホテルを含む）'].astype(float)
                     / df19['総人口'].astype(float) * 10000)   # 旅館・ホテル数/万人
df19['宿泊密度']  = (df19['延べ宿泊者数'].astype(float)
                     / df19['総人口'].astype(float))            # 延べ宿泊者数/人
df19['高齢化率']  = (df19['65歳以上人口'].astype(float)
                     / df19['総人口'].astype(float) * 100)      # 65歳以上%
df19['外国人比率'] = (df19['外国人延べ宿泊者数'].astype(float)
                      / df19['延べ宿泊者数'].astype(float) * 100)  # インバウンド%
df19['消費支出']  = df19['消費支出（二人以上の世帯）'].astype(float)  # 円

# 地域ラベル
df19['地域'] = df19['都道府県'].map(region_map)
print("地域マッピング未対応:", df19[df19['地域'].isna()]['都道府県'].tolist())

print("=== 2019年度データ（47都道府県）基本統計 ===")
summary_cols = ['旅館密度', '宿泊密度', '外国人比率', '高齢化率', '消費支出']
print(df19[summary_cols].describe().round(2))
print()

# ── 重回帰分析（OLS） ────────────────────────────────────────────────────────
X_vars = ['旅館密度', '宿泊密度', '高齢化率', '外国人比率']
y_var  = '消費支出'

df_reg = df19[X_vars + [y_var, '都道府県', '地域']].dropna().reset_index(drop=True)
print(f"回帰分析サンプル数: {len(df_reg)}")

X_raw = df_reg[X_vars].values.astype(float)
y_raw = df_reg[y_var].values.astype(float)

# OLS（非標準化）
X_ols = sm.add_constant(X_raw)
model = sm.OLS(y_raw, X_ols).fit()
print("=== 重回帰分析結果（非標準化） ===")
print(f"  R² = {model.rsquared:.4f}, Adj.R² = {model.rsquared_adj:.4f}")
print(f"  F-stat = {model.fvalue:.3f}, F p-val = {model.f_pvalue:.6f}")
for i, v in enumerate(X_vars):
    stars = ('***' if model.pvalues[i+1] < 0.001 else
             '**' if model.pvalues[i+1] < 0.01 else
             '*' if model.pvalues[i+1] < 0.05 else 'n.s.')
    print(f"  {v}: coef={model.params[i+1]:.2f}, p={model.pvalues[i+1]:.4f} {stars}")
print()

# 標準化（z-score）
X_mean = X_raw.mean(axis=0)
X_sd   = X_raw.std(axis=0, ddof=1)   # ddof=1 for sample std
y_mean = y_raw.mean()
y_sd   = y_raw.std(ddof=1)
X_std  = (X_raw - X_mean) / X_sd
y_std  = (y_raw - y_mean) / y_sd

X_ols_std = sm.add_constant(X_std)
model_std = sm.OLS(y_std, X_ols_std).fit()

beta_std  = model_std.params[1:]
se_std    = model_std.bse[1:]
pvals_std = model_std.pvalues[1:]

print("=== 標準化偏回帰係数 ===")
for i, v in enumerate(X_vars):
    stars = ('***' if pvals_std[i] < 0.001 else
             '**' if pvals_std[i] < 0.01 else
             '*' if pvals_std[i] < 0.05 else 'n.s.')
    print(f"  {v}: β={beta_std[i]:.4f}, SE={se_std[i]:.4f}, p={pvals_std[i]:.4f} {stars}")
print()

# VIF
vif_vals = [variance_inflation_factor(X_raw, i) for i in range(len(X_vars))]
print("=== VIF（多重共線性診断） ===")
for v, vf in zip(X_vars, vif_vals):
    flag = " ← 要注意(VIF>5)" if vf > 5 else ""
    print(f"  {v}: VIF = {vf:.2f}{flag}")
print()

# Pearson相関
r_val, p_val = stats.pearsonr(df19['旅館密度'], df19['消費支出'])
print(f"=== Pearson相関：旅館密度 vs 消費支出 ===")
print(f"  r = {r_val:.4f}, p = {p_val:.6f}")
print()

# ═══════════════════════════════════════════════════════════════════
# Figure 1: 都道府県別旅館密度ランキング（棒グラフ）
# ═══════════════════════════════════════════════════════════════════
df19_sorted = df19.sort_values('旅館密度', ascending=False).reset_index(drop=True)
bar_colors  = [region_colors[df19_sorted.loc[i, '地域']] for i in range(len(df19_sorted))]
national_avg = df19['旅館密度'].mean()

fig, ax = plt.subplots(figsize=(14, 6))
ax.bar(range(len(df19_sorted)), df19_sorted['旅館密度'],
       color=bar_colors, alpha=0.85, edgecolor='white', linewidth=0.5)
ax.axhline(national_avg, color='black', linewidth=1.8, linestyle='--')

ax.set_xticks(range(len(df19_sorted)))
ax.set_xticklabels(df19_sorted['都道府県'], rotation=90, fontsize=8)
ax.set_ylabel('旅館・ホテル数（施設数/万人）', fontsize=12)
ax.set_title('図1：都道府県別宿泊施設密度（旅館・ホテル数/万人）ランキング〈2019年〉',
             fontsize=13, fontweight='bold', pad=14)

legend_patches = [Patch(color=v, label=k, alpha=0.85) for k, v in region_colors.items()]
legend_patches.append(
    plt.Line2D([0], [0], color='black', linewidth=1.8, linestyle='--',
               label=f'全国平均 {national_avg:.1f}施設/万人')
)
ax.legend(handles=legend_patches, fontsize=9, loc='upper right', ncol=2)

ax.set_xlim(-0.6, len(df19_sorted) - 0.4)
ax.set_ylim(0, df19_sorted['旅館密度'].max() * 1.15)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_3_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("Figure 1 saved.")

# ═══════════════════════════════════════════════════════════════════
# Figure 2: 宿泊施設密度 vs 消費支出の散布図
# ═══════════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(10, 7))

for region, grp in df19.groupby('地域'):
    ax.scatter(grp['旅館密度'], grp['消費支出'] / 10000,
               color=region_colors[region], s=60, alpha=0.85,
               label=region, zorder=3)

# 都道府県ラベル
for _, row in df19.iterrows():
    ax.annotate(row['都道府県'], (row['旅館密度'], row['消費支出'] / 10000),
                fontsize=6.5, ha='left', va='bottom', alpha=0.75,
                xytext=(2, 2), textcoords='offset points')

# 回帰直線
slope, intercept, r_lin, p_lin, _ = stats.linregress(
    df19['旅館密度'], df19['消費支出'] / 10000)
x_plot = np.linspace(df19['旅館密度'].min() - 0.5, df19['旅館密度'].max() + 0.5, 200)
ax.plot(x_plot, intercept + slope * x_plot, color='#333333',
        linewidth=1.8, linestyle='-', zorder=2)

p_str = f'p = {p_lin:.4f}' if p_lin >= 0.0001 else 'p < 0.0001'
ax.text(0.97, 0.05,
        f'r = {r_lin:.3f}\n{p_str}\nR² = {r_lin**2:.3f}',
        transform=ax.transAxes, ha='right', va='bottom', fontsize=11,
        bbox=dict(boxstyle='round,pad=0.4', facecolor='white', alpha=0.85))

ax.set_xlabel('宿泊施設密度（旅館・ホテル数/万人）', fontsize=12)
ax.set_ylabel('消費支出（万円/世帯）', fontsize=12)
ax.set_title('図2：宿泊施設密度と消費支出の関係〈2019年・47都道府県〉',
             fontsize=13, fontweight='bold', pad=12)

legend_patches = [Patch(color=v, label=k, alpha=0.85) for k, v in region_colors.items()]
ax.legend(handles=legend_patches, fontsize=9, loc='upper left', framealpha=0.85)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_3_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("Figure 2 saved.")

# ═══════════════════════════════════════════════════════════════════
# Figure 3: 地域別消費支出の時系列推移（2012–2022）
# ═══════════════════════════════════════════════════════════════════
df_ts = df_b[df_b['年度'].between(2012, 2022)].copy()
df_ts['地域'] = df_ts['都道府県'].map(region_map)
df_ts = df_ts.dropna(subset=['消費支出（二人以上の世帯）', '地域'])

national_ts = (df_ts.groupby('年度')['消費支出（二人以上の世帯）'].mean() / 10000)
region_ts   = (df_ts.groupby(['年度', '地域'])['消費支出（二人以上の世帯）'].mean() / 10000)

fig, ax = plt.subplots(figsize=(11, 6))

for region, color in region_colors.items():
    lvl = region_ts.index.get_level_values('地域')
    if region in lvl:
        data = region_ts.xs(region, level='地域')
        ax.plot(data.index, data.values, color=color, linewidth=2.0,
                marker='o', markersize=4.5, label=region, alpha=0.9)

ax.plot(national_ts.index, national_ts.values, color='black',
        linewidth=3.0, marker='D', markersize=5, label='全国平均', zorder=5)

# COVID-19 シェーディング
ax.axvspan(2019.7, 2021.3, alpha=0.07, color='red')
ax.text(2020.5, national_ts.max() * 1.01, 'COVID-19', ha='center',
        fontsize=9, color='red', alpha=0.7)

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('消費支出（万円/世帯）', fontsize=12)
ax.set_title('図3：地域別消費支出の時系列推移〈2012–2022年〉',
             fontsize=13, fontweight='bold', pad=12)
ax.set_xticks(range(2012, 2023))
ax.set_xticklabels(range(2012, 2023), rotation=45, fontsize=9)
ax.legend(fontsize=9, loc='lower left', framealpha=0.85)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_3_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("Figure 3 saved.")

# ═══════════════════════════════════════════════════════════════════
# Figure 4: 標準化偏回帰係数プロット（横棒 + 95%CI）
# ═══════════════════════════════════════════════════════════════════
var_labels = [
    '宿泊施設密度\n（旅館・ホテル数/万人）',
    '延べ宿泊者数\n（人口あたり）',
    '高齢化率\n（65歳以上人口比）',
    'インバウンド比率\n（外国人宿泊者割合）',
]

ci_lower = beta_std - 1.96 * se_std
ci_upper = beta_std + 1.96 * se_std

def bar_color(b, p):
    if p >= 0.05:
        return '#aaaaaa'
    return '#e05c5c' if b > 0 else '#4e9af1'

def sig_label(p):
    if p < 0.001:   return '***'
    elif p < 0.01:  return '**'
    elif p < 0.05:  return '*'
    else:           return 'n.s.'

colors4 = [bar_color(b, p) for b, p in zip(beta_std, pvals_std)]
y_pos   = np.arange(len(X_vars))

fig, ax = plt.subplots(figsize=(9, 5.5))
ax.barh(y_pos, beta_std, color=colors4, alpha=0.85, height=0.55, zorder=3)
ax.errorbar(beta_std, y_pos,
            xerr=[beta_std - ci_lower, ci_upper - beta_std],
            fmt='none', color='#333333', capsize=5, linewidth=1.5, zorder=4)

for i, (b, ci_u, ci_l, p) in enumerate(zip(beta_std, ci_upper, ci_lower, pvals_std)):
    x_ann = ci_u + 0.015 if b >= 0 else ci_l - 0.015
    ha_ann = 'left' if b >= 0 else 'right'
    ax.text(x_ann, i, sig_label(p), ha=ha_ann, va='center',
            fontsize=11, color='#333333', fontweight='bold')

ax.axvline(0, color='black', linewidth=0.9)
ax.set_yticks(y_pos)
ax.set_yticklabels(var_labels, fontsize=10)
ax.set_xlabel('標準化偏回帰係数 (β)', fontsize=12)
ax.set_title('図4：重回帰分析の標準化偏回帰係数〈2019年・47都道府県〉',
             fontsize=13, fontweight='bold', pad=12)

legend_items = [
    Patch(color='#e05c5c', alpha=0.85, label='正の有意効果'),
    Patch(color='#4e9af1', alpha=0.85, label='負の有意効果'),
    Patch(color='#aaaaaa', alpha=0.85, label='非有意 (n.s.)'),
]
ax.legend(handles=legend_items, fontsize=9, loc='lower right')
ax.text(0.02, 0.03,
        f'R² = {model.rsquared:.3f}, Adj.R² = {model.rsquared_adj:.3f}  (N={len(df_reg)})',
        transform=ax.transAxes, fontsize=9, color='#555555')

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
x_min = min(ci_lower.min() - 0.20, -0.40)
x_max = max(ci_upper.max() + 0.20, 0.60)
ax.set_xlim(x_min, x_max)

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_3_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("Figure 4 saved.")

# ── 最終確認 ────────────────────────────────────────────────────────────────
print()
print("=== 最終確認 ===")
for i in range(1, 5):
    path = os.path.join(FIG_DIR, f'2019_U5_3_fig{i}.png')
    status = 'OK' if os.path.exists(path) else 'MISSING'
    size   = os.path.getsize(path) // 1024 if os.path.exists(path) else 0
    print(f"  fig{i}: {status} ({size}KB)")

print()
print("DONE: 2019_U5_3_shorei")