面板数据回归

遗漏变量偏差

遗漏变量偏差(Omitted variable bias，简称OV bias)表示还有模型遗漏的变量和XY同时相关。OV bias会导致回归假设\(E(u|X)=0\)不成立，从而使\(\beta\)的估计结果不一致(not consistent)，也就是说，\(\beta\)的估计结果是不准确的。

OV bias非常常见，因为现实中的问题涉及的因素非常多，很难把所有因素都考虑在内。为了使回归系数的估计结果更加准确，需要使用一些方法解决OV bias的问题，面板数据回归就是其中一种方法。

面板数据回归(panel data regression)专门处理一个个体被多次观测的情况。比如现在有336个observations，分别是美国48个州7年的数据，每个州一年只有一个观测，但是一个州共被观测到7次，类似下面这样

1
2
3
4
5
6
7
8
city  year  x1   x2   y
A     2000  0.5  0.4  0.3
A     2001
...             .....
A     2006
B     2000
B     2001
...(共336行)

如果我们不关注这些观测值属于哪个州，哪一年，只把它看做336*3的数据,用y ~ x1 + x2进行线性回归，就容易产生OV bias

举一个例子

  • y是交通事故发生率
  • x1是啤酒税
  • x2是法定驾驶年龄

这里没有考虑“交通拥堵情况”这个指标。而根据经验

  • 交通越拥堵，事故发生率越高
  • 西部的州，交通更不拥堵，且啤酒税更低
  • 交通拥堵情况不随时间变化（每一年都差不多，但是在不同州有差异）

那么“交通拥堵情况”这个指标就会产生OV bias。类似的变量可能还有“对驾车和喝酒的态度”、“道路质量”等等因素。

对于这种不随时间变化但又同时影响XY的变量，我们没有办法全部找出来，但是可以用一个指标来表示每一个州的“特质”，这个指标就可以控制住这类的变量。控制时间因素的指标也是一样。这就是面板数据回归的思想。

固定效应

这里考虑固定个体效应和固定时间效应。在只考虑一个X的情况下，固定个体效应的回归方程有两种形式（公式字母含义大同小异，请读者自行参考计量教材）

\[Y_{it} = \beta_1X_{it}+\alpha_i + u_{it}\] 和

\[Y_{it} = \beta_0 + \beta_1X_{it}+\gamma_2D2_i+...+\gamma_nDn_i + u_{it}\] 首先要说明，上面两个公式的本质是一样的。

  • 两式需要估计的参数数量是一样的
  • 两个公式对每一个州来说，都是\(X_{it}\)加对应州的一个截距（相同州的样本截距项相同，不同州不同，即用每个州自己的截距项固定住州的特征）

求解出这些系数有两种方法

第一种直接把固定效应当做k-1个0-1变量进行OLS回归

\[Y_{it} = \beta_0 + \beta_1X_{it}+\gamma_2D2_i+...+\gamma_nDn_i + u_{it}\]

第二种是先对所有样本按照州分组，每组内计算Y和X的均值，每个样本减去对应组均值，拿新数据进行回归。

\[Y_{it} = \alpha_i +\beta_1X_{it}+u_{it}\] 回归结果必满足

\[\frac1T\sum_{t=1}^TY_{it}=\alpha_i + \beta_1\frac1T\sum_{t=1}^TX_{it}+\frac1T\sum_{t=1}^Tu_{it}\] 经过变换

\[Y_{it}-\frac1T\sum_{t=1}^TY_{it}= \beta_1(X_{it}-\frac1T\sum_{t=1}^TX_{it})+(u_{it}-\frac1T\sum_{t=1}^Tu_{it})\] 得到新数据

\[\tilde Y_{it}=\beta_1\tilde X_{it}+\tilde u_{it}\] 第二种形式的理解：将各自州的特性减掉，得到“纯净”的数据放在一起进行回归。

时间效应同理，这里不再赘述。

R语言面板数据回归

R语言中进行面板数据回归使用plm包，安装后导入

1
2
library(plm)                      # 导入包
data("Grunfeld", package = "plm") # 导入数据

数据形式如下 

1
2
3
4
5
6
7
8
> head(Grunfeld)
  firm year   inv  value capital
1    1 1935 317.6 3078.5     2.8
2    1 1936 391.8 4661.7    52.6
3    1 1937 410.6 5387.1   156.9
4    1 1938 257.7 2792.2   209.2
5    1 1939 330.8 4313.2   203.4
6    1 1940 461.2 4643.9   207.2

不考虑个体效应和时间效应的情况下，进行OLS回归（结果和lm函数回归相同）

1
2
fit_lr <- plm(inv ~ value + capital,
              data = Grunfeld, model = "pooling")

结果

1
2
3
> coef(fit_lr)
(Intercept)       value     capital 
-42.7143694   0.1155622   0.2306785

固定个体效应(两种方式)

1
2
3
4
5
fit_individual1 <- plm(inv ~ value + capital, data = Grunfeld, 
                      model = "within", effect = "individual")
                      
fit_individual2 <- plm(inv ~ value + capital, data = Grunfeld, 
                      model = "within", index = "firm")

结果

1
2
3
4
5
6
7
> coef(fit_individual1)
    value   capital 
0.1101238 0.3100653 

> coef(fit_individual2)
    value   capital 
0.1101238 0.3100653

固定时间效应类似。时间和个体同时固定代码如下

1
2
fit_both <- plm(inv ~ value + capital, data = Grunfeld,       
                model = "within", effect = "twoways")

结果

1
2
3
> coef(fit_both)
    value   capital 
0.1177159 0.3579163

随机效应

