コロナ禍の地域経済
第一波による雇用・消費への影響分析

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)

# 分析期間: 2012〜2023
df_b = df_b[df_b['年度'].between(2012, 2023)].copy()

df_b['COVID期'] = df_b['年度'].isin([2020, 2021]).astype(int)

# 消費支出（目的変数）
df_b['消費支出'] = pd.to_numeric(df_b['消費支出（二人以上の世帯）'], errors='coerce')

# 求人倍率（月間有効求人数 / 月間有効求職者数）
df_b['求人倍率'] = (
    pd.to_numeric(df_b['月間有効求人数（一般）'], errors='coerce') /
    pd.to_numeric(df_b['月間有効求職者数（一般）'], errors='coerce')
)

# 旅館密度（旅館営業施設数 / 総人口 × 10,000）
df_b['旅館密度'] = (
    pd.to_numeric(df_b['旅館営業施設数（ホテルを含む）'], errors='coerce') /
    pd.to_numeric(df_b['総人口'], errors='coerce') * 10000
)

# 高齢化率（65歳以上人口 / 総人口 × 100）
df_b['高齢化率'] = (
    pd.to_numeric(df_b['65歳以上人口'], errors='coerce') /
    pd.to_numeric(df_b['総人口'], errors='coerce') * 100
)

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

df_b['地域'] = df_b['都道府県'].map(region_map)

df_2019_density = df_b[df_b['年度'] == 2019][['都道府県', '旅館密度']].dropna()
top10_prefs  = df_2019_density.nlargest(10, '旅館密度')['都道府県'].tolist()
bot10_prefs  = df_2019_density.nsmallest(10, '旅館密度')['都道府県'].tolist()

df_b['観光グループ'] = df_b['都道府県'].apply(
    lambda x: '高観光依存（上位10県）' if x in top10_prefs
    else ('低観光依存（下位10県）' if x in bot10_prefs else None)
)

print("\n=== 旅館密度 上位10県 ===")
print(df_2019_density.nlargest(10, '旅館密度')[['都道府県', '旅館密度']].to_string(index=False))
print("\n=== 旅館密度 下位10県 ===")
print(df_2019_density.nsmallest(10, '旅館密度')[['都道府県', '旅館密度']].to_string(index=False))

df_2019c = df_b[df_b['年度'] == 2019][['都道府県', '消費支出', '地域']].rename(columns={'消費支出': '消費2019'})
df_2020c = df_b[df_b['年度'] == 2020][['都道府県', '消費支出']].rename(columns={'消費支出': '消費2020'})
df_change = df_2019c.merge(df_2020c, on='都道府県').dropna()
df_change['変化率'] = (df_change['消費2020'] - df_change['消費2019']) / df_change['消費2019'] * 100
df_change = df_change.sort_values('変化率')

print("\n=== 2019→2020年 消費支出変化率（下位5県） ===")
print(df_change.head(5)[['都道府県', '変化率']].to_string(index=False))
print("\n=== 2019→2020年 消費支出変化率（上位5県） ===")
print(df_change.tail(5)[['都道府県', '変化率']].to_string(index=False))

print("\n" + "="*60)
print("=== HTML転記用 統計サマリ ===")
print("="*60)
print(f"分析対象: 47都道府県 × 12年度 = {len(df_panel_all)}観測値（欠損除外後）")
print(f"モデル: {model_type}")
print(f"R² = {r2:.4f}")
print(f"N = {n_obs}")
print()
for var, (coef, se, pval) in coef_results.items():
    sig = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.'))
    ci_lo, ci_hi = coef - 1.96*se, coef + 1.96*se
    print(f"{var}: β={coef:+.1f} (95%CI: {ci_lo:+.1f}〜{ci_hi:+.1f}), SE={se:.1f}, p={pval:.4f} {sig}")
