# -*- coding: utf-8 -*-
"""보고서1 6장(진단과 보완) — 외생 6종 진단 계산 함수 모음.

노트북 01(입수)·02(분석)·03(결과)이 이 모듈의 함수만 호출한다. 경로·상수는 config.py에 둔다.

원본 스크립트와의 대응(정리하여 옮긴 것)
  모형 추정·전망   industry_pipeline/run_grid.py, run_forecast.py  (차분·수준 명세만, VECM 제외)
  외생 반사실      scratchpad/analyze_all_exog.py, oil_effect_compare.py
  선형 SHAP        scratchpad/shap_decompose.py, shap_oil_all.py
  두바이유 시산    scratchpad/oil_dubai_sim.py

계산 구조
  1. 모형     : 방정식별 베이지언 VARX(정규 사전, IO 가중 미네소타). 산업 패널은 y=100·Δlog VA(차분),
                총량 모형은 y=100·log VA(수준). 내생 4시차, 외생 동시+2시차, 상수.
  2. 반사실   : 외생 하나를 고정한 채 재전망한 2035년 총 실질VA와 기준 전망(S1)의 차이.
                차분·로그차분 외생은 0(수준 불변), 수준 외생은 2025Q4 값으로 고정한다.
  3. 선형 SHAP: 전망 40분기 동안 각 산업 Δ(100·log VA)를 상수(drift)·내생 시차(dynamics)·외생 6종의
                가법 기여로 나눈다. 외생 기여 = 계수 × 외생 경로(기준점 0). 합계는 10년 로그성장과 같다.
  4. 유가     : LOIL_us = 100·ln(FRED POILWTIUSDM 분기 평균). 시나리오 유가 충격(oil_g 편차)은
                누적합(cumsum)으로 수준에 더해져 일시 충격이 영구 수준차로 남는다.
  5. 두바이유 : 외생부 축약 OLS로 WTI·두바이 유가 계수와 10년 기여를 비교한다(원본 시산과 같은 방식).
  6. 정비판   : 2026-10 정비 외생(실질금리 1차 차분)으로 다시 추정하고, 차분 명세 산업 모형의 상수항을
                5장과 같은 드리프트 재설정(anchor_drift, 5장 패키지 varx_engine.anchor_drift와 같은 식)으로 바꾼 뒤
                선형 SHAP·유가 진단을 다시 계산한다(02 10절).
"""
from __future__ import annotations

import hashlib
import io
import json
import urllib.error
import urllib.parse
import urllib.request
from dataclasses import dataclass

import numpy as np
import pandas as pd

EXOG = ["유가", "세계수요", "실질금리", "실질환율", "고용_취업자", "우리경기"]
TF = {"유가": "diff", "세계수요": "level", "실질금리": "level",
      "실질환율": "diff", "고용_취업자": "dlog", "우리경기": "level"}   # VARX 투입 형태(적분차수 정합)
TF_POST = {**TF, "실질금리": "diff"}   # 2026-10 정비판 투입 형태(실질금리 1차 차분) — 정비 후 반사실(02 9절)에만 쓴다
P, Q, LAM = 4, 2, 0.2              # 내생 시차, 외생 시차(동시+Q), 미네소타 전체 수축


# ── 공통 유틸 ──────────────────────────────────────────────────────────
def quarters(start: str, end: str) -> list[str]:
    """'2000Q1'~'2025Q4' 형식의 분기 목록."""
    y, q = int(start[:4]), int(start[5])
    y1, q1 = int(end[:4]), int(end[5])
    out = []
    while (y, q) <= (y1, q1):
        out.append(f"{y}Q{q}")
        q += 1
        if q > 4:
            q, y = 1, y + 1
    return out


SAMPLE = quarters("2000Q1", "2025Q4")   # 추정 표본 104분기
FUT = quarters("2026Q1", "2035Q4")      # 전망 40분기


def sha256(path: str) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(1 << 20), b""):
            h.update(chunk)
    return h.hexdigest()


def cagr(x0: float, x1: float, years: int = 10) -> float:
    """연평균 증가율(%)."""
    return 100.0 * ((x1 / x0) ** (1.0 / years) - 1.0)


