合成控制法进阶：矩阵补全与因子模型扩展

合成控制法（Synthetic Control Method, SCM）作为过去二十年间因果推断领域最具影响力的创新之一，通过构建反事实框架解决了政策评估中的"不可观测结果"难题。传统SCM在处理低维面板数据时表现出色，但当面对高维协变量、时变混杂因子或稀疏干预场景时，其线性组合假设显得力不从心。本文将揭示**矩阵补全（Matrix Completion）与因子模型（Factor Models）**如何为SCM注入新的数学活力，使其在数字化实验时代突破维数灾难的桎梏。

I. 经典合成控制法的数学框架与局限

1.1 潜在结果框架与反事实合成

设我们有 $J+1$ 个单位，单位 $i=1$ 为被干预的处理组，$i=2,\ldots,J+1$ 为未受干预的对照组。观测期为 $t=1,\ldots,T_0,T_0+1,\ldots,T$，其中 $T_0$ 为干预发生时间点。

对于处理组单位，我们观测到的结果变量为：

$$Y_{1t}^{\text{observed}} = \begin{cases} Y_{1t}^0, & t \leq T_0 \\ Y_{1t}^1, & t > T_0 \end{cases}$$

SCM的核心目标是估计 干预效应 $\tau_{1t} = Y_{1t}^1 - Y_{1t}^0$，其中 $Y_{1t}^0$ 是反事实结果。

经典SCM通过凸组合构建合成控制：

$$\hat{Y}_{1t}^0 = \sum_{i=2}^{J+1} w_i Y_{it}, \quad \text{s.t.} \quad w_i \geq 0, \sum_{i=2}^{J+1} w_i = 1$$

权重 $w$ 通过最小化干预前拟合误差确定：

$$w^* = \arg\min_{w \in \mathcal{W}} \sum_{t=1}^{T_0} (Y_{1t} - \sum_{i=2}^{J+1} w_i Y_{it})^2 + \lambda \sum_{k=1}^K (X_{1k} - \sum_{i=2}^{J+1} w_i X_{ik})^2$$

其中 $X_{ik}$ 为预测变量。

实例分析：加州控烟法案评估

Abadie et al. (2010) 研究加州1988年99号控烟法案对人均香烟销量的影响。数据包含加州和38个控制州1980-2000年的年度数据。关键挑战在于：** smoker-related variables ** 与 ** 州特异趋势**的混淆。

原始数据的典型结构如下表：

# 经典SCM实现：classic_scm.py
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from typing import List, Dict, Tuple

class ClassicSCM:
    def __init__(self, data: pd.DataFrame, 
                 unit_col: str, time_col: str, 
                 outcome_col: str, treat_unit: str, 
                 treat_period: int):
        """
        经典合成控制法实现
        
        参数：
        - data: 面板数据，需包含unit_id, time, outcome
        - treat_unit: 处理单位名称
        - treat_period: 干预开始时间点
        """
        self.data = data.set_index([unit_col, time_col])
        self.unit_col = unit_col
        self.time_col = time_col
        self.outcome_col = outcome_col
        self.treat_unit = treat_unit
        self.treat_period = treat_period
        
        # 分离处理组与对照组
        self.Y_treat = self.data.xs(treat_unit, level=unit_col)[outcome_col]
        self.Y_control = self.data.drop(treat_unit, level=unit_col)[outcome_col].unstack(level=0)
        
        # 时间索引
        self.time_pre = self.Y_treat.index[self.Y_treat.index < treat_period]
        self.time_post = self.Y_treat.index[self.Y_treat.index >= treat_period]
        
        self.weights = None
        self.Y_synth = None
    
    def fit(self, X_predictors: List[str] = None, 
            lambda_val: float = 0.0) -> np.ndarray:
        """
        估计合成控制权重
        
        X_predictors: 用于匹配的预测变量名
        lambda_val: 正则化参数
        """
        Y_pre_treat = self.Y_treat.loc[self.time_pre].values
        Y_pre_control = self.Y_control.loc[self.time_pre].values
        
        n_controls = Y_pre_control.shape[1]
        
        # 定义损失函数
        def loss(w):
            # 主损失：干预前拟合误差
            pred_error = np.sum((Y_pre_treat - Y_pre_control @ w)**2)
            
            # 预测变量损失
            pred_var_error = 0.0
            if X_predictors:
                X_treat = self.data.xs(self.treat_unit)[X_predictors].mean().values
                X_control = self.data.drop(self.treat_unit)[X_predictors].groupby(level=1).mean().values
                pred_var_error = np.sum((X_treat - X_control @ w)**2)
            
            return pred_error + lambda_val * pred_var_error
        
        # 约束条件（非负且和为1）
        constraints = [
            {'type': 'eq', 'fun': lambda w: np.sum(w) - 1}
        ]
        bounds = [(0, 1) for _ in range(n_controls)]
        
        # 初始值
        w0 = np.ones(n_controls) / n_controls
        
        # 优化
        result = minimize(loss, w0, method='SLSQP', 
                         bounds=bounds, constraints=constraints)
        
        self.weights = result.x
        return self.weights
    
    def predict_counterfactual(self) -> pd.Series:
        """预测反事实结果"""
        if self.weights is None:
            raise ValueError("必须先调用fit()")
        
        # 使用对照组的加权平均作为反事实
        Y_post_control = self.Y_control.loc[self.time_post].values
        Y_synth_post = Y_post_control @ self.weights
        
        # 合并前后期
        Y_synth_pre = self.Y_treat.loc[self.time_pre].values  # 实际观测
        Y_synth = np.concatenate([Y_synth_pre, Y_synth_post])
        
        all_time = np.concatenate([self.time_pre, self.time_post])
        
        self.Y_synth = pd.Series(Y_synth, index=all_time, name='synthetic')
        return self.Y_synth
    
    def treatment_effect(self) -> pd.Series:
        """计算干预效应"""
        Y_obs = self.Y_treat.reindex(self.Y_synth.index)
        return Y_obs - self.Y_synth
    
    def mape_pre_treatment(self) -> float:
        """干预前MAPE拟合优度"""
        if self.Y_synth is None:
            self.predict_counterfactual()
        
        mape = np.mean(np.abs(
            self.Y_treat.loc[self.time_pre] - self.Y_synth.loc[self.time_pre]
        ) / self.Y_treat.loc[self.time_pre])
        
        return mape

# 数据加载与预处理
def load_california_data() -> pd.DataFrame:
    """
    加载加州控烟数据并进行预处理
    
    原始数据通常格式：
    - state: 州名
    - year: 年份
    - cigsale: 人均香烟销量（包/人/年）
    - other predictors: 收入、年龄结构等
    """
    # 此处为模拟生成符合特征的数据
    np.random.seed(42)
    states = ['California'] + [f'State_{i}' for i in range(1, 39)]
    years = range(1980, 2001)
    
    data = []
    for state in states:
        for year in years:
            # 生成带有趋势和州特异性的数据
            base_trend = 120 - (year - 1980) * 2 + np.random.normal(0, 2)
            state_effect = np.random.normal(0, 5)
            
            # 加州在1988年后的真实干预效应
            if state == 'California' and year >= 1989:
                base_trend -= 15  # 政策效应使销量下降
                
            data.append({
                'state': state,
                'year': year,
                'cigsale': max(0, base_trend + state_effect)
            })
    
    df = pd.DataFrame(data)
    return df

# 使用示例
if __name__ == "__main__":
    df = load_california_data()
    scm = ClassicSCM(
        data=df,
        unit_col='state',
        time_col='year',
        outcome_col='cigsale',
        treat_unit='California',
        treat_period=1989
    )
    
    # 拟合权重
    weights = scm.fit()
    print("合成控制权重:", {f'State_{i}': w for i, w in enumerate(weights, 1)})
    
    # 预测反事实
    Y_synth = scm.predict_counterfactual()
    te = scm.treatment_effect()
    
    print(f"干预前MAPE: {scm.mape_pre_treatment():.2%}")
    print(f"1990年干预效应: {te.loc[1990]:.2f} 包/人")