在接触panel data的时候，经常会在软件中见到fixed effect和random effect，前者就是我们所说的固定效应，而后者称为随机效应。

随机效应不是本文讨论的重点，但是这是一个曾经困扰我很久的点，所以我希望在这里留下一些关于随机效应的学习资料。

随机效应是多阶段回归模型，只在特定条件下才需要使用。随机效应的形式不是唯一的，在plm包的plm函数中，当参数为model = "random"时，指定了随机效应。同时，还需要指定另一个参数

1
random.method = c("swar", "walhus", "amemiya", "nerlove", "kinla")

这是不同理论下的不同随机效应形式，用不同的方法做出来的结果会有很大差异，详见这里。

这里的每一种方法都有一篇论文做支撑，如果想要真正了解每一种随机效应是怎么做的，需要阅读这些论文。论文可以到plm帮助文档中的Details部分查看。

了解随机效应理论的另外一种好办法是直接看实现它们的源代码，通过查看plm函数源代码，看参数是如何调用的，抽丝剥茧即可找到函数定义位置。最后找到函数定义位置只要在R语言控制台输入getAnywhere(ercomp.formula)即可查看。

除此之外，plm包的学习除了看帮助文档中的参数解释和代码示例，还可以看plm包的论文，这个其实在包的User guides, package vignettes and other documentation.部分可以找到。

固定效应代码实现验证

上文提到，面板数据回归系数估计有两种形式，一种是多个0-1变量的OLS回归，一种是通过减去均值变换后进行回归。下面我们来分别使用这两种方法与plm包结果进行对比。以固定个体效应为例

首先贴出plm包的结果

1
2
3
4
5
> fit_individual1 <- plm(inv ~ value + capital, data = Grunfeld, 
+                        model = "within", effect = "individual")
> coef(fit_individual1)
    value   capital 
0.1101238 0.3100653

1.用多个0-1变量的OLS回归

\[Y_{it} = \beta_0 + \beta_1X_{it}+\gamma_2D2_i+...+\gamma_nDn_i + u_{it}\]

其实相当于在最简单OLS回归的基础上，加了一个离散变量firm，离散变量加入回归中，R语言会自动转化为多个0-1变量

1
2
fit_01 <- lm(inv ~ value + capital + factor(firm),
             data = Grunfeld)

得到结果如下 

1
2
3
4
5
6
7
> coef(fit_01)
   (Intercept)          value        capital  factor(firm)2  factor(firm)3
   -70.2967175      0.1101238      0.3100653    172.2025312   -165.2751236
 factor(firm)4  factor(firm)5  factor(firm)6  factor(firm)7  factor(firm)8
    42.4874229    -44.3200953     47.1354223      3.7432439     12.7510602
 factor(firm)9 factor(firm)10 
   -16.9255550     63.7288739

我们只要关心value和capital两个变量的系数即可，发现它们和plm包给出的结果是一样的。

2.通过减去均值变换后进行回归

\[Y_{it}-\frac1T\sum_{t=1}^TY_{it}= \beta_1(X_{it}-\frac1T\sum_{t=1}^TX_{it})+(u_{it}-\frac1T\sum_{t=1}^Tu_{it})\]

1
2
3
4
5
6
7
8
9
library(dplyr)
# 定义变换函数
get_new <- function(x) x - mean(x) 
# 数据变换（分组减去个体均值）
new_Grunfeld <- Grunfeld %>% group_by(firm) %>% 
  mutate_at(vars(inv:capital), get_new)
  
fit_mean <- lm(inv ~ 0 + value + capital, # 加0表示不含有截距项
               data = new_Grunfeld)

得到结果如下

1
2
3
> coef(fit_mean)
    value   capital 
0.1101238 0.3100653

结果也和plm包相同。

面板回归的统计理解

从上面的实现代码可以看出，从统计的角度来看，固定个体效应，就是把表示个体的那个变量加入回归模型之中，如果考虑时间效应，则再把离散的时间变量加入进来。

这个过程非常简单，也很符合直觉。拿到一个5列的数据，当然是把5列数据都加进去跑模型。离散变量和连续变量没有本质区别，都是变量，都是信息，没有理由放着不用。只不过离散型变量需要转化为多个变量进行回归而已。

但是这个问题到了计量经济学里面，却变复杂了许多。不仅新定义了“面板数据回归”、“固定效应”这些概念；增加了公式的推导；而且非常注重解释这样做的经济学含义，力图说清楚这样做是如何消除了OV bias的影响，如何使回归结果更可靠。统计人学了半天，拨开云雾，才发现原来这就是加了一个离散型变量。

这个现象让我开始思考统计学与计量经济学的区别。

计量经济学会使用大量的统计方法，但是又不只局限于使用，他们更在乎的是模型的理解，更关注因果关系而不是相关关系。这是计量经济学的目的所决定的，他们要探究一些经济现象，和现象背后的原因。他们只能谨慎地关注少数几个变量，考虑各种因素，使用各种方法，尽量消除它们的内生性，使结果更可靠。

我们经常说统计学是万金油，跟业务百搭。统计学研究创造各种方法，为业务问题提供工具。当业务问题不满足于使用工具，延伸出自己的一套思路的时候，它们就独立成一门新的学科了，计量经济学是如此，机器学习也是如此。

我们可以看到这样一个有趣的现象：统计中的离散型变量回归转化为dummy variable(0-1变量)的思想，在计量经济学中叫panel data（其实DID也有用到这个思想），在机器学习中叫one hot编码。计量经济学拿着统计学工具更加谨慎地探索，机器学习拿着统计学方法更加暴力地计算！