feat: Lifetime PD (50yr) - Belkin & Suchower + Vasicek model

- Belkin & Suchower (1998) credit cycle index (Zt) estimation via WLS
- Vasicek single-factor conditional PD/TM model
- Macro-Zt OLS regression with stepwise variable selection
- 3-scenario (boom/neutral/recession) 50yr PD projection
- Statistical validation suite (ADF, Ljung-Box, R2, ARCH)
- BOK ECOS API integration with fallback data
- Visualization module (7 chart types)
- Detailed theoretical methodology docs/methodology.md

models/__init__.py

@@ -0,0 +1 @@
+# Core models: 신용사이클, Vasicek, 거시연계 모형

models/credit_cycle.py

@@ -0,0 +1,279 @@
+"""
+Belkin & Suchower (1998) 신용사이클 인덱스 Zt 추정 모듈
+
+핵심 방법론:
+  X_i = √ρ · Z + √(1-ρ) · Y_i
+  
+여기서:
+  X_i: 차입자 i의 신용도 변화 (표준정규)
+  Z:   체계적 요인 (credit cycle index, 표준정규)
+  Y_i: 개별적 요인 (표준정규, Z와 독립)
+  ρ:   자산상관계수
+
+TTC 전이행렬의 누적확률 임계값을 Φ⁻¹로 변환한 후,
+관측 연도별 전이행렬과 모델 전이행렬 사이의 WLS를 최소화하여 Zt 추정.
+
+참고문헌:
+- Belkin, B., Suchower, S., & Forest, L.R. (1998). 
+  "A One-Parameter Representation of Credit Risk and Transition Matrices"
+- Basel Committee on Banking Supervision (2005). 
+  "An Explanatory Note on the Basel II IRB Risk Weight Functions"
+"""
+
+import numpy as np
+from scipy.stats import norm
+from scipy.optimize import minimize_scalar, minimize
+from typing import Dict, Tuple, Optional
+import logging
+
+logger = logging.getLogger(__name__)
+
+
+def compute_thresholds(ttc_matrix: np.ndarray) -> np.ndarray:
+    """
+    TTC 전이행렬에서 등급 경계 임계값(thresholds) 산출
+    
+    각 시작등급 i에 대해, 누적 전이확률의 역정규분포로 임계값 산출:
+      d_{i,j} = Φ⁻¹(Σ_{k≤j} p̄_{i,k})
+    
+    Parameters
+    ----------
+    ttc_matrix : np.ndarray
+        N×N TTC 전이행렬 (행 합 = 1)
+    
+    Returns
+    -------
+    np.ndarray
+        N×N 임계값 행렬 (마지막 열은 항상 +∞)
+    """
+    n = ttc_matrix.shape[0]
+    thresholds = np.full((n, n), np.inf)
+    
+    for i in range(n):
+        cum_prob = 0.0
+        for j in range(n - 1):
+            cum_prob += ttc_matrix[i, j]
+            # 누적확률을 [1e-10, 1-1e-10] 범위로 클리핑 (Φ⁻¹ 발산 방지)
+            cum_prob_clipped = np.clip(cum_prob, 1e-10, 1.0 - 1e-10)
+            thresholds[i, j] = norm.ppf(cum_prob_clipped)
+    
+    return thresholds
+
+
+def model_transition_prob(
+    thresholds: np.ndarray,
+    z: float,
+    rho: float,
+    i: int,
+    j: int
+) -> float:
+    """
+    Z 조건부 전이확률 계산
+    
+    p_{ij}(Z) = Φ((d_{i,j} - √ρ·Z) / √(1-ρ)) - Φ((d_{i,j-1} - √ρ·Z) / √(1-ρ))
+    
+    Parameters
+    ----------
+    thresholds : np.ndarray - 임계값 행렬
+    z : float              - 신용사이클 인덱스
+    rho : float            - 자산상관계수
+    i : int                - 시작 등급 인덱스
+    j : int                - 목표 등급 인덱스
+    
+    Returns
+    -------
+    float : 조건부 전이확률
+    """
+    sqrt_rho = np.sqrt(rho)
+    sqrt_1_rho = np.sqrt(1.0 - rho)
+    
+    # 상한 임계값
+    d_upper = thresholds[i, j]
+    upper = norm.cdf((d_upper - sqrt_rho * z) / sqrt_1_rho)
+    
+    # 하한 임계값 (j=0이면 -∞)
+    if j == 0:
+        lower = 0.0
+    else:
+        d_lower = thresholds[i, j - 1]
+        lower = norm.cdf((d_lower - sqrt_rho * z) / sqrt_1_rho)
+    
+    return max(upper - lower, 0.0)
+
+
+def model_transition_matrix(
+    thresholds: np.ndarray,
+    z: float,
+    rho: float
+) -> np.ndarray:
+    """
+    Z 조건부 전체 전이행렬 산출
+    """
+    n = thresholds.shape[0]
+    tm = np.zeros((n, n))
+    
+    for i in range(n - 1):  # D행은 흡수상태
+        for j in range(n):
+            tm[i, j] = model_transition_prob(thresholds, z, rho, i, j)
+        # 행 합 정규화 (수치오차 보정)
+        row_sum = tm[i].sum()
+        if row_sum > 0:
+            tm[i] /= row_sum
+    
+    # D행: 흡수상태
+    tm[-1, -1] = 1.0
+    
+    return tm
+
+
+def zt_objective(
+    z: float,
+    observed_tm: np.ndarray,
+    thresholds: np.ndarray,
+    rho: float,
+    weights: Optional[np.ndarray] = None
+) -> float:
+    """
+    Zt 추정을 위한 WLS 목적함수
+    
+    minimize_Z Σ_{i,j} w_{ij} * (p_{ij}^{obs} - p_{ij}^{model}(Z))²
+    
+    Parameters
+    ----------
+    z : float              - 신용사이클 인덱스 후보값
+    observed_tm : np.ndarray - 관측된 전이행렬
+    thresholds : np.ndarray  - TTC 임계값
+    rho : float              - 자산상관계수
+    weights : np.ndarray     - 가중치 행렬 (기본: 부도열에 높은 가중치)
+    """
+    n = observed_tm.shape[0]
+    
+    if weights is None:
+        # 가중치: 부도열(D)에 10배 가중, 대각에 5배, 나머지 1배
+        weights = np.ones((n, n))
+        weights[:, -1] = 10.0     # 부도 전이확률에 높은 가중
+        for i in range(n):
+            weights[i, i] = 5.0   # 잔류 확률에도 가중
+    
+    wss = 0.0
+    for i in range(n - 1):  # D행 제외
+        for j in range(n):
+            p_obs = observed_tm[i, j]
+            p_model = model_transition_prob(thresholds, z, rho, i, j)
+            wss += weights[i, j] * (p_obs - p_model) ** 2
+    
+    return wss
+
+
+def estimate_zt(
+    observed_tm: np.ndarray,
+    thresholds: np.ndarray,
+    rho: float,
+    z_bounds: Tuple[float, float] = (-4.0, 4.0)
+) -> float:
+    """
+    단일 연도의 Zt 추정
+    
+    scipy.optimize.minimize_scalar로 WLS 목적함수 최소화
+    
+    Parameters
+    ----------
+    observed_tm : np.ndarray  - 해당 연도 관측 전이행렬
+    thresholds : np.ndarray   - TTC 임계값
+    rho : float               - 자산상관계수
+    z_bounds : tuple           - Z 탐색 범위
+    
+    Returns
+    -------
+    float : 추정된 Zt 값
+    """
+    result = minimize_scalar(
+        zt_objective,
+        bounds=z_bounds,
+        method="bounded",
+        args=(observed_tm, thresholds, rho)
+    )
+    
+    return result.x
+
+
+def estimate_zt_series(
+    transition_matrices: Dict[int, np.ndarray],
+    ttc_matrix: np.ndarray,
+    rho: float = 0.20
+) -> Dict[int, float]:
+    """
+    전체 기간에 대한 Zt 시계열 추정
+    
+    Parameters
+    ----------
+    transition_matrices : Dict[int, np.ndarray]
+        연도별 관측 전이행렬
+    ttc_matrix : np.ndarray
+        TTC 전이행렬
+    rho : float
+        자산상관계수
+    
+    Returns
+    -------
+    Dict[int, float]
+        {연도: Zt값} 딕셔너리
+    """
+    logger.info("TTC 전이행렬에서 임계값 산출 중...")
+    thresholds = compute_thresholds(ttc_matrix)
+    
+    zt_series = {}
+    years = sorted(transition_matrices.keys())
+    
+    logger.info(f"Zt 시계열 추정 중 ({years[0]}-{years[-1]}, rho={rho})...")
+    
+    for year in years:
+        observed_tm = transition_matrices[year]
+        z_hat = estimate_zt(observed_tm, thresholds, rho)
+        zt_series[year] = z_hat
+        logger.debug(f"  {year}: Zt = {z_hat:+.4f}")
+    
+    logger.info(f"Zt 추정 완료. 범위: [{min(zt_series.values()):.3f}, {max(zt_series.values()):.3f}]")
+    
+    return zt_series
+
+
+def estimate_rho_and_zt(
+    transition_matrices: Dict[int, np.ndarray],
+    ttc_matrix: np.ndarray,
+    rho_bounds: Tuple[float, float] = (0.05, 0.50)
+) -> Tuple[float, Dict[int, float]]:
+    """
+    자산상관계수 ρ와 Zt 시계열 동시 추정 (NLS)
+    
+    총 목적함수 = Σ_t Σ_{i,j} w_{ij} * (p_{ij,t}^{obs} - p_{ij,t}^{model}(Z_t(ρ), ρ))²
+    
+    외부 루프: ρ 탐색
+    내부 루프: 각 연도별 Zt 추정 (ρ 고정)
+    
+    Returns
+    -------
+    Tuple[float, Dict[int, float]]
+        (최적 ρ, Zt 시계열)
+    """
+    years = sorted(transition_matrices.keys())
+    
+    def total_objective(rho):
+        thresholds = compute_thresholds(ttc_matrix)
+        total_wss = 0.0
+        for year in years:
+            observed_tm = transition_matrices[year]
+            z_hat = estimate_zt(observed_tm, thresholds, rho)
+            total_wss += zt_objective(z_hat, observed_tm, thresholds, rho)
+        return total_wss
+    
+    logger.info(f"ρ 동시 추정 중 (범위: {rho_bounds})...")
+    result = minimize_scalar(total_objective, bounds=rho_bounds, method="bounded")
+    optimal_rho = result.x
+    
+    logger.info(f"최적 ρ = {optimal_rho:.4f}")
+    
+    # 최적 ρ로 Zt 재추정
+    zt_series = estimate_zt_series(transition_matrices, ttc_matrix, optimal_rho)
+    
+    return optimal_rho, zt_series

