Python算法
面板数据模型原理与Python实现
本文详解面板数据模型（Panel Data Model），涵盖定义、特征及与传统单维度模型的对比。深入推导固定效应（FE）与随机效应（RE）的数学原理，包括组内离差变换与广义最小二乘法（GLS）。提供统计检验方法（F 检验、Hausman 检验）指导模型选择，并给出基于 Python linearmodels 库的完整代码实现，包含数据模拟、预处理、模型估计及可视化全流程，适用于数学建模与计量经济学分析。
本文详解面板数据模型（Panel Data Model），涵盖定义、特征及与传统单维度模型的对比。深入推导固定效应（FE）与随机效应（RE）的数学原理，包括组内离差变换与广义最小二乘法（GLS）。提供统计检验方法（F 检验、Hausman 检验）指导模型选择，并给出基于 Python linearmodels 库的完整代码实现，包含数据模拟、预处理、模型估计及可视化全流程，适用于数学建模与计量经济学分析。
面板数据（Panel Data）是同时包含横截面维度（Cross-Sectional Dimension）和时间维度（Time Dimension）的二维结构化数据，其数学表示为 ${y_{it}, X_{it}}_{i=1,t=1}^{N,T}$，其中：
面板数据的核心特征是'跟踪同一批个体的动态变化'，例如：
传统计量模型仅利用单一维度的信息，无法处理未观测异质性（Unobserved Heterogeneity）——即个体层面或时间层面无法直接测量的特征/冲击，最终导致估计结果有偏且不一致。具体局限如下：
仅使用同一时间点的个体数据（如 2023 年 1000 个家庭的收入与消费），模型形式为： $$ y_i = \alpha + X_i\beta + \epsilon_i $$ 其中 $i$ 为个体，无时间维度。
局限：忽略了个体未观测异质性（如家庭理财能力、遗传禀赋）。若未观测异质性与解释变量相关（如理财能力强的家庭收入更高），则违反 OLS'解释变量与误差项不相关'的假设，导致遗漏变量偏差（Omitted Variable Bias）。数学上，若真实模型为 $y_i = \alpha + X_i\beta + \alpha_i + \epsilon_i$（$\alpha_i$ 为未观测异质性），且 $Cov(X_i, \alpha_i) \neq 0$，则横截面回归的估计量为： $$ \hat{\beta}^{CS} = \beta + (X'X)^{-1}X'\alpha $$ 由于 $Cov(X_i, \alpha_i) \neq 0$，$\hat{\beta}^{CS}$ 是有偏且不一致的。
仅使用单个个体的时间序列数据（如某企业 2010-2020 年的投资与利润），模型形式为： $$ y_t = \alpha + X_t\beta + \epsilon_t $$ 其中 $t$ 为时间，无个体维度。
局限：
面板数据通过结合横截面与时间维度的双重信息，从根本上解决了单维度模型的局限，具体优势包括：
面板数据模型的核心目标是将被解释变量的总变异分解为三个部分： $$ \text{总变异} = \text{可观测因素变异} + \text{未观测异质性变异} + \text{随机误差变异} $$
数学上，双向固定效应模型（最常用的面板模型）的变异分解为： $$ y_{it} = \underbrace{\alpha_i}{\text{个体固定效应}} + \underbrace{\lambda_t}{\text{时间固定效应}} + \underbrace{X_{it}\beta}{\text{可观测因素}} + \underbrace{\epsilon{it}}_{\text{随机误差}} $$
各部分的含义如下：
模型的关键是如何处理未观测异质性。根据 $\alpha_i$（或 $\lambda_t$）的性质（固定参数 vs 随机变量），面板数据模型分为两类：
面板数据模型的基础分类需同时考虑个体异质性与时间异质性，具体如下：
| 模型类型 | 假设条件 | 模型形式 |
|---|---|---|
| 混合回归（POLS） | 无个体异质性（$\alpha_1=\alpha_2=...=\alpha_N$）、无时间异质性（$\lambda_1=...=\lambda_T$） | $y_{it} = \alpha + X_{it}\beta + \epsilon_{it}$ |
| 单向个体固定效应（FE） | 个体异质性为固定参数，与 $X_{it}$ 相关；无时间异质性 | $y_{it} = \alpha_i + X_{it}\beta + \epsilon_{it}$ |
| 单向时间固定效应（TE） | 时间异质性为固定参数，与 $X_{it}$ 相关；无个体异质性 | $y_{it} = \lambda_t + X_{it}\beta + \epsilon_{it}$ |
| 双向固定效应（Two-Way FE） | 个体异质性与时间异质性均为固定参数，且至少一个与 $X_{it}$ 相关 | $y_{it} = \alpha_i + \lambda_t + X_{it}\beta + \epsilon_{it}$ |
| 随机效应（RE） | 个体异质性为随机变量（$\alpha_i \sim iid(0,\sigma_\alpha^2)$），与 $X_{it}$ 不相关；无时间异质性 | $y_{it} = \alpha + X_{it}\beta + (\alpha_i + \epsilon_{it})$ |
固定效应模型的核心是消除与解释变量相关的未观测异质性，最常用的方法是组内离差变换（Within Transformation）。以下以单向个体固定效应模型为例，详细推导其原理。
单向个体固定效应模型的设定为： $$ y_{it} = \alpha_i + X_{it}\beta + \epsilon_{it}, \quad i=1,...,N; \ t=1,...,T $$
关键假设：
由于 $\alpha_i$ 是不随时间变化的固定参数，我们可以通过个体时间均值来消除它。具体步骤如下：
对每个个体 $i$，计算被解释变量与解释变量的时间均值（即该个体在所有时间点的平均值）： $$ \bar{y}i = \frac{1}{T} \sum{t=1}^T y_{it}, \quad \bar{X}i = \frac{1}{T} \sum{t=1}^T X_{it}, \quad \bar{\epsilon}i = \frac{1}{T} \sum{t=1}^T \epsilon_{it} $$
将原模型对时间 $t$ 求和并除以 $T$，得到均值模型： $$ \bar{y}_i = \alpha_i + \bar{X}_i\beta + \bar{\epsilon}_i $$
将原模型与均值模型相减，得到离差模型（Within Model）： $$ y_{it} - \bar{y}i = (X{it} - \bar{X}i)\beta + (\epsilon{it} - \bar{\epsilon}_i) $$
记离差变量为： $$ \tilde{y}{it} = y{it} - \bar{y}i, \quad \tilde{X}{it} = X_{it} - \bar{X}i, \quad \tilde{\epsilon}{it} = \epsilon_{it} - \bar{\epsilon}_i $$
则离差模型简化为： $$ \tilde{y}{it} = \tilde{X}{it}\beta + \tilde{\epsilon}_{it} $$
离差变换完全消除了个体固定效应 $\alpha_i$——因为 $\alpha_i$ 不随时间变化，其时间均值仍为 $\alpha_i$，故 $\alpha_i - \alpha_i = 0$。
更重要的是，离差变换仅保留了解释变量的'个体内变异'（即同一个体随时间的变化），而剔除了'个体间变异'（即不同个体的差异）。这意味着固定效应模型的估计结果仅依赖于个体内部的动态变化，而非个体之间的静态差异——这正是固定效应模型解决'未观测异质性与解释变量相关'问题的关键！
离差模型满足 OLS 的基本假设（解释变量与误差项不相关），因此我们可以用普通最小二乘法（OLS）估计离差模型，得到固定效应估计量：
$$ \hat{\beta}^{FE} = \left( \sum_{i=1}^N \sum_{t=1}^T \tilde{X}{it}'\tilde{X}{it} \right)^{-1} \sum_{i=1}^N \sum_{t=1}^T \tilde{X}{it}'\tilde{y}{it} $$
在得到 $\hat{\beta}^{FE}$ 后，个体固定效应 $\alpha_i$ 的估计值为： $$ \hat{\alpha}_i = \bar{y}_i - \bar{X}_i\hat{\beta}^{FE} $$
其含义是：个体 $i$ 的未观测特征对被解释变量的平均影响（控制了解释变量的个体均值后）。
若同时存在个体固定效应与时间固定效应，则模型为双向固定效应模型： $$ y_{it} = \alpha_i + \lambda_t + X_{it}\beta + \epsilon_{it} $$
此时需使用双重离差变换（同时消除 $\alpha_i$ 与 $\lambda_t$），具体步骤如下：
首先计算所有观测值的总均值（即所有个体、所有时间的平均值）： $$ \bar{y} = \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T y_{it}, \quad \bar{X} = \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T X_{it} $$
将原模型减去个体时间均值与时间个体均值，再加上总均值（避免重复消除）： $$ y_{it} - \bar{y}_i - \bar{y}t + \bar{y} = (X{it} - \bar{X}_i - \bar{X}t + \bar{X})\beta + (\epsilon{it} - \bar{\epsilon}_i - \bar{\epsilon}_t + \bar{\epsilon}) $$
其中 $\bar{y}t = \frac{1}{N} \sum{i=1}^N y_{it}$（时间 $t$ 的个体均值），$\bar{X}t = \frac{1}{N} \sum{i=1}^N X_{it}$（时间 $t$ 的解释变量均值）。
记双重离差变量为： $$ \ddot{y}{it} = y{it} - \bar{y}i - \bar{y}t + \bar{y}, \quad \ddot{X}{it} = X{it} - \bar{X}_i - \bar{X}_t + \bar{X} $$
则双向固定效应模型的离差形式为： $$ \ddot{y}{it} = \ddot{X}{it}\beta + \ddot{\epsilon}_{it} $$
同样使用 OLS 估计该模型，得到双向固定效应估计量 $\hat{\beta}^{Two-Way FE}$。
固定效应估计量的一致性（Consistency）要求：当 $N \to \infty$（或 $T \to \infty$）时，$\hat{\beta}^{FE} \to \beta$（真实参数）。
一致性的核心前提是离差变量与离差误差项不相关，即： $$ plim_{N\to\infty} \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T \tilde{X}{it}'\tilde{\epsilon}{it} = 0 $$
根据离差变量的定义，$\tilde{X}{it} = X{it} - \bar{X}i$，$\tilde{\epsilon}{it} = \epsilon_{it} - \bar{\epsilon}i$。因此： $$ \sum{i=1}^N \sum_{t=1}^T \tilde{X}{it}'\tilde{\epsilon}{it} = \sum_{i=1}^N \sum_{t=1}^T (X_{it} - \bar{X}i)'(\epsilon{it} - \bar{\epsilon}_i) $$
展开后： $$ = \sum_{i=1}^N \left[ \sum_{t=1}^T X_{it}'\epsilon_{it} - \bar{X}i'\sum{t=1}^T \epsilon_{it} - \sum_{t=1}^T X_{it}'\bar{\epsilon}_i + T\bar{X}_i'\bar{\epsilon}_i \right] $$
由于 $\bar{X}i = \frac{1}{T}\sum{t=1}^T X_{it}$，$\bar{\epsilon}i = \frac{1}{T}\sum{t=1}^T \epsilon_{it}$，因此： $$ \bar{X}i'\sum{t=1}^T \epsilon_{it} = T\bar{X}i'\bar{\epsilon}i, \quad \sum{t=1}^T X{it}'\bar{\epsilon}_i = T\bar{X}_i'\bar{\epsilon}_i $$
代入后，中间两项抵消，最终简化为： $$ \sum_{i=1}^N \sum_{t=1}^T \tilde{X}{it}'\tilde{\epsilon}{it} = \sum_{i=1}^N \sum_{t=1}^T X_{it}'\epsilon_{it} - T\sum_{i=1}^N \bar{X}_i'\bar{\epsilon}_i $$
根据严格外生性假设（$E(\epsilon_{it} | X_{i1},..., X_{iT}) = 0$），有： $$ E\left( \sum_{t=1}^T X_{it}'\epsilon_{it} \right) = 0, \quad E\left( T\bar{X}_i'\bar{\epsilon}_i \right) = 0 $$
因此，当 $N \to \infty$ 时，样本均值收敛于总体均值： $$ plim_{N\to\infty} \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T \tilde{X}{it}'\tilde{\epsilon}{it} = 0 $$
即固定效应估计量是一致的。
随机效应模型的核心是将个体异质性视为随机变量，并利用广义最小二乘法（GLS）处理复合误差项的相关性。
随机效应模型的设定为： $$ y_{it} = \alpha + X_{it}\beta + u_{it} $$
其中复合误差项为： $$ u_{it} = \alpha_i + \epsilon_{it} $$
关键假设：
由于 $\alpha_i$ 是个体共有的随机项，同一个体的不同时间观测值的误差项存在序列相关性。具体来说：
复合误差项的方差为： $$ Var(u_{it}) = Var(\alpha_i + \epsilon_{it}) = Var(\alpha_i) + Var(\epsilon_{it}) = \sigma_\alpha^2 + \sigma_\epsilon^2 $$
同一个体 $i$、不同时间 $t$ 与 $s$ 的协方差为： $$ Cov(u_{it}, u_{is}) = Cov(\alpha_i + \epsilon_{it}, \alpha_i + \epsilon_{is}) = Var(\alpha_i) + Cov(\alpha_i, \epsilon_{is}) + Cov(\epsilon_{it}, \alpha_i) + Cov(\epsilon_{it}, \epsilon_{is}) $$
根据假设，$\alpha_i$ 与 $\epsilon_{it}$ 独立，$\epsilon_{it}$ 与 $\epsilon_{is}$ 独立，因此： $$ Cov(u_{it}, u_{is}) = \sigma_\alpha^2 $$
对于个体 $i$，其 $T$ 个观测值的误差项向量为 $u_i = (u_{i1}, u_{i2},..., u_{iT})'$，则其方差 - 协方差矩阵为： $$ \Omega_i = E(u_i u_i') = \sigma_\alpha^2 \iota_T \iota_T' + \sigma_\epsilon^2 I_T $$
其中 $\iota_T$ 是 $T$ 维全 1 向量（$\iota_T = (1,1,...,1)'$），$I_T$ 是 $T$ 维单位矩阵。
同一个体内不同时间观测值的误差相关性用**个体内相关系数（Intraclass Correlation）**衡量： $$ \rho = \frac{Cov(u_{it}, u_{is})}{Var(u_{it})} = \frac{\sigma_\alpha^2}{\sigma_\alpha^2 + \sigma_\epsilon^2} $$
$\rho$ 的取值范围为 $[0,1]$：
由于复合误差项存在序列相关性，OLS 估计虽然无偏，但不有效（即标准误偏大，无法利用相关性信息）。因此需使用广义最小二乘法（GLS）——通过变换矩阵将复合误差转换为同方差且不相关的误差项。
GLS 的核心是找到变换矩阵 $P$，使得 $P u_i$ 的方差 - 协方差矩阵为单位矩阵 $I_T$。根据方差 - 协方差矩阵 $\Omega_i$ 的结构，我们可以构造如下变换矩阵： $$ P = I_T - \theta \cdot \frac{1}{T}\iota_T \iota_T' $$
其中 $\theta$ 为变换参数，计算公式为： $$ \theta = 1 - \frac{\sigma_\epsilon}{\sqrt{\sigma_\epsilon^2 + T\sigma_\alpha^2}} $$
我们希望变换后的误差项 $P u_i$ 的方差为 1，即： $$ Var(P u_{it}) = 1 $$
代入 $P$ 的定义： $$ Var(u_{it} - \theta \bar{u}i) = Var(u{it}) + \theta^2 Var(\bar{u}i) - 2\theta Cov(u{it}, \bar{u}_i) = 1 $$
其中 $\bar{u}i = \frac{1}{T}\sum{t=1}^T u_{it}$，且：
将上述结果代入方差公式： $$ (\sigma_\alpha^2 + \sigma_\epsilon^2) + \theta^2 \cdot \frac{\sigma_\alpha^2 + \sigma_\epsilon^2}{T} - 2\theta \cdot \frac{\sigma_\alpha^2 + \sigma_\epsilon^2}{T} = 1 $$
提取公因子 $(\sigma_\alpha^2 + \sigma_\epsilon^2)$： $$ (\sigma_\alpha^2 + \sigma_\epsilon^2) \left[ 1 + \frac{\theta^2 - 2\theta}{T} \right] = 1 $$
为了简化，令 $\theta = 1 - \frac{\sigma_\epsilon}{\sqrt{\sigma_\epsilon^2 + T\sigma_\alpha^2}}$，则： $$ 1 + \frac{\theta^2 - 2\theta}{T} = 1 - \frac{(1 - \theta)^2}{T} = \frac{\sigma_\epsilon^2}{\sigma_\epsilon^2 + T\sigma_\alpha^2} $$
代入后： $$ (\sigma_\alpha^2 + \sigma_\epsilon^2) \cdot \frac{\sigma_\epsilon^2}{\sigma_\epsilon^2 + T\sigma_\alpha^2} = 1 $$
显然成立，因此变换参数 $\theta$ 的形式正确。
使用变换矩阵 $P$ 对原模型进行 GLS 变换，得到： $$ P y_i = P \alpha \iota_T + P X_i \beta + P u_i $$
其中 $y_i = (y_{i1},..., y_{iT})'$，$X_i = (X_{i1},..., X_{iT})'$。
由于 $P \iota_T = \iota_T - \theta \cdot \frac{1}{T}\iota_T \iota_T' \iota_T = \iota_T - \theta \iota_T = (1 - \theta)\iota_T$，因此变换后的模型可写为： $$ y_{it}^* = \alpha^* + X_{it}^\beta + \epsilon_{it}^ $$
其中：
对变换后的模型使用 OLS，得到随机效应估计量： $$ \hat{\beta}^{RE} = \left( \sum_{i=1}^N \sum_{t=1}^T (X_{it} - \theta \bar{X}i)'(X{it} - \theta \bar{X}i) \right)^{-1} \sum{i=1}^N \sum_{t=1}^T (X_{it} - \theta \bar{X}i)'(y{it} - \theta \bar{y}_i) $$
变换参数 $\theta$ 依赖于 $\sigma_\alpha^2$ 与 $\sigma_\epsilon^2$，而这两个方差是未知的，因此需用方差组分估计（Variance Components Estimation）从样本中估计，具体步骤如下：
组内方差是固定效应模型的残差方差，计算公式为： $$ \hat{\sigma}\epsilon^2 = \frac{1}{NT - N - K} \sum{i=1}^N \sum_{t=1}^T (\tilde{y}{it} - \tilde{X}{it}\hat{\beta}^{FE})^2 $$
其中：
组间方差是混合回归模型的个体均值残差方差减去组内方差的均值，计算公式为： $$ \hat{\sigma}\alpha^2 = \frac{1}{N - K - 1} \sum{i=1}^N (\bar{y}_i - \bar{X}i\hat{\beta}^{POLS})^2 - \frac{\hat{\sigma}\epsilon^2}{T} $$
其中：
将估计的 $\hat{\sigma}\alpha^2$ 与 $\hat{\sigma}\epsilon^2$ 代入变换参数公式，得到： $$ \hat{\theta} = 1 - \frac{\hat{\sigma}\epsilon}{\sqrt{\hat{\sigma}\epsilon^2 + T\hat{\sigma}_\alpha^2}} $$
此时的 GLS 称为可行广义最小二乘法（FGLS）——因为我们用样本估计的方差组分代替了真实方差。
原假设 $H_0$：所有个体固定效应相同（$\alpha_1=\alpha_2=...=\alpha_N=\alpha$），即混合回归模型成立；备择假设 $H_1$：至少有一个个体固定效应不同（$\exists i \neq j, \alpha_i \neq \alpha_j$），即固定效应模型成立。
F 检验的核心是比较混合回归与固定效应模型的拟合优度——若固定效应模型的拟合优度显著高于混合回归，则拒绝原假设。
F 统计量的公式为： $$ F = \frac{(RSSR - RSSFE)/ (N - 1)}{RSSFE / (NT - N - K)} $$
其中：
因为固定效应模型有 $N$ 个个体固定效应，而混合回归模型只有 1 个截距项，因此约束数为 $N-1$（即要求 $N$ 个固定效应等于 1 个共同值）。
若 $F > F_{\alpha}(N-1, NT-N-K)$（$\alpha$ 为显著性水平，通常取 0.05），则拒绝原假设，选择固定效应模型；否则选择混合回归模型。
原假设 $H_0$：$Cov(X_{it}, \alpha_i) = 0$（随机效应模型的核心假设成立，即个体异质性与解释变量不相关）；备择假设 $H_1$：$Cov(X_{it}, \alpha_i) \neq 0$（固定效应模型更合适）。
Hausman 检验的本质是比较固定效应（FE）与随机效应（RE）估计量的差异：
因此，若两个估计量的差异显著不为零，则拒绝原假设，选择固定效应模型。
Hausman 统计量的公式为： $$ H = (\hat{\beta}^{FE} - \hat{\beta}^{RE})' \left( Var(\hat{\beta}^{FE}) - Var(\hat{\beta}^{RE}) \right)^{-1} (\hat{\beta}^{FE} - \hat{\beta}^{RE}) $$
令 $\hat{\delta} = \hat{\beta}^{FE} - \hat{\beta}^{RE}$，则 Hausman 统计量可写为： $$ H = \hat{\delta}' Var(\hat{\delta})^{-1} \hat{\delta} $$
其中 $Var(\hat{\delta}) = Var(\hat{\beta}^{FE}) - Var(\hat{\beta}^{RE})$——因为 FE 与 RE 估计量在 $H_0$ 成立时是渐近独立的（证明略）。
当 $N \to \infty$ 时，Hausman 统计量服从卡方分布： $$ H \sim \chi^2(K) $$
其中 $K$ 为解释变量的个数（即估计量 $\hat{\beta}$ 的维度）。
若 $H > \chi^2_{\alpha}(K)$（$\alpha$ 为显著性水平，通常取 0.05），则拒绝原假设，选择固定效应模型；否则选择随机效应模型。
面板数据的预处理是模型估计的基础，直接影响结果的可靠性。具体步骤如下：
计算变量的个体内变异与个体间变异，判断固定效应模型是否适用（若个体内变异很小，则固定效应模型可能无法估计出显著结果）：
根据研究问题与数据特征，选择合适的模型并进行初步估计：
通过统计检验选择最优模型，具体步骤如下：
模型估计后，需进行诊断检验（判断模型假设是否成立）与稳健性检验（判断结果是否可靠）。
聚类稳健标准误的计算公式为： $$ Var(\hat{\beta}^{FE}){cluster} = \left( \sum{i=1}^N \tilde{X}_i'\tilde{X}i \right)^{-1} \left( \sum{i=1}^N \tilde{X}_i'\tilde{\epsilon}_i \tilde{\epsilon}_i'\tilde{X}i \right) \left( \sum{i=1}^N \tilde{X}_i'\tilde{X}_i \right)^{-1} $$
其中 $\tilde{X}i = (\tilde{X}{i1},..., \tilde{X}{iT})'$，$\tilde{\epsilon}i = (\tilde{\epsilon}{i1},..., \tilde{\epsilon}{iT})'$。
稳健性检验是判断结果是否可靠的关键，常用方法如下：
模型估计与检验完成后，需清晰、准确地解释结果，并给出政策建议。
| 模型类型 | 适用场景 | 局限性 |
|---|---|---|
| 混合回归（POLS） | 未观测异质性与解释变量完全不相关 | 忽略异质性会导致严重偏差 |
| 固定效应（FE） | 未观测异质性与解释变量相关（如家庭能力与收入） | 无法估计不随时间变化的变量；依赖个体内变异 |
| 随机效应（RE） | 未观测异质性与解释变量不相关（如企业地理位置与广告支出） | 若异质性与解释变量相关，估计量有偏；依赖随机抽样假设 |
| 双向固定效应（Two-Way FE） | 同时存在个体与时间未观测异质性，且至少一个与解释变量相关 | 无法估计不随时间/个体变化的变量；计算复杂度高 |
面板数据的维度（$N$ 与 $T$ 的大小）直接影响模型的选择与估计结果：
内生性是面板数据模型中最常见的问题之一，主要来源包括双向因果、测量误差与遗漏变量。以下是常用的处理方法：
面板数据模型是处理多维度数据的核心工具，其优势在于：
模型选择的核心依据是未观测异质性与解释变量的相关性：
在实际应用中，需严格遵循数据预处理→模型估计→统计检验→稳健性检验→结果解释的流程，确保结果的可靠性与政策含义的准确性。
参考文献（供深入学习）：
模拟企业面板数据（2015-2019 年，200 家企业），研究研发投入对企业利润的影响，使用双向固定效应模型，包含个体（企业）固定效应和时间固定效应。
import numpy as np
import pandas as pd
import statsmodels.api as sm
from linearmodels.panel import PanelOLS
from linearmodels.panel import PooledOLS
from linearmodels.iv import IV2SLS
import matplotlib.pyplot as plt
# 设置 matplotlib 字体以支持中文显示
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
# 1. 数据模拟函数
def simulate_panel_data(n_firms=200, n_years=5, seed=42):
"""
模拟企业面板数据
:param n_firms: 企业数量
:param n_years: 时间跨度（年数）
:param seed: 随机种子，保证结果可复现
:return: 面板数据 DataFrame
"""
# 设置随机种子
np.random.seed(seed)
# 企业 ID 和年份
firms = np.repeat(range(n_firms), n_years)
years = np.tile(range(2015, 2015 + n_years), n_firms)
# 个体固定效应（企业固有特征）
firm_effects = np.repeat(np.random.normal(10, 2, n_firms), n_years)
# 时间固定效应（宏观冲击）
time_effects = np.tile(np.random.normal(5, 1, n_years), n_firms)
# 研发投入（解释变量，随时间和企业变化）
R_D = np.random.normal(20, 5, n_firms * n_years) + firm_effects * 0.2 + time_effects * 0.1
# 控制变量：资产规模
assets = np.random.normal(100, 20, n_firms * n_years) + firm_effects * 0.5
# 随机误差项
epsilon = np.random.normal(0, 3, n_firms * n_years)
# 被解释变量：利润（真实关系：利润 = 0.8*研发投入 + 0.2*资产规模 + 个体效应 + 时间效应 + 随机误差）
profit = 0.8 * R_D + 0.2 * assets + firm_effects + time_effects + epsilon
# 构建 DataFrame
data = pd.DataFrame({
'firm_id': firms,
'year': years,
'R_D': R_D,
'assets': assets,
'profit': profit
})
# 设置面板数据索引（个体，时间）
data = data.set_index(['firm_id', 'year'])
return data
# 2. 面板数据预处理函数
def preprocess_panel_data(data):
"""
面板数据预处理（缺失值、异常值处理）
:param data: 原始面板数据
:return: 预处理后面板数据
"""
# 检查缺失值并删除
data = data.dropna()
# 异常值处理：缩尾处理（前后 1%）
for col in ['R_D', 'assets', 'profit']:
q1 = data[col].quantile(0.01)
q99 = data[col].quantile(0.99)
data[col] = np.clip(data[col], q1, q99)
return data
# 3. 模型估计与检验函数
def estimate_panel_models(data):
"""
估计面板数据模型并进行检验
:param data: 预处理后面板数据
:return: 各模型估计结果
"""
# 准备解释变量和被解释变量
y = data['profit']
X = data[['R_D', 'assets']]
# 3.1 混合回归模型（POLS）
pooled = PooledOLS(y, sm.add_constant(X)).fit()
# 3.2 单向个体固定效应模型（FE）
fe_individual = PanelOLS(y, X, entity_effects=True).fit()
# 3.3 双向固定效应模型（Two-Way FE）
two_way_fe = PanelOLS(y, X, entity_effects=True, time_effects=True).fit()
# 3.4 随机效应模型（RE）
# 注意：Linearmodels 库中 RE 需要使用 RandomEffects，但需先安装最新版本
try:
from linearmodels.panel import RandomEffects
re_model = RandomEffects(y, sm.add_constant(X)).fit()
except ImportError:
re_model = None
return pooled, fe_individual, two_way_fe, re_model
# 4. Hausman 检验函数
def hausman_test(fe_results, re_results):
"""
执行 Hausman 检验（固定效应 vs 随机效应）
:param fe_results: 固定效应模型结果
:param re_results: 随机效应模型结果
:return: Hausman 检验统计量和 p 值
"""
if re_results is None:
return "RandomEffects 模型未安装"
# 选取公共变量进行比较（去除 RE 中的常数项等 FE 中没有的变量）
common_params = list(set(fe_results.params.index) & set(re_results.params.index))
if not common_params:
return "没有公共参数，无法进行 Hausman 检验"
fe_beta = fe_results.params[common_params]
re_beta = re_results.params[common_params]
fe_cov = fe_results.cov.loc[common_params, common_params]
re_cov = re_results.cov.loc[common_params, common_params]
diff = fe_beta - re_beta
var_diff = fe_cov - re_cov
# 计算 Hausman 统计量
# 注意：如果 var_diff 不可逆，可能需要使用广义逆，这里使用 pinv 以防万一
hausman_stat = diff.T @ np.linalg.pinv(var_diff) @ diff
# 计算 p 值（卡方分布，自由度=解释变量个数）
from scipy.stats import chi2
df = len(diff)
p_value = 1 - chi2.cdf(hausman_stat, df)
return f"Hausman 统计量：{hausman_stat:.4f}, p 值：{p_value:.4f}"
# 主程序
if __name__ == "__main__":
# 生成模拟数据
panel_data = simulate_panel_data()
# 数据预处理
processed_data = preprocess_panel_data(panel_data)
# 模型估计
pooled_result, fe_individual_result, two_way_fe_result, re_result = estimate_panel_models(processed_data)
# 输出模型结果
print("混合回归模型结果：")
print(pooled_result)
print("\n单向个体固定效应模型结果：")
print(fe_individual_result)
print("\n双向固定效应模型结果：")
print(two_way_fe_result)
# Hausman 检验
print("\nHausman 检验结果：")
hausman_result = hausman_test(fe_individual_result, re_result)
print(hausman_result)
# 可视化结果
plt.figure(figsize=(12, 8))
# 绘制研发投入与利润的关系
plt.scatter(panel_data['R_D'], panel_data['profit'], alpha=0.3, label='原始数据')
# 绘制双向固定效应模型的拟合线（仅展示核心解释变量 R_D 的影响）
# 注意：固定效应模型的拟合需控制个体和时间效应，此处简化为展示系数方向
plt.axline((20, two_way_fe_result.params['assets']*100 + two_way_fe_result.params['R_D']*20),
slope=two_way_fe_result.params['R_D'],
color='red',
label=f'双向 FE 拟合线（斜率：{two_way_fe_result.params["R_D"]:.2f}）')
plt.xlabel('研发投入（R_D）')
plt.ylabel('企业利润')
plt.title('研发投入与企业利润的关系')
plt.legend()
plt.grid(True)
# 保存图像
plt.savefig('R_D_profit_relation.png', dpi=300)
plt.show()
这份代码是数学建模比赛中面板数据双向固定效应模型的完整研究流程实现，涵盖「模拟数据→预处理→多模型估计→Hausman 检验→可视化」全环节，核心目标是验证研发投入对企业利润的影响。接下来从数学原理、代码逻辑、建模细节三方面逐块解析：
import numpy as np
import pandas as pd
import statsmodels.api as sm
from linearmodels.panel import PanelOLS, PooledOLS, RandomEffects
import matplotlib.pyplot as plt
# 中文显示设置
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
numpy/pandas：用于数据模拟与结构化处理（面板数据的核心载体是 pd.DataFrame）。linearmodels：专门为面板数据设计的库，支持直接指定个体/时间固定效应（无需手动生成哑变量），比 statsmodels 更高效。statsmodels：仅用于给混合回归模型添加常数项。simulate_panel_data模拟的是平衡面板数据（200 企业×5 年=1000 条观测），真实关系符合双向固定效应模型： $$ \text{profit}{it} = 0.8 \cdot R{D,it} + 0.2 \cdot \text{assets}{it} + \alpha_i + \lambda_t + \varepsilon{it} $$ 其中：
def simulate_panel_data(n_firms=200, n_years=5, seed=42):
np.random.seed(seed) # 保证模拟结果可复现（比赛必填）
# 1. 生成面板维度：[firm_id, year] 双索引（linearmodels 面板数据的必要格式）
firms = np.repeat(range(n_firms), n_years)
years = np.tile(range(2015, 2015 + n_years), n_firms)
# 2. 个体固定效应 α_i：每个企业 1 个值，正态分布 N(10, 2²)
firm_effects = np.repeat(np.random.normal(10, 2, n_firms), n_years)
# 3. 时间固定效应 λ_t：每年 1 个值，正态分布 N(5, 1²)
time_effects = np.tile(np.random.normal(5, 1, n_years), n_firms)
# 4. 解释变量：研发投入 R_D（与 α_i、λ_t 相关，更接近现实）
R_D = np.random.normal(20, 5, n_firms * n_years) + firm_effects * 0.2 + time_effects * 0.1
# 5. 控制变量：资产规模（与个体固定效应 α_i 相关，合理假设）
assets = np.random.normal(100, 20, n_firms * n_years) + firm_effects * 0.5
# 6. 随机误差项
epsilon = np.random.normal(0, 3, n_firms * n_years)
# 7. 被解释变量：利润（严格遵循真实双向 FE 模型）
profit = 0.8 * R_D + 0.2 * assets + firm_effects + time_effects + epsilon
# 8. 构建面板 DataFrame 并设置双索引
data = pd.DataFrame({
'firm_id': firms,
'year': years,
'R_D': R_D,
'assets': assets,
'profit': profit
})
data = data.set_index(['firm_id', 'year']) # linearmodels 要求：(个体，时间) 双索引
return data
linearmodels 识别「个体维度（firm_id）」和「时间维度（year）」的唯一方式，必须放在预处理前。preprocess_panel_datadef preprocess_panel_data(data):
data = data.dropna() # 缺失值处理：删除（模拟数据无缺失，实际数据需补充）
# 异常值处理：缩尾法（Winsorize），前后 1% 分位数截断
for col in ['R_D', 'assets', 'profit']:
q1 = data[col].quantile(0.01)
q99 = data[col].quantile(0.99)
data[col] = np.clip(data[col], q1, q99) # 把小于 q1 的设为 q1，大于 q99 的设为 q99
return data
estimate_panel_models为验证「双向固定效应模型的有效性」，同时估计 4 种面板模型做对比：
def estimate_panel_models(data):
y = data['profit'] # 被解释变量
X = data[['R_D', 'assets']] # 核心解释变量 + 控制变量
# 1. 混合回归：需加常数项（sm.add_constant）
pooled = PooledOLS(y, sm.add_constant(X)).fit()
# 2. 单向个体固定效应：entity_effects=True（无需加常数项，α_i 已包含个体异质性）
fe_individual = PanelOLS(y, X, entity_effects=True).fit()
# 3. 双向固定效应：同时指定 entity_effects=True（个体）和 time_effects=True（时间）
two_way_fe = PanelOLS(y, X, entity_effects=True, time_effects=True).fit()
# 4. 随机效应：需加常数项，若库未安装则返回 None
try:
re_model = RandomEffects(y, sm.add_constant(X)).fit()
except ImportError:
re_model = None
return pooled, fe_individual, two_way_fe, re_model
PanelOLS 的常数项：默认不包含常数项，因为固定效应 $\alpha_i$ 已涵盖个体层面的截距差异；若需全局常数项，需手动添加 sm.add_constant(X)。time_effects=True 自动生成时间哑变量，无需手动创建（如 2015-2019 年生成 4 个哑变量）。hausman_testHausman 检验用于模型选择：
def hausman_test(fe_results, re_results):
if re_results is None:
return "RandomEffects 模型未安装"
# 1. 筛选公共参数（FE 无常数项，RE 有，需统一）
common_params = list(set(fe_results.params.index) & set(re_results.params.index))
if not common_params:
return "无公共参数，无法检验"
# 2. 提取系数与协方差矩阵
fe_beta = fe_results.params[common_params]
re_beta = re_results.params[common_params]
fe_cov = fe_results.cov.loc[common_params, common_params] # FE 的协方差
re_cov = re_results.cov.loc[common_params, common_params] # RE 的协方差
# 3. 计算 Hausman 统计量（用广义逆 pinv 避免方差矩阵不可逆）
diff = fe_beta - re_beta
var_diff = fe_cov - re_cov
hausman_stat = diff.T @ np.linalg.pinv(var_diff) @ diff # np.pinv 是广义逆
# 4. 计算 p 值
from scipy.stats import chi2
df = len(diff)
p_value = 1 - chi2.cdf(hausman_stat, df)
return f"Hausman 统计量：{hausman_stat:.4f}, p 值：{p_value:.4f}"
np.linalg.pinv 可替代逆矩阵，保证检验可执行。R_D 和 assets 的系数。if __name__ == "__main__":
# 1. 生成模拟数据
panel_data = simulate_panel_data()
# 2. 数据预处理
processed_data = preprocess_panel_data(panel_data)
# 3. 模型估计
pooled_result, fe_individual_result, two_way_fe_result, re_result = estimate_panel_models(processed_data)
# 4. 输出结果（核心看双向 FE 的 R_D 系数是否接近真实值 0.8）
print("混合回归结果：", pooled_result, sep=" ")
print(" 单向个体 FE 结果：", fe_individual_result, sep=" ")
print(" 双向 FE 结果：", two_way_fe_result, sep=" ")
# 5. Hausman 检验（若 RE 安装，结果应拒绝原假设，选 FE）
print(" Hausman 检验：", hausman_test(fe_individual_result, re_result), sep=" ")
# 6. 可视化：研发投入与利润的关系（控制 assets=100 的简化拟合）
plt.figure(figsize=(12, 8))
plt.scatter(panel_data['R_D'], panel_data['profit'], alpha=0.3, label='原始数据')
# 双向 FE 拟合线：固定 assets=100（均值），斜率为 R_D 的边际效应
plt.axline((20, two_way_fe_result.params['assets']*100 + two_way_fe_result.params['R_D']*20),
slope=two_way_fe_result.params['R_D'],
color='red',
label=f'双向 FE 拟合线（斜率：{two_way_fe_result.params["R_D"]:.2f}）')
plt.xlabel('研发投入（R_D）')
plt.ylabel('企业利润')
plt.title('研发投入与企业利润的关系')
plt.legend()
plt.savefig('R_D_profit_relation.png', dpi=300) # 保存高清图（比赛要求）
plt.show()
0.8（误差由随机误差项 $\varepsilon_{it}$ 导致）。p 值<0.05，拒绝原假设，说明需用固定效应模型。simulate_panel_data 替换为 pd.read_csv 读取真实面板数据，注意必须保证双索引格式。IV2SLS（工具变量法）。这份代码是面板数据双向固定效应模型的标准化建模模板，从数据生成到结果输出完全符合数学建模比赛的要求。其核心优势在于：
linearmodels，代码简洁且不易出错。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
解析常见 curl 参数并生成 fetch、axios、PHP curl 或 Python requests 示例代码。 在线工具，curl 转代码在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML 转 Markdown 互为补充。 在线工具，Markdown 转 HTML在线工具，online
将 HTML 片段转为 GitHub Flavored Markdown，支持标题、列表、链接、代码块与表格等；浏览器内处理，可链接预填。 在线工具，HTML 转 Markdown在线工具，online