def q_to_annual(df_q: pd.DataFrame) -> pd.DataFrame:
    """분기 흐름(십억원) → 연도 합. index 'YYYYQq' → 연도."""
    yr = pd.Series([int(q[:4]) for q in df_q.index], index=df_q.index)
    return df_q.groupby(yr).sum()


# ── 1. 입수 ────────────────────────────────────────────────────────────
def read_exog_hist(path: str) -> pd.DataFrame:
    """외생 6종 정상형 과거값(2000Q1~2025Q4). 행=분기, 열=EXOG."""
    return pd.read_csv(path, encoding="utf-8-sig").set_index("분기")[EXOG]


def read_exog_scen(path: str) -> dict[str, pd.DataFrame]:
    """외생 6종 정상형 시나리오 경로(S1~S5, 2026Q1~2035Q4) → {시나리오: DataFrame}."""
    df = pd.read_csv(path, encoding="utf-8-sig")
    return {s: g.set_index("분기")[EXOG] for s, g in df.groupby("시나리오", sort=True)}


def read_exog_levels(path: str) -> pd.DataFrame:
    """외생 6종 수준값(시나리오, 분기, 2000Q1~2035Q4). 유가는 100·ln(달러/배럴)."""
    return pd.read_csv(path, encoding="utf-8-sig")


def to_stationary(levels: pd.DataFrame, tf: dict | None = None) -> pd.DataFrame:
    """수준 외생 → VARX 투입 정상형(수준·차분·100·Δlog). 원본 build_varx_exog.py와 같은 변환. tf는 투입 형태표(기본 TF)."""
    TF_ = tf or TF
    rows = []
    for s, g in levels.groupby("시나리오", sort=True):
        g = g.reset_index(drop=True)
        out = pd.DataFrame({"시나리오": s, "분기": g["분기"]})
        for c in EXOG:
            x = g[c].astype(float)
            if TF_[c] == "level":
                v = x
            elif TF_[c] == "diff":
                v = x.diff()
            else:
                v = 100 * np.log(x).diff()
            out[c] = v.bfill().round(4)
        rows.append(out)
    return pd.concat(rows, ignore_index=True)


def read_fred_csv(path: str) -> pd.Series:
    """FRED 월별 CSV(observation_date, 값) → Series(index=월초 Timestamp)."""
    d = pd.read_csv(path)
    d.columns = ["date", "value"]
    d["date"] = pd.to_datetime(d["date"])
    d["value"] = pd.to_numeric(d["value"], errors="coerce")
    return d.dropna().set_index("date")["value"]


def fetch_fred(series_id: str, api_key: str | None = None, timeout: int = 30) -> pd.Series:
    """FRED에서 월별 계열을 받는다. 키가 있으면 API, 없으면 공개 CSV(fredgraph)를 쓴다.

    오류가 나도 요청 주소(키 포함)는 메시지에 넣지 않는다.
    """
    try:
        if api_key:
            q = urllib.parse.urlencode({"series_id": series_id, "api_key": api_key, "file_type": "json"})
            with urllib.request.urlopen(f"https://api.stlouisfed.org/fred/series/observations?{q}", timeout=timeout) as r:
                obs = json.loads(r.read().decode("utf-8"))["observations"]
            s = pd.Series({pd.Timestamp(o["date"]): o["value"] for o in obs})
        else:
            with urllib.request.urlopen(f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={series_id}",
                                        timeout=timeout) as r:
                d = pd.read_csv(io.StringIO(r.read().decode("utf-8")))
            s = pd.Series(d.iloc[:, 1].values, index=pd.to_datetime(d.iloc[:, 0]))
    except urllib.error.HTTPError as e:
        raise RuntimeError(f"FRED 응답 오류(HTTP {e.code}): {series_id}") from None
    except Exception as e:  # noqa: BLE001  (주소·키가 섞이지 않도록 유형만 표시)
        raise RuntimeError(f"FRED 연결 실패({type(e).__name__}): {series_id}") from None
    s = pd.to_numeric(s, errors="coerce").dropna().sort_index()
    s.index.name = "date"
    return s.rename("value")