print()
# 消費支出2020変化
covid_mean = df_change['変化率'].mean()
covid_worst = df_change.iloc[0]
covid_best  = df_change.iloc[-1]
print(f"2019→2020消費変化（全国平均）: {covid_mean:+.2f}%")
print(f"最大減少: {covid_worst['都道府県']} {covid_worst['変化率']:+.2f}%")
print(f"最大増加: {covid_best['都道府県']} {covid_best['変化率']:+.2f}%")
print()
# 求人倍率2020
job_2019 = df_b[df_b['年度']==2019]['求人倍率'].mean()
job_2020 = df_b[df_b['年度']==2020]['求人倍率'].mean()
print(f"求人倍率（全国平均）: 2019年={job_2019:.2f} → 2020年={job_2020:.2f} (変化: {job_2020-job_2019:+.2f})")
print()
# 観光依存グループ比較
print("=== 観光依存グループ別 2019→2020消費変化率 ===")
df_change_group = df_change.copy()
df_change_group['観光グループ'] = df_change_group['都道府県'].apply(
    lambda x: '高観光依存（上位10）' if x in top10_prefs
    else ('低観光依存（下位10）' if x in bot10_prefs else None)
)
for grp in ['高観光依存（上位10）', '低観光依存（下位10）']:
    sub = df_change_group[df_change_group['観光グループ']==grp]['変化率']
    print(f"{grp}: 平均 {sub.mean():+.2f}% (N={len(sub)})")
print("="*60)

fig, ax = plt.subplots(figsize=(10, 5.5))

region_yearly = df_b.groupby(['年度', '地域'])['消費支出'].mean().reset_index()

for region, color in region_colors.items():
    sub = region_yearly[region_yearly['地域'] == region].sort_values('年度')
    if len(sub) > 0:
        ax.plot(sub['年度'], sub['消費支出'] / 1000, marker='o', markersize=4,
                label=region, color=color, linewidth=1.8)

# COVID期シェーディング
ax.axvspan(2019.5, 2021.5, alpha=0.15, color='#e74c3c', zorder=0)
ax.axvline(x=2020, color='#e74c3c', linestyle='--', linewidth=1.2, alpha=0.7)
ax.text(2020.05, ax.get_ylim()[0] if ax.get_ylim()[0] != 0 else 250,
        'COVID-19\n第一波', fontsize=9, color='#c0392b', va='bottom')

ax.set_xlabel('年度', fontsize=11)
ax.set_ylabel('消費支出（千円/月）', fontsize=11)
ax.set_title('地域別 消費支出の推移（2012〜2023年）', fontsize=13, fontweight='bold')
ax.legend(loc='upper left', fontsize=9, ncol=2, framealpha=0.8)
ax.set_xticks(range(2012, 2024))
ax.tick_params(axis='x', rotation=45)
ax.grid(axis='y', alpha=0.3, linestyle='--')

# y軸範囲を確定してからテキスト位置を修正
ymin, ymax = ax.get_ylim()
for txt in ax.texts:
    txt.set_y(ymin + (ymax - ymin) * 0.02)

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2020_U1_fig1.png')
plt.savefig(fig1_path, bbox_inches='tight')
plt.close()
print(f"\n図1保存: {fig1_path}")

fig, ax = plt.subplots(figsize=(10, 5.5))

region_job = df_b.groupby(['年度', '地域'])['求人倍率'].mean().reset_index()

for region, color in region_colors.items():
    sub = region_job[region_job['地域'] == region].sort_values('年度')
    if len(sub) > 0:
        ax.plot(sub['年度'], sub['求人倍率'], marker='o', markersize=4,
                label=region, color=color, linewidth=1.8)

# COVID期シェーディング
ax.axvspan(2019.5, 2021.5, alpha=0.15, color='#e74c3c', zorder=0)
ax.axvline(x=2020, color='#e74c3c', linestyle='--', linewidth=1.2, alpha=0.7)
ax.axhline(y=1.0, color='gray', linestyle=':', linewidth=1.2, alpha=0.8, label='倍率=1.0')

ymin2, ymax2 = ax.get_ylim()
ax.text(2020.05, ymin2 + (ymax2 - ymin2) * 0.02,
        'COVID-19\n第一波', fontsize=9, color='#c0392b', va='bottom')

