# -*- coding: utf-8 -*- """산업 실질부가가치 베이지언 VARX 12옵션 비교 — 계산 함수 모음. 노트북 01(입수)·02(분석)·03(결과)이 이 모듈의 함수만 호출한다. 경로·상수는 config.py에 둔다. 출처 industry_pipeline/run_grid.py(2026-08-28 판)의 추정·평가 엔진을 정리한 사본이다. 계산식, 기본값(P=4, Q=2, λ=0.2), 사전분산 산식, 평가 구간은 원본과 같다. 바꾼 점은 다음과 같다. - 경로와 상수를 인수로 받는다(원본은 config_env.py의 절대경로를 읽는다). - RMSE를 반올림하지 않고 반환한다(원본은 소수 셋째 자리로 반올림해 grid_summary_*.csv에 저장했다). - 파일 쓰기를 노트북으로 옮겼다. industry_pipeline/industry_models.py(정규-역위샤트 켤레 사후평균, 다단계 예측)는 grid_summary 산출에 쓰이지 않은 이전 엔진이므로 옮기지 않았다. 분해 함수 decompose_recursive는 scratchpad/shap_decompose.py (통합 번들 industry/code/shap_decompose.py)의 누적 기여 계산을 일반화한 것이다. 모형(산업 i, 분기 t) y_t = Σ_{l=1..P} A_l y_{t-l} [+ α·ect_{t-1}] [+ Σ_{s=0..Q} Γ_s x_{t-s}] + c + e_t 명세 차분 : y = Δ(100·log VA) 수준 : y = 100·log VA, 자기 1차 시차 계수의 사전평균 1(임의보행) VECM : y = Δ(100·log VA), 오차수정항 ect_i = 100·log VA_i − 100·log 총량(1기 시차) 사전 계수별 독립 정규분포, 평균 0(수준의 자기 1차항만 1) 자기 시차 λ²/l² 교차 시차 λ²/l² · σ_i²/σ_j² (미네소타) λ²/l² · σ_i²/σ_j² · (|W_ij| + 0.05) (산업연관표, W=후방연쇄 강도) 오차수정항 λ², 외생 4λ², 상수 10⁶ 추정 방정식별 사후평균 b_i = (V_i⁻¹ + X'X/σ_i²)⁻¹ (V_i⁻¹ b0_i + X'y_i/σ_i²), σ_i = 단변량 AR(P) 잔차 표준편차(고정). 표본추출(MCMC)이 없으므로 결과는 결정적이다. 평가 모든 명세를 수준(100·log VA)의 1기-앞 예측으로 환산해 구간별 RMSE를 계산한다. train : 2000–2021 적합값, valid : train 계수로 2022–2023 예측, test : train+valid 재적합 계수로 2024–2025 예측. """ from __future__ import annotations import os import numpy as np import pandas as pd P_DEFAULT, Q_DEFAULT, LAM_DEFAULT = 4, 2, 0.2 IO_FLOOR = 0.05 # 산업연관 가중 (|W_ij| + 0.05): 연관이 없어도 계수를 0으로 고정하지 않는다 EXOG_SCALE = 4.0 # 외생 계수 사전분산 = 4λ² CONST_VAR = 1e6 # 상수항 사전분산(사실상 무정보) SPECS = ["차분", "수준", "VECM"] PRIORS = ["minnesota", "io"] EXOG_MODES = ["varx", "var"] # ── 0. 표본·분할 ───────────────────────────────────────────────────── 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 def split_tags(qs: list[str], train_end: str, valid_end: str) -> pd.Series: """분기 → 'train'/'valid'/'test'.""" return pd.Series({q: "train" if q <= train_end else ("valid" if q <= valid_end else "test") for q in qs}) def options() -> list[dict]: """12개 옵션(명세 3 × 사전 2 × 외생 2). 순서는 원본 run_grid.py와 같다.""" return [{"옵션": f"{s}_{pr}_{em}", "명세": s, "사전": pr, "외생": em} for s in SPECS for pr in PRIORS for em in EXOG_MODES] # ── 1. 입수: 원자료 읽기 ───────────────────────────────────────────── def read_levels(path: str, gdp_col: str, qs: list[str]) -> tuple[pd.DataFrame, pd.Series]: """수준 패널(십억 원) → (100·log 산업 VA, 100·log 총량). 표본 밖 분기는 버린다.""" lev = pd.read_csv(path, index_col=0, encoding="utf-8-sig") inds = [c for c in lev.columns if c != gdp_col] L = 100.0 * np.log(lev[inds].reindex(qs)) gdp = 100.0 * np.log(lev[gdp_col].reindex(qs)) return L, gdp def read_exog(path: str, cols: list[str], qs: list[str]) -> pd.DataFrame: """외생 6종(정상형) 과거값을 표본 분기에 맞춘다.""" return pd.read_csv(path, encoding="utf-8-sig").set_index("분기").reindex(qs)[cols] def read_io(path: str | None, inds: list[str]) -> pd.DataFrame | None: """산업연관 강도 행렬을 패널 산업 순서로 정렬한다(원본 run_grid.load와 같은 규칙). 모든 산업이 행렬에 있으면 그대로 정렬하고, 절반 이상이면 빠진 칸을 0.05로 채우며, 그보다 적거나 파일이 없으면 None(= 산업연관표 사전이 미네소타와 같아짐)을 돌려준다. """ if path is None or not os.path.exists(path): return None w = pd.read_csv(path, index_col=0, encoding="utf-8-sig") common = [c for c in inds if c in w.columns and c in w.index] if len(common) == len(inds): return w.loc[inds, inds] if len(common) >= 0.5 * len(inds): return w.reindex(index=inds, columns=inds).fillna(IO_FLOOR) return None def normalize_io(W: pd.DataFrame, floor: float = IO_FLOOR) -> pd.DataFrame: """62·69 행렬과 같은 정규화: 비대각 양수 평균을 1로, 하한 0.05, 대각 1(industry_pipeline/build_io_linkage.py).""" A = W.astype(float).values.copy() n = A.shape[0] off = A[~np.eye(n, dtype=bool)] m = off[off > 0].mean() if (off > 0).any() else 1.0 A = np.clip(A / m, floor, None) np.fill_diagonal(A, 1.0) return pd.DataFrame(A, index=W.index, columns=W.columns) def offdiag_summary(W: pd.DataFrame) -> dict: """비대각 원소의 분포 요약(사전 가중의 크기 비교용).""" A = W.values.astype(float) off = A[~np.eye(A.shape[0], dtype=bool)] return {"평균": off.mean(), "중앙값": np.median(off), "최소": off.min(), "최대": off.max(), "평균 가중(|W|+0.05)": (np.abs(off) + IO_FLOOR).mean()} # ── 2. 추정 ────────────────────────────────────────────────────────── def ar_sigma(Y: pd.DataFrame, p: int = P_DEFAULT) -> np.ndarray: """변수별 단변량 AR(p)+상수 OLS 잔차 표준편차(자유도 보정). 0에 가까우면 1로 둔다.""" 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, ect: pd.DataFrame | None = None, p: int = P_DEFAULT, q: int = Q_DEFAULT): """설계행렬. 열 순서: [y_{t-1}…y_{t-p}] [ect_{t-1}] [x_t…x_{t-q}] [1]. groups는 열마다 (종류, 시차, 변수).""" n = Y.shape[1] rows = list(range(max(p, q), len(Y.index))) 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 ect is not None: cols.append(ect.values[[r - 1 for r in rows]]) groups += [("ect", 1, 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: list, sigma: np.ndarray, n: int, prior: str, W: pd.DataFrame | None, lam: float = LAM_DEFAULT) -> np.ndarray: """계수 사전분산 V(k×n). 열 i가 방정식 i(피설명 산업)이다.""" k = len(groups) V = np.ones((k, n)) Wm = W.values if W 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 prior == "io" and Wm is not None: base *= (abs(Wm[i, j]) + IO_FLOOR) # 연관이 강할수록 느슨하게, 약할수록 0 쪽으로 수축 V[c, i] = base elif kind == "ect": V[c, i] = lam ** 2 elif kind == "exog": V[c, i] = (lam ** 2) * EXOG_SCALE else: V[c, i] = CONST_VAR return V def fit_eqwise(Ym: np.ndarray, X: np.ndarray, V: np.ndarray, sigma: np.ndarray, delta1: bool) -> np.ndarray: """방정식별 정규 사전 × 정규 우도(σ_i 고정)의 사후평균 B(k×n).""" 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 # 수준 명세: 자기 1차 시차 사전평균 1 iV0 = np.diag(1.0 / V[:, i]) B[:, i] = np.linalg.solve(iV0 + XtX / s2, iV0 @ b0 + (X.T @ Ym[:, i]) / s2) return B # ── 3. 평가 ────────────────────────────────────────────────────────── def predict_1step(B, spec, L, G, exo, use_x, gdp, seg_idx, p=P_DEFAULT, q=Q_DEFAULT): """구간 seg_idx의 각 분기를 실제 과거값으로 1기 앞 예측하고 수준(100·log VA)으로 환산한다.""" idx = list(L.index) preds, actual = [], [] src = L if spec == "수준" else G for qq in seg_idx: pos = idx.index(qq) if pos < p: continue row = [src.values[pos - l] for l in range(1, p + 1)] if spec == "VECM": row.append(L.values[pos - 1] - (gdp.values[pos - 1] if gdp is not None else 0.0)) if use_x and exo is not None: for s in range(0, q + 1): row.append(exo.values[pos - s]) row.append(np.array([1.0])) xrow = np.concatenate([np.atleast_1d(r) for r in row]) if np.isnan(xrow).any() or np.isnan(L.values[pos]).any(): continue # 차분의 첫 시차가 없는 분기는 건너뛴다 yhat = xrow @ B preds.append(yhat if spec == "수준" else L.values[pos - 1] + yhat) actual.append(L.values[pos]) return np.array(preds), np.array(actual) def rmse(a: np.ndarray, b: np.ndarray) -> float: """산업·분기를 모두 모은 RMSE(단위: 100·log, 근사적으로 % 수준 오차).""" return float(np.sqrt(((a - b) ** 2).mean())) def fit_option(L, gdp, exo, W, qset, spec, prior, exog_mode, sigma, p=P_DEFAULT, q=Q_DEFAULT, lam=LAM_DEFAULT): """한 옵션을 분기 집합 qset으로 적합해 (B, groups)를 돌려준다.""" G = L.diff() Y = L if spec == "수준" else G use_x = exog_mode == "varx" ect = L.subtract(gdp, axis=0) if (spec == "VECM" and gdp is not None) else None Ysub = Y.loc[qset].dropna() ecs = ect.loc[Ysub.index] if ect is not None else None exs = exo.loc[Ysub.index] if use_x else None Ym, X, groups = design(Ysub, exs, use_x, ecs, p, q) V = prior_var(groups, sigma, L.shape[1], prior, W, lam) return fit_eqwise(Ym, X, V, sigma, delta1=(spec == "수준")), groups def run_grid(L, gdp, exo, W, split, p=P_DEFAULT, q=Q_DEFAULT, lam=LAM_DEFAULT) -> pd.DataFrame: """12개 옵션을 추정·평가한다. 반환: 옵션별 RMSE_train·valid·test(반올림 전), valid 오름차순.""" G = L.diff() tr = [x for x in L.index if split[x] == "train"] va = [x for x in L.index if split[x] == "valid"] te = [x for x in L.index if split[x] == "test"] sig_cache = {} rows = [] for o in options(): spec, prior, em = o["명세"], o["사전"], o["외생"] use_x = em == "varx" Y = L if spec == "수준" else G if spec not in sig_cache: sig_cache[spec] = ar_sigma(Y.dropna(), p) # σ_i: 전체 표본(2000–2025)의 AR(p) 잔차(원본과 같음) sig = sig_cache[spec] B, _ = fit_option(L, gdp, exo, W, tr, spec, prior, em, sig, p, q, lam) pin, ain = predict_1step(B, spec, L, G, exo, use_x, gdp, tr[p:], p, q) pv, av = predict_1step(B, spec, L, G, exo, use_x, gdp, va, p, q) B2, _ = fit_option(L, gdp, exo, W, tr + va, spec, prior, em, sig, p, q, lam) pt, at = predict_1step(B2, spec, L, G, exo, use_x, gdp, te, p, q) rows.append({**o, "RMSE_train": rmse(pin, ain), "RMSE_valid": rmse(pv, av), "RMSE_test": rmse(pt, at), "n_valid": int(av.size), "n_test": int(at.size)}) return pd.DataFrame(rows).sort_values("RMSE_valid", kind="stable").reset_index(drop=True) def exog_coef_table(B: np.ndarray, groups: list, inds: list[str], exog_names: list[str], name: str) -> pd.DataFrame: """외생 변수 name의 산업별 계수(동시·1기·2기 시차와 합계). B는 fit_option의 반환값(열 = 방정식).""" j = exog_names.index(name) rows = {f"시차{s}": B[[c for c, g in enumerate(groups) if g == ("exog", s, j)][0], :] for s in range(Q_DEFAULT + 1)} out = pd.DataFrame(rows, index=inds) out["합계"] = out.sum(axis=1) return out def sign_summary(coef: pd.DataFrame, weights: pd.Series) -> dict: """계수표 → 음(−) 비율과 부가가치 가중 평균(동시항·합계).""" w = weights.reindex(coef.index) / weights.reindex(coef.index).sum() return {"산업 수": len(coef), "동시항 음수 비율(%)": 100 * (coef["시차0"] < 0).mean(), "합계 음수 비율(%)": 100 * (coef["합계"] < 0).mean(), "가중평균 동시항": float((w * coef["시차0"]).sum()), "가중평균 합계": float((w * coef["합계"]).sum())} def select(df: pd.DataFrame, decimals: int = 3) -> list[str]: """valid RMSE(소수 decimals 자리 기준)가 가장 작은 옵션 목록. 동률이면 여러 개를 돌려준다.""" r = df["RMSE_valid"].round(decimals) return df.loc[r == r.min(), "옵션"].tolist() def axis_means(grids: dict[str, pd.DataFrame], col: str = "RMSE_valid") -> pd.DataFrame: """명세·사전·외생 축별 평균 RMSE(분류별·전체).""" out = [] for axis in ["명세", "사전", "외생"]: for g, d in grids.items(): for lvl, v in d.groupby(axis)[col].mean().items(): out.append({"축": axis, "수준": lvl, "분류": g, col: v}) return pd.DataFrame(out).pivot_table(index=["축", "수준"], columns="분류", values=col) # ── 4. 기여 분해(SHAP 방식, 선형 VARX) ───────────────────────────────── def decompose_recursive(B, L, exo, horizon, exog_names, p=P_DEFAULT, q=Q_DEFAULT): """차분 VARX의 재귀 예측을 상수(드리프트)·내생 시차(동학)·외생 변수별 기여로 가법 분해한다. 선형 모형에서 SHAP 값은 '계수 × 입력값(기준 0)'이므로 각 분기 예측 Δ(100·log VA)를 열 블록별 곱으로 나눌 수 있다. horizon 분기의 내생 시차는 실적 대신 직전 예측값을 쓴다(재귀). 외생은 주어진 exo 경로를 쓴다. 반환: 산업별 누적 기여(합계 = 누적 예측 로그성장) DataFrame. """ allq = list(L.index) G = L.diff() gv = {t: G.loc[t].values.astype(float) for t in allq} inds = list(L.columns) acc = {k: np.zeros(len(inds)) for k in ["드리프트", "동학", *exog_names]} total = np.zeros(len(inds)) for t in horizon: ti = allq.index(t) blocks = [("endog", gv[allq[ti - l]]) for l in range(1, p + 1)] blocks += [("exog", exo.loc[allq[ti - s]].values.astype(float)) for s in range(0, q + 1)] blocks += [("const", np.array([1.0]))] xrow = np.concatenate([b for _, b in blocks]) yhat = xrow @ B off = 0 for kind, v in blocks: slab = B[off:off + len(v), :] if kind == "const": acc["드리프트"] += v @ slab elif kind == "endog": acc["동학"] += v @ slab else: for jx, name in enumerate(exog_names): acc[name] += v[jx] * slab[jx, :] off += len(v) gv[t] = yhat # 다음 분기의 시차로 예측값을 쓴다 total += yhat out = pd.DataFrame(acc, index=inds) out.insert(0, "누적예측", total) return out