02 분석 — 외생 반사실, 선형 SHAP, 유가 진단

모형. 산업 패널은 $y_t = 100\,\Delta\log VA_t$(차분 명세), 총량 모형은 $y_t = 100\log GDP_t$(수준 명세)이다. 두 모형 모두 방정식별 베이지언 VARX이다.

$$y_t = c + \sum_{l=1}^{4} A_l\, y_{t-l} + \sum_{s=0}^{2} \Gamma_s\, x_{t-s} + u_t$$

계수의 사전분포는 미네소타형 정규분포이며, 교차 시차의 사전분산은 산업연관 강도 $W_{ij}$가 클수록 커진다(IO 사전). 사후평균은 닫힌 해로 구한다. 전망은 S1 외생 경로 $x_{2026Q1:2035Q4}$로 1기-앞 재귀 계산한다.

외생 반사실. 외생 $k$의 2035년 순효과는 기준 전망과 $k$를 고정한 재전망의 2035년 총 실질부가가치 차이이다. 차분·로그차분 외생(유가·실질환율·고용)은 0(수준 불변), 수준 외생(세계수요·실질금리·우리경기)은 2025Q4 값으로 고정한다. CAGR 기여는 $\mathrm{CAGR}(V_{2025}, V_{2035}) - \mathrm{CAGR}(V_{2025}, V_{2035} - \text{효과})$이다.

선형 SHAP. 전망 40분기의 $\hat y_{i,t}$를 상수(drift)·내생 시차(dynamics)·외생 6종으로 나눈다. 외생 $k$의 기여는 $\sum_t \sum_{s} \Gamma_{s,ik}\, x_{k,t-s}$, 곧 계수 × 외생 경로(기준점 0)이다. 분해의 합은 산업 $i$의 2025Q4→2035Q4 로그성장(×100)과 같다. 외생이 내생 시차를 거쳐 다음 분기로 전파되는 몫은 dynamics에 들어가므로, 외생 SHAP은 직접 효과이다(5절에서 반사실과 비교한다).

산출물. output/ 아래 model_fit_check.csv, model_totals.csv, exog_effects.csv, exog_paths_summary.csv, exog_persistence.csv, shap_36.csv·shap_62.csv·shap_69.csv, shap_check.csv, oil_rank.csv, exog_top_industries_62.csv, oil_shap_vs_cf_62.csv, battery_parent_check.csv, oil_source_check.csv, oil_source_variants.csv, oil_level_gap.csv, oil_dollar_paths.csv, oil_effect_compare.csv, oil_comovement.csv, oil_coef_controls.csv, oil_benchmark.csv, dubai_sim.csv, exog_effects_postreform.csv·exog_effects_before_after.csv(9절), 정비판(10절) drift_anchor_jeongbi.csv, shap_36/62/69_jeongbi.csv·.json, shap_check_jeongbi.csv, shap_summary_jeongbi.csv, oil_rank_jeongbi.csv, oil_factor_jeongbi.csv, exog_top_industries_62_jeongbi.csv, oil_shap_vs_cf_62_jeongbi.csv, exog_effects_jeongbi.csv, oil_effect_compare_jeongbi.csv, exog_persistence_jeongbi.csv, exog_paths_summary_jeongbi.csv, oil_dollar_paths_jeongbi.csv, (FULL=1) full_effects_by_scenario.csv

import os, json
import config as C            # BLAS 스레드 설정이 적용되도록 numpy보다 먼저 불러온다
import numpy as np, pandas as pd
import exog_diagnosis as F
pd.set_option("display.float_format", lambda v: f"{v:,.3f}")
pd.set_option("display.width", 200)
C.setup()
import time
T0 = time.time()
eh = F.read_exog_hist(C.EXOG_HIST)
es = F.read_exog_scen(C.EXOG_SCEN)
lv = F.read_exog_levels(C.EXOG_LEVELS)
dsge = pd.read_csv(C.DSGE_EXOG, encoding="utf-8-sig")
obs = pd.read_csv(C.BOKDPM_OBS)
wti_m, dub_m = F.read_fred_csv(C.WTI_CSV), F.read_fred_csv(C.DUB_CSV)
ex_s = es[C.SCEN]
out = lambda name: os.path.join(C.OUTPUT, name)   # noqa: E731

Python 3.12.10 · pandas 2.2.3 · numpy 2.4.6
  ONLINE=0 (FRED 유가 재수집)  FULL=0 (확장 점검)  FRED 키 없음(공개 CSV 사용)
  외생 과거(정상형)     있음  data\inputs\exog_hist_v2.csv
  외생 시나리오(정상형)   있음  data\inputs\exog_scen_v2.csv
  외생 수준          있음  data\inputs\exog_scenarios_v2.csv
  시나리오 충격        있음  data\inputs\dsge_exog_paths.csv
  BOKDPM 관측      있음  data\inputs\bokdpm_obs_2000_2025.csv
  FRED WTI       있음  data\inputs\fred_POILWTIUSDM.csv
  FRED 두바이       있음  data\inputs\fred_POILDUBUSDM.csv
  동결 계수          있음  data\inputs\varx_coef_frozen.npz
  패널 총량          있음  data\inputs\levels_agg.csv
  패널 36          있음  data\inputs\levels_sa.csv
  패널 62          있음  data\inputs\levels_62.csv
  패널 69          있음  data\inputs\levels_69.csv

1. 모형 재추정과 동결 계수 대조

네 모형을 전체 표본(2000Q1–2025Q4)으로 다시 추정하고, 원 파이프라인(industry_pipeline/run_forecast.py)이 같은 입력으로 추정한 계수(varx_coef_frozen.npz)와 비교한다. 차이는 부동소수점 연산 순서에서 생기는 수준(1e-12 이하)이어야 한다.

frozen = np.load(C.COEF_FROZEN, allow_pickle=False)
M, rows = {}, []
for name, spec in C.MODELS.items():
    m = F.fit_varx(name, spec["levels"], spec["gdp"], spec["spec"], eh, spec["w"])
    M[name] = m
    B0, i0 = frozen[f"B_{spec['key']}"], list(frozen[f"inds_{spec['key']}"])
    rows.append({"모형": name, "명세": spec["spec"], "산업 수": len(m.industries), "계수 행렬": f"{m.B.shape[0]}×{m.B.shape[1]}",
                 "산업 순서 일치": i0 == m.industries, "동결 계수와 최대 |차이|": float(np.abs(m.B - B0).max())})
fit = pd.DataFrame(rows).set_index("모형")
fit.to_csv(out("model_fit_check.csv"), encoding="utf-8-sig")
assert fit["산업 순서 일치"].all() and (fit["동결 계수와 최대 |차이|"] < 1e-8).all(), "재추정 계수가 동결 계수와 다르다"
fit.style.format({"동결 계수와 최대 |차이|": "{:.1e}"})

	명세	산업 수	계수 행렬	산업 순서 일치	동결 계수와 최대 |차이|
모형
총량	수준	1	23×1	True	0.0e+00
36	차분	36	163×36	True	1.6e-14
62	차분	62	267×62	True	5.6e-14
69	차분	69	295×69	True	7.2e-14

2. 외생 6종 반사실 분해 — 총량 모형과 62 바텀업

같은 외생 고정 규칙을 총량 모형(실질 GDP, 수준 명세)과 62분류 산업 모형(산업 합계, 차분 명세)에 적용한다. 두 모형의 2025년 수준이 다른 것은 총량 모형이 실질 GDP를, 바텀업이 산업 실질부가가치 합계(기초가격)를 전망하기 때문이다.

R = {"총량": F.counterfactual_effects(M["총량"], ex_s), "62": F.counterfactual_effects(M["62"], ex_s)}
eff = F.effects_table(R["총량"], R["62"])
eff.to_csv(out("exog_effects.csv"), encoding="utf-8-sig")
tot = pd.DataFrame([{"모형": k, "2025(조원)": r["v25"], "2035_S1(조원)": r["v35"], "CAGR(%)": F.cagr(r["v25"], r["v35"])}
                    for k, r in R.items()]).set_index("모형")
tot.to_csv(out("model_totals.csv"), encoding="utf-8-sig")
display(tot)
eff.sort_values("총량_효과_조원", ascending=False)

	2025(조원)	2035_S1(조원)	CAGR(%)