def quarterly_100ln(monthly: pd.Series, start: str = "1999-10-01", end: str | None = None) -> pd.Series:
    """월별 유가 → 분기 평균 → 100·ln. index=PeriodIndex(Q)."""
    m = monthly[monthly.index >= pd.Timestamp(start)]
    if end:
        m = m[m.index <= pd.Timestamp(end)]
    q = m.groupby(m.index.to_period("Q")).mean()
    return 100.0 * np.log(q)


# ── 2. 모형(원본 run_grid.py·run_forecast.py에서 옮김) ─────────────────────
@dataclass
class Model:
    """추정된 VARX 한 개. L=100·log VA(분기×산업), exo=외생 정상형 과거값, B=계수(k×n)."""
    name: str
    spec: str
    L: pd.DataFrame
    exo: pd.DataFrame
    B: np.ndarray
    groups: list

    @property
    def industries(self) -> list[str]:
        return list(self.L.columns)


def load_panel(levels_csv: str, gdp_col: str, exog_hist: pd.DataFrame, w_csv: str | None):
    """산업 실질VA 패널(십억원)을 100·log로, IO 연관행렬을 산업 순서에 맞춰 읽는다."""
    lev = pd.read_csv(levels_csv, index_col=0)
    inds = [c for c in lev.columns if c != gdp_col]
    L = 100.0 * np.log(lev[inds].reindex(SAMPLE))
    exo = exog_hist.reindex(SAMPLE)[EXOG]
    Wio = None
    if w_csv:
        w = pd.read_csv(w_csv, index_col=0)
        common = [c for c in inds if c in w.columns and c in w.index]
        if len(common) == len(inds):
            Wio = w.loc[inds, inds]
        elif len(common) >= 0.5 * len(inds):
            Wio = w.reindex(index=inds, columns=inds).fillna(0.05)   # 일부 미부여 산업은 최소 수축
    return L, exo, Wio


def ar_sigma(Y: pd.DataFrame) -> np.ndarray:
    """산업별 단변량 AR(P) 잔차 표준편차(미네소타 척도)."""
    s = np.zeros(Y.shape[1])
    V = Y.values
    for i in range(Y.shape[1]):
        y = V[:, i]
        yy = y[P:]
        Z = np.column_stack([y[P - l:-l] for l in range(1, P + 1)] + [np.ones(len(yy))])
        b, *_ = np.linalg.lstsq(Z, yy, rcond=None)
        s[i] = (yy - Z @ b).std(ddof=Z.shape[1]) or 1.0
    return np.where(s > 1e-6, s, 1.0)


def design(Y: pd.DataFrame, exo: pd.DataFrame | None, use_x: bool = True):
    """설계행렬: [내생 1..P 시차 | 외생 0..Q 시차 | 상수]."""
    n = Y.shape[1]
    rows = list(range(max(P, Q), len(Y)))
    Ym = Y.values[rows]
    cols, groups = [], []
    for l in range(1, P + 1):
        cols.append(Y.values[[r - l for r in rows]])
        groups += [("endog", l, j) for j in range(n)]
    if use_x and exo is not None:
        for s in range(0, Q + 1):
            cols.append(exo.values[[r - s for r in rows]])
            groups += [("exog", s, j) for j in range(exo.shape[1])]
    cols.append(np.ones((len(rows), 1)))
    groups += [("const", 0, 0)]
    return Ym, np.hstack(cols), groups


def prior_var(groups, sigma: np.ndarray, n: int, Wio: pd.DataFrame | None) -> np.ndarray:
    """IO 가중 미네소타 사전분산. 교차 시차는 산업연관이 강할수록 느슨하다."""
    V = np.ones((len(groups), n))
    Wm = Wio.values if Wio is not None else None
    for c, (kind, l, j) in enumerate(groups):
        for i in range(n):
            if kind == "endog":
                if i == j:
                    base = LAM ** 2 / l ** 2
                else:
                    base = LAM ** 2 / l ** 2 * (sigma[i] ** 2 / sigma[j] ** 2)
                    if Wm is not None:
                        base *= abs(Wm[i, j]) + 0.05
                V[c, i] = base
            elif kind == "exog":
                V[c, i] = LAM ** 2 * 4.0
            else:
                V[c, i] = 1e6
    return V


