外国人技能実習生の実習地選択における
就労産業・国籍・都道府県別クラスター効果の導出

TABLE4 = [
    # (表示名, 係数, 標準誤差, 有意記号)
    ('x1 就労産業クラスター',          2.299,  0.8178, '**'),
    ('x1^2 就労産業クラスター(2乗)',  -2.495,  0.9774, '*'),
    ('x2 国籍クラスター',              0.3785, 0.1587, '*'),
    ('x3 有効求人倍率',              -0.2047, 0.1625, ''),
    ('x4 製造業事業所数/万人',        -0.1001, 0.0172, '***'),
    ('x5 着工床面積/万人',           -0.0187, 0.0582, ''),
    ('x6 住宅地価格',               -0.0433, 0.0071, '***'),
    ('x3×x4 交差項',                 0.0103, 0.0068, ''),
]
N_OBS, ADJ_R2 = 188, 0.98284
print(f"  {'独立変数':<26}{'係数':>10}{'標準誤差':>10}{'有意':>6}")
print("  " + "-" * 52)
for name, b, se, sig in TABLE4:
    print(f"  {name:<26}{b:>10.4f}{se:>10.4f}{sig:>6}")
print("  " + "-" * 52)
print(f"  サンプルサイズ N = {N_OBS} ／ 自由度調整済み決定係数 = {ADJ_R2}")

# 原論文 表5（2016年を基準=0 とした年度固定効果）
YEARS      = [2016, 2017, 2018, 2019]
TIME_COEF  = [0.0000, 0.0667, 0.1287, 0.2254]
TIME_SE    = [None,   0.0197, 0.0334, 0.0345]
TIME_SIG   = ['***',  '**',   '***',  '***']
print(f"  {'年度':<8}{'時間効果係数':>12}{'標準誤差':>10}{'有意':>6}")
print("  " + "-" * 36)
for y, c, se, sig in zip(YEARS, TIME_COEF, TIME_SE, TIME_SIG):
    se_s = '―' if se is None else f'{se:.4f}'
    print(f"  {y:<8}{c:>12.4f}{se_s:>10}{sig:>6}")
print("  " + "-" * 36)

# 原論文 表6 の生の個体効果 αi（PDFから全47都道府県を転記）。愛知は基準県で αi=0。
ALPHA = {
    '北海道':-1.676,'青森':-1.728,'岩手':-0.9746,'宮城':-1.463,'秋田':-1.176,'山形':-0.5994,
    '福島':-0.8719,'茨城':-0.6030,'栃木':-0.4443,'群馬':-0.1008,'埼玉':-0.8639,'千葉':-1.442,
    '東京':-0.5908,'神奈川':-1.043,'新潟':-0.5530,'富山':0.2044,'石川':0.0059,'福井':0.3584,
    '山梨':-0.7100,'長野':-0.2901,'岐阜':0.5550,'静岡':-0.0541,'愛知':0.0000,'三重':-0.3649,
    '滋賀':-0.6049,'京都':-0.7951,'大阪':-0.1931,'兵庫':-0.9517,'奈良':-1.280,'和歌山':-1.033,
    '鳥取':-1.128,'島根':-0.9412,'岡山':-0.7145,'広島':-0.3725,'山口':-1.352,'徳島':-0.8800,
    '香川':-0.2273,'愛媛':-0.5517,'高知':-1.064,'福岡':-1.443,'佐賀':-1.061,'長崎':-1.432,
    '熊本':-1.212,'大分':-1.232,'宮崎':-1.337,'鹿児島':-1.402,'沖縄':-1.962,
}
# 原論文 表7 が標準化に用いた定数（平均・標準偏差＝報告値）
REP_MEAN, REP_SD = -0.916, 0.526
prefs  = list(ALPHA.keys())
alpha  = np.array([ALPHA[p] for p in prefs])
# 表7 の標準化（原論文の報告パラメータを使用）: CV_i = (αi - 平均) / 標準偏差
alpha_std = (alpha - REP_MEAN) / REP_SD
print(f"  標準化定数（報告値）: 平均 {REP_MEAN} / 標準偏差 {REP_SD} を使用")
print(f"  （表6 αi 実測の平均 {alpha.mean():.3f}・SD {alpha.std():.3f} とは僅差だが、")
print(f"    個々の標準化値は表7の記載値を正しく再現する＝報告値の転記）")
order = np.argsort(-alpha_std)
print("\n  標準化個体効果（表7の再現）上位10:")
for r in order[:10]:
    print(f"    {r+1:>2}位  {prefs[r]:<4} 標準化={alpha_std[r]:+.4f}  (生の αi={alpha[r]:+.4f})")

df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)
# 都道府県行のみ（地域コード = R + 5桁、全国 R00000 を除外）
df = df_raw[df_raw['地域コード'].astype(str).str.match(r'^R\d{5}$', na=False)].copy()
df['年度'] = df['年度'].astype(int)
# 論文の分析期間 2016〜2019 に限定（47都道府県 × 4年 = 188 と一致するはず）
d = df[df['年度'].between(2016, 2019)].copy()

pop   = pd.to_numeric(d['総人口'], errors='coerce')
req   = pd.to_numeric(d['月間有効求人数（一般）'], errors='coerce')
seek  = pd.to_numeric(d['月間有効求職者数（一般）'], errors='coerce')
floor = pd.to_numeric(d['着工建築物床面積'], errors='coerce')
land  = pd.to_numeric(d['標準価格（平均価格）（住宅地）'], errors='coerce')

d['x3_有効求人倍率'] = req / seek                 # = 有効求人数 / 有効求職者数
d['x5_着工床面積対人口'] = floor / pop            # = 着工建築物床面積 / 総人口（m2/人）
d['x6_住宅地価格'] = land / 10000                 # = 標準価格（住宅地）/ 10000（万円/m2）

print(f"  対象: {d['都道府県'].nunique()} 都道府県 × {d['年度'].nunique()} 年度 = {len(d)} 観測（論文 N=188）")
print()
# 原論文 表2 の記述統計（報告値）
REPORTED = {
    'x3_有効求人倍率':   dict(mean=1.3955, sd=0.2578, mn=0.8846, mx=1.9998),
    'x5_着工床面積対人口': dict(mean=1.0362, sd=0.1594, mn=0.6783, mx=1.5569),
    'x6_住宅地価格':     dict(mean=5.0888, sd=5.6117, mn=1.3300, mx=37.4300),
}
LABEL = {'x3_有効求人倍率':'x3 有効求人倍率',
         'x5_着工床面積対人口':'x5 着工床面積/人口',
         'x6_住宅地価格':'x6 住宅地価格(万円/m2)'}
print(f"  {'変数':<20}{'統計量':<6}{'再計算値':>12}{'報告値(表2)':>12}{'一致度':>9}")
print("  " + "-" * 62)
recomp = {}
for col, rep in REPORTED.items():
    s = d[col]
    got = dict(mean=s.mean(), sd=s.std(), mn=s.min(), mx=s.max())
    recomp[col] = got
    for stat, jp in [('mean','平均'), ('sd','標準偏差'), ('mn','最小'), ('mx','最大')]:
        ratio = got[stat] / rep[stat] if rep[stat] else float('nan')
        head = LABEL[col] if stat == 'mean' else ''
        print(f"  {head:<20}{jp:<6}{got[stat]:>12.4f}{rep[stat]:>12.4f}{ratio:>8.1%}")