모형
총량	2,326.267	2,970.379	2.474
62	2,178.489	3,506.996	4.876

	변환	총량_효과_조원	총량_CAGR기여_%p	바텀업_효과_조원	바텀업_CAGR기여_%p
외생
실질금리	level	162.495	0.575	354.211	1.111
고용_취업자	dlog	62.309	0.217	31.793	0.095
유가	diff	14.596	0.050	65.473	0.197
실질환율	diff	4.797	0.017	9.949	0.030
세계수요	level	1.598	0.006	-12.023	-0.036
우리경기	level	-0.816	-0.003	7.624	0.023

실질금리 효과가 가장 크고 부호가 양(+)이다. 효과의 크기는 계수와 전망 기간 외생 이동폭의 곱이므로, 전망에서 실질금리가 얼마나 움직이는지를 함께 본다. 수준 외생은 2025Q4 값으로 고정되므로, 효과는 2025Q4 이후 경로의 변화분에서 나온다.

s1 = lv[lv["시나리오"] == C.SCEN]
yr = s1["분기"].str[:4].astype(int)
paths = pd.DataFrame({"2025 연평균": s1[yr == 2025][F.EXOG].mean(), "2035 연평균(S1)": s1[yr == 2035][F.EXOG].mean(),
                      "2025Q4(고정값)": s1[s1["분기"] == "2025Q4"][F.EXOG].iloc[0]})
paths["변화(2035−2025)"] = paths["2035 연평균(S1)"] - paths["2025 연평균"]
paths.to_csv(out("exog_paths_summary.csv"), encoding="utf-8-sig")
paths.loc[["실질금리", "세계수요", "우리경기"]]

	2025 연평균	2035 연평균(S1)	2025Q4(고정값)	변화(2035−2025)
실질금리	0.994	2.776	0.939	1.782
세계수요	0.996	-0.001	1.022	-0.997
우리경기	-0.313	-0.073	0.102	0.239

3. 시나리오 폭 — 일시 충격과 영구 수준차

다섯 시나리오의 연평균 외생 수준에서 최댓값과 최솟값의 차이(폭)를 2026년과 2035년에 비교한다. 2035년 폭이 2026년 폭의 5%를 넘고 0.01보다 크면 '영구', 아니면 '일시'로 판정한다(원 스크립트 기준).

per = F.scenario_widths(lv)
per.to_csv(out("exog_persistence.csv"), encoding="utf-8-sig")
per

	폭_2026	폭_2035	판정
외생
유가	25.156	51.750	영구
세계수요	1.790	0.000	일시
실질금리	4.803	0.000	일시
실질환율	11.869	0.004	일시
고용_취업자	531.605	8.630	일시
우리경기	2.943	0.000	일시

4. 선형 SHAP — 36·62·69분류

세 해상도의 산업 모형 각각에 대해 S1 전망을 분해한다. 분해 합계가 실제 전망 로그성장과 같은지(가법성), 공표 산출물(data/published/shap_*.json, 소수 둘째 자리)과 같은지를 확인한다.

SH, rows = {}, []
for r in ("36", "62", "69"):
    m = M[r]
    sh = F.shap_decompose(m, ex_s)
    fc = F.forecast(m, ex_s)
    growth = 100 * np.log(fc.loc["2035Q4"] / np.exp(m.L.loc["2025Q4"] / 100))
    cmp = F.compare_published_shap(sh, C.SHAP_PUBLISHED[r])
    rows.append({"분류": r, "산업 수": len(sh), "가법성 오차(최대)": float((growth - sh["total"]).abs().max()),
                 "공표 대조 값 수": cmp["n_values"], "반올림 후 최대 |차이|": cmp["max_diff_rounded"],
                 "반올림 전 최대 |차이|": cmp["max_diff_raw"], "누락 산업": cmp["missing"]})
    sh.to_csv(out(f"shap_{r}.csv"), encoding="utf-8-sig")
    SH[r] = sh
shap_chk = pd.DataFrame(rows).set_index("분류")
shap_chk.to_csv(out("shap_check.csv"), encoding="utf-8-sig")
assert (shap_chk["가법성 오차(최대)"] < 1e-8).all(), "SHAP 분해 합계가 전망 로그성장과 다르다"
shap_chk.style.format({"가법성 오차(최대)": "{:.1e}", "반올림 후 최대 |차이|": "{:.4f}", "반올림 전 최대 |차이|": "{:.4f}"})

	산업 수	가법성 오차(최대)	공표 대조 값 수	반올림 후 최대 |차이|	반올림 전 최대 |차이|	누락 산업
분류
36	36	9.7e-13	324	0.0000	0.0050	0
62	62	1.1e-12	558	0.0000	0.0050	0
69	69	5.3e-13	621	0.0000	0.0050	0

유가 기여(10년 로그성장, %p)의 상위·하위 5개 산업과 |기여|가 가장 큰 산업을 해상도별로 정리한다.

rows = []
for r, sh in SH.items():
    rk = F.oil_rank(sh)
    for i, (ind, v) in enumerate(rk["top"].items(), 1):
        rows.append({"분류": r, "구분": "상위(+)", "순위": i, "산업": ind, "유가 기여(%p)": v, "|기여| 최대": ind == rk["max_abs"]})
    for i, (ind, v) in enumerate(rk["bottom"].items(), 1):
        rows.append({"분류": r, "구분": "하위(−)", "순위": i, "산업": ind, "유가 기여(%p)": v, "|기여| 최대": ind == rk["max_abs"]})
orank = pd.DataFrame(rows)
orank.to_csv(out("oil_rank.csv"), index=False, encoding="utf-8-sig")
orank.pivot_table(index=["구분", "순위"], columns="분류", values="산업", aggfunc="first")

	분류	36	62	69
구분	순위
상위(+)	1	코크스 및 석유정제품 제조업	14.석탄ㆍ석유제품	14.석탄ㆍ석유제품
	2	기계 및 장비 제조업	25.기타기계ㆍ장비	이차전지
	3	예술, 스포츠 및 여가관련 서비스업	59.문화서비스업	25.기타기계ㆍ장비
	4	운송장비 제조업	27.기타 운송장비	27.기타 운송장비(잔여)
	5	전기장비 제조업	60.스포츠ㆍ오락서비스업	59.문화서비스업
하위(−)	1	주거용 건물 건설업	43.통신업	토목건설
	2	토목건설업	29.전기ㆍ가스ㆍ증기ㆍ공기조절공급업	건물건설건축보수
	3	통신업	33.건설업	43.통신업
	4	건축보수업	46.금융업	29.전기ㆍ가스ㆍ증기ㆍ공기조절공급업(잔여)
	5	금융 및 보험업	48.금융ㆍ보험관련서비스업	46.금융업

rows = []
for k in F.EXOG:
    s = SH["62"][k].sort_values(ascending=False)
    rows.append({"외생": k, "상위1": f"{s.index[0]} {s.iloc[0]:+.1f}", "상위2": f"{s.index[1]} {s.iloc[1]:+.1f}",
                 "하위1": f"{s.index[-1]} {s.iloc[-1]:+.1f}", "하위2": f"{s.index[-2]} {s.iloc[-2]:+.1f}"})
top62 = pd.DataFrame(rows).set_index("외생")
top62.to_csv(out("exog_top_industries_62.csv"), encoding="utf-8-sig")
top62

	상위1	상위2	하위1	하위2
외생
유가	14.석탄ㆍ석유제품 +31.4	25.기타기계ㆍ장비 +11.3	43.통신업 -4.4	29.전기ㆍ가스ㆍ증기ㆍ공기조절공급업 -4.0
세계수요	29.전기ㆍ가스ㆍ증기ㆍ공기조절공급업 +2.0	27.기타 운송장비 +1.9	31.하수ㆍ폐수ㆍ분뇨처리업 -2.6	32.폐기물ㆍ자원재활용서비스업 -2.4
실질금리	23.의료ㆍ정밀ㆍ광학기기ㆍ시계 +53.8	32.폐기물ㆍ자원재활용서비스업 +50.6	7.음료품 -21.2	29.전기ㆍ가스ㆍ증기ㆍ공기조절공급업 -20.0
실질환율	20.1차 금속 +1.1	21.금속가공제품 +0.9	7.음료품 -0.9	6.식료품 -0.9
고용_취업자	51.연구개발업 +5.0	53.과학기술관련전문서비스업 +4.8	49.부동산업 -2.7	43.통신업 -2.0
우리경기	2.임업 +1.3	1.농업 +1.3	50.장비ㆍ용품임대업 -1.0	54.사업지원서비스업 -1.0