models/macro_model.py

@@ -0,0 +1,307 @@
+"""
+거시경제 변수 ↔ Zt 연계 통계모형
+
+Zt(신용사이클 인덱스)를 거시경제변수로 설명하는 회귀모형을 구축하고,
+미래 거시 시나리오에 따른 Zt 전망을 생성합니다.
+
+모형:
+  Z_t = β₀ + β₁·GDP_growth + β₂·Unemployment + β₃·Base_Rate 
+      + β₄·CD_Rate + β₅·CPI_growth + β₆·Leading_Index + ε_t
+
+방법론 참고:
+- IMF (2021). "IFRS 9 and CECL Compatible Estimation for Top-Down Solvency Stress Testing"
+- ECB (2019). "Scenario Design for IFRS 9 Expected Credit Loss Estimation"
+- Fed (2022). "Dodd-Frank Act Stress Test Methodology"
+"""
+
+import numpy as np
+import pandas as pd
+import statsmodels.api as sm
+from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox
+from statsmodels.stats.stattools import durbin_watson
+from statsmodels.stats.outliers_influence import variance_inflation_factor
+from scipy import stats
+from typing import Dict, List, Optional, Tuple
+import logging
+import warnings
+
+logger = logging.getLogger(__name__)
+warnings.filterwarnings("ignore", category=FutureWarning)
+
+
+class MacroZtModel:
+    """
+    거시경제변수 → Zt 회귀모형
+    
+    Features:
+    - OLS 다중회귀
+    - 변수 선택 (Stepwise AIC/BIC)
+    - 잔차 진단 (ADF, Ljung-Box, Breusch-Pagan, DW)
+    - VIF 다중공선성 체크
+    - 시나리오별 Zt 예측
+    """
+    
+    def __init__(self):
+        self.model = None
+        self.result = None
+        self.selected_vars = None
+        self.scaler_params = {}  # 정규화 파라미터
+    
+    def fit(
+        self,
+        zt_series: pd.Series,
+        macro_data: pd.DataFrame,
+        method: str = "stepwise_aic",
+        standardize: bool = True
+    ) -> "MacroZtModel":
+        """
+        Zt ~ 거시변수 회귀모형 적합
+        
+        Parameters
+        ----------
+        zt_series : pd.Series
+            index=연도, values=Zt 추정값
+        macro_data : pd.DataFrame
+            index=연도, columns=거시변수
+        method : str
+            변수 선택 방법:
+            - "all": 모든 변수 사용
+            - "stepwise_aic": Forward stepwise (AIC 기준)
+            - "stepwise_bic": Forward stepwise (BIC 기준)
+        standardize : bool
+            거시변수 표준화 여부
+        
+        Returns
+        -------
+        self
+        """
+        # 인덱스 정렬 및 교집합
+        common_years = sorted(set(zt_series.index) & set(macro_data.index))
+        if len(common_years) < 5:
+            raise ValueError(f"공통 데이터 포인트가 부족합니다: {len(common_years)}개")
+        
+        y = zt_series.loc[common_years].values.astype(float)
+        X = macro_data.loc[common_years].copy()
+        
+        # 결측치 처리
+        X = X.ffill().bfill().dropna(axis=1)
+        
+        # 표준화
+        if standardize:
+            for col in X.columns:
+                mean = X[col].mean()
+                std = X[col].std()
+                if std > 0:
+                    self.scaler_params[col] = {"mean": mean, "std": std}
+                    X[col] = (X[col] - mean) / std
+                else:
+                    X = X.drop(columns=[col])
+        
+        # 변수 선택
+        if method == "all":
+            self.selected_vars = list(X.columns)
+        elif method.startswith("stepwise"):
+            criterion = "aic" if "aic" in method else "bic"
+            self.selected_vars = self._stepwise_selection(y, X, criterion)
+        else:
+            self.selected_vars = list(X.columns)
+        
+        if not self.selected_vars:
+            logger.warning("변수 선택 결과 선택된 변수가 없습니다. 전체 변수 사용.")
+            self.selected_vars = list(X.columns)
+        
+        # 최종 모형 적합
+        X_selected = sm.add_constant(X[self.selected_vars].values)
+        self.model = sm.OLS(y, X_selected)
+        self.result = self.model.fit()
+        
+        logger.info(f"회귀모형 적합 완료: 선택변수 = {self.selected_vars}")
+        logger.info(f"  R² = {self.result.rsquared:.4f}, "
+                    f"Adj.R² = {self.result.rsquared_adj:.4f}, "
+                    f"AIC = {self.result.aic:.2f}")
+        
+        return self
+    
+    def _stepwise_selection(
+        self,
+        y: np.ndarray,
+        X: pd.DataFrame,
+        criterion: str = "aic"
+    ) -> List[str]:
+        """Forward Stepwise 변수 선택"""
+        remaining = list(X.columns)
+        selected = []
+        current_score = np.inf
+        
+        while remaining:
+            scores = {}
+            for var in remaining:
+                trial_vars = selected + [var]
+                X_trial = sm.add_constant(X[trial_vars].values)
+                try:
+                    model = sm.OLS(y, X_trial).fit()
+                    score = model.aic if criterion == "aic" else model.bic
+                    scores[var] = score
+                except Exception:
+                    continue
+            
+            if not scores:
+                break
+            
+            best_var = min(scores, key=scores.get)
+            best_score = scores[best_var]
+            
+            if best_score < current_score:
+                selected.append(best_var)
+                remaining.remove(best_var)
+                current_score = best_score
+                logger.debug(f"  + {best_var} ({criterion.upper()} = {best_score:.2f})")
+            else:
+                break
+        
+        return selected
+    
+    def predict(self, macro_scenario: pd.DataFrame) -> np.ndarray:
+        """
+        거시 시나리오로 Zt 예측
+        
+        Parameters
+        ----------
+        macro_scenario : pd.DataFrame
+            columns에 selected_vars가 포함되어야 함
+        
+        Returns
+        -------
+        np.ndarray : Zt 예측값 배열
+        """
+        if self.result is None:
+            raise ValueError("모형이 적합되지 않았습니다. fit()을 먼저 실행하세요.")
+        
+        X = macro_scenario[self.selected_vars].copy()
+        
+        # 학습 데이터와 동일한 표준화 적용
+        for col in X.columns:
+            if col in self.scaler_params:
+                mean = self.scaler_params[col]["mean"]
+                std = self.scaler_params[col]["std"]
+                X[col] = (X[col] - mean) / std
+        
+        X_const = sm.add_constant(X.values, has_constant="add")
+        return self.result.predict(X_const)
+    
+    def diagnostics(self) -> Dict[str, any]:
+        """
+        회귀 모형 진단 결과 반환
+        
+        Returns
+        -------
+        dict with keys:
+            - r_squared, adj_r_squared
+            - f_stat, f_pvalue
+            - aic, bic
+            - durbin_watson
+            - ljung_box (p-value)
+            - breusch_pagan (p-value)
+            - vif (각 변수별)
+            - coefficients (DataFrame)
+        """
+        if self.result is None:
+            return {}
+        
+        diag = {
+            "r_squared": self.result.rsquared,
+            "adj_r_squared": self.result.rsquared_adj,
+            "f_stat": self.result.fvalue,
+            "f_pvalue": self.result.f_pvalue,
+            "aic": self.result.aic,
+            "bic": self.result.bic,
+            "n_obs": int(self.result.nobs),
+            "selected_vars": self.selected_vars,
+        }
+        
+        # Durbin-Watson
+        residuals = self.result.resid
+        diag["durbin_watson"] = durbin_watson(residuals)
+        
+        # Ljung-Box (자기상관 검정)
+        try:
+            lb_result = acorr_ljungbox(residuals, lags=[5], return_df=True)
+            diag["ljung_box_stat"] = lb_result["lb_stat"].values[0]
+            diag["ljung_box_pvalue"] = lb_result["lb_pvalue"].values[0]
+        except Exception:
+            diag["ljung_box_pvalue"] = np.nan
+        
+        # Breusch-Pagan (이분산 검정)
+        try:
+            bp_stat, bp_pvalue, _, _ = het_breuschpagan(
+                residuals, self.result.model.exog
+            )
+            diag["breusch_pagan_stat"] = bp_stat
+            diag["breusch_pagan_pvalue"] = bp_pvalue
+        except Exception:
+            diag["breusch_pagan_pvalue"] = np.nan
+        
+        # VIF (다중공선성)
+        try:
+            X = self.result.model.exog
+            vif_values = {}
+            var_names = ["const"] + self.selected_vars
+            for i in range(X.shape[1]):
+                vif_values[var_names[i]] = variance_inflation_factor(X, i)
+            diag["vif"] = vif_values
+        except Exception:
+            diag["vif"] = {}
+        
+        # 계수 요약
+        coef_df = pd.DataFrame({
+            "변수": ["const"] + self.selected_vars,
+            "계수": self.result.params,
+            "표준오차": self.result.bse,
+            "t값": self.result.tvalues,
+            "p값": self.result.pvalues,
+        })
+        diag["coefficients"] = coef_df
+        
+        return diag
+    
+    def summary(self) -> str:
+        """모형 요약 출력"""
+        if self.result is None:
+            return "모형이 적합되지 않았습니다."
+        return str(self.result.summary())
+    
+    def residual_series(self) -> np.ndarray:
+        """잔차 시계열 반환"""
+        if self.result is None:
+            return np.array([])
+        return self.result.resid
+
+
+def build_macro_zt_model(
+    zt_dict: Dict[int, float],
+    macro_df: pd.DataFrame,
+    method: str = "stepwise_aic"
+) -> MacroZtModel:
+    """
+    편의 함수: Zt 딕셔너리 + 거시 DataFrame → 회귀모형 구축
+    
+    Parameters
+    ----------
+    zt_dict : Dict[int, float]
+        {연도: Zt값}
+    macro_df : pd.DataFrame
+        index=연도, columns=거시변수
+    method : str
+        변수 선택 방법
+    
+    Returns
+    -------
+    MacroZtModel : 적합된 모형
+    """
+    zt_series = pd.Series(zt_dict, name="Zt")
+    zt_series.index.name = "YEAR"
+    
+    model = MacroZtModel()
+    model.fit(zt_series, macro_df, method=method)
+    
+    return model