1.2 经典SCM的三大理论局限

这些局限在数字时代被放大：互联网平台的A/B测试涉及数千个用户细分群、上百个协变量，经典SCM的优化问题变得病态。

II. 矩阵补全视角：重构反事实的核范数魔法

2.1 从面板数据到低秩矩阵

关键洞察：面板数据$Y_{it}$可视为一个$N \times T$矩阵$\mathbf{Y}$，其中$N=J+1$为单位数，$T$为时期数。干预后处理组单元的数据是结构缺失的。

假设真实数据生成过程（DGP）为：

$$\mathbf{Y} = \mathbf{L}^* + \mathbf{E}, \quad \mathbf{L}^* = \mathbf{U}\mathbf{V}^\top$$

其中$\mathbf{L}^*$是低秩成分（秩$r \ll \min(N,T)$），$\mathbf{E}$为噪声。干预效应对应于$\mathbf{L}^*$的缺失元素估计误差。

实例分析：电商平台促销活动评估

某电商在2023年"618"期间对100个高端用户群进行专属促销，需评估GMV提升。数据维度为100群×365天，但处理组的180天干预后数据缺失。传统SCM需从2000个对照群中选权重，计算不可行。而矩阵补全将问题转化为：利用全平台用户群的低维潜因子结构，直接补全缺失块。

# matrix_completion_scm.py
import numpy as np
from typing import Optional, Tuple
from scipy.sparse.linalg import svds

class MatrixCompletionSCM:
    """
    基于核范数最小化的合成控制
    
    数学形式：
    min_L ||P_Ω(Y - L)||_F^2 + λ||L||_*
    其中P_Ω是观测投影算子，||·||_*是核范数
    """
    
    def __init__(self, rank: int = 5, lambda_reg: float = 0.01, 
                 max_iter: int = 1000, tol: float = 1e-6):
        """
        参数：
        - rank: 潜在因子秩的先验估计
        - lambda_reg: 核范数正则化强度
        - max_iter: ADMM迭代次数
        - tol: 收敛阈值
        """
        self.rank = rank
        self.lambda_reg = lambda_reg
        self.max_iter = max_iter
        self.tol = tol
        
        self.L_est = None  # 低秩估计
        self.W = None      # 权重矩阵（自动学习）
    
    def fit(self, Y: np.ndarray, 
            treated_unit_idx: int,
            treat_start: int) -> np.ndarray:
        """
        Y: 完整的面板数据矩阵 (N x T)
        treated_unit_idx: 处理组在矩阵中的行索引
        treat_start: 干预开始列索引
        """
        self.N, self.T = Y.shape
        self.treated_idx = treated_unit_idx
        self.treat_start = treat_start
        
        # 构造掩码矩阵：1为观测，0为缺失
        self.M = np.ones_like(Y)
        self.M[treated_unit_idx, treat_start:] = 0  # 干预后数据缺失
        
        # 使用ADMM求解核范数最小化
        self.L_est = self._admm_solver(Y, self.M)
        
        # 从L_est中提取合成权重（解析解）
        self._extract_synthetic_weights()
        
        return self.W
    
    def _admm_solver(self, Y: np.ndarray, M: np.ndarray) -> np.ndarray:
        """ADMM算法求解核范数正则化问题"""
        # ADMM变量初始化
        L = np.zeros_like(Y)  # 低秩成分
        Z = np.zeros_like(Y)  # 辅助变量
        U = np.zeros_like(Y)  # 对偶变量
        
        rho = 1.0  # 增广拉格朗日参数
        
        for iteration in range(self.max_iter):
            # L-step: 近端算子（奇异值阈值）
            C = Z - U
            # 对C进行SVD
            U_svd, s, Vt = np.linalg.svd(C, full_matrices=False)
            # 奇异值软阈值
            s_thresh = np.maximum(s - self.lambda_reg / rho, 0)
            # 重建L
            L_new = U_svd @ np.diag(s_thresh) @ Vt
            
            # Z-step: 最小二乘拟合观测值
            # Z = argmin ||M*(Y - Z)||_F^2 + (ρ/2)||Z - (L + U)||_F^2
            numerator = M * Y + rho * (L_new + U)
            denominator = M + rho
            Z_new = numerator / denominator
            
            # U-step: 更新对偶变量
            U_new = U + (L_new - Z_new)
            
            # 收敛检查
            primal_residual = np.linalg.norm(L_new - Z_new, 'fro')
            dual_residual = rho * np.linalg.norm(Z_new - Z, 'fro')
            
            if primal_residual < self.tol and dual_residual < self.tol:
                print(f"ADMM converged at iteration {iteration}")
                break
            
            L, Z, U = L_new, Z_new, U_new
            
            # 动态调整rho
            if primal_residual > 10 * dual_residual:
                rho *= 2
            elif dual_residual > 10 * primal_residual:
                rho /= 2
        
        return L
    
    def _extract_synthetic_weights(self) -> None:
        """
        从低秩估计中提取合成权重
        
        关键洞察：L_est[treated] ≈ Σ w_i * L_est[control_i]
        通过最小化干预前重构误差求解w
        """
        # 使用干预前的低秩成分
        L_pre = self.L_est[:, :self.treat_start]
        
        # 分离处理组与对照组
        L_treated_pre = L_pre[self.treated_idx, :]  # (T_pre,)
        L_control_pre = L_pre[np.arange(self.N) != self.treated_idx, :]  # (N-1, T_pre)
        
        # 岭回归求解权重（比单纯形约束更灵活）
        from sklearn.linear_model import Ridge
        
        # 转置：将时间作为特征，单位作为样本
        reg = Ridge(alpha=1e-6, fit_intercept=False)
        reg.fit(L_control_pre.T, L_treated_pre)
        
        self.weights = reg.coef_
        
        # 投影到单纯形（可选，保持经济解释性）
        self.weights = self._project_to_simplex(self.weights)
        
    def _project_to_simplex(self, v: np.ndarray) -> np.ndarray:
        """将权重投影到概率单纯形"""
        n = len(v)
        if v.sum() == 1 and (v >= 0).all():
            return v
        
        # 使用O(n log n)算法
        u = np.sort(v)[::-1]
        cssv = np.cumsum(u)
        rho = np.nonzero(u * np.arange(1, n+1) > (cssv - 1))[0][-1]
        theta = (cssv[rho] - 1) / (rho + 1)
        
        return np.maximum(v - theta, 0)
    
    def predict_counterfactual(self) -> np.ndarray:
        """生成反事实预测"""
        if self.L_est is None:
            raise ValueError("必须先调用fit()")
        
        # 处理组的反事实 = 低秩估计的对应部分
        return self.L_est[self.treated_idx, self.treat_start:]
    
    def treatment_effect(self, Y_obs_post: np.ndarray) -> np.ndarray:
        """计算干预效应"""
        Y_cf = self.predict_counterfactual()
        return Y_obs_post - Y_cf

# 实例：高维用户分群评估
def demo_matrix_completion():
    """
    模拟100个用户群、365天的GMV数据
    其中10个群在180天后接受促销干预
    """
    np.random.seed(123)
    N, T = 100, 365
    treat_start = 180
    n_treated = 10
    
    # 生成低秩矩阵（秩=5）
    rank = 5
    U = np.random.randn(N, rank)
    V = np.random.randn(T, rank)
    Y_true = U @ V.T
    
    # 添加噪声
    Y_obs = Y_true + np.random.randn(N, T) * 0.1
    
    # 模拟干预效应：处理组后期增加
    treatment_effect = np.random.randn(n_treated, T - treat_start) * 5 + 10
    Y_obs[:n_treated, treat_start:] += treatment_effect
    
    # 应用矩阵补全SCM
    mc_scm = MatrixCompletionSCM(rank=5, lambda_reg=0.1)
    
    # 评估第一个处理单位
    mc_scm.fit(Y_obs, treated_unit_idx=0, treat_start=treat_start)
    
    Y_cf = mc_scm.predict_counterfactual()
    te_est = Y_cf - Y_true[0, treat_start:]  # 使用真实反事实对比
    
    print(f"矩阵补全SCM权重稀疏度: {(mc_scm.weights > 0.01).sum()}/{len(mc_scm.weights)}")
    print(f"干预效应估计RMSE: {np.sqrt(np.mean(te_est**2)):.3f}")
    
    return mc_scm, Y_obs, Y_true