5. SHAP(직접 효과)과 반사실(전파 포함)의 차이 — 62분류 유가

선형 SHAP의 외생 기여는 해당 분기의 직접 효과만 담는다. 유가를 고정한 재전망과 기준 전망의 산업별 10년 로그성장 차이는 산업 간 시차 전파까지 포함한 총효과이다. 두 값의 차이가 dynamics로 넘어간 전파분이다.

cf = F.industry_counterfactual(M["62"], ex_s, "유가")
sv = pd.DataFrame({"SHAP_유가(직접)": SH["62"]["유가"], "반사실_유가(총효과)": cf})
sv["전파분"] = sv["반사실_유가(총효과)"] - sv["SHAP_유가(직접)"]
sv.to_csv(out("oil_shap_vs_cf_62.csv"), encoding="utf-8-sig")
print(f"산업 간 상관 {sv.iloc[:, 0].corr(sv.iloc[:, 1]):.3f} · |전파분| 평균 {sv['전파분'].abs().mean():.2f}%p")
sv.sort_values("반사실_유가(총효과)", ascending=False).head(6)

산업 간 상관 0.903 · |전파분| 평균 1.92%p

	SHAP_유가(직접)	반사실_유가(총효과)	전파분
14.석탄ㆍ석유제품	31.359	32.188	0.829
9.섬유ㆍ의복	9.088	11.237	2.148
10.가죽제품	9.264	11.216	1.951
25.기타기계ㆍ장비	11.306	11.185	-0.121
59.문화서비스업	11.150	10.558	-0.592
60.스포츠ㆍ오락서비스업	10.361	10.554	0.193

69분류에서 이차전지의 유가 기여가 높게 나타난다. 이차전지 계열은 전기장비 실질부가가치에 시변 몫(KOSIS 생산지수 기반)을 곱해 만들고 별도 방정식으로 추정하므로, 모산업 잔여 부문과의 성장률 상관과 기여를 함께 확인한다.

L69 = M["69"].L.diff().dropna()
bat = pd.DataFrame({"유가 기여(%p)": SH["69"].loc[["이차전지", "24.전기장비(잔여)"], "유가"],
                    "10년 로그성장(%p)": SH["69"].loc[["이차전지", "24.전기장비(잔여)"], "total"]})
bat["분기 성장률 상관(이차전지와)"] = [1.0, L69["이차전지"].corr(L69["24.전기장비(잔여)"])]
bat.to_csv(out("battery_parent_check.csv"), encoding="utf-8-sig")
bat

	유가 기여(%p)	10년 로그성장(%p)	분기 성장률 상관(이차전지와)
ind
이차전지	15.807	115.900	1.000
24.전기장비(잔여)	10.019	42.602	0.117

6. 유가 진단

6.1 원천 — VARX 유가 과거값은 두바이유이다

01_입수.ipynb 4절에서 확인한 대로, VARX 외생 유가의 과거값은 100·ln(두바이유 분기 평균)과 일치하고 WTI와는 다르다. 같은 모형을 FRED 두바이유·WTI의 100·Δln 분기값으로 각각 다시 추정해, 유가 원천이 공표 수치에 주는 영향을 계산한다. 전망 기간의 유가 경로(S1)는 세 경우 모두 같게 둔다.

s1o = lv[lv["시나리오"] == "S1"].set_index("분기")["유가"].reindex(F.SAMPLE)
wq = F.quarterly_100ln(wti_m, start="2000-01-01"); wq.index = wq.index.astype(str)
dq = F.quarterly_100ln(dub_m, start="2000-01-01"); dq.index = dq.index.astype(str)
gap_w = s1o - wq.reindex(F.SAMPLE)
src_chk = pd.DataFrame([
    {"비교": "BOKDPM 관측 LOIL_us − 100·ln(WTI)", "최대 |차이|": F.loil_observation_check(obs, wti_m)["차이"].abs().max(), "최소": np.nan, "최대": np.nan},
    {"비교": "VARX 유가 − 100·ln(두바이)", "최대 |차이|": (s1o - dq.reindex(F.SAMPLE)).abs().max(), "최소": np.nan, "최대": np.nan},
    {"비교": "VARX 유가 − 100·ln(WTI)", "최대 |차이|": gap_w.abs().max(), "최소": gap_w.min(), "최대": gap_w.max()},
    {"비교": "분기 변화 상관: VARX 유가 vs WTI", "최대 |차이|": np.nan, "최소": np.nan, "최대": s1o.diff().corr(wq.reindex(F.SAMPLE).diff())},
]).set_index("비교")
src_chk.to_csv(out("oil_source_check.csv"), encoding="utf-8-sig")
display(src_chk.style.format("{:.6f}", na_rep="—"))

rows = []
for tag, ex in [("공표 입력(exog_hist)", eh), ("FRED 두바이 100·Δln", F.exog_with_oil(eh, dub_m)), ("FRED WTI 100·Δln", F.exog_with_oil(eh, wti_m))]:
    ma = F.fit_varx("총량", C.MODELS["총량"]["levels"], "총량GDP", "수준", ex)
    mb = F.fit_varx("62", C.MODELS["62"]["levels"], "총량GDP", "차분", ex, C.MODELS["62"]["w"])
    ra, rb = F.counterfactual_effects(ma, ex_s), F.counterfactual_effects(mb, ex_s)
    sh = F.shap_decompose(mb, ex_s)
    rows.append({"유가 과거값": tag, "총량 2035(조원)": ra["v35"], "총량 유가효과(조원)": ra["effects"]["유가"],
                 "총량 유가효과(%)": 100 * ra["effects"]["유가"] / ra["v35"], "62 2035(조원)": rb["v35"],
                 "62 유가효과(조원)": rb["effects"]["유가"], "62 유가효과(%)": 100 * rb["effects"]["유가"] / rb["v35"],
                 "62 실질금리효과(조원)": rb["effects"]["실질금리"], "석유제품 유가 SHAP(%p)": sh.loc["14.석탄ㆍ석유제품", "유가"]})
var = pd.DataFrame(rows).set_index("유가 과거값")
var.to_csv(out("oil_source_variants.csv"), encoding="utf-8-sig")
var

	최대 |차이|	최소	최대
비교
BOKDPM 관측 LOIL_us − 100·ln(WTI)	0.000000	—	—
VARX 유가 − 100·ln(두바이)	0.000050	—	—
VARX 유가 − 100·ln(WTI)	30.309548	-30.309548	19.618954
분기 변화 상관: VARX 유가 vs WTI	—	—	0.962280

	총량 2035(조원)	총량 유가효과(조원)	총량 유가효과(%)	62 2035(조원)	62 유가효과(조원)	62 유가효과(%)	62 실질금리효과(조원)	석유제품 유가 SHAP(%p)
유가 과거값
공표 입력(exog_hist)	2,970.379	14.596	0.491	3,506.996	65.473	1.867	354.211	31.359
FRED 두바이 100·Δln	2,970.379	14.596	0.491	3,506.996	65.473	1.867	354.211	31.359
FRED WTI 100·Δln	2,923.696	5.731	0.196	3,511.477	57.114	1.627	342.823	34.708

FRED 두바이유로 다시 추정하면 공표 수치가 그대로 재현된다(차이는 입력 파일의 소수 넷째 자리 반올림). 문서상 정의인 WTI로 추정하면 총량 모형의 유가 효과가 크게 줄고, 총량 기준 전망 자체도 낮아진다. 유가 원천의 선택은 공표된 '두바이유 교체 시산'(7절)이 시사하는 것보다 결과에 더 큰 영향을 준다.

6.2 시나리오 유가 충격의 누적 — 일시 충격이 영구 수준차가 된다

시나리오 유가 수준은 기준 경로에 충격 성장률 편차(oil_g − S1)의 누적합을 더해 만든다. 충격은 8분기 뒤 0이 되지만 누적합은 그대로 남는다.