models/vasicek.py

@@ -0,0 +1,218 @@
+"""
+Vasicek 단일팩터 모델 기반 조건부 PD 및 전이행렬 모듈
+
+핵심 공식:
+  PD_PIT(Z) = Φ( (Φ⁻¹(PD_TTC) - √ρ · Z) / √(1-ρ) )
+
+이 모듈은 Belkin & Suchower의 임계값 방식 대신, 
+Vasicek 공식을 직접 적용하는 간편 버전도 제공합니다.
+
+참고문헌:
+- Vasicek, O. (2002). "The Distribution of Loan Portfolio Value"
+- Basel Committee (2005). "An Explanatory Note on the Basel II IRB Risk Weight Functions"
+- Merton, R.C. (1974). "On the Pricing of Corporate Debt"
+"""
+
+import numpy as np
+from scipy.stats import norm
+from typing import Optional
+import logging
+
+logger = logging.getLogger(__name__)
+
+
+def conditional_pd(pd_ttc: float, z: float, rho: float) -> float:
+    """
+    Vasicek 공식으로 PIT PD 계산
+    
+    PD_PIT(Z) = Φ( (Φ⁻¹(PD_TTC) - √ρ · Z) / √(1-ρ) )
+    
+    Parameters
+    ----------
+    pd_ttc : float  - TTC (Through-the-Cycle) 부도확률
+    z : float       - 체계적 요인 (Z > 0: 호황, Z < 0: 불황)
+    rho : float     - 자산상관계수 (0 < ρ < 1)
+    
+    Returns
+    -------
+    float : PIT (Point-in-Time) 부도확률
+    
+    Examples
+    --------
+    >>> conditional_pd(0.02, 0, 0.20)   # Z=0이면 PD_PIT = PD_TTC
+    0.02
+    >>> conditional_pd(0.02, -2, 0.20)  # 불황시 PD 상승
+    0.1016...
+    >>> conditional_pd(0.02, 2, 0.20)   # 호황시 PD 하락
+    0.0024...
+    """
+    if pd_ttc <= 0:
+        return 0.0
+    if pd_ttc >= 1:
+        return 1.0
+    
+    sqrt_rho = np.sqrt(rho)
+    sqrt_1_rho = np.sqrt(1.0 - rho)
+    
+    numerator = norm.ppf(pd_ttc) - sqrt_rho * z
+    pd_pit = norm.cdf(numerator / sqrt_1_rho)
+    
+    return float(np.clip(pd_pit, 0.0, 1.0))
+
+
+def conditional_pd_array(pd_ttc_array: np.ndarray, z: float, rho: float) -> np.ndarray:
+    """
+    벡터화된 Vasicek 공식 (등급별 TTC PD 배열 → PIT PD 배열)
+    """
+    pd_ttc_clipped = np.clip(pd_ttc_array, 1e-10, 1.0 - 1e-10)
+    
+    sqrt_rho = np.sqrt(rho)
+    sqrt_1_rho = np.sqrt(1.0 - rho)
+    
+    numerator = norm.ppf(pd_ttc_clipped) - sqrt_rho * z
+    pd_pit = norm.cdf(numerator / sqrt_1_rho)
+    
+    return np.clip(pd_pit, 0.0, 1.0)
+
+
+def conditional_transition_matrix(
+    ttc_tm: np.ndarray,
+    z: float,
+    rho: float
+) -> np.ndarray:
+    """
+    임계값 기반 Z-조건부 전이행렬 산출
+    
+    TTC 전이행렬로부터 누적확률 임계값을 산출하고,
+    Z 값에 따라 조건부 전이확률을 계산합니다.
+    
+    Parameters
+    ----------
+    ttc_tm : np.ndarray  - N×N TTC 전이행렬
+    z : float            - 체계적 요인
+    rho : float          - 자산상관계수
+    
+    Returns
+    -------
+    np.ndarray : N×N 조건부 전이행렬
+    """
+    n = ttc_tm.shape[0]
+    sqrt_rho = np.sqrt(rho)
+    sqrt_1_rho = np.sqrt(1.0 - rho)
+    
+    # 임계값 산출 (누적확률 → Φ⁻¹)
+    thresholds = np.full((n, n), np.inf)
+    for i in range(n):
+        cum_prob = 0.0
+        for j in range(n - 1):
+            cum_prob += ttc_tm[i, j]
+            cum_prob_clipped = np.clip(cum_prob, 1e-10, 1.0 - 1e-10)
+            thresholds[i, j] = norm.ppf(cum_prob_clipped)
+    
+    # 조건부 전이행렬 계산
+    cond_tm = np.zeros((n, n))
+    
+    for i in range(n - 1):
+        for j in range(n):
+            d_upper = thresholds[i, j]
+            upper = norm.cdf((d_upper - sqrt_rho * z) / sqrt_1_rho)
+            
+            if j == 0:
+                lower = 0.0
+            else:
+                d_lower = thresholds[i, j - 1]
+                lower = norm.cdf((d_lower - sqrt_rho * z) / sqrt_1_rho)
+            
+            cond_tm[i, j] = max(upper - lower, 0.0)
+        
+        # 행 합 정규화
+        row_sum = cond_tm[i].sum()
+        if row_sum > 0:
+            cond_tm[i] /= row_sum
+    
+    # D행: 흡수상태
+    cond_tm[-1, -1] = 1.0
+    
+    return cond_tm
+
+
+def multi_period_pd(
+    annual_tm: np.ndarray,
+    horizon: int,
+    initial_grade_idx: Optional[int] = None
+) -> np.ndarray:
+    """
+    전이행렬 거듭제곱으로 다기간 누적/한계 PD 계산
+    
+    Parameters
+    ----------
+    annual_tm : np.ndarray   - 1년 전이행렬
+    horizon : int             - 예측 기간 (년)
+    initial_grade_idx : int   - 특정 등급만 계산 (None이면 전체)
+    
+    Returns
+    -------
+    np.ndarray
+        shape (horizon, N-1): 연도별 각 등급의 누적 PD
+        또는 shape (horizon,): 특정 등급의 누적 PD
+    """
+    n = annual_tm.shape[0]
+    cumulative_tm = np.eye(n)
+    
+    cumulative_pds = []
+    for t in range(1, horizon + 1):
+        cumulative_tm = cumulative_tm @ annual_tm
+        # 부도열(마지막 열)이 누적 PD
+        if initial_grade_idx is not None:
+            cumulative_pds.append(cumulative_tm[initial_grade_idx, -1])
+        else:
+            cumulative_pds.append(cumulative_tm[:-1, -1].copy())
+    
+    return np.array(cumulative_pds)
+
+
+def marginal_pd_from_cumulative(cumulative_pds: np.ndarray) -> np.ndarray:
+    """
+    누적 PD에서 한계 PD(Marginal PD) 계산
+    
+    Marginal PD(t) = Cumulative PD(t) - Cumulative PD(t-1)
+    """
+    if cumulative_pds.ndim == 1:
+        marginal = np.diff(cumulative_pds, prepend=0.0)
+    else:
+        first_row = np.zeros((1, cumulative_pds.shape[1]))
+        marginal = np.diff(cumulative_pds, axis=0, prepend=first_row)
+    
+    return np.maximum(marginal, 0.0)
+
+
+def survival_probability(cumulative_pds: np.ndarray) -> np.ndarray:
+    """생존확률 = 1 - 누적 PD"""
+    return 1.0 - cumulative_pds
+
+
+def annualized_pd(cumulative_pd: float, horizon: int) -> float:
+    """
+    누적 PD를 연환산 PD로 변환
+    
+    AnnualizedPD = 1 - (1 - CumulativePD)^(1/horizon)
+    """
+    if cumulative_pd >= 1.0:
+        return 1.0
+    return 1.0 - (1.0 - cumulative_pd) ** (1.0 / horizon)
+
+
+def worst_case_pd(pd_ttc: float, rho: float, confidence: float = 0.999) -> float:
+    """
+    Basel II IRB 방식 Worst-Case PD (99.9% 신뢰수준)
+    
+    WCPD = Φ( (Φ⁻¹(PD) + √ρ · Φ⁻¹(confidence)) / √(1-ρ) )
+    """
+    if pd_ttc <= 0:
+        return 0.0
+    
+    sqrt_rho = np.sqrt(rho)
+    sqrt_1_rho = np.sqrt(1.0 - rho)
+    
+    numerator = norm.ppf(pd_ttc) + sqrt_rho * norm.ppf(confidence)
+    return float(norm.cdf(numerator / sqrt_1_rho))