def fit_eqwise(Ym: np.ndarray, X: np.ndarray, V: np.ndarray, sigma: np.ndarray, delta1: bool) -> np.ndarray:
    """방정식별 정규 사전 사후평균(닫힌 해). 수준 명세는 자기 1시차 사전평균 1."""
    k, n = X.shape[1], Ym.shape[1]
    XtX = X.T @ X
    B = np.zeros((k, n))
    for i in range(n):
        s2 = sigma[i] ** 2
        b0 = np.zeros(k)
        if delta1:
            b0[i] = 1.0
        iV0 = np.diag(1.0 / V[:, i])
        B[:, i] = np.linalg.solve(iV0 + XtX / s2, iV0 @ b0 + (X.T @ Ym[:, i]) / s2)
    return B


def fit_varx(name: str, levels_csv: str, gdp_col: str, spec: str, exog_hist: pd.DataFrame,
             w_csv: str | None = None) -> Model:
    """VARX를 전체 표본(2000Q1~2025Q4)으로 추정한다. spec='차분'(산업 패널) 또는 '수준'(총량)."""
    L, exo, Wio = load_panel(levels_csv, gdp_col, exog_hist, w_csv)
    Y = L if spec == "수준" else L.diff()
    sig = ar_sigma(Y.dropna())
    Ysub = Y.dropna()
    Ym, X, groups = design(Ysub, exo.loc[Ysub.index], True)
    V = prior_var(groups, sig, L.shape[1], Wio)
    B = fit_eqwise(Ym, X, V, sig, delta1=(spec == "수준"))
    return Model(name, spec, L, exo, B, groups)


def forecast(m: Model, exo_future: pd.DataFrame) -> pd.DataFrame:
    """외생 경로로 2026Q1~2035Q4를 1기-앞 재귀 전망한다. 반환: 분기×산업 실질VA(십억원)."""
    efull = pd.concat([m.exo[EXOG], exo_future[EXOG]])
    Ld = m.L.diff()
    Lv = {q: m.L.loc[q].values.astype(float) for q in m.L.index}
    Gv = {q: Ld.loc[q].values.astype(float) for q in m.L.index}
    allq = list(m.L.index) + FUT
    for t in FUT:
        ti = allq.index(t)
        src = Lv if m.spec == "수준" else Gv
        cols = [src[allq[ti - l]] for l in range(1, P + 1)]
        cols += [efull.loc[allq[ti - s]].values.astype(float) for s in range(0, Q + 1)]
        cols.append(np.array([1.0]))
        yhat = np.concatenate([np.atleast_1d(c) for c in cols]) @ m.B
        if m.spec == "수준":
            Lv[t] = yhat
            Gv[t] = yhat - Lv[allq[ti - 1]]
        else:
            Lv[t] = Lv[allq[ti - 1]] + yhat
            Gv[t] = yhat
    return pd.DataFrame({t: np.exp(Lv[t] / 100.0) for t in FUT}, index=m.industries).T