ax.set_xlabel('年度', fontsize=11)
ax.set_ylabel('求人倍率（求人数/求職者数）', fontsize=11)
ax.set_title('地域別 求人倍率の推移（2012〜2023年）', fontsize=13, fontweight='bold')
ax.legend(loc='upper left', fontsize=9, ncol=2, framealpha=0.8)
ax.set_xticks(range(2012, 2024))
ax.tick_params(axis='x', rotation=45)
ax.grid(axis='y', alpha=0.3, linestyle='--')

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2020_U1_fig2.png')
plt.savefig(fig2_path, bbox_inches='tight')
plt.close()
print(f"図2保存: {fig2_path}")

var_labels = {
    'COVID期':        'COVID期ダミー\n（2020・2021年）',
    '求人倍率':       '求人倍率',
    'COVID×旅館密度': 'COVID×旅館密度\n（観光依存×コロナ）',
    '高齢化率':       '高齢化率（%）',
}

fig, ax = plt.subplots(figsize=(8, 5))

vars_order = [v for v in ['COVID期', '求人倍率', 'COVID×旅館密度', '高齢化率'] if v in coef_results]
y_pos = list(range(len(vars_order)))
coefs  = [coef_results[v][0] for v in vars_order]
ses    = [coef_results[v][1] for v in vars_order]
pvals  = [coef_results[v][2] for v in vars_order]
ci95   = [1.96 * se for se in ses]

colors_bar = []
for p in pvals:
    if p < 0.01:
        colors_bar.append('#c0392b')
    elif p < 0.05:
        colors_bar.append('#e67e22')
    else:
        colors_bar.append('#7f8c8d')

bars = ax.barh(y_pos, coefs, xerr=ci95, color=colors_bar, alpha=0.85,
               capsize=5, error_kw={'linewidth': 1.5}, height=0.55)

ax.axvline(x=0, color='black', linewidth=1.2, linestyle='-')

for i, (coef, se, pval) in enumerate(zip(coefs, ses, pvals)):
    sig_str = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.'))
    x_text = coef + 1.96 * se + abs(max(coefs + [ci95[j] for j in range(len(ci95))])) * 0.03
    ax.text(coef, i, f' {coef:+.0f}  {sig_str}',
            va='center', ha='left' if coef >= 0 else 'right',
            fontsize=10, fontweight='bold', color='#2c3e50')

ax.set_yticks(y_pos)
ax.set_yticklabels([var_labels[v] for v in vars_order], fontsize=10)
ax.set_xlabel('推定係数（消費支出への効果：円/月）', fontsize=10)
ax.set_title(f'固定効果パネル推定：主要係数と95%信頼区間\n（{model_type}）',
             fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3, linestyle='--')

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#c0392b', alpha=0.85, label='p < 0.01 (***)'),
    Patch(facecolor='#e67e22', alpha=0.85, label='p < 0.05 (*)'),
    Patch(facecolor='#7f8c8d', alpha=0.85, label='n.s.'),
]
ax.legend(handles=legend_elements, loc='lower right', fontsize=9)

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2020_U1_fig3.png')
plt.savefig(fig3_path, bbox_inches='tight')
plt.close()
print(f"図3保存: {fig3_path}")

vars_needed = ['都道府県', '年度', '消費支出', 'COVID期', '求人倍率', '旅館密度', '高齢化率']
df_panel_all = df_b[vars_needed].dropna().copy()

coef_results = {}
res_summary = None
model_type = 'OLS（フォールバック）'

try:
    from linearmodels.panel import PanelOLS
    # COVID期ダミーは時間固定効果に完全吸収されるため、
    # 「COVID期 × 観光依存度」の交互作用項を用いてEntity FEのみで推定
    df_panel_fe = df_panel_all.copy()
    df_panel_fe['COVID_旅館'] = df_panel_fe['COVID期'] * df_panel_fe['旅館密度']
    df_panel_fe = df_panel_fe.set_index(['都道府県', '年度'])
    mod = PanelOLS.from_formula(
        '消費支出 ~ COVID期 + 求人倍率 + COVID_旅館 + EntityEffects',
        data=df_panel_fe
    )
    res = mod.fit(cov_type='clustered', cluster_entity=True)
    print("=== PanelOLS（Entity FE）推定結果 ===")
    print(res.summary)
    coef_results = {
        'COVID期':       (res.params['COVID期'],     res.std_errors['COVID期'],     res.pvalues['COVID期']),
        '求人倍率':      (res.params['求人倍率'],    res.std_errors['求人倍率'],    res.pvalues['求人倍率']),
        'COVID×旅館密度': (res.params['COVID_旅館'], res.std_errors['COVID_旅館'],  res.pvalues['COVID_旅館']),
    }
    model_type = 'Entity FE パネル推定（PanelOLS）'
    r2 = res.rsquared
    n_obs = res.nobs
    print(f"\nR²={r2:.4f}, N={n_obs}")