gap = F.oil_level_gap_check(lv, dsge)
gap.to_csv(out("oil_level_gap.csv"), encoding="utf-8-sig")
assert gap["최대|수준차−누적합|"].max() < 1e-6, "시나리오 유가 수준차가 충격 누적합과 다르다"
dollars = F.oil_dollar_paths(lv)
dollars.to_csv(out("oil_dollar_paths.csv"), encoding="utf-8-sig")
display(gap)
dollars.loc[[2000, 2008, 2012, 2020, 2025, 2030, 2035]]

	충격_분기수	첫분기_oil_g편차	2035Q4_oil_g편차	2035Q4_수준차(100ln)	최대|수준차−누적합|
시나리오
S1	0	0.000	0.000	0.000	0.000
S2	8	-4.000	0.000	-18.000	0.000
S3	8	-1.500	0.000	-6.750	0.000
S4	0	0.000	0.000	0.000	0.000
S5	8	7.500	0.000	33.750	0.000

시나리오	S1	S2	S3	S4	S5
연도
2000	26.185	26.185	26.185	26.185	26.185
2008	89.135	89.135	89.135	89.135	89.135
2012	109.027	109.027	109.027	109.027	109.027
2020	42.123	42.123	42.123	42.123	42.123
2025	69.227	69.227	69.227	69.227	69.227
2030	99.706	83.282	93.198	99.706	139.732
2035	163.364	136.453	152.701	163.364	228.945

S2·S3·S5의 2035Q4 수준차(−18.0, −6.75, +33.75, 100·ln)가 그대로 2035년 시나리오 폭(51.75)이 된다. 과거 유가는 큰 폭으로 오르내리지만, 기준 경로(S1)는 2025년 이후 단조 상승해 2035년 연평균 $163에 이른다.

6.3 유가 순효과 — 총량 모형과 62 바텀업

rows = []
for name, r in R.items():
    e = r["effects"]["유가"]
    rows.append({"모형": "총량 예측모형" if name == "총량" else "62 바텀업", "2025(조원)": r["v25"], "2035(조원)": r["v35"],
                 "유가 고정 2035(조원)": r["v35"] - e, "유가 효과(조원)": e, "유가 효과(%)": 100 * e / r["v35"],
                 "CAGR(%)": F.cagr(r["v25"], r["v35"]), "유가 고정 CAGR(%)": F.cagr(r["v25"], r["v35"] - e),
                 "CAGR 기여(%p)": F.cagr(r["v25"], r["v35"]) - F.cagr(r["v25"], r["v35"] - e)})
oil_cmp = pd.DataFrame(rows).set_index("모형")
oil_cmp.to_csv(out("oil_effect_compare.csv"), encoding="utf-8-sig")
print(f"바텀업/총량 유가 효과 배율: {oil_cmp.iloc[1]['유가 효과(조원)'] / oil_cmp.iloc[0]['유가 효과(조원)']:.2f}배")
oil_cmp

바텀업/총량 유가 효과 배율: 4.49배

	2025(조원)	2035(조원)	유가 고정 2035(조원)	유가 효과(조원)	유가 효과(%)	CAGR(%)	유가 고정 CAGR(%)	CAGR 기여(%p)
모형
총량 예측모형	2,326.267	2,970.379	2,955.783	14.596	0.491	2.474	2.424	0.050
62 바텀업	2,178.489	3,506.996	3,441.522	65.473	1.867	4.876	4.679	0.197

6.4 양(+)의 유가 계수 — 공행과 인과의 구분

한국은 원유 순수입국이므로 순수한 공급 충격이라면 유가 상승은 비용 요인이어야 한다. 그런데 추정 계수는 양(+)이다. 아래 표는 이 부호가 세계 경기와의 공행에서 비롯되었을 가능성을 점검한다. (1) 유가 변화와 세계수요(미국 산출갭)의 상관은 2000년대에 높고 이후 낮다. (2) 총량 실질 GDP 성장률을 유가 변화에 회귀하면 다른 외생을 통제할 때 계수가 줄지만 양(+)으로 남는다. 통제 변수가 미국 밖 수요와 투기 수요를 담지 못하므로 공행 성분이 일부만 제거된다. 공급·수요 충격을 구조적으로 분리하는 작업(Kilian, 2009)은 후속 과제로 남긴다.

com = F.oil_comovement(eh, C.COMOVE_PERIODS)
com.to_csv(out("oil_comovement.csv"), encoding="utf-8-sig")
display(com)
y = M["총량"].L["총량"].diff()   # 실질 GDP 분기 성장률(100·Δlog)
rows = []
for label, cols in [("유가만", ["유가"]), ("+세계수요", ["유가", "세계수요"]), ("+외생 5종", ["유가"] + F.DUBAI_CTRL)]:
    for a, b in [("2000Q1", "2025Q4"), ("2000Q1", "2008Q4"), ("2009Q1", "2025Q4")]:
        bb, se, n = F.ols_coef(y.loc[a:b], eh.loc[a:b, cols], "유가")
        rows.append({"통제": label, "기간": f"{a}–{b}", "유가 계수": bb, "표준오차": se, "t값": bb / se, "관측치": n})
coef = pd.DataFrame(rows)
coef.to_csv(out("oil_coef_controls.csv"), index=False, encoding="utf-8-sig")
coef.pivot(index="통제", columns="기간", values="유가 계수").loc[["유가만", "+세계수요", "+외생 5종"]]

	분기수	corr(유가, 세계수요)	corr(유가, 우리경기)
기간
2000Q1–2025Q4	104	0.170	0.136
2000Q1–2008Q4	36	0.433	0.244
2009Q1–2025Q4	68	0.059	0.020

기간	2000Q1–2008Q4	2000Q1–2025Q4	2009Q1–2025Q4
통제
유가만	0.039	0.034	0.029
+세계수요	0.036	0.033	0.029
+외생 5종	0.020	0.023	0.020

세계수요(미국 산출갭) 하나만 통제하면 유가 계수는 거의 변하지 않고(전 기간 0.034 → 0.033), 외생 5종을 모두 통제하면 약 3분의 1 줄어든다(0.023). 유가와 미국 산출갭의 상관이 0에 가까운 2009년 이후에도 계수는 양(+)이다. 따라서 미국 산출갭과의 공행만으로는 양(+)의 부호를 설명하지 못하며, 중국 등 미국 밖 수요와 원유시장 고유의 수요 충격을 구분하는 지표가 필요하다. 이는 Kilian(2009) 분해(원유 공급·세계 실물수요·원유시장 고유 수요)를 후속 과제로 두는 근거이다. 이 절의 회귀는 부호의 출처를 가늠하는 보조 분석이며 인과 효과의 추정치가 아니다.

7. 두바이유 교체 시산 — 원 시산 재현

원 시산(oil_dubai_sim.py)은 FRED WTI·두바이유를 분기 평균 100·ln으로 바꾸고, 산업 6개의 ΔlogVA(×100)를 유가 차분과 통제 5종에 OLS로 회귀해 유가 계수를 비교한다. 10년 기여는 유가 계수 × 전망 유가 누적차분(S1, 2025→2035 연평균)이다. 전체 베이지언 VARX가 아닌 외생부 축약 회귀이다.

원 시산은 동결 FRED 파일 전체(1999-10–2026-07)로 상관을 계산했으므로 표본 끝은 2026Q3(7월 한 달)이다. 원 출력의 기간 표기(2000Q1~2025Q4)와 맞춘 계산을 함께 둔다. 산업 회귀는 산업 패널이 2025Q4에서 끝나므로 두 경우 모두 2000Q2–2025Q4이다.

bm = pd.DataFrame({"원 시산 표본(파일 전체)": F.oil_benchmark(wti_m, dub_m),
                   "2000Q1–2025Q4": F.oil_benchmark(wti_m, dub_m, end="2025-12-31")}).T
bm.to_csv(out("oil_benchmark.csv"), encoding="utf-8-sig")
o25, o35, sumdiff = F.oil_sumdiff(lv)
print(f"전망 유가 누적차분(S1): {o35:.1f} − {o25:.1f} = {sumdiff:.2f} (100·ln, ${np.exp(o25/100):.0f} → ${np.exp(o35/100):.0f})")
dub = F.dubai_reduced_ols(C.MODELS["62"]["levels"], eh, wti_m, dub_m, sumdiff, C.DUBAI_INDUSTRIES)
dub.to_csv(out("dubai_sim.csv"), encoding="utf-8-sig")
display(bm[["첫분기", "끝분기", "분기수", "corr", "기울기", "WTI_$", "DUB_$", "스프레드_100ln"]])
dub