if __name__ == "__main__":
    mc_scm, Y_obs, Y_true = demo_matrix_completion()

2.2 核范数正则化的统计意义

核范数$\|\mathbf{L}\|_* = \sum_i \sigma_i(\mathbf{L})$是矩阵秩的凸松弛。相比经典SCM的硬约束，矩阵补全提供软正则化：

关键优势：当$J \gg T_0$（对照组远多于时间维度）时，经典SCM的优化问题病态，而矩阵补全通过低秩约束保持良态。

III. 因子模型扩展：捕捉时变结构的动态合成

3.1 从静态权重到动态载荷

经典SCM假设合成权重$w$在干预前后时间不变，这在许多场景中不现实。因子模型捕捉了时变因子载荷：

数据生成过程：

$$Y_{it} = \lambda_i' F_t + \epsilon_{it}, \quad i=1,\ldots,N; t=1,\ldots,T$$

其中$F_t \in \mathbb{R}^r$为时变因子，$\lambda_i$为单位载荷。处理组的反事实为：

$$Y_{1t}^0 = \lambda_1' F_t$$

合成控制的因子视角：通过对照组估计$(\hat{F}_t, \hat{\lambda}_i)$后，用处理组的载荷$\lambda_1$生成反事实。

实例分析：网约车平台动态定价效应评估

某平台在2023年7月对5个城市实施动态定价（最高溢价30%），需评估对司机在线时长的影响。传统SCM假设合成城市权重固定，但忽略了：① 节假日效应因子 ② 天气冲击因子 ③ 竞品活动因子的时变影响。因子模型SCM可动态调整各因子权重。

# factor_model_scm.py
import numpy as np
from scipy.linalg import svd
from sklearn.preprocessing import StandardScaler
import pandas as pd

class FactorModelSCM:
    """
    因子模型扩展的合成控制
    
    模型设定：
    Y_it = λ_i'F_t + ε_it
    其中F_t通过对照组数据估计，λ_1通过干预前数据估计
    """
    
    def __init__(self, n_factors: int = 3, 
                 factor_model: str = "PCA",
                 intercept: bool = True):
        """
        参数：
        - n_factors: 潜在因子数r
        - factor_model: "PCA"或"IPCA"(迭代PCA)
        - intercept: 是否包含单位固定效应
        """
        self.n_factors = n_factors
        self.factor_model = factor_model
        self.intercept = intercept
        
        # 估计结果
        self.F = None          # 因子矩阵 (T x r)
        self.lambda_loadings = None  # 载荷矩阵 (N x r)
        self.error_component = None  # 误差项
    
    def fit(self, Y: np.ndarray, treated_idx: int, 
            treat_start: int, **kwargs) -> 'FactorModelSCM':
        """
        两阶段估计
        
        阶段1: 使用对照组估计因子F_t
        阶段2: 使用处理组干预前数据估计λ_1
        """
        self.N, self.T = Y.shape
        self.treated_idx = treated_idx
        self.treat_start = treat_start
        
        # 阶段1: 因子估计（仅使用对照组数据）
        control_mask = np.ones(self.N, dtype=bool)
        control_mask[treated_idx] = False
        
        # 对照组数据 (N-1, T)
        Y_control = Y[control_mask, :]
        
        if self.factor_model == "PCA":
            self.F = self._estimate_pca_factors(Y_control)
        elif self.factor_model == "IPCA":
            self.F = self._estimate_ipca_factors(Y_control, **kwargs)
        else:
            raise ValueError("Unsupported factor model")
        
        # 阶段2: 载荷估计（使用全样本干预前）
        Y_pre = Y[:, :treat_start]
        self.lambda_loadings = self._estimate_loadings(Y_pre, self.F[:treat_start, :])
        
        # 计算误差成分
        self.error_component = Y - self.lambda_loadings @ self.F.T
        
        return self
    
    def _estimate_pca_factors(self, Y_control: np.ndarray) -> np.ndarray:
        """PCA估计因子"""
        # 时间维度作为主成分方向
        # 对Y_control的转置进行SVD: (T x N_ctrl)
        T = Y_control.shape[1]
        
        # 中心化
        Y_centered = Y_control.T - Y_control.mean(axis=1)  # 时间维度中心化
        
        # SVD
        U, s, Vt = np.linalg.svd(Y_centered, full_matrices=False)
        
        # 前r个主成分作为因子
        F = U[:, :self.n_factors] * s[:self.n_factors]  # (T x r)
        
        return F
    
    def _estimate_ipca_factors(self, Y_control: np.ndarray, 
                               max_iter: int = 100) -> np.ndarray:
        """迭代PCA（处理缺失数据更稳健）"""
        T = Y_control.shape[1]
        N_ctrl = Y_control.shape[0]
        
        # 初始化
        F = np.random.randn(T, self.n_factors)
        
        for iteration in range(max_iter):
            # 固定F，估计载荷
            Lambda = Y_control @ F @ np.linalg.inv(F.T @ F + 0.01 * np.eye(self.n_factors))
            
            # 固定载荷，估计F
            F_new = Y_control.T @ Lambda @ np.linalg.inv(Lambda.T @ Lambda + 0.01 * np.eye(self.n_factors))
            
            # 标准化
            F_new = StandardScaler().fit_transform(F_new)
            
            # 收敛检查
            if np.linalg.norm(F_new - F, 'fro') / np.linalg.norm(F, 'fro') < 1e-6:
                break
            
            F = F_new
        
        return F
    
    def _estimate_loadings(self, Y_pre: np.ndarray, 
                           F_pre: np.ndarray) -> np.ndarray:
        """通过时间序列回归估计载荷"""
        # Y_pre: (N, T_pre), F_pre: (T_pre, r)
        N = Y_pre.shape[0]
        
        lambda_est = np.zeros((N, self.n_factors))
        
        for i in range(N):
            # 每个单位的载荷回归
            # λ_i = (F'F)^{-1}F'y_i
            lambda_est[i, :] = np.linalg.lstsq(F_pre, Y_pre[i, :], rcond=None)[0]
        
        return lambda_est
    
    def predict_counterfactual(self) -> np.ndarray:
        """生成反事实路径"""
        if self.F is None or self.lambda_loadings is None:
            raise ValueError("必须先调用fit()")
        
        # 处理组的反事实: λ_1'F_t
        lambda_treated = self.lambda_loadings[self.treated_idx, :]  # (r,)
        
        # 全时段预测
        Y_cf_all = lambda_treated @ self.F.T  # (T,)
        
        # 仅返回干预后
        return Y_cf_all[self.treat_start:]
    
    def factor_contribution(self) -> pd.DataFrame:
        """返回各因子对处理组的贡献度"""
        lambda_treated = self.lambda_loadings[self.treated_idx, :]
        
        contributions = pd.DataFrame({
            'factor_id': range(self.n_factors),
            'loading': lambda_treated,
            'explained_variance': lambda_treated**2
        })
        
        contributions['contribution_pct'] = (
            contributions['explained_variance'] / 
            contributions['explained_variance'].sum() * 100
        )
        
        return contributions.sort_values('contribution_pct', ascending=False)
    
    def dynamic_weights(self, t: int) -> np.ndarray:
        """计算时间t的动态合成权重"""
        # 在时刻t，权重为 λ_i'F_t 的归一化
        contributions = self.lambda_loadings @ self.F[t, :]  # (N,)
        
        # 排除处理组自身
        weights = contributions.copy()
        weights[self.treated_idx] = 0
        
        # 归一化到非负（投影）
        weights = np.maximum(weights, 0)
        if weights.sum() > 0:
            weights /= weights.sum()
        
        return weights