def anchor_drift(m: Model, target: float) -> tuple[Model, float, float]:
    """드리프트 재설정(2026-10 정비, 5장 패키지 varx_engine.anchor_drift와 같은 식)을 적용한 모형을 돌려준다.

    차분 VARX Δy = c + Σ_l A_l Δy_{t−l} + Γx 의 장기 평균은 μ = (I − ΣA)⁻¹(c + Γx̄)이다.
    δ = 표본(2000Q1–2025Q4) 산업 합계 평균 성장률(연율) − target 으로 두고 c ← c − (I − ΣA)·(δ/4)·1 을 쓰면
    μ가 모든 산업에서 δ/4(분기 100·log)씩 낮아진다. 외생·동학 계수와 산업 간 성장 격차는 그대로이다.
    반환: (상수항을 바꾼 모형, δ(연 %p), 표본 산업 합계 평균 성장률(연 %))
    """
    assert m.spec == "차분", "드리프트 재설정은 차분 명세에만 정의한다."
    n = len(m.industries)
    g_sample = float(400.0 * np.log(np.exp(m.L / 100.0).sum(axis=1)).diff().mean())
    d_ann = g_sample - target
    lagsum = m.B[:P * n, :].reshape(P, n, n).sum(axis=(0, 1))       # 방정식 i: Σ_l Σ_j A_l[i, j]
    B2 = m.B.copy()
    B2[-1, :] = m.B[-1, :] - (1.0 - lagsum) * d_ann / 4.0
    return Model(m.name, m.spec, m.L, m.exo, B2, m.groups), d_ann, g_sample


def history_annual(m: Model) -> pd.DataFrame:
    """과거 분기 실질VA(십억원)의 연도 합."""
    return q_to_annual(np.exp(m.L / 100.0))


# ── 3. 외생 반사실 ─────────────────────────────────────────────────────
def freeze_exog(m: Model, ex_s: pd.DataFrame, k: str, tf: dict | None = None) -> pd.DataFrame:
    """외생 k를 고정한 경로: 차분·로그차분 → 0, 수준 → 2025Q4 값. tf는 투입 형태표(기본 TF)."""
    tf = tf or TF
    exf = ex_s.copy()
    exf[k] = 0.0 if tf[k] != "level" else float(m.exo[k].iloc[-1])
    return exf


def counterfactual_effects(m: Model, ex_s: pd.DataFrame, tf: dict | None = None) -> dict:
    """외생별 2035년 총 실질VA 순효과(조원) = 기준 전망 − 해당 외생 고정 전망."""
    tot = lambda exf: q_to_annual(forecast(m, exf)).sum(axis=1).loc[2035] / 1e3   # noqa: E731
    base = tot(ex_s)
    v25 = history_annual(m).sum(axis=1).loc[2025] / 1e3
    eff = {k: base - tot(freeze_exog(m, ex_s, k, tf)) for k in EXOG}
    return {"v25": v25, "v35": base, "effects": eff}


def effects_table(agg: dict, bu: dict, tf: dict | None = None) -> pd.DataFrame:
    """총량·바텀업 반사실 결과 → 외생별 순효과(조원)·CAGR 기여(%p) 표."""
    tf = tf or TF
    rows = []
    for k in EXOG:
        r = {"외생": k, "변환": tf[k]}
        for tag, R in (("총량", agg), ("바텀업", bu)):
            e = R["effects"][k]
            r[f"{tag}_효과_조원"] = e
            r[f"{tag}_CAGR기여_%p"] = cagr(R["v25"], R["v35"]) - cagr(R["v25"], R["v35"] - e)
        rows.append(r)
    return pd.DataFrame(rows).set_index("외생")


def industry_counterfactual(m: Model, ex_s: pd.DataFrame, k: str, tf: dict | None = None) -> pd.Series:
    """외생 k 고정 시 산업별 10년 로그성장 차이(100·Δlog, %p) — 내생 전파를 포함한 총효과."""
    b = forecast(m, ex_s).loc["2035Q4"]
    c = forecast(m, freeze_exog(m, ex_s, k, tf)).loc["2035Q4"]
    return 100.0 * np.log(b / c)


def scenario_widths(levels: pd.DataFrame, years=(2026, 2035)) -> pd.DataFrame:
    """시나리오 폭(연평균 기준 max−min)과 일시/영구 판정(원본 analyze_all_exog.py 기준)."""
    df = levels.copy()
    df["연도"] = df["분기"].astype(str).str[:4].astype(int)
    rows = []
    for k in EXOG:
        w = {}
        for y in years:
            g = df[df["연도"] == y].groupby("시나리오")[k].mean()
            w[y] = g.max() - g.min()
        a, b = w[years[0]], w[years[1]]
        kind = "영구" if abs(b) > 0.05 * abs(a) + 1e-9 and abs(b) > 0.01 else "일시"
        rows.append({"외생": k, f"폭_{years[0]}": a, f"폭_{years[1]}": b, "판정": kind})
    return pd.DataFrame(rows).set_index("외생")