except Exception as e:
    print(f'PanelOLS failed: {e}, falling back to OLS')
    # エンティティダミー + 年度ダミーをOLSで実装（2Way FE相当）
    df_ols = df_panel_all.copy()
    df_ols['COVID_旅館'] = df_ols['COVID期'] * df_ols['旅館密度']
    df_ols = pd.get_dummies(df_ols, columns=['都道府県', '年度'], drop_first=True)
    bool_cols = df_ols.select_dtypes(include='bool').columns
    df_ols[bool_cols] = df_ols[bool_cols].astype(int)
    y_ols = df_ols['消費支出']
    X_cols = ['COVID期', '求人倍率', 'COVID_旅館', '高齢化率'] + \
             [c for c in df_ols.columns if c.startswith('都道府県_') or c.startswith('年度_')]
    X_ols = sm.add_constant(df_ols[X_cols])
    res_ols = sm.OLS(y_ols, X_ols).fit(cov_type='HC1')
    print("=== OLS（FEフォールバック）推定結果 ===")
    print(res_ols.summary())
    coef_results = {
        'COVID期':        (res_ols.params['COVID期'],     res_ols.bse['COVID期'],     res_ols.pvalues['COVID期']),
        '求人倍率':       (res_ols.params['求人倍率'],    res_ols.bse['求人倍率'],    res_ols.pvalues['求人倍率']),
        'COVID×旅館密度': (res_ols.params['COVID_旅館'],  res_ols.bse['COVID_旅館'],  res_ols.pvalues['COVID_旅館']),
        '高齢化率':       (res_ols.params['高齢化率'],    res_ols.bse['高齢化率'],    res_ols.pvalues['高齢化率']),
    }
    r2 = res_ols.rsquared
    n_obs = int(res_ols.nobs)
    print(f"\nR²={r2:.4f}, N={n_obs}")

# 推定係数サマリ表示
print("\n=== 主要係数サマリ ===")
for var, (coef, se, pval) in coef_results.items():
    sig = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.'))
    print(f"{var}: coef={coef:+.1f}, SE={se:.1f}, p={pval:.4f} {sig}")

fig, ax = plt.subplots(figsize=(12, 7))

bar_colors = [region_colors.get(r, '#aaaaaa') for r in df_change['地域']]
bars = ax.bar(range(len(df_change)), df_change['変化率'], color=bar_colors, alpha=0.85, edgecolor='white', linewidth=0.5)

ax.axhline(y=0, color='black', linewidth=1.2)
ax.axhline(y=df_change['変化率'].mean(), color='#2c3e50', linewidth=1.5,
           linestyle='--', alpha=0.8, label=f"全国平均: {df_change['変化率'].mean():+.1f}%")

ax.set_xticks(range(len(df_change)))
ax.set_xticklabels(df_change['都道府県'].tolist(), rotation=90, fontsize=8)
ax.set_ylabel('消費支出変化率（%）', fontsize=11)
ax.set_title('2019年→2020年 消費支出変化率（都道府県別）\nCOVID-19 第一波の地域格差', fontsize=13, fontweight='bold')
ax.grid(axis='y', alpha=0.3, linestyle='--')

# 地域凡例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, alpha=0.85, label=r) for r, c in region_colors.items()]
ax.legend(handles=legend_elements, loc='upper left', fontsize=9, ncol=2, framealpha=0.85)
ax.legend(handles=legend_elements + [
    plt.Line2D([0], [0], color='#2c3e50', linestyle='--', linewidth=1.5, label=f"全国平均: {df_change['変化率'].mean():+.1f}%")
], loc='upper left', fontsize=9, ncol=2, framealpha=0.85)

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2020_U1_fig4.png')
plt.savefig(fig4_path, bbox_inches='tight')
plt.close()
print(f"図4保存: {fig4_path}")