# 网约车动态定价案例
def demo_ride_hailing():
    """
    模拟50个城市、180天的司机在线时长数据
    其中5个城市在第90天开始动态定价
    数据包含3个潜在因子：周末效应、天气冲击、竞品活动
    """
    np.random.seed(42)
    N, T = 50, 180
    treat_start = 90
    n_factors = 3
    
    # 生成因子
    # F1: 周末效应（周期7）
    F1 = np.sin(np.arange(T) * 2 * np.pi / 7) + np.random.randn(T) * 0.1
    
    # F2: 天气冲击（AR1过程）
    F2 = np.zeros(T)
    for t in range(1, T):
        F2[t] = 0.8 * F2[t-1] + np.random.randn() * 0.5
    
    # F3: 竞品活动（泊松跳跃）
    F3 = np.random.poisson(0.1, T) * 5
    
    F = np.column_stack([F1, F2, F3])
    
    # 生成载荷
    Lambda = np.random.randn(N, n_factors)
    
    # 基础数据
    Y = Lambda @ F.T + np.random.randn(N, T) * 0.5
    
    # 模拟干预效应：动态定价提升司机在线时长
    # 但对不同因子载荷的城市效果异质
    treated_cities = np.arange(5)
    for i in treated_cities:
        # 效应强度与λ_i1（周末敏感性）正相关
        effect = 5 + Lambda[i, 0] * 3
        Y[i, treat_start:] += effect
    
    # 因子模型SCM
    fm_scm = FactorModelSCM(n_factors=3, factor_model="PCA")
    fm_scm.fit(Y, treated_idx=0, treat_start=treat_start)
    
    # 预测与评估
    Y_cf = fm_scm.predict_counterfactual()
    Y_true_cf = Lambda[0, :] @ F[treat_start:, :].T
    
    print("因子贡献度:")
    print(fm_scm.factor_contribution())
    
    print(f"\n动态定价效应RMSE: {np.sqrt(np.mean((Y_cf - Y_true_cf)**2)):.3f}")
    
    # 查看权重动态性
    weights_t90 = fm_scm.dynamic_weights(t=90)
    weights_t120 = fm_scm.dynamic_weights(t=120)
    
    print(f"\n第90天权重离散度: {np.std(weights_t90):.4f}")
    print(f"第120天权重离散度: {np.std(weights_t120):.4f}")
    
    return fm_scm, Y, F

if __name__ == "__main__":
    fm_scm, Y, F = demo_ride_hailing()

3.2 IPCA因子估计：处理非平衡面板

在真实场景中，对照组可能存在缺失数据（如某些州统计延迟）。Iterative PCA通过EM算法思想迭代估计：