# ── 4. 선형 SHAP ───────────────────────────────────────────────────────
def shap_decompose(m: Model, ex_s: pd.DataFrame) -> pd.DataFrame:
    """전망 40분기 Δ(100·log VA)의 가법 분해: drift(상수)·dynamics(내생 시차)·외생 6종.

    외생 기여는 해당 시점 계수 × 외생값(기준점 0)이다. 외생이 내생 시차를 거쳐 다음 분기로
    전파되는 몫은 dynamics에 들어간다(직접 효과만 외생에 귀속). 합계 = 2025Q4→2035Q4 로그성장(×100).
    """
    assert m.spec == "차분", "선형 SHAP은 차분 명세 산업 모형에 적용한다."
    inds = m.industries
    efull = pd.concat([m.exo[EXOG], ex_s[EXOG]])
    Ld = m.L.diff()
    Gv = {q: Ld.loc[q].values.astype(float) for q in m.L.index}
    allq = list(m.L.index) + FUT
    n = len(inds)
    drift, dyn = np.zeros(n), np.zeros(n)
    ex = {k: np.zeros(n) for k in EXOG}
    for t in FUT:
        ti = allq.index(t)
        blocks = [("endog", Gv[allq[ti - l]]) for l in range(1, P + 1)]
        blocks += [("exog", efull.loc[allq[ti - s]].values.astype(float)) for s in range(0, Q + 1)]
        blocks.append(("const", np.array([1.0])))
        Gv[t] = np.concatenate([b for _, b in blocks]) @ m.B
        off = 0
        for kind, v in blocks:
            w = len(v)
            slab = m.B[off:off + w, :]
            if kind == "const":
                drift += v @ slab
            elif kind == "endog":
                dyn += v @ slab
            else:
                for jx, k in enumerate(EXOG):
                    ex[k] += v[jx] * slab[jx, :]
            off += w
    out = pd.DataFrame({"drift": drift, "dynamics": dyn, **ex}, index=inds)
    out.insert(0, "total", out.sum(axis=1))
    out.index.name = "ind"
    return out


def exog_coef_sum(m: Model, k: str) -> pd.Series:
    """산업별 외생 k 계수의 합 Σ_s Γ_{s,ik}(동시+Q시차). 외생 SHAP ≈ 이 합 × 전망 기간 외생 누적값."""
    jx = EXOG.index(k)
    rows = [c for c, (kind, _, j) in enumerate(m.groups) if kind == "exog" and j == jx]
    return pd.Series(m.B[rows, :].sum(axis=0), index=m.industries)


def compare_published_shap(shap: pd.DataFrame, json_path: str) -> dict:
    """공표 shap_{36,62,69}.json(소수 둘째 자리)과 계산값의 최대 차이."""
    pub = json.load(open(json_path, encoding="utf-8"))
    cols = ["total", "drift", "dynamics"] + EXOG
    pdf = pd.DataFrame(pub["rows"]).set_index("ind")[cols]
    diff = (shap[cols].reindex(pdf.index).round(2) - pdf).abs()
    raw = (shap[cols].reindex(pdf.index) - pdf).abs()
    return {"n_values": int(pdf.size), "n_ind": len(pdf), "max_diff_rounded": float(diff.max().max()),
            "max_diff_raw": float(raw.max().max()), "missing": int(shap.reindex(pdf.index)["total"].isna().sum())}


def oil_rank(shap: pd.DataFrame, k: int = 5) -> dict:
    """유가 기여 상위(+)·하위(−) k개와 |기여| 최대 산업."""
    s = shap["유가"]
    return {"top": s.sort_values(ascending=False).head(k), "bottom": s.sort_values().head(k),
            "max_abs": s.abs().idxmax()}