전망 유가 누적차분(S1): 509.6 − 423.7 = 85.86 (100·ln, $69 → $163)

	첫분기	끝분기	분기수	corr	기울기	WTI_$	DUB_$	스프레드_100ln
원 시산 표본(파일 전체)	2000Q1	2026Q3	107	0.954	0.925	79.710	77.132	3.287
2000Q1–2025Q4	2000Q1	2025Q4	104	0.961	0.931	59.646	63.595	-6.411

	WTI계수	DUB계수	WTI기여	DUB기여	변화
산업
14.석탄ㆍ석유제품	0.290	0.308	24.879	26.464	1.585
25.기타기계ㆍ장비	0.104	0.109	8.888	9.323	0.435
27.기타 운송장비	0.085	0.086	7.297	7.407	0.110
33.건설업	-0.012	-0.003	-0.988	-0.298	0.690
43.통신업	-0.001	-0.010	-0.078	-0.864	-0.785
46.금융업	0.016	0.011	1.407	0.910	-0.497

9. 정비판 외생으로 본 반사실(2026-10, 참고)

3장 기준선 정비(정상상태 교체, 실질금리 1차 차분, 인력수급 전망 노동공급, 레버 연결) 뒤의 외생 경로(data/inputs/postreform/)로 같은 반사실 분해를 다시 한다. 실질금리는 1차 차분이므로 고정 규칙이 '변화 0(2025Q4 수준 유지)'이 되어 이전 판의 '2025Q4 수준 고정'과 같은 뜻이다. 62 바텀업은 드리프트 재설정 전 모형이다. 이 절은 정비의 효과를 보이는 참고 자료이며 검증 대상이 아니다(검증은 이전 판 공표 수치를 대상으로 한다). 드리프트 재설정 모형으로 다시 계산한 선형 SHAP과 유가 진단은 10절에 있다.

eh_p = F.read_exog_hist(C.POST_EXOG_HIST)
ex_p = F.read_exog_scen(C.POST_EXOG_SCEN)[C.SCEN]
Mp = {k: F.fit_varx(k, C.MODELS[k]["levels"], C.MODELS[k]["gdp"], C.MODELS[k]["spec"], eh_p, C.MODELS[k]["w"]) for k in ("총량", "62")}
Rp = {k: F.counterfactual_effects(Mp[k], ex_p, F.TF_POST) for k in Mp}
eff_p = F.effects_table(Rp["총량"], Rp["62"], F.TF_POST)
eff_p.to_csv(out("exog_effects_postreform.csv"), encoding="utf-8-sig")
cmp = pd.concat({"이전 판": eff[["총량_효과_조원", "바텀업_효과_조원"]], "정비판": eff_p[["총량_효과_조원", "바텀업_효과_조원"]]}, axis=1)
cmp.to_csv(out("exog_effects_before_after.csv"), encoding="utf-8-sig")
cmp.style.format("{:+.1f}")

	이전 판		정비판
	총량_효과_조원	바텀업_효과_조원	총량_효과_조원	바텀업_효과_조원
외생
유가	+14.6	+65.5	+1.1	+4.7
세계수요	+1.6	-12.0	+3.2	-1.3
실질금리	+162.5	+354.2	+0.4	+1.4
실질환율	+4.8	+9.9	-4.5	-13.2
고용_취업자	+62.3	+31.8	-0.9	-0.1
우리경기	-0.8	+7.6	-3.5	+0.6

10. 정비판 선형 SHAP·유가 진단(2026-10)

9절의 정비판 외생(data/inputs/postreform/)으로 총량·36·62·69 모형을 다시 추정하고, 산업 모형(차분 명세)의 상수항을 5장과 같은 드리프트 재설정으로 바꾼 뒤 4–6절의 선형 SHAP, 유가 최민감 산업, SHAP 대 반사실, 유가 순효과를 다시 계산한다. 이 절의 결과는 03의 표 8–11·그림 4와 검증 셀(정비판 기록 수치)에 쓰인다.

• 드리프트 재설정. 표본(2000Q1–2025Q4) 산업 합계의 평균 성장률을 3장 BOKDPM 잠재성장률 정상상태(연 1.8%)로 낮추는 공통 절편 보정이다. $\delta$ = 표본 평균 성장률 − 1.8(연 %p)로 두고 산업 $i$의 상수항을 $(1-\sum_{l,j}A_{l,ij})\,\delta/4$만큼 낮춘다. 외생·동학 계수는 그대로이므로 선형 SHAP에서 달라지는 것은 과거추세(drift)뿐이고, 그 차이는 $-10\,(1-\sum_{l,j}A_{l,ij})\,\delta$ %p이다. 총량 모형(수준 명세)은 재설정하지 않는다.
• 5장과의 관계. 62분류 S1 분해는 5장 재현 패키지 02 10절(output/shap_62.csv, 전망 갱신 카드의 요인 분해)과 같은 계산이다. 5장 기록 사본(data/published/jeongbi/shap_62_s5.csv)과 값 단위로 대조하고, 이웃 폴더에 5장 패키지가 있으면 그 output/shap_62.csv와도 대조한다.
• 대상. 분해와 반사실은 정합화 이전 바텀업 경로를 대상으로 한다. 최종 전망은 같은 연도의 모든 산업에 같은 비례계수를 곱하므로(5장), 통제 총량 맞춤의 효과는 산업 간에 같다.

10.1 정비판 모형과 드리프트 재설정

for rel, h in C.JB_SHA256.items():
    assert F.sha256(os.path.join(C.PUB_DIR, rel)) == h, f"5장 기록 사본 {rel}의 SHA-256이 다르다"
da5 = pd.read_csv(C.S5_DRIFT, index_col=0, encoding="utf-8-sig")
da5.index = da5.index.astype(str)
s1_cagr = lambda mm: F.cagr(F.history_annual(mm).sum(axis=1).loc[2025], F.q_to_annual(F.forecast(mm, ex_p)).sum(axis=1).loc[2035])   # noqa: E731
Mj, Mraw, rows = {}, {}, []
for name, spec in C.MODELS.items():
    m = F.fit_varx(name, spec["levels"], spec["gdp"], spec["spec"], eh_p, spec["w"])
    r = {"모형": name, "명세": spec["spec"], "산업 수": len(m.industries)}
    if spec["spec"] == "차분":
        Mraw[name] = m
        m, d_ann, g_s = F.anchor_drift(m, C.ANCHOR_GROWTH)
        r |= {"표본 합계 평균 성장률(%)": g_s, "δ(연 %p)": d_ann, "S1 연평균 재설정 전(%)": s1_cagr(Mraw[name]), "S1 연평균 재설정 후(%)": s1_cagr(m),
              "5장 δ": da5.loc[name, "δ(연 %p)"], "5장 S1 연평균 재설정 후(%)": da5.loc[name, "바텀업 S1 연평균 — 재설정 후(%)"]}
    else:
        r |= {"S1 연평균 재설정 후(%)": s1_cagr(m)}
    Mj[name] = m
    rows.append(r)
anc = pd.DataFrame(rows).set_index("모형")
anc["5장과 최대 |차이|"] = np.maximum((anc["δ(연 %p)"] - anc["5장 δ"]).abs(), (anc["S1 연평균 재설정 후(%)"] - anc["5장 S1 연평균 재설정 후(%)"]).abs())
anc.to_csv(out("drift_anchor_jeongbi.csv"), encoding="utf-8-sig")
assert (anc["5장과 최대 |차이|"].dropna() < 1e-6).all(), "드리프트 재설정 δ·바텀업 연평균이 5장 기록과 다르다"
anc.style.format("{:.4f}", subset=[c for c in anc.columns if c not in ("명세", "산업 수", "5장과 최대 |차이|")], na_rep="—") \
   .format("{:.1e}", subset=["5장과 최대 |차이|"], na_rep="—")

	명세	산업 수	S1 연평균 재설정 후(%)	표본 합계 평균 성장률(%)	δ(연 %p)	S1 연평균 재설정 전(%)	5장 δ	5장 S1 연평균 재설정 후(%)	5장과 최대 |차이|
모형
총량	수준	1	1.3931	—	—	—	—	—	—
36	차분	36	2.4570	3.3240	1.5240	3.9652	1.5240	2.4570	5.8e-09
62	차분	62	2.5187	3.3240	1.5240	3.9958	1.5240	2.5187	1.1e-08
69	차분	69	2.5250	3.3240	1.5240	4.0021	1.5240	2.5250	1.3e-08

