python面板数据分析代码_【译】用python做计量之面板数据模型

【译】用python做计量之面板数据模型

雁陎 • 2020 年 06 月 25 日

Loading…

原文：[https://bashtage.github.io/linearmodels/panel/index.html](https://bashtage.github.io/linearmodels/panel/index.html)

## 介绍

面板数据是具有时间序列(T)和截面(N)两个维度组织的数据。在面板数据的大多数经典应用中，观测的数量N大，而时间段的数量T小(通常在2到5之间)。这些估计量的渐近理论多是在N固定而T会发散的前提下得到的的。

大多数面板模型旨在估计下面方程的参数：

$$

y_{i t}=x_{i t} \beta+\alpha_{i}+\epsilon_{i t}

$$

其中i为观测，t为时间，$\beta$包含了要估计的参数。$\alpha_i$是个体的异质性，无法观测。$\epsilon_{i t}$是特异性误差(idiosyncratic errors)，与$\alpha_i$和$x_{it}$不相关。

所有的模型都需要两个输入：

* `dependent`，即因变量$y_{it}$

* `exog`，即自变量$x_{it}$

并使用不同的技术处理$\alpha_i$

* [`PanelOLS`](https://bashtage.github.io/linearmodels/panel/panel/linearmodels.panel.model.PanelOLS.html#linearmodels.panel.model.PanelOLS “linearmodels.panel.model.PanelOLS”)固定效应，即面板数据中随个体变化但不随时间变化的一类变量方法。这在数学上等价于为每个观测包含一个虚拟变量，尽管出于性能原因，实际上并不是这样实现的。

* [`BetweenOLS`](https://bashtage.github.io/linearmodels/panel/panel/linearmodels.panel.model.BetweenOLS.html#linearmodels.panel.model.BetweenOLS “linearmodels.panel.model.BetweenOLS”)平均数估计。即对面板数据的每个个体求平均数，然后利用自变量和因变量的N组观测值估计参数。

* [`FirstDifferenceOLS`](https://bashtage.github.io/linearmodels/panel/panel/linearmodels.panel.model.FirstDifferenceOLS.html#linearmodels.panel.model.FirstDifferenceOLS “linearmodels.panel.model.FirstDifferenceOLS”)一阶差分估计。即对个体固定效应模型中的回归量与被回归量的差分变量构成的模型的参数进行OLS估计。

* [`RandomEffects`](https://bashtage.github.io/linearmodels/panel/panel/linearmodels.panel.model.RandomEffects.html#linearmodels.panel.model.RandomEffects “linearmodels.panel.model.RandomEffects”)随机效应。当固定效应独立于回归变量时，RandomEffects 使用一个拟差分来有效地估计。 然而，当固定效应和回归变量之间存在依赖关系时，估计量就不再具有一致性。

* [`PooledOLS`](https://bashtage.github.io/linearmodels/panel/panel/linearmodels.panel.model.PooledOLS.html#linearmodels.panel.model.PooledOLS “linearmodels.panel.model.PooledOLS”)混合最小二乘估计。Pooledols 忽略了固定效应，当效应与回归因子无关时，它是一致的但是低效的。

Panelols 比其他的估计方法更通用，可以用于两类效应的模型(例如，实体(即固定)和时间效应)

> 模型中的$\alpha_i$，与回归自变量相关，称之为固定效应模型；与回归自变量不相关，称之为随机效应模型。

>

> 固定效应将α消掉，随机效应则将其放入误差项，然后探索方差结构。

模型规范类似于 statsmodels。 下面例子估计了一个固定效应回归面板上的劳动者的工资，模型中劳动者的工资是经验平方的函数，如果一个男人结婚了虚拟变量设为1，如果这个男人是工会成员同样设定虚拟变量。

“`python

from linearmodels.panel import PanelOLS

from linearmodels.datasets import wage_panel

import statsmodels.api as sm

data = wage_panel.load()

data = data.set_index([‘nr’,’year’])

dependent = data.lwage

exog = sm.add_constant(data[[‘expersq’,’married’,’union’]])

mod = PanelOLS(dependent, exog, entity_effects=True)

res = mod.fit(cov_type=’unadjusted’)

res

“`

结果：

“`python

PanelOLS Estimation Summary

================================================================================

Dep. Variable: lwage R-squared: 0.1365

Estimator: PanelOLS R-squared (Between): -0.0674

No. Observations: 4360 R-squared (Within): 0.1365

Date: Wed, Jun 24 2020 R-squared (Overall): 0.0270

Time: 23:03:18 Log-likelihood -1439.0

Cov. Estimator: Unadjusted

F-statistic: 200.87

Entities: 545 P-value 0.0000

Avg Obs: 8.0000 Distribution: F(3,3812)

Min Obs: 8.0000

Max Obs: 8.0000 F-statistic (robust): 200.87

P-value 0.0000

Time periods: 8 Distribution: F(3,3812)

Avg Obs: 545.00

Min Obs: 545.00

Max Obs: 545.00

Parameter Estimates

==============================================================================

Parameter Std. Err. T-stat P-value Lower CI Upper CI

——————————————————————————

const 1.3953 0.0123 113.50 0.0000 1.3712 1.4194

expersq 0.0037 0.0002 19.560 0.0000 0.0033 0.0041

married 0.1073 0.0182 5.8992 0.0000 0.0717 0.1430

union 0.0828 0.0198 4.1864 0.0000 0.0440 0.1215

==============================================================================

F-test for Poolability: 9.3360

P-value: 0.0000

Distribution: F(544,3812)

Included effects: Entity

“`

与statamodel类似，面板模型也可以使用类似 r 的公式来指定。 下面的模型和前面的一模一样。 注意使用特殊变量 EntityEffects 来包含固定效果。

“`python

mod = PanelOLS.from_formula(‘lwage ~ 1 + expersq + union + married + EntityEffects’,data)

“`

## 数据格式

使用面板数据回归模型，需要以下的数据格式

* 二维数据结构：MultiIndex DataFrames，其中外部索引是个体，内部索引是时间，可以使用pandas

* 三维数据结构：3D structures，维度0是变量，维度1是个体索引，维度2是时间索引，还可以使用具有维度(t，n)的2D 数据结构，该数据结构被视为具有维度(1，t，n)的3D 数据结构。 这些3D数据结构可以是 pandas、 NumPy 或 xarray。

> 因为本人习惯用pandas，所以就以pandas为例，其他numpy或者xarry自行去原文档查看。

最精确的数据格式是MultiIndex DataFrame，因为只有单个列可以保留面板中的所有类型。

下面示例使用培训数据，使用 set_index 命令构造 MultiIndex DataFrame。 个体索引为 fcode，时间索引为year。

“`python

from linearmodels.datasets import jobtraining

data = jobtraining.load()

print(data.head())

mi_data = data.set_index([‘fcode’, ‘year’])

print(mi_data.head())

“`

“`python

employ sales avgsal scrap rework tothrs union grant \

fcode year

410032 1987 100.0 47000000.0 35000.0 NaN NaN 12.0 0 0

1988 131.0 43000000.0 37000.0 NaN NaN 8.0 0 0

1989 123.0 49000000.0 39000.0 NaN NaN 8.0 0 0

410440 1987 12.0 1560000.0 10500.0 NaN NaN 12.0 0 0

1988 13.0 1970000.0 11000.0 NaN NaN 12.0 0 0

d89 d88 … grant_1 clscrap cgrant clemploy clsales \

fcode year …

410032 1987 0 0 … 0 NaN 0 NaN NaN

1988 0 1 … 0 NaN 0 0.270027 -0.088949

1989 1 0 … 0 NaN 0 -0.063013 0.130621

410440 1987 0 0 … 0 NaN 0 NaN NaN

1988 0 1 … 0 NaN 0 0.080043 0.233347

lavgsal clavgsal cgrant_1 chrsemp clhrsemp

fcode year

410032 1987 10.463100 NaN NaN NaN NaN

1988 10.518670 0.055570 0.0 -8.946565 -1.165385

1989 10.571320 0.052644 0.0 0.198597 0.047832

410440 1987 9.259130 NaN NaN NaN NaN

1988 9.305651 0.046520 0.0 0.000000 0.000000

[5 rows x 28 columns]

“`

当只引用单个序列时，可以使用 MultiIndex Series 表示形式。

“`python

from linearmodels import PanelOLS

mod = PanelOLS(mi_data.lscrap, mi_data.hrsemp, entity_effects=True)

print(mod.fit())

“`

## 案例

这些例子涵盖了可用于估计面板模型的模型。 最初的示例都忽略协方差选项，因此使用默认的经典协方差，该协方差适用于 homsedastic (同方差)数据。 备选的协方差选项将在本文档的末尾描述。

### 加载数据

这些例子都利用了工资面板数据，参考文献

> Vella and M. Verbeek (1998), “Whose Wages Do Unions Raise? A Dynamic Model of Unionism and Wage Rate Determination for Young Men,” Journal of Applied Econometrics 13, 163-183.

这组数据包括20世纪80年代男性的工资和特征。 个体标识符为`nr`，标识的时间为`year`。 这些数据在 Jeffrey Wooldridge 的《计量经济学导论》第14章中得到了广泛的应用。同时，再添加一列`year`，方便添加时间虚拟变量。

“`python

from linearmodels.datasets import wage_panel

import pandas as pd

data = wage_panel.load()

year = pd.Categorical(data.year)

data = data.set_index([‘nr’, ‘year’])

data[‘year’] = year

print(wage_panel.DESCR)

print(data.head())

“`

“`python

nr person identifier

year 1980 to 1987

black =1 if black

exper labor market experience

hisp =1 if Hispanic

hours annual hours worked

married =1 if married

educ years of schooling

union =1 if in union

lwage log(wage)

expersq exper^2

occupation Occupation code

“`

### 基本回归：混合回归

仅仅是一个普通的 OLS，方便理解各种面板数据结构。 它作为一个基本模型是很有用的。 在这里，工资的对数被设定为因变量，该模型使用所有的变量和时间虚拟变量。

“`python

from linearmodels.panel import PooledOLS

import statsmodels.api as sm

exog_vars = [‘black’,’hisp’,’exper’,’expersq’,’married’, ‘educ’, ‘union’, ‘year’]

exog = sm.add_constant(data[exog_vars])

mod = PooledOLS(data.lwage, exog)

pooled_res = mod.fit()

print(pooled_res)

“`

### 具有不相关效应的参数估计

在对面板数据建模时，通常会考虑超出 OLS 有效估计范围的模型。最常见的是误差成分模型，它在标准OLS模型中增加了一个附加项，其中[Math Processing Error(数学处理误差)]影响个体i的所有数值，当[数学处理误差]与[数学处理误差]中的回归因子不相关时，可以使用随机效应模型来有效估计该模型的参数。

#### 随机效应RE

随机效应模型与混合OLS模型几乎相同，但它考虑到了模型的结构，因此效率更高。随机效应采用了一种quasi-demeaning策略，即减去个体内数值的时间平均数，以考虑到共同的冲击。

“`python

from linearmodels.panel import RandomEffects

mod = RandomEffects(data.lwage, exog)

re_res = mod.fit()

print(re_res)

“`

模型的拟合度相当相似，尽管experience的回归发生了很大变化，其重要性也发生了变化。这在一定程度上是可以解释的，因为加入了年份虚拟变量，会拟合experience的趋势，所以只有横截面的差异才是重要的。随机效应估计器中的准差异取决于一个量，这个量取决于特异性冲击和共同冲击的相对方差。这可以使用`variance_decomposition`方法来访问。

“`python

re_res.variance_decomposition

“`

“`python

Effects 0.106946

Residual 0.123324

Percent due to Effects 0.464438

Name: Variance Decomposition, dtype: float64

“`

系数[数学处理误差]决定了demeaning发生的程度。当这个值为1时，RE模型就会还原为pooled模型，因为当effects没有方差时，就会出现这种情况。当面板不平衡时，它会在不同的实体之间变化，但在这个平衡的面板中，所有的值都是一样的。

#### 平均数估计BE

平均数估计是一种替代性的，通常效率较低的方法，可以用来估计模型参数。它特别简单，因为它首先计算[数学处理误差]和[数学处理误差]的时间平均数，然后用这些平均数进行简单的回归。

剔除了年份虚拟变量，因为平均化消除了因年份而产生的差异。 expersq也被剔除，因为它与exper相当共线性。这些结果与前面的模型大致相似。

“`python

from linearmodels.panel import BetweenOLS

exog_vars = [‘black’,’hisp’,’exper’,’married’, ‘educ’, ‘union’]

exog = sm.add_constant(data[exog_vars])

mod = BetweenOLS(data.lwage, exog)

be_res = mod.fit()

print(be_res)

“`

### 具有相关效应的参数估计FE

当效应与回归因子相关时，RE和BE估计量不一致。通常的解决方案是使用PanelOLS中的固定效应。固定效应应用于个体维度时称为实体效应，应用于时间维度时称为时间效应。

#### 实体效应

通过设置`entity_effects=True`参数，实体效应被包括在内。这相当于为每个个体加入了虚拟变量。在本面板中，这将增加545个虚拟变量，模型的估计速度将大大降低。PanelOLS实际上并没有使用虚拟变量，而是使用分组demeaning来达到同样的效果。

时间不变变量：在使用实体效应时，不能包含时间变量，因为一旦demeaned，这些变量将全部为0。这里也不包括exper，因为一旦纳入实体效应和时间虚拟值，exper是完全共线性的。

“`python

from linearmodels.panel import PanelOLS

exog_vars = [‘expersq’, ‘union’, ‘married’, ‘year’]

exog = sm.add_constant(data[exog_vars])

mod = PanelOLS(data.lwage, exog, entity_effects=True)

fe_res = mod.fit()

print(fe_res)

“`

#### 时间效应

可以使用`time_effects=True`添加时间效应。这里去掉了时间虚拟变量。请注意，核心系数是相同的。唯一的变化是对集合性的检验统计，因为 “效应 “不同时包括实体和时间，而之前只包括实体。

“`python

exog_vars = [‘expersq’,’union’,’married’]

exog = sm.add_constant(data[exog_vars])

mod = PanelOLS(data.lwage, exog, entity_effects=True, time_effects=True)

fe_te_res = mod.fit()

print(fe_te_res)

“`

**效应VS虚拟变量**

对于可以一致估计的变量，如通常的大N、固定T面板中的时间效应，没有必要将其作为效应列入。它们可以作为虚拟变量来实施。

**其他选项**

当使用PanelOLS和效应时，还有一些额外的选项可以使用。`use_lsdv`可以用来强制使用虚拟变量来估计模型。这不是必须的，而且在大多数情况下速度会比较慢。`auto_df`指示PanelOLS确定要进行的自由度调整。与大多数其他估计不同，效果的处理取决于协方差估计器。例如，当使用经典协方差估计器时，效应算作损失的自由度。当使用实体效应和按实体聚类时，则不会。设置`auto_df=False`将允许实体被计算或不计算，这取决于`count_effects`的值。

#### 一阶差分

在可能存在相关性的情况下，一阶差分是固定效应的一种替代方法。在使用一阶差分时，必须排除时间不变的变量。此外，只能包含一个线性时间趋势变量，因为它看起来像一个常数。这个变量将吸收数据中所有的时间趋势，因此对这些变量的解释可能是具有挑战性的。

“`python

from linearmodels.panel import FirstDifferenceOLS

exog_vars = [‘exper’,’expersq’, ‘union’, ‘married’]

exog = data[exog_vars]

mod = FirstDifferenceOLS(data.lwage, exog)

fd_res = mod.fit()

print(fd_res)

“`

### 比较模型

模型结果可以使用`compare`进行比较。compare接受结果列表、结果字典，其中键被解释为模型名称。

“`python

from linearmodels.panel import compare

print(compare({‘BE’:be_res,’RE’:re_res,’Pooled’:pooled_res}))

“`

“`python

Model Comparison

===============================================================

BE RE Pooled

—————————————————————

Dep. Variable lwage lwage lwage

Estimator BetweenOLS RandomEffects PooledOLS

No. Observations 545 4360 4360

Cov. Est. Unadjusted Unadjusted Unadjusted

R-squared 0.2155 0.1806 0.1893

R-Squared (Within) 0.1141 0.1799 0.1692

R-Squared (Between) 0.2155 0.1853 0.2066

R-Squared (Overall) 0.1686 0.1828 0.1893

F-statistic 24.633 68.409 72.459

P-value (F-stat) 0.0000 0.0000 0.0000

===================== ============ =============== ============

const 0.2836 0.0234 0.0921

(1.5897) (0.1546) (1.1761)

black -0.1414 -0.1394 -0.1392

(-2.8915) (-2.9054) (-5.9049)

hisp 0.0100 0.0217 0.0160

(0.2355) (0.5078) (0.7703)

exper 0.0278 0.1058 0.0672

(2.4538) (6.8706) (4.9095)

married 0.1416 0.0638 0.1083

(3.4346) (3.8035) (6.8997)

educ 0.0913 0.0919 0.0913

(8.5159) (8.5744) (17.442)

union 0.2587 0.1059 0.1825

(5.6214) (5.9289) (10.635)

expersq -0.0047 -0.0024

(-6.8623) (-2.9413)

year.1981 0.0404 0.0583

(1.6362) (1.9214)

year.1982 0.0309 0.0628

(0.9519) (1.8900)

year.1983 0.0202 0.0620

(0.4840) (1.6915)

year.1984 0.0430 0.0905

(0.8350) (2.2566)

year.1985 0.0577 0.1092

(0.9383) (2.5200)

year.1986 0.0918 0.1420

(1.2834) (3.0580)

year.1987 0.1348 0.1738

(1.6504) (3.5165)

—————————————————————

T-stats reported in parentheses

“`

### 协方差选项

#### 异方差性稳健共方差

Heteroskedasticity Robust Covariance

White的稳健协方差可以通过设置`cov_type=’robust’`来使用。这个估计器对某些类型的规范问题增加了一些稳健性，但在使用固定效应(实体效应)时不应使用，因为它不再是稳健的。取而代之的是需要一个聚类协方差。

“`python

exog_vars = [‘black’,’hisp’,’exper’,’expersq’,’married’, ‘educ’, ‘union’]

exog = sm.add_constant(data[exog_vars])

mod = PooledOLS(data.lwage, exog)

robust = mod.fit(cov_type=’robust’)

“`

#### 按实体聚类

通常要聚类的变量是实体或实体和时间。可以使用`cov_type=’clusterered’`和附加关键字参数`cluster_entity=True`和/或`cluster_time=True`来实现。

“`python

clust_entity = mod.fit(cov_type=’clustered’, cluster_entity=True)

clust_entity_time = mod.fit(cov_type=’clustered’, cluster_entity=True, cluster_time=True)

“`

`OrderedDict` 用于保存比较模型的结果。这允许对模型进行命名，也允许指定模型的顺序。一个标准的 dict 将产生有效的随机顺序。

“`python

from collections import OrderedDict

res = OrderedDict()

res[‘Robust’] = robust

res[‘Entity’] = clust_entity

res[‘Entity-Time’] = clust_entity_time

print(compare(res))

“`

对实体的聚类降低了整个T统计量。这表明每个实体的残差有重要的相关性。按两者进行聚类也会降低t统计量，这表明数据中存在横截面依赖性。注：按实体聚类解决的是跨时间的相关性，按时间聚类控制的是一个时间段内实体之间的相关性。

#### 其他聚类

其他聚类可以用输入的聚类直接传递整数组(1列或2列，或1-d数组)。本例按职业进行聚类，由于只有9组，而且通常聚类标准误差的理论要求聚类的数量要大，所以可能不是一个可靠的变量。

“`python

clust_entity = mod.fit(cov_type=’clustered’, clusters=data.occupation)

print(data.occupation.value_counts())

print(clust_entity)

“`

## 使用公式指定模型

所有的模型都可以使用公式来指定，和statsmodels中的公式类似。单变量回归的基础公式语法为

“`python

y ~ 1 + x

“`

用于BetweenOLS、PooledOLS和RandomEffects的公式是完全标准的，与statsmodels相同。FirstDifferenceOLS几乎完全相同，但需要注意的是，该模型不能包含截距。

实现效应(实体、时间或其他)的PanelOLS对公式有一个小的扩展，它允许将实体效应或时间效应(或两者)指定为公式的一部分。虽然不能使用公式接口指定其他效应，但在使用公式时，这些效应可以作为一个可选的参数包含在内。

### 加载数据

使用公式时，需要一个MultiIndex pandas数据框架，其中的指数是实体时间。这里用statsmodels提供的 “企业投资的决定因素 “中的Grunfeld数据来说明公式的使用。这个数据集包含了企业投资、市场价值和工厂资本存量的数据。

“`python

from statsmodels.datasets import grunfeld

data = grunfeld.load_pandas().data

data = data.set_index([‘firm’,’year’])

print(data.head())

“`

### 实体效应

实体效果使用特殊命令EntityEffects来指定。默认情况下不包含常数，所以如果需要常数，应该在公式中包含1+。当包含效应时，无论是否包含常数，模型和拟合都是相同的。

当一个常数被包含在模型中时，并施加额外的约束条件，即效应数为0，这样就可以使用因变量和回归变量的大平均数来确定常数。

“`python

from linearmodels import PanelOLS

mod = PanelOLS.from_formula(‘invest ~ value + capital + EntityEffects’, data=data)

print(mod.fit())

mod = PanelOLS.from_formula(‘invest ~ 1 + value + capital + EntityEffects’, data=data)

print(mod.fit())

“`

### 实体和时间效应

时间效应同样可以使用TimeEffect纳入。在许多模型中，时间效应可以被一致估计，因此它们可以等价地包含在使用分类变量的回归变量集中。

“`python

mod = PanelOLS.from_formula(‘invest ~ 1 + value + capital + EntityEffects + TimeEffects’, data=data)

print(mod.fit())

“`

### 平均数估计

其他的面板模型都是直截了当的，为了完整起见，将其包括在内。

“`python

from linearmodels import BetweenOLS, FirstDifferenceOLS, PooledOLS

mod = BetweenOLS.from_formula(‘invest ~ 1 + value + capital’, data=data)

print(mod.fit())

“`

### 一阶差分

第一个差分模型绝对不能包含一个常数，因为这在差分后是无法确定的。

“`python

mod = FirstDifferenceOLS.from_formula(‘invest ~ value + capital’, data=data)

print(mod.fit())

“`

### 混合回归

混合OLS估计是PanelOLS在没有效应时的一种特殊情况。它与statsmodels(或WLS)中的OLS实际上是相同的，但为了完整起见，将其包括在内。

“`python

mod = PooledOLS.from_formula(‘invest ~ 1 + value + capital’, data=data)

print(mod.fit())

“`

---

版权属于：雁陎

本文链接：https://www.sitstars.com/archives/92/

转载时须注明出处及本声明