# ── 5. 유가 특이사항 ───────────────────────────────────────────────────
def loil_observation_check(obs: pd.DataFrame, wti_m: pd.Series) -> pd.DataFrame:
    """BOKDPM 관측 LOIL_us와 100·ln(WTI 분기 평균) 비교."""
    o = obs.set_index(pd.PeriodIndex(obs["quarter"], freq="Q"))["LOIL_us"]
    w = quarterly_100ln(wti_m, start="2000-01-01").reindex(o.index)
    return pd.DataFrame({"LOIL_us_관측": o, "100ln_WTI": w, "차이": o - w})


def smoothed_vs_observed(levels: pd.DataFrame, obs: pd.DataFrame) -> pd.DataFrame:
    """VARX 유가(=BOKDPM 평활 LOIL_us, 시나리오 공통 과거)와 관측 LOIL_us(=100·ln WTI) 비교."""
    s1 = levels[levels["시나리오"] == "S1"].set_index("분기")["유가"]
    o = obs.set_index("quarter")["LOIL_us"]
    d = pd.DataFrame({"평활_LOIL": s1.reindex(o.index), "관측_LOIL": o})
    d["관측−평활"] = d["관측_LOIL"] - d["평활_LOIL"]
    return d


def oil_level_gap_check(levels: pd.DataFrame, dsge_exog: pd.DataFrame, n: int = 40) -> pd.DataFrame:
    """시나리오 유가 수준차(S−S1)가 oil_g 편차의 누적합과 같은지, 충격 지속기간과 2035 수준차."""
    fut = levels[levels["분기"] >= "2026Q1"]
    base = fut[fut["시나리오"] == "S1"]["유가"].to_numpy()
    g1 = dsge_exog[dsge_exog["시나리오"] == "S1"]["oil_g"].to_numpy()[:n]
    rows = []
    for s in sorted(levels["시나리오"].unique()):
        gap = fut[fut["시나리오"] == s]["유가"].to_numpy() - base
        dev = dsge_exog[dsge_exog["시나리오"] == s]["oil_g"].to_numpy()[:n] - g1
        rows.append({"시나리오": s, "충격_분기수": int((np.abs(dev) > 1e-12).sum()),
                     "첫분기_oil_g편차": dev[0], "2035Q4_oil_g편차": dev[-1],
                     "2035Q4_수준차(100ln)": gap[-1], "최대|수준차−누적합|": float(np.abs(gap - np.cumsum(dev)).max())})
    return pd.DataFrame(rows).set_index("시나리오")


def oil_dollar_paths(levels: pd.DataFrame) -> pd.DataFrame:
    """연도별 유가($/배럴) = exp(분기 100·ln 연평균/100). 행=연도, 열=시나리오."""
    df = levels.copy()
    df["연도"] = df["분기"].astype(str).str[:4].astype(int)
    return np.exp(df.groupby(["연도", "시나리오"])["유가"].mean().unstack() / 100.0)


def oil_comovement(exog_hist: pd.DataFrame, periods) -> pd.DataFrame:
    """유가 변화와 세계수요(미국 산출갭)·우리경기의 상관(기간별)."""
    rows = []
    for a, b in periods:
        s = exog_hist.loc[a:b]
        rows.append({"기간": f"{a}–{b}", "분기수": len(s), "corr(유가, 세계수요)": s["유가"].corr(s["세계수요"]),
                     "corr(유가, 우리경기)": s["유가"].corr(s["우리경기"])})
    return pd.DataFrame(rows).set_index("기간")


def ols_coef(y: pd.Series, X: pd.DataFrame, col: str) -> tuple[float, float, int]:
    """상수 포함 OLS에서 col 계수와 표준오차."""
    df = pd.concat([y.rename("y"), X], axis=1).dropna()
    Z = np.column_stack([np.ones(len(df))] + [df[c].values for c in X.columns])
    b, *_ = np.linalg.lstsq(Z, df["y"].values, rcond=None)
    e = df["y"].values - Z @ b
    s2 = e @ e / (len(df) - Z.shape[1])
    se = np.sqrt(np.diag(s2 * np.linalg.inv(Z.T @ Z)))
    j = 1 + list(X.columns).index(col)
    return float(b[j]), float(se[j]), len(df)