세 분류의 산업 합계는 같은 국민계정 총량이므로 δ는 세 분류에서 같다(연 약 1.52%p). 바텀업 S1 연평균은 5장 표(재설정 후 2.46·2.52·2.53%)와 같다.

10.2 선형 SHAP — 36·62·69분류(정비판)

점검: (가) 성분 합 = 합계, (나) 합계 = 전망 로그성장(2025Q4→2035Q4), (다) 과거추세의 재설정분 = $-10\,(1-\sum A)\,\delta$, (라) 62분류는 5장 분해와 같은 값.

s5_62 = pd.read_csv(C.S5_SHAP62, index_col=0, encoding="utf-8-sig")
SHJ, rows = {}, []
for r in ("36", "62", "69"):
    m, m0 = Mj[r], Mraw[r]
    sh = F.shap_decompose(m, ex_p)
    sh0 = F.shap_decompose(m0, ex_p)
    fc = F.forecast(m, ex_p)
    growth = 100 * np.log(fc.loc["2035Q4"] / np.exp(m.L.loc["2025Q4"] / 100))
    n_ = len(m.industries)
    lagsum = m0.B[:F.P * n_, :].reshape(F.P, n_, n_).sum(axis=(0, 1))
    reset = (sh["drift"] - sh0["drift"]) + 10 * (1 - lagsum) * anc.loc[r, "δ(연 %p)"]
    row = {"분류": r, "산업 수": n_, "(가) 성분 합 − 합계": float((sh.drop(columns="total").sum(axis=1) - sh["total"]).abs().max()),
           "(나) 합계 − 전망 로그성장": float((growth - sh["total"]).abs().max()), "(다) 재설정분 − (−10(1−ΣA)δ)": float(reset.abs().max()),
           "(라) 5장 기록과 최대 |차이|": np.nan, "(라) 5장 패키지 산출물과 최대 |차이|": np.nan}
    if r == "62":
        assert list(s5_62.index) == list(sh.index), "5장 분해와 산업 순서가 다르다"
        row["(라) 5장 기록과 최대 |차이|"] = float((sh[s5_62.columns] - s5_62.values).abs().max().max())
        if os.path.exists(C.S5_SHAP62_LIVE):   # 개발 폴더·통합 번들처럼 5장 패키지가 이웃에 있을 때
            live = pd.read_csv(C.S5_SHAP62_LIVE, index_col=0, encoding="utf-8-sig")
            row["(라) 5장 패키지 산출물과 최대 |차이|"] = float((sh[live.columns] - live.reindex(sh.index).values).abs().max().max())
    rows.append(row)
    sh.to_csv(out(f"shap_{r}_jeongbi.csv"), encoding="utf-8-sig")
    meta = {"version": "2026-10 정비판", "scenario": C.SCEN, "period": "2025Q4–2035Q4", "unit": "10년 로그 성장 기여(100·Δln, %p)",
            "target": "정합화 이전 바텀업(드리프트 재설정 후)", "delta": round(float(anc.loc[r, "δ(연 %p)"]), 3)}
    with open(out(f"shap_{r}_jeongbi.json"), "w", encoding="utf-8") as f:
        json.dump({"exog": F.EXOG, "meta": meta, "rows": [{"ind": i, **{c: round(float(v), 2) for c, v in x.items()}} for i, x in sh.iterrows()]},
                  f, ensure_ascii=False, indent=1)
    SHJ[r] = sh
shap_chk_j = pd.DataFrame(rows).set_index("분류")
shap_chk_j.to_csv(out("shap_check_jeongbi.csv"), encoding="utf-8-sig")
assert (shap_chk_j.iloc[:, 1:4] < 1e-8).all().all(), "정비판 SHAP 점검 (가)–(다)를 통과하지 못했다"
assert shap_chk_j.loc["62", "(라) 5장 기록과 최대 |차이|"] < 1e-8, "정비판 62분류 SHAP이 5장 분해와 다르다"
if not np.isnan(shap_chk_j.loc["62", "(라) 5장 패키지 산출물과 최대 |차이|"]):
    assert shap_chk_j.loc["62", "(라) 5장 패키지 산출물과 최대 |차이|"] < 1e-8, "정비판 62분류 SHAP이 5장 패키지 output/shap_62.csv와 다르다"
else:
    print("5장 패키지가 이웃 폴더에 없어 기록 사본(data/published/jeongbi/)과만 대조하였다.")
shap_chk_j.style.format("{:.1e}", subset=shap_chk_j.columns[1:], na_rep="—")

	산업 수	(가) 성분 합 − 합계	(나) 합계 − 전망 로그성장	(다) 재설정분 − (−10(1−ΣA)δ)	(라) 5장 기록과 최대 |차이|	(라) 5장 패키지 산출물과 최대 |차이|
분류
36	36	0.0e+00	4.6e-13	6.4e-14	—	—
62	62	0.0e+00	3.9e-13	8.2e-14	1.4e-14	1.4e-14
69	69	0.0e+00	1.2e-12	8.7e-14	—	—

요인별 크기를 이전 판(4절)과 비교한다. 평균 절대 기여는 산업별 기여의 절댓값을 산업 평균한 것이며, 외생은 산업마다 여섯 외생 기여의 절댓값을 더한 뒤 평균한다. 이전 판에서 외생 기여의 대부분을 차지하던 실질금리(수준 투입)가 1차 차분 투입에서 사라지고, 과거추세는 재설정만큼 작아진다.

rows = []
for r in ("36", "62", "69"):
    for tag, sh in (("이전 판", SH[r]), ("정비판", SHJ[r])):
        rows.append({"분류": r, "판": tag, "합계 산업 평균": sh["total"].mean(), "과거추세": sh["drift"].abs().mean(), "산업 동학": sh["dynamics"].abs().mean(),
                     "외생(6종 절댓값 합)": sh[F.EXOG].abs().sum(axis=1).mean(), "실질금리 최대 |기여|": sh["실질금리"].abs().max(),
                     "유가 최대 |기여|": sh["유가"].abs().max()})
shap_sum_j = pd.DataFrame(rows).set_index(["분류", "판"])
shap_sum_j.to_csv(out("shap_summary_jeongbi.csv"), encoding="utf-8-sig")
shap_sum_j.style.format("{:.2f}")

		합계 산업 평균	과거추세	산업 동학	외생(6종 절댓값 합)	실질금리 최대 |기여|	유가 최대 |기여|
분류	판
36	이전 판	29.08	24.61	8.27	34.62	56.84	41.06
	정비판	4.14	21.04	5.00	3.88	0.25	8.47
62	이전 판	34.79	29.02	18.61	33.42	53.81	31.36
	정비판	10.37	23.53	14.21	3.30	0.16	6.21
69	이전 판	37.15	30.76	22.48	33.23	54.02	30.83
	정비판	12.57	25.75	15.74	3.39	0.16	6.08

10.3 유가 최민감 산업(정비판)

4절과 같은 방식으로 유가 기여의 상위·하위 5개 산업과 |기여|가 가장 큰 산업을 정리한다. 유가 SHAP은 산업별 유가 계수 합($\sum_s \Gamma_s$)과 전망 기간 유가 누적 변화(S1, 2025Q4→2035Q4의 100·Δln 합)의 곱에 가깝다. 아래 둘째 표는 이 두 요인을 이전 판과 나눠 보인다.

rows = []
for r, sh in SHJ.items():
    rk = F.oil_rank(sh)
    for i, (ind, v) in enumerate(rk["top"].items(), 1):
        rows.append({"분류": r, "구분": "상위(+)", "순위": i, "산업": ind, "유가 기여(%p)": v, "|기여| 최대": ind == rk["max_abs"]})
    for i, (ind, v) in enumerate(rk["bottom"].items(), 1):
        rows.append({"분류": r, "구분": "하위(−)", "순위": i, "산업": ind, "유가 기여(%p)": v, "|기여| 최대": ind == rk["max_abs"]})
orank_j = pd.DataFrame(rows)
orank_j.to_csv(out("oil_rank_jeongbi.csv"), index=False, encoding="utf-8-sig")
display(orank_j.pivot_table(index=["구분", "순위"], columns="분류", values="산업", aggfunc="first"))