** IPCA算法步骤**：
I. 初始化：用列均值填充缺失，PCA得$F^{(0)}$
II. E步：固定$F$，通过回归估计$\lambda_i$（仅使用观测时段）
III. M步：固定$\Lambda$，更新$F = (\Lambda'\Lambda)^{-1}\Lambda'Y$
IV. 重复：直至$\|F^{(k+1)}-F^{(k)}\|<\epsilon$

IV. 完整代码实现与加州控烟法案重分析

4.1 端到端实验框架

现在整合三种方法（经典SCM、矩阵补全、因子模型），对加州控烟法案进行系统性重分析，并构建可重复实验的Pipeline。

# pipeline.py
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit
from classic_scm import ClassicSCM
from matrix_completion_scm import MatrixCompletionSCM
from factor_model_scm import FactorModelSCM
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

class SCMExperimentPipeline:
    """
    合成控制法实验流水线
    
    功能：
    1. 数据加载与清洗
    2. 模型拟合与交叉验证
    3. 多方法对比评估
    4. 效应可视化与稳健性检验
    """
    
    def __init__(self, data_path: str, config: Dict):
        self.data_path = data_path
        self.config = config
        
        # 实验配置
        self.treat_unit = config['treat_unit']
        self.treat_period = config['treat_period']
        self.outcome_var = config['outcome_var']
        self.unit_var = config['unit_var']
        self.time_var = config['time_var']
        
        # 加载数据
        self.raw_data = self._load_data()
        self.panel_data = self._create_panel()
        
        # 结果存储
        self.results = {}
        self.models = {}
    
    def _load_data(self) -> pd.DataFrame:
        """加载并清洗原始数据"""
        df = pd.read_csv(self.data_path)
        
        # 数据质量检查
        print(f"数据加载完成: {df.shape}")
        print(f"缺失值统计:\n{df.isnull().sum().head()}")
        
        # 处理缺失
        df = df.dropna(subset=[self.outcome_var])
        
        # 时间格式化
        df[self.time_var] = pd.to_numeric(df[self.time_var])
        
        return df
    
    def _create_panel(self) -> pd.DataFrame:
        """构造平衡面板"""
        # 创建完整索引
        units = self.raw_data[self.unit_var].unique()
        times = range(
            self.raw_data[self.time_var].min(),
            self.raw_data[self.time_var].max() + 1
        )
        
        # 检查平衡性
        panel_counts = self.raw_data.groupby(self.unit_var)[self.time_var].nunique()
        print(f"面板平衡性:\n{panel_counts.describe()}")
        
        # 填充缺失（使用前值填充）
        panel = self.raw_data.set_index([self.unit_var, self.time_var])
        panel = panel.reindex(
            pd.MultiIndex.from_product([units, times], names=[self.unit_var, self.time_var])
        )
        panel[self.outcome_var] = panel.groupby(level=0)[self.outcome_var].fillna(method='ffill')
        
        return panel.reset_index()
    
    def prepare_features(self, feature_cols: List[str]) -> np.ndarray:
        """提取预测变量矩阵"""
        return self.panel_data.pivot(
            index=self.unit_var, columns=self.time_var, values=feature_cols
        ).values
    
    def run_classic_scm(self) -> Dict[str, float]:
        """运行经典SCM"""
        print("\n=== 运行经典SCM ===")
        
        # 数据准备
        outcome_matrix = self.panel_data.pivot(
            index=self.unit_var, columns=self.time_var, 
            values=self.outcome_var
        ).values
        
        # 获取单位索引
        units = self.panel_data[self.unit_var].unique()
        treat_idx = np.where(units == self.treat_unit)[0][0]
        
        # 实例化模型
        scm = ClassicSCM(
            data=self.panel_data,
            unit_col=self.unit_var,
            time_col=self.time_var,
            outcome_col=self.outcome_var,
            treat_unit=self.treat_unit,
            treat_period=self.treat_period
        )
        
        # 拟合（无预测变量）
        weights = scm.fit()
        
        # 预测与评估
        Y_cf = scm.predict_counterfactual()
        te = scm.treatment_effect()
        
        # 计算指标
        mape = scm.mape_pre_treatment()
        post_rmspe = np.sqrt(
            np.mean(te.loc[self.config['post_period']]**2)
        )
        
        result = {
            'method': 'Classic SCM',
            'weights': {units[i]: w for i, w in enumerate(weights) if w > 0.01},
            'pre_mape': mape,
            'post_rmspe': post_rmspe,
            'avg_te': te.loc[self.config['post_period']].mean(),
            'model': scm
        }
        
        self.results['classic'] = result
        self.models['classic'] = scm
        
        print(f"权重非零数: {(weights > 0.01).sum()}")
        print(f"干预前MAPE: {mape:.2%}")
        print(f"平均干预效应: {result['avg_te']:.2f}")
        
        return result
    
    def run_matrix_completion_scm(self, rank: int = 5, lambda_reg: float = 0.1) -> Dict:
        """运行矩阵补全SCM"""
        print("\n=== 运行矩阵补全SCM ===")
        
        # 准备矩阵
        Y_matrix = self.panel_data.pivot(
            index=self.unit_var, columns=self.time_var, 
            values=self.outcome_var
        ).values
        
        units = self.panel_data[self.unit_var].unique()
        treat_idx = np.where(units == self.treat_unit)[0][0]
        
        # 时间索引转换
        times = self.panel_data[self.time_var].unique()
        treat_time_idx = np.where(times == self.treat_period)[0][0]
        
        # 运行模型
        mc_scm = MatrixCompletionSCM(rank=rank, lambda_reg=lambda_reg)
        mc_scm.fit(Y_matrix, treated_idx=treat_idx, treat_start=treat_time_idx)
        
        # 评估
        Y_obs_post = Y_matrix[treat_idx, treat_time_idx:]
        Y_cf = mc_scm.predict_counterfactual()
        
        # 计算干预前拟合优度
        Y_pre_treat = Y_matrix[treat_idx, :treat_time_idx]
        Y_synth_pre = mc_scm.L_est[treat_idx, :treat_time_idx]
        mape = np.mean(np.abs(Y_pre_treat - Y_synth_pre) / Y_pre_treat)
        
        te_est = Y_obs_post - Y_cf
        post_rmspe = np.sqrt(np.mean(te_est**2))
        
        result = {
            'method': 'Matrix Completion SCM',
            'weights': {units[i]: w for i, w in enumerate(mc_scm.weights) if w > 0.01},
            'pre_mape': mape,
            'post_rmspe': post_rmspe,
            'avg_te': te_est.mean(),
            'model': mc_scm,
            'L_est': mc_scm.L_est
        }
        
        self.results['matrix_completion'] = result
        self.models['matrix_completion'] = mc_scm
        
        print(f"核范数: {np.linalg.norm(mc_scm.L_est, 'nuc'):.2f}")
        print(f"平均干预效应: {result['avg_te']:.2f}")
        
        return result
    
    def run_factor_scm(self, n_factors: int = 3) -> Dict:
        """运行因子模型SCM"""
        print("\n=== 运行因子模型SCM ===")
        
        # 准备矩阵
        Y_matrix = self.panel_data.pivot(
            index=self.unit_var, columns=self.time_var, 
            values=self.outcome_var
        ).values
        
        units = self.panel_data[self.unit_var].unique()
        treat_idx = np.where(units == self.treat_unit)[0][0]
        
        # 时间索引
        times = self.panel_data[self.time_var].unique()
        treat_time_idx = np.where(times == self.treat_period)[0][0]
        
        # 运行模型
        fm_scm = FactorModelSCM(n_factors=n_factors, factor_model="PCA")
        fm_scm.fit(Y_matrix, treated_idx=treat_idx, treat_start=treat_time_idx)
        
        # 评估
        Y_obs_post = Y_matrix[treat_idx, treat_time_idx:]
        Y_cf = fm_scm.predict_counterfactual()
        
        # 计算干预前拟合优度
        Y_pre_treat = Y_matrix[treat_idx, :treat_time_idx]
        Y_synth_pre = fm_scm.L_est[treat_idx, :treat_time_idx] if fm_scm.L_est is not None else fm_scm.lambda_loadings[treat_idx] @ fm_scm.F[:treat_time_idx].T
        mape = np.mean(np.abs(Y_pre_treat - Y_synth_pre) / Y_pre_treat)
        
        te_est = Y_obs_post - Y_cf
        post_rmspe = np.sqrt(np.mean(te_est**2))
        
        result = {
            'method': 'Factor Model SCM',
            'factor_contributions': fm_scm.factor_contribution(),
            'pre_mape': mape,
            'post_rmspe': post_rmspe,
            'avg_te': te_est.mean(),
            'model': fm_scm,
            'F': fm_scm.F,
            'lambda': fm_scm.lambda_loadings
        }
        
        self.results['factor_model'] = result
        self.models['factor_model'] = fm_scm
        
        print(f"因子贡献度:\n{result['factor_contributions']}")
        print(f"平均干预效应: {result['avg_te']:.2f}")
        
        return result
    
    def run_robustness_checks(self) -> Dict:
        """稳健性检验"""
        print("\n=== 运行稳健性检验 ===")
        
        checks = {}
        
        # 1. 安慰剂检验（随机干预点）
        placebo_effects = []
        for placebo_time in [1985, 1986, 1987]:  # 真实干预前
            try:
                scm_placebo = ClassicSCM(
                    data=self.panel_data,
                    unit_col=self.unit_col,
                    time_col=self.time_col,
                    outcome_col=self.outcome_var,
                    treat_unit=self.treat_unit,
                    treat_period=placebo_time
                )
                scm_placebo.fit()
                te_placebo = scm_placebo.treatment_effect()
                placebo_effects.append(te_placebo.loc[placebo_time:].mean())
            except:
                continue
        
        checks['placebo_avg_effect'] = np.mean(placebo_effects) if placebo_effects else 0
        checks['placebo_std'] = np.std(placebo_effects) if placebo_effects else 1
        
        # 2. 留出法交叉验证（时间序列）
        Y_matrix = self.panel_data.pivot(
            index=self.unit_var, columns=self.time_var, 
            values=self.outcome_var
        ).values
        
        tscv = TimeSeriesSplit(n_splits=3)
        cv_scores = []
        
        for train_idx, val_idx in tscv.split(Y_matrix.T):  # 时间维度切分
            # 使用train段作为干预前
            train_start = train_idx[0]
            train_end = train_idx[-1]
            
            # 模拟SCM拟合
            # 计算在验证集的预测误差
            pass  # 简化，实际需实现
        
        checks['cv_rmspe'] = np.mean(cv_scores) if cv_scores else 0
        
        self.results['robustness'] = checks
        return checks
    
    def visualize_results(self, save_path: str = None):
        """可视化对比结果"""
        fig, axes = plt.subplots(2, 2, figsize=(16, 12))
        
        # 1. 时间序列对比
        ax1 = axes[0, 0]
        
        # 获取真实观测
        units = self.panel_data[self.unit_var].unique()
        treat_idx = np.where(units == self.treat_unit)[0][0]
        times = self.panel_data[self.time_var].unique()
        
        Y_matrix = self.panel_data.pivot(
            index=self.unit_var, columns=self.time_var, 
            values=self.outcome_var
        ).values
        
        # 真实观测
        Y_obs = Y_matrix[treat_idx, :]
        ax1.plot(times, Y_obs, 'k-', linewidth=2, label='Observed (California)')
        
        # 各模型反事实
        for method, result in self.results.items():
            if method == 'robustness':
                continue
            
            model = result['model']
            if hasattr(model, 'L_est'):
                Y_cf = model.L_est[treat_idx, :]
            elif hasattr(model, 'predict_counterfactual'):
                Y_cf_pre = Y_matrix[treat_idx, :model.treat_start]
                Y_cf_post = model.predict_counterfactual()
                Y_cf = np.concatenate([Y_cf_pre, Y_cf_post])
            else:
                continue
            
            ax1.plot(times, Y_cf, '--', label=f'{result["method"]}', alpha=0.8)
        
        ax1.axvline(x=self.treat_period, color='r', linestyle=':')
        ax1.set_title('Observed vs Synthetic Control')
        ax1.legend()
        ax1.set_xlabel('Year')
        ax1.set_ylabel(self.outcome_var)
        
        # 2. 干预效应对比
        ax2 = axes[0, 1]
        
        for method, result in self.results.items():
            if method == 'robustness':
                continue
            
            model = result['model']
            if isinstance(model, ClassicSCM):
                te = model.treatment_effect()
                ax2.plot(te.loc[self.config['post_period']], 
                        label=f'{method}', marker='o')
        
        ax2.axhline(y=0, color='k', linestyle='-')
        ax2.set_title('Treatment Effect Over Time')
        ax2.legend()
        ax2.set_xlabel('Year')
        ax2.set_ylabel('Effect')
        
        # 3. 权重稀疏性对比
        ax3 = axes[1, 0]
        
        methods = []
        sparsity = []
        for method, result in self.results.items():
            if method == 'robustness' or 'weights' not in result:
                continue
            
            methods.append(result['method'])
            sparsity.append((len(result['weights'])))
        
        ax3.bar(methods, sparsity)
        ax3.set_title('Number of Non-zero Weights')
        ax3.set_ylabel('Count')
        plt.setp(ax3.get_xticklabels(), rotation=45)
        
        # 4. 因子贡献（因子模型特有）
        ax4 = axes[1, 1]
        if 'factor_model' in self.results:
            factor_contrib = self.results['factor_model']['factor_contributions']
            ax4.bar(factor_contrib['factor_id'], 
                   factor_contrib['contribution_pct'])
            ax4.set_title('Factor Contributions')
            ax4.set_xlabel('Factor ID')
            ax4.set_ylabel('Contribution (%)')
        
        plt.tight_layout()
        
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
        
        return fig
    
    def generate_report(self) -> str:
        """生成实验报告"""
        report = f"""
        合成控制法实验报告
        生成时间: {datetime.now().isoformat()}
        
        实验配置:
        - 处理单位: {self.treat_unit}
        - 干预时间: {self.treat_period}
        - 样本量: {len(self.panel_data[self.unit_var].unique())} 单位 × 
                 {len(self.panel_data[self.time_var].unique())} 时期
        
        模型对比结果:
        """
        
        for method, result in self.results.items():
            if method == 'robustness':
                continue
            
            report += f"""
            {result['method']}:
            - 干预前MAPE: {result['pre_mape']:.2%}
            - 干预后RMSPE: {result['post_rmspe']:.3f}
            - 平均干预效应: {result['avg_te']:.2f}
            """
        
        if 'robustness' in self.results:
            robust = self.results['robustness']
            report += f"""
            
            稳健性检验:
            - 安慰剂效应均值: {robust['placebo_avg_effect']:.3f}
            - 安慰剂标准差: {robust['placebo_std']:.3f}
            """
        
        return report

# 加州控烟法案实验配置
config = {
    'treat_unit': 'California',
    'treat_period': 1989,
    'outcome_var': 'cigsale',
    'unit_var': 'state',
    'time_var': 'year',
    'post_period': range(1989, 2001)
}

# 使用模拟数据运行（实际应加载真实数据）
if __name__ == "__main__":
    # 生成模拟数据
    from classic_scm import load_california_data
    df_simulated = load_california_data()
    
    # 保存为CSV供流水线加载
    df_simulated.to_csv('synthetic_california_data.csv', index=False)
    
    # 运行实验
    pipeline = SCMExperimentPipeline(
        data_path='synthetic_california_data.csv',
        config=config
    )
    
    # 执行三种方法
    pipeline.run_classic_scm()
    pipeline.run_matrix_completion_scm(rank=3, lambda_reg=0.05)
    pipeline.run_factor_scm(n_factors=2)
    
    # 稳健性检验
    pipeline.run_robustness_checks()
    
    # 可视化
    fig = pipeline.visualize_results(save_path='scm_comparison.png')
    
    # 生成报告
    report = pipeline.generate_report()
    print(report)
    
    # 保存结果
    with open('scm_results.json', 'w') as f:
        import json
        # 处理numpy数组
        def default(o):
            if isinstance(o, np.ndarray):
                return o.tolist()
            return str(o)
        
        json.dump(pipeline.results, f, indent=2, default=default)

V. 高维数据与正则化：从Ridge到Elastic Net

5.1 高维预测变量的挑战

当预测变量数$K$超过时间维度$T_0$时，经典SCM的优化问题病态。引入正则化：

扩展目标函数：

$$w^* = \arg\min_{w \in \mathcal{W}} \|Y_{1,\text{pre}} - Y_{\text{ctrl,pre}} w\|_2^2 + \lambda \left( \alpha \|X_1 - X_{\text{ctrl}} w\|_2^2 + (1-\alpha) \|w\|_1 \right)$$

其中$\|w\|_1$引入稀疏性，自动选择相关对照组。

# regularized_scm.py
from sklearn.linear_model import ElasticNet, LassoCV, RidgeCV
from typing import Union

class RegularizedSCM:
    """
    正则化合成控制法
    
    支持L1、L2、Elastic Net三种正则化
    自动通过交叉验证选择λ
    """
    
    def __init__(self, penalty: str = "elasticnet", 
                 alpha: float = 1.0,
                 l1_ratio: float = 0.5,
                 cv_folds: int = 3):
        """
        参数：
        - penalty: "l1"|"l2"|"elasticnet"
        - alpha: 正则化强度
        - l1_ratio: Elastic Net中L1比例
        - cv_folds: 交叉验证折数
        """
        self.penalty = penalty
        self.alpha = alpha
        self.l1_ratio = l1_ratio
        self.cv_folds = cv_folds
        
        self.weights = None
        self.model = None
        self.coef_path = None
    
    def fit(self, Y_pre_treat: np.ndarray, 
            Y_pre_control: np.ndarray,
            X_treat: np.ndarray = None,
            X_control: np.ndarray = None) -> np.ndarray:
        """
        拟合正则化SCM
        
        参数：
        - Y_pre_treat: (T_pre,) 处理组干预前结果
        - Y_pre_control: (T_pre, N_ctrl) 对照组干预前结果
        - X_treat: (K,) 处理组预测变量
        - X_control: (K, N_ctrl) 对照组预测变量
        """
        T_pre, N_ctrl = Y_pre_control.shape
        
        # 特征矩阵构建：Y + X
        features_y = Y_pre_control.T  # (N_ctrl, T_pre)
        
        if X_control is not None:
            features = np.hstack([features_y, X_control.T])
            target = np.concatenate([Y_pre_treat, X_treat])
        else:
            features = features_y
            target = Y_pre_treat
        
        # 标准化（重要！避免不同量纲影响）
        from sklearn.preprocessing import StandardScaler
        scaler = StandardScaler()
        features_scaled = scaler.fit_transform(features)
        target_scaled = scaler.transform(target.reshape(1, -1)).flatten()
        
        # 选择模型
        if self.penalty == "l1":
            model = LassoCV(cv=self.cv_folds, fit_intercept=False, 
                           positive=True)  # L1支持非负约束
        elif self.penalty == "l2":
            model = RidgeCV(cv=self.cv_folds, fit_intercept=False)
        else:  # elasticnet
            # 需要手动搜索
            from sklearn.model_selection import GridSearchCV
            en = ElasticNet(fit_intercept=False, max_iter=10000)
            param_grid = {'alpha': np.logspace(-3, 3, 20), 
                         'l1_ratio': [0.1, 0.3, 0.5, 0.7, 0.9]}
            model = GridSearchCV(en, param_grid, cv=self.cv_folds)
        
        # 拟合
        model.fit(features_scaled, target_scaled)
        
        if hasattr(model, 'best_params_'):
            print(f"最优参数: {model.best_params_}")
        
        # 获取权重（注意：仅取对照组部分）
        weights_scaled = model.coef_[:N_ctrl]
        
        # 反标准化权重（需保证非负）
        std_factor = np.sqrt(scaler.var_[:N_ctrl])
        self.weights = np.maximum(weights_scaled / std_factor, 0)
        
        # 投影到单纯形
        self.weights = self._project_to_simplex(self.weights)
        
        self.model = model
        return self.weights
    
    def _project_to_simplex(self, v: np.ndarray) -> np.ndarray:
        """投影到概率单纯形（同上）"""
        n = len(v)
        if v.sum() == 1 and (v >= 0).all():
            return v
        
        u = np.sort(v)[::-1]
        cssv = np.cumsum(u)
        rho = np.nonzero(u * np.arange(1, n+1) > (cssv - 1))[0][-1]
        theta = (cssv[rho] - 1) / (rho + 1)
        
        return np.maximum(v - theta, 0)
    
    def get_sparse_weights(self, threshold: float = 1e-4) -> Dict[int, float]:
        """获取大于阈值的非零权重"""
        return {i: w for i, w in enumerate(self.weights) if w > threshold}
    
    def stability_selection(self, Y_pre_treat: np.ndarray,
                           Y_pre_control: np.ndarray,
                           n_bootstraps: int = 100) -> Dict:
        """稳定性选择：评估权重的鲁棒性"""
        stability_scores = np.zeros(Y_pre_control.shape[1])
        
        for b in range(n_bootstraps):
            # 时间维度Bootstrap
            bootstrap_idx = np.random.choice(
                len(Y_pre_treat), 
                size=len(Y_pre_treat), 
                replace=True
            )
            
            Y_pre_treat_boot = Y_pre_treat[bootstrap_idx]
            Y_pre_control_boot = Y_pre_control[bootstrap_idx, :]
            
            # 拟合
            self.fit(Y_pre_treat_boot, Y_pre_control_boot)
            
            # 记录选择
            stability_scores += (self.weights > 1e-4).astype(float)
        
        stability_scores /= n_bootstraps
        
        return {
            'scores': stability_scores,
            'stable_weights': {i: s for i, s in enumerate(stability_scores) if s > 0.8}
        }

# 高维协变量案例
def demo_high_dimensional():
    """
    模拟50个州，每个州有100个社会经济指标
    仅30个时期，构成高维问题
    """
    np.random.seed(42)
    N, T, K = 50, 30, 100
    
    # 生成低秩数据
    rank = 5
    U = np.random.randn(N, rank)
    V = np.random.randn(T, rank)
    L = U @ V.T
    
    # 高维协变量
    X = np.random.randn(N, K)
    
    # 真实权重：仅与X前10个维度相关
    true_weights = np.random.randn(K, 1)
    true_weights[10:] = 0  # 稀疏
    
    # 结果变量
    Y = L + X @ true_weights + np.random.randn(N, T) * 0.1
    
    # 处理组
    treat_idx = 0
    treat_start = 20
    
        # 干预效应
    Y[treat_idx, treat_start:] += 2.0
    
    # 测试不同正则化
    models = {
        'Classic': RegularizedSCM(penalty="l2", alpha=0.0),  # 无正则化
        'Ridge': RegularizedSCM(penalty="l2", alpha=1.0),
        'Lasso': RegularizedSCM(penalty="l1"),
        'ElasticNet': RegularizedSCM(penalty="elasticnet", l1_ratio=0.5)
    }
    
    results = {}
    Y_pre_treat = Y[treat_idx, :treat_start]
    Y_pre_control = Y[1:, :treat_start].T
    X_treat = X[treat_idx, :]
    X_control = X[1:, :].T
    
    for name, model in models.items():
        try:
            w = model.fit(Y_pre_treat, Y_pre_control, X_treat, X_control)
            
            # 预测
            Y_cf = Y[1:, treat_start:].T @ w
            te_est = Y[treat_idx, treat_start:] - Y_cf
            
            results[name] = {
                'weights': w,
                'rmse': np.sqrt(np.mean(te_est**2)),
                'n_nonzero': (w > 0.01).sum(),
                'stability': model.stability_selection(Y_pre_treat, Y_pre_control)
            }
            
            print(f"{name}: RMSE={results[name]['rmse']:.3f}, "
                  f"Nonzero={results[name]['n_nonzero']}")
        except Exception as e:
            print(f"{name} failed: {e}")
    
    # 可视化权重分布
    import matplotlib.pyplot as plt
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    for idx, (name, res) in enumerate(results.items()):
        ax = axes.flat[idx]
        ax.hist(res['weights'], bins=20, alpha=0.7)
        ax.set_title(f'{name} Weight Distribution')
        ax.set_xlabel('Weight Value')
        ax.set_ylabel('Frequency')
        
        # 标注稀疏性
        ax.text(0.05, 0.95, f"Non-zero: {res['n_nonzero']}", 
                transform=ax.transAxes, verticalalignment='top',
                bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    plt.tight_layout()
    plt.savefig('weight_distribution_comparison.png')
    
    return results, Y, X

if __name__ == "__main__":
    results, Y, X = demo_high_dimensional()

5.2 稳定性选择与推断

在计算因果推断中，稳定性选择（Stability Selection）评估权重对数据扰动的鲁棒性：

** Meinshausen准则 **：若变量在至少$\pi$比例的子样本中被选中，则视为稳定。推荐$\pi=0.8$。

# 稳定性选择代码已集成在RegularizedSCM.stability_selection中

VI. 生产部署与性能优化

6.1 实时反事实预测API

将SCM部署为微服务，支持在线政策模拟与What-if分析。

# deployment/api_server.py
from fastapi import FastAPI, HTTPException, Depends
from pydantic import BaseModel, Field
from typing import List, Dict, Optional
import redis
import hashlib
import json
import numpy as np
import pickle
from concurrent.futures import ThreadPoolExecutor

from factor_model_scm import FactorModelSCM
from matrix_completion_scm import MatrixCompletionSCM

# 数据模型
class FitRequest(BaseModel):
    panel_data: List[Dict[str, float]] = Field(..., 
        description="面板数据列表，每项包含unit, time, outcome")
    treat_unit: str
    treat_period: int
    method: str = Field("factor", regex="^(classic|matrix|factor)$")
    hyperparams: Optional[Dict] = None

class PredictRequest(BaseModel):
    model_id: str
    post_periods: List[int]
    confidence_interval: bool = False

class WhatIfRequest(BaseModel):
    model_id: str
    alternative_treat_period: int
    magnitude_multiplier: float = 1.0

app = FastAPI(title="Synthetic Control API", version="1.1.0")

# Redis缓存（生产环境使用集群）
redis_client = redis.Redis(host='localhost', port=6379, db=0)

# 模型存储
class ModelStore:
    def __init__(self, redis_client: redis.Redis):
        self.redis = redis_client
        self.local_cache = {}
    
    def save_model(self, model, model_id: str):
        """序列化并存储模型"""
        # 使用pickle序列化
        pickled = pickle.dumps(model)
        
        # 存储到Redis
        self.redis.set(f"scm:model:{model_id}", pickled)
        
        # 元数据
        meta = {
            'model_id': model_id,
            'type': type(model).__name__,
            'saved_at': datetime.now().isoformat()
        }
        self.redis.set(f"scm:meta:{model_id}", json.dumps(meta))
        
        return model_id
    
    def load_model(self, model_id: str):
        """加载模型"""
        # 检查本地缓存
        if model_id in self.local_cache:
            return self.local_cache[model_id]
        
        # 从Redis加载
        pickled = self.redis.get(f"scm:model:{model_id}")
        if not pickled:
            raise HTTPException(status_code=404, detail="Model not found")
        
        model = pickle.loads(pickled)
        self.local_cache[model_id] = model
        
        return model
    
    def generate_id(self, data_hash: str) -> str:
        """基于数据哈希生成模型ID"""
        return f"scm_{hashlib.sha256(data_hash.encode()).hexdigest()[:12]}"

model_store = ModelStore(redis_client)

# 后台任务执行器
executor = ThreadPoolExecutor(max_workers=4)

@app.post("/model/fit", response_model=Dict)
async def fit_model(request: FitRequest):
    """
    异步拟合模型
    
    工作流程：
    1. 验证数据完整性
    2. 根据哈希检查是否已存在模型
    3. 提交到线程池异步训练
    4. 返回任务ID供查询
    """
    # 数据转换
    df = pd.DataFrame(request.panel_data)
    
    # 生成数据哈希
    data_str = df.to_string() + request.treat_unit + str(request.treat_period)
    model_id = model_store.generate_id(data_str)
    
    # 检查缓存
    if redis_client.exists(f"scm:model:{model_id}"):
        return {
            "status": "cached",
            "model_id": model_id,
            "message": "Model already exists"
        }
    
    # 异步训练
    future = executor.submit(_fit_model_task, request, model_id)
    
    # 存储任务状态
    task_id = f"task_{model_id}"
    redis_client.setex(task_id, 3600, json.dumps({
        "status": "running",
        "model_id": model_id,
        "submitted_at": datetime.now().isoformat()
    }))
    
    return {
        "status": "submitted",
        "task_id": task_id,
        "model_id": model_id
    }

def _fit_model_task(request: FitRequest, model_id: str):
    """后台训练任务"""
    try:
        df = pd.DataFrame(request.panel_data)
        
        # 根据方法选择模型
        if request.method == "classic":
            # 经典SCM需要重构ClassicSCM类以支持dict输入
            pass
        elif request.method == "matrix":
            rank = request.hyperparams.get("rank", 3)
            lambda_reg = request.hyperparams.get("lambda_reg", 0.1)
            
            # 构造矩阵
            Y_matrix = df.pivot(index='unit', columns='time', values='outcome').values
            
            units = df['unit'].unique()
            treat_idx = np.where(units == request.treat_unit)[0][0]
            times = df['time'].unique()
            treat_time_idx = np.where(times == request.treat_period)[0][0]
            
            # 拟合
            model = MatrixCompletionSCM(rank=rank, lambda_reg=lambda_reg)
            model.fit(Y_matrix, treated_idx=treat_idx, treat_start=treat_time_idx)
        
        elif request.method == "factor":
            n_factors = request.hyperparams.get("n_factors", 3)
            
            Y_matrix = df.pivot(index='unit', columns='time', values='outcome').values
            
            units = df['unit'].unique()
            treat_idx = np.where(units == request.treat_unit)[0][0]
            times = df['time'].unique()
            treat_time_idx = np.where(times == request.treat_period)[0][0]
            
            model = FactorModelSCM(n_factors=n_factors)
            model.fit(Y_matrix, treated_idx=treat_idx, treat_start=treat_time_idx)
        
        # 保存
        model_store.save_model(model, model_id)
        
        # 更新任务状态
        redis_client.setex(f"task_{model_id}", 3600, json.dumps({
            "status": "completed",
            "model_id": model_id,
            "completed_at": datetime.now().isoformat()
        }))
        
    except Exception as e:
        redis_client.setex(f"task_{model_id}", 3600, json.dumps({
            "status": "failed",
            "error": str(e)
        }))

@app.get("/task/{task_id}")
async def get_task_status(task_id: str):
    """查询任务状态"""
    status = redis_client.get(task_id)
    if not status:
        raise HTTPException(status_code=404, detail="Task not found")
    
    return json.loads(status)

@app.post("/model/predict")
async def predict_counterfactual(request: PredictRequest):
    """预测反事实结果"""
    model = model_store.load_model(request.model_id)
    
    if isinstance(model, (MatrixCompletionSCM, FactorModelSCM)):
        Y_cf = model.predict_counterfactual()
    else:
        raise HTTPException(status_code=400, detail="Unsupported model type")
    
    result = {
        "model_id": request.model_id,
        "counterfactual": Y_cf.tolist(),
        "periods": request.post_periods
    }
    
    if request.confidence_interval:
        # 通过Bootstrap计算CI
        ci_lower, ci_upper = _bootstrap_ci(model, Y_cf, n_bootstraps=200)
        result['ci_lower'] = ci_lower.tolist()
        result['ci_upper'] = ci_upper.tolist()
    
    return result

def _bootstrap_ci(model, Y_cf_point, n_bootstraps: int = 200):
    """Bootstrap置信区间"""
    bootstrap_samples = []
    
    for b in range(n_bootstraps):
        # 对因子或低秩成分进行Bootstrap
        if hasattr(model, 'F'):
            # 因子Bootstrap（时间 block bootstrap）
            F_boot = _block_bootstrap(model.F, block_size=5)
            lambda_treated = model.lambda_loadings[model.treated_idx]
            Y_cf_boot = lambda_treated @ F_boot[model.treat_start:, :].T
        elif hasattr(model, 'L_est'):
            # 残差Bootstrap
            residuals = model.error_component[model.treated_idx, model.treat_start:]
            resid_boot = np.random.choice(residuals, size=len(residuals))
            Y_cf_boot = Y_cf_point + resid_boot
        
        bootstrap_samples.append(Y_cf_boot)
    
    bootstrap_samples = np.array(bootstrap_samples)
    ci_lower = np.percentile(bootstrap_samples, 2.5, axis=0)
    ci_upper = np.percentile(bootstrap_samples, 97.5, axis=0)
    
    return ci_lower, ci_upper

def _block_bootstrap(F: np.ndarray, block_size: int):
    """块Bootstrap保留时间依赖性"""
    T, r = F.shape
    n_blocks = T // block_size
    
    # 随机打乱块
    block_indices = np.random.choice(n_blocks, n_blocks, replace=True)
    
    F_boot = []
    for idx in block_indices:
        start = idx * block_size
        end = min((idx + 1) * block_size, T)
        F_boot.append(F[start:end, :])
    
    return np.vstack(F_boot)

@app.post("/model/whatif")
async def what_if_analysis(request: WhatIfRequest):
    """What-if分析：模拟不同干预情景"""
    model = model_store.load_model(request.model_id)
    
    if not isinstance(model, FactorModelSCM):
        raise HTTPException(status_code=400, detail="What-if only supported for factor models")
    
    # 修改干预时间点
    original_treat_start = model.treat_start
    model.treat_start = request.alternative_treat_period
    
    # 预测
    Y_cf_alt = model.predict_counterfactual()
    
    # 效应缩放
    Y_cf_alt *= request.magnitude_multiplier
    
    # 恢复原模型
    model.treat_start = original_treat_start
    
    return {
        "model_id": request.model_id,
        "alternative_cf": Y_cf_alt.tolist(),
        "assumptions": "Factor structure unchanged; treatment effect is multiplicative"
    }

@app.get("/health")
async def health_check():
    """健康检查"""
    redis_status = redis_client.ping()
    return {
        "status": "healthy" if redis_status else "degraded",
        "redis_connected": redis_status,
        "active_models": len(redis_client.keys("scm:model:*"))
    }

# Docker化配置
"""
FROM python:3.9-slim

WORKDIR /app

COPY requirements.txt .
RUN pip install --no-cache-dir -r requirements.txt

COPY . .

EXPOSE 8000

CMD ["uvicorn", "api_server:app", 
     "--host", "0.0.0.0", 
     "--port", "8000",
     "--workers", "4"]
"""

# 部署脚本
"""
docker build -t scm-api:v1.0 .
docker run -d -p 8000:8000 \
  -e REDIS_HOST=redis-cluster.internal \
  scm-api:v1.0
"""

# Kubernetes部署
"""
apiVersion: apps/v1
kind: Deployment
metadata:
  name: scm-api
spec:
  replicas: 3
  selector:
    matchLabels:
      app: scm-api
  template:
    metadata:
      labels:
        app: scm-api
    spec:
      containers:
      - name: api
        image: scm-api:v1.0
        env:
        - name: REDIS_HOST
          valueFrom:
            configMapKeyRef:
              name: redis-config
              key: host
        resources:
          requests:
            memory: "4Gi"
            cpu: "1000m"
          limits:
            memory: "8Gi"
            cpu: "2000m"
"""

if __name__ == "__main__":
    import uvicorn
    uvicorn.run(app, host="0.0.0.0", port=8000)