# ── 6. 두바이유 교체 시산(원본 oil_dubai_sim.py) ──────────────────────────
DUBAI_CTRL = ["세계수요", "실질금리", "실질환율", "고용_취업자", "우리경기"]


def oil_benchmark(wti_m: pd.Series, dub_m: pd.Series, end: str | None = None) -> dict:
    """WTI·두바이 분기 100·ln의 로그차분 상관·회귀와 마지막 분기 수준(스프레드)."""
    oil = pd.DataFrame({"WTI": quarterly_100ln(wti_m, end=end), "DUB": quarterly_100ln(dub_m, end=end)}).dropna()
    d = oil.diff().dropna()
    slope, icpt = np.polyfit(d["WTI"], d["DUB"], 1)
    last = oil.iloc[-1]
    return {"첫분기": str(d.index[0]), "끝분기": str(d.index[-1]), "분기수": len(d),
            "corr": float(np.corrcoef(d["WTI"], d["DUB"])[0, 1]), "기울기": float(slope), "절편": float(icpt),
            "WTI_100ln": float(last["WTI"]), "DUB_100ln": float(last["DUB"]),
            "WTI_$": float(np.exp(last["WTI"] / 100)), "DUB_$": float(np.exp(last["DUB"] / 100)),
            "스프레드_100ln": float(last["WTI"] - last["DUB"])}


def oil_sumdiff(levels: pd.DataFrame, scen: str = "S1") -> tuple[float, float, float]:
    """전망 유가 누적차분 = 연평균 100·ln(2035) − 연평균 100·ln(2025). (2025, 2035, 차이)."""
    df = levels[levels["시나리오"] == scen].copy()
    yr = df["분기"].astype(str).str[:4].astype(int)
    a, b = df[yr == 2025]["유가"].mean(), df[yr == 2035]["유가"].mean()
    return float(a), float(b), float(b - a)


def dubai_reduced_ols(levels62_csv: str, exog_hist: pd.DataFrame, wti_m: pd.Series, dub_m: pd.Series,
                      sumdiff: float, industries: list[str]) -> pd.DataFrame:
    """산업 ΔlogVA(×100)를 유가 차분 + 통제 5종에 OLS 회귀해 WTI·두바이 유가 계수와 10년 기여를 비교한다.

    두바이 전망 경로는 WTI와 같은 로그성장으로 근사하므로 10년 기여 = 유가 계수 × 누적차분(공통).
    전체 베이지언 VARX가 아닌 외생부 축약 회귀이며, 벤치마크 교체의 민감도를 가늠하는 시산이다.
    """
    lev = pd.read_csv(levels62_csv, index_col=0)
    lev.index = pd.PeriodIndex(lev.index, freq="Q")
    ex = exog_hist.copy()
    ex.index = pd.PeriodIndex(ex.index, freq="Q")
    dW = quarterly_100ln(wti_m).diff().dropna()
    dD = quarterly_100ln(dub_m).diff().dropna()
    rows = []
    for ind in industries:
        y = 100 * np.log(lev[ind]).diff()
        bw = ols_coef(y, pd.concat([dW.rename("oil"), ex[DUBAI_CTRL]], axis=1), "oil")[0]
        bd = ols_coef(y, pd.concat([dD.rename("oil"), ex[DUBAI_CTRL]], axis=1), "oil")[0]
        rows.append({"산업": ind, "WTI계수": bw, "DUB계수": bd, "WTI기여": bw * sumdiff, "DUB기여": bd * sumdiff,
                     "변화": (bd - bw) * sumdiff})
    return pd.DataFrame(rows).set_index("산업")


def exog_with_oil(exog_hist: pd.DataFrame, monthly: pd.Series) -> pd.DataFrame:
    """외생 과거값의 유가 열을 FRED 월별 계열의 분기 100·Δln으로 바꾼다(FULL 확장 점검용)."""
    d = quarterly_100ln(monthly).diff()
    d.index = [str(p) for p in d.index]
    out = exog_hist.copy()
    out["유가"] = d.reindex(out.index).values
    return out