cum_prev, cum_jb = float(ex_s["유가"].sum()), float(ex_p["유가"].sum())
rows = []
for r in ("36", "62", "69"):
    top = orank_j[(orank_j["분류"] == r) & orank_j["|기여| 최대"]]["산업"].iloc[0]
    rows.append({"분류": r, "산업": top, "이전 판 유가 기여(%p)": SH[r].loc[top, "유가"], "정비판 유가 기여(%p)": SHJ[r].loc[top, "유가"],
                 "이전 판 유가 계수 합": F.exog_coef_sum(M[r], "유가")[top], "정비판 유가 계수 합": F.exog_coef_sum(Mj[r], "유가")[top],
                 "이전 판 유가 누적(100·ln)": cum_prev, "정비판 유가 누적(100·ln)": cum_jb,
                 "산업 순위 상관(이전 판 대비)": SH[r]["유가"].rank().corr(SHJ[r]["유가"].rank())})
oil_fac_j = pd.DataFrame(rows).set_index("분류")
oil_fac_j.to_csv(out("oil_factor_jeongbi.csv"), encoding="utf-8-sig")
print(f"S1 유가 누적 변화(2025Q4→2035Q4, 100·ln): 이전 판 {cum_prev:.2f} → 정비판 {cum_jb:.2f} "
      f"(10년간 유가 {np.exp(cum_prev / 100):.2f}배 → {np.exp(cum_jb / 100):.2f}배)")
oil_fac_j.style.format({c: "{:.3f}" if "계수" in c else "{:.2f}" for c in oil_fac_j.columns if c != "산업"})

	분류	36	62	69
구분	순위
상위(+)	1	코크스 및 석유정제품 제조업	14.석탄ㆍ석유제품	14.석탄ㆍ석유제품
	2	기계 및 장비 제조업	25.기타기계ㆍ장비	25.기타기계ㆍ장비
	3	예술, 스포츠 및 여가관련 서비스업	27.기타 운송장비	이차전지
	4	운송장비 제조업	26.자동차	27.기타 운송장비(잔여)
	5	섬유 및 가죽제품 제조업	59.문화서비스업	조선업
하위(−)	1	토목건설업	43.통신업	43.통신업
	2	통신업	52.사업관련전문서비스업	52.사업관련전문서비스업
	3	주거용 건물 건설업	53.과학기술관련전문서비스업	토목건설
	4	건축보수업	51.연구개발업	53.과학기술관련전문서비스업
	5	금융 및 보험업	46.금융업	건물건설건축보수

S1 유가 누적 변화(2025Q4→2035Q4, 100·ln): 이전 판 98.05 → 정비판 22.43 (10년간 유가 2.67배 → 1.25배)

	산업	이전 판 유가 기여(%p)	정비판 유가 기여(%p)	이전 판 유가 계수 합	정비판 유가 계수 합	이전 판 유가 누적(100·ln)	정비판 유가 누적(100·ln)	산업 순위 상관(이전 판 대비)
분류
36	코크스 및 석유정제품 제조업	41.06	8.47	0.439	0.440	98.05	22.43	0.98
62	14.석탄ㆍ석유제품	31.36	6.21	0.337	0.339	98.05	22.43	0.97
69	14.석탄ㆍ석유제품	30.83	6.08	0.332	0.334	98.05	22.43	0.98

세 분류 모두 유가에 가장 민감한 산업은 이전 판과 같이 석유정제(석탄·석유제품)이며, 산업별 유가 기여의 순위도 거의 같다(순위 상관 0.97 이상). 유가 기여의 크기는 세 분류에서 약 5분의 1로 줄었다. 석유정제의 유가 계수 합은 정비 전후에 거의 같으므로(소수 셋째 자리 차이), 감소는 계수가 아니라 S1 유가 경로에서 온다. 3장 기준선 정비에서 유가 정상상태 상승률을 연 9.88%에서 2.3%(실질 불변)로 바꾸어 10년간 유가 누적 상승폭이 줄었기 때문이다. 유가 계수의 양(+) 부호(6.4절)는 정비 후에도 그대로이다.

rows = []
for k in F.EXOG:
    s = SHJ["62"][k].sort_values(ascending=False)
    rows.append({"외생": k, "상위1": f"{s.index[0]} {s.iloc[0]:+.2f}", "상위2": f"{s.index[1]} {s.iloc[1]:+.2f}",
                 "하위1": f"{s.index[-1]} {s.iloc[-1]:+.2f}", "하위2": f"{s.index[-2]} {s.iloc[-2]:+.2f}"})
top62_j = pd.DataFrame(rows).set_index("외생")
top62_j.to_csv(out("exog_top_industries_62_jeongbi.csv"), encoding="utf-8-sig")
top62_j

	상위1	상위2	하위1	하위2
외생
유가	14.석탄ㆍ석유제품 +6.21	25.기타기계ㆍ장비 +2.28	43.통신업 -0.90	52.사업관련전문서비스업 -0.83
세계수요	29.전기ㆍ가스ㆍ증기ㆍ공기조절공급업 +1.87	27.기타 운송장비 +1.65	31.하수ㆍ폐수ㆍ분뇨처리업 -2.30	32.폐기물ㆍ자원재활용서비스업 -2.16
실질금리	48.금융ㆍ보험관련서비스업 +0.16	62.수리ㆍ개인서비스업 +0.15	33.건설업 -0.12	6.식료품 -0.11
실질환율	31.하수ㆍ폐수ㆍ분뇨처리업 +1.85	3.어업 +1.75	20.1차 금속 -3.11	54.사업지원서비스업 -2.02
고용_취업자	56.교육서비스업 +0.04	46.금융업 +0.03	2.임업 -0.03	1.농업 -0.03
우리경기	2.임업 +1.72	4.석탄ㆍ원유ㆍ천연가스 +1.68	54.사업지원서비스업 -1.69	50.장비ㆍ용품임대업 -1.61

10.4 SHAP(직접 효과)과 반사실(전파 포함) — 62분류 유가(정비판)

cf_j = F.industry_counterfactual(Mj["62"], ex_p, "유가", F.TF_POST)
sv_j = pd.DataFrame({"SHAP_유가(직접)": SHJ["62"]["유가"], "반사실_유가(총효과)": cf_j})
sv_j["전파분"] = sv_j["반사실_유가(총효과)"] - sv_j["SHAP_유가(직접)"]
sv_j.to_csv(out("oil_shap_vs_cf_62_jeongbi.csv"), encoding="utf-8-sig")
print(f"산업 간 상관 {sv_j.iloc[:, 0].corr(sv_j.iloc[:, 1]):.3f}(이전 판 {sv.iloc[:, 0].corr(sv.iloc[:, 1]):.3f}) · "
      f"|전파분| 평균 {sv_j['전파분'].abs().mean():.2f}%p(이전 판 {sv['전파분'].abs().mean():.2f}%p)")
sv_j.sort_values("반사실_유가(총효과)", ascending=False).head(6)

산업 간 상관 0.898(이전 판 0.903) · |전파분| 평균 0.56%p(이전 판 1.92%p)

	SHAP_유가(직접)	반사실_유가(총효과)	전파분
14.석탄ㆍ석유제품	6.208	7.099	0.891
37.항공운송업	1.622	2.581	0.959
38.창고ㆍ운송보조서비스업	1.682	2.544	0.863
35.육상운송업	1.478	2.537	1.059
36.수상운송업	1.640	2.488	0.849
59.문화서비스업	1.866	2.295	0.429

10.5 유가 순효과 — 총량 모형과 62 바텀업(정비판)

6.3절과 같은 반사실(유가 고정 재전망)이다. 62 바텀업은 드리프트 재설정 모형이며, 9절의 재설정 전 모형 결과를 함께 둔다. 재설정은 상수항만 바꾸므로 산업별 로그 효과는 같고, 조 원 단위 효과는 2035년 수준이 낮아진 만큼 작아진다. 정비판 5장에서 통제 총량은 BOKDPM 실질GDP이고 총량 모형은 시나리오 비율에만 쓰이므로, 총량 모형의 값은 두 모형의 외생 반응을 비교하는 기준이다.

Rj = {"총량": F.counterfactual_effects(Mj["총량"], ex_p, F.TF_POST), "62": F.counterfactual_effects(Mj["62"], ex_p, F.TF_POST)}
eff_j = F.effects_table(Rj["총량"], Rj["62"], F.TF_POST)
eff_j.to_csv(out("exog_effects_jeongbi.csv"), encoding="utf-8-sig")
rows = []
for label, r in (("총량 예측모형", Rj["총량"]), ("62 바텀업(드리프트 재설정)", Rj["62"]), ("62 바텀업(재설정 전, 9절)", Rp["62"])):
    e = r["effects"]["유가"]
    rows.append({"모형": label, "2025(조원)": r["v25"], "2035(조원)": r["v35"], "유가 고정 2035(조원)": r["v35"] - e, "유가 효과(조원)": e,
                 "유가 효과(%)": 100 * e / r["v35"], "CAGR(%)": F.cagr(r["v25"], r["v35"]), "유가 고정 CAGR(%)": F.cagr(r["v25"], r["v35"] - e),
                 "CAGR 기여(%p)": F.cagr(r["v25"], r["v35"]) - F.cagr(r["v25"], r["v35"] - e)})
oil_cmp_j = pd.DataFrame(rows).set_index("모형")
oil_cmp_j.to_csv(out("oil_effect_compare_jeongbi.csv"), encoding="utf-8-sig")
print(f"바텀업(재설정)/총량 유가 효과 배율: {oil_cmp_j.iloc[1]['유가 효과(조원)'] / oil_cmp_j.iloc[0]['유가 효과(조원)']:.2f}배 "
      f"(이전 판 {oil_cmp.iloc[1]['유가 효과(조원)'] / oil_cmp.iloc[0]['유가 효과(조원)']:.2f}배)")
oil_cmp_j

바텀업(재설정)/총량 유가 효과 배율: 3.69배 (이전 판 4.49배)

	2025(조원)	2035(조원)	유가 고정 2035(조원)	유가 효과(조원)	유가 효과(%)	CAGR(%)	유가 고정 CAGR(%)	CAGR 기여(%p)
모형
총량 예측모형	2,326.267	2,671.423	2,670.349	1.074	0.040	1.393	1.389	0.004
62 바텀업(드리프트 재설정)	2,178.489	2,793.745	2,789.779	3.966	0.142	2.519	2.504	0.015
62 바텀업(재설정 전, 9절)	2,178.489	3,223.383	3,218.716	4.667	0.145	3.996	3.981	0.015

유가 순효과는 두 모형 모두 2035년 수준의 0.2% 미만으로 작아졌고(이전 판 총량 0.49%, 바텀업 1.87%), 방향은 여전히 양(+)이다. 바텀업이 총량보다 크게 반응하는 차이(6.3절)는 줄었지만 남아 있다. SHAP(직접 효과)과 반사실(총효과)의 산업 간 상관은 이전 판과 비슷하게 높아, 유가 기여의 산업 순위는 두 방법에서 같은 그림을 준다.

10.6 정비판 외생 경로와 시나리오 폭

3장 정비판 외생 수준값(postreform/exog_levels_2000_2035.csv)에 정비판 변환(실질금리 1차 차분)을 적용해 9·10절의 정상형 입력과 같은지 확인하고, 2절·3절·6.2절과 같은 방식으로 S1 외생 경로, 시나리오 폭, 유가 경로(달러/배럴)를 계산한다. 정비판은 인구 레버를 경제활동인구에 연결하였으므로 취업자 수도 2035년에 시나리오 간 수준 차이가 남는다(인구 추계 가정의 차이이며 충격의 누적이 아니다).

lv_p = F.read_exog_levels(C.POST_EXOG_LEVELS)
st_p = F.to_stationary(lv_p, F.TF_POST)
es_p = F.read_exog_scen(C.POST_EXOG_SCEN)
d_hist_p = (st_p[(st_p["시나리오"] == "S1") & (st_p["분기"] < "2026Q1")].set_index("분기")[F.EXOG] - eh_p).abs().max().max()
d_scen_p = max((g.set_index("분기")[F.EXOG] - es_p[s]).abs().max().max() for s, g in st_p[st_p["분기"] >= "2026Q1"].groupby("시나리오"))
print(f"정비판 수준값 재변환 최대 차이: 과거 {d_hist_p:.1e} · 시나리오 {d_scen_p:.1e}")
assert d_hist_p < 1e-4 and d_scen_p < 1e-4, "정비판 수준값과 정상형 파일의 변환 관계가 맞지 않는다"
per_j = F.scenario_widths(lv_p)
per_j.to_csv(out("exog_persistence_jeongbi.csv"), encoding="utf-8-sig")
s1p = lv_p[lv_p["시나리오"] == C.SCEN]
yrp = s1p["분기"].str[:4].astype(int)
paths_j = pd.DataFrame({"2025 연평균": s1p[yrp == 2025][F.EXOG].mean(), "2035 연평균(S1)": s1p[yrp == 2035][F.EXOG].mean()})
paths_j["변화(2035−2025)"] = paths_j["2035 연평균(S1)"] - paths_j["2025 연평균"]
paths_j["이전 판 2035 연평균(S1)"] = paths["2035 연평균(S1)"]
paths_j.to_csv(out("exog_paths_summary_jeongbi.csv"), encoding="utf-8-sig")
dollars_j = F.oil_dollar_paths(lv_p)
dollars_j.to_csv(out("oil_dollar_paths_jeongbi.csv"), encoding="utf-8-sig")
display(per_j.style.format({"폭_2026": "{:,.2f}", "폭_2035": "{:,.3f}"}))
display(paths_j.style.format("{:,.3f}"))
dollars_j.loc[[2025, 2026, 2030, 2035]].round(1)

정비판 수준값 재변환 최대 차이: 과거 0.0e+00 · 시나리오 0.0e+00

	폭_2026	폭_2035	판정
외생
유가	25.16	51.750	영구
세계수요	1.79	0.000	일시
실질금리	5.63	0.000	일시
실질환율	13.65	0.004	일시
고용_취업자	681.87	943.421	영구
우리경기	3.44	0.000	일시

	2025 연평균	2035 연평균(S1)	변화(2035−2025)	이전 판 2035 연평균(S1)
유가	423.740	436.821	13.081	509.598
세계수요	0.996	-0.002	-0.998	-0.001
실질금리	0.994	1.312	0.317	2.776
실질환율	827.317	833.318	6.002	825.095
고용_취업자	28,771.771	28,778.034	6.263	31,362.543
우리경기	-0.313	-0.089	0.224	-0.073

시나리오	S1	S2	S3	S4	S5
연도
2025	69.200	69.200	69.200	69.200	69.200
2026	63.900	58.500	61.800	63.900	75.200
2030	70.400	58.800	65.800	70.400	98.600
2035	78.900	65.900	73.800	78.900	110.600

정비판에서도 2035년에 시나리오 간 수준 차이가 남는 외생은 유가(충격 누적)와 취업자 수(인구 추계 가정)이다. 기준 시나리오의 두바이유는 유가 정상상태 상승률을 실질 불변(연 2.3%)으로 바꾸어 2035년 연평균이 이전 판보다 크게 낮다. 실질금리·세계수요·국내 경기의 2035년 수준은 시나리오 간에 같다.

8. 확장 점검(FULL=1)

FULL=1이면 공표 수치 검증에 쓰지 않는 확장 점검을 더한다. 다섯 시나리오 각각에서 총량·62 반사실 순효과를 계산하고, S1에서 36·69분류 바텀업 반사실을 계산한다. 기본 실행은 이 절을 건너뛴다.

if C.FULL:
    rows = []
    for s in sorted(es):
        for name in ("총량", "36", "62", "69"):
            if name in ("36", "69") and s != C.SCEN:
                continue
            r = F.counterfactual_effects(M[name], es[s])
            for k in F.EXOG:
                rows.append({"시나리오": s, "모형": name, "외생": k, "2035(조원)": r["v35"], "효과(조원)": r["effects"][k],
                             "CAGR 기여(%p)": F.cagr(r["v25"], r["v35"]) - F.cagr(r["v25"], r["v35"] - r["effects"][k])})
    full = pd.DataFrame(rows)
    full.to_csv(out("full_effects_by_scenario.csv"), index=False, encoding="utf-8-sig")
    display(full.pivot_table(index="외생", columns=["모형", "시나리오"], values="효과(조원)").round(1))
else:
    print("FULL=0: 확장 점검을 건너뛴다.")
print(f"02 분석 소요 {time.time() - T0:.1f}초")

FULL=0: 확장 점검을 건너뛴다.
02 분석 소요 1.4초