第七讲 计量回归基础与应用（II）

樊潇彦 复旦大学经济学院 金融数据

---

setwd("D:\\...\\Ch07")
rm(list=ls()) 
install.packages(c("multcomp","HH","MASS"))  # 方差分析
install.packages(c("AER","plm")）            # 面板数据模型
load("Ch06_Data.RData")

1. 广义线性模型（GLM）

1.1 定义与分类

广义线性模型（Generalized linear model）是线性模型在以下两方面的推广：

$$E\left(y|X\right)=G\left(X\beta\right)$$

• 连接函数 $G$ 可以有多种形式，当 $G\left(X\beta\right)=X\beta$ 时就还原为线性模型；
• 误差项 $u$ 可以有多种分布，当 $y$ 连续时，$u$ 可以假定服从高斯分布、伽玛分布等，而当 $y$ 离散时，则可能服从二项分布。

模型类别	连接函数 $G$	误差项 $u$ 的分布
线性模型	线性函数 $X \beta$	高斯分布（正态分布）
Logit 模型	Logist函数 $\frac{e^{X\beta}}{1+e^{X\beta}}$	二项分布
Probit 模型	正态累积分布函数 $\Phi\left(X\beta\right)$	二项分布

1.2 程序实现

# Logit模型
fit.logit <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, 
                     data=glmdata, family=binomial())
summary(fit.logit)
# Probit模型
fit.probit <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, 
                     data=glmdata, family = binomial(link = "probit"))
summary(fit.probit)
# 模型预测与比较：其他变量取均值，看随着婚姻质量的变化，发生婚外情的概率如何变化
testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
                       age = mean(glmdata$age),
                       yearsmarried = mean(glmdata$yearsmarried),
                       religiousness = mean(glmdata$religiousness))
testdata$pre_logit <- predict(fit.logit, newdata=testdata, type="response") 
testdata$pre_probit <- predict(fit.probit, newdata=testdata, type="response")
# 比较两种模型的预测结果，相差很小
round(data.frame(testdata$pre_logit,testdata$pre_probit),3)

2. 方差分析（AOV）

2.1 定义

当影响因素 $A$ 是离散变量时，按照 $A$ 将数据 $y$ 分组，将总方差 $S_T$ 分解为组间方差（因素效应 $S_A$）和组内方差（误差效应 $S_E$），进而分析各部分对总变异的贡献。

$$S_T=S_A+S_E$$

以单因素一元分析为例，推导过程如下：

$$\begin{array}{*{20}{l}} {{S_T} = Var\left( {{y_{ia}}} \right) = E{{\left( {{y_{ia}} - \bar y} \right)}^2} = \sum\limits_{a = 1}^A {\sum\limits_{i = 1}^{{N_a}} {{{\left( {{y_{ia}} - \bar y} \right)}^2}} } }\\ { ~~~~~= \sum\limits_{a = 1}^A {\sum\limits_{i = 1}^{{N_a}} {{{\left( {{y_{ia}} - {{\bar y}_a} + {{\bar y}_a} - \bar y} \right)}^2}} }}\\ {~~~~~= \underbrace {\sum\limits_{a = 1}^A {\sum\limits_{i = 1}^{{N_a}} {{{\left( {{y_{ia}} - {{\bar y}_a}} \right)}^2}} } }_{{S_E}} + \underbrace {\sum\limits_{a = 1}^A {{N_a}{{\left( {{{\bar y}_a} - \bar y} \right)}^2}} }_{{S_A}}} \end{array}$$

其中 ${{\bar y}_a} = \frac{1}{{{N_a}}}\sum\limits_{i = 1}^{{N_a}} {{y_{ia}}}$ 为各组均值，$\bar y = \frac{1}{N}\sum\limits_{a = 1}^A {\sum\limits_{i = 1}^{{N_a}} {{y_{ia}}}}$ 为总平均值。

2.2 分类

• 单因素一元分析：例如把全国的收入差距分解为省间和省内的差异
  $$y_{i,a}=A_a+\varepsilon_{i,a} \Rightarrow S_T=S_A+S_E$$
• 多因素一元分析：例如把全国的收入差距按地区和行业两个因素，分解为地区间差异、行业间差异、地区和行业间的交互差异，以及其他差异。
  $${y_{i,a,b}} = {A_a} + {B_b} + {A_a}{B_b} + {\varepsilon _{i,a,b}}\\ \Rightarrow S = {S_A} + {S_B} + {S_{AB}} + {S_E}$$
• 单因素多元分析：例如我们想分析经济和社会发展程度的差异，用人均GDP、平均受教育程度等多个指标 $K>1$ 作为被解释变量，将其分解为省间和省内两部分差异。
  $$y_{i,k,a}=A_a + \varepsilon_{i,k,a}\Rightarrow S_T=S_A+S_E$$

2.3 程序实现

# 单因素一元分析
library(multcomp)
data(litter)
attach(litter)
aggregate(weight, by=list(gesttime), FUN=mean) # 分组均值
aov1 <- aov(weight ~ gesttime)                 # 服药时间单因素
summary(aov1)
## 从F-statistic和p-value来看，等同于线性回归：
ols1 <- lm(weight ~ gesttime, litter)           
summary(ols1)
# 双因素一元分析
aov2 <- aov(weight ~ gesttime + dose)          # 服药时间和剂量双因素，无交互项
summary(aov2)
library(HH)
ancova(weight ~ gesttime + dose, data=litter)  # 做图查看
aov2c <- aov(weight ~ gesttime * dose)         # 服药时间和剂量双因素，有交互项 
summary(aov2c)
ancova(weight ~ gesttime * dose, data=litter)
detach(litter)
# 单因素多元分析
library(MASS)
attach(UScereal)
y <- cbind(calories, fat, sugars)              # 3个被解释变量
aggregate(y, by=list(shelf), FUN=mean)         # 分组计算均值
fit <- manova(y ~ shelf)                       # 一个因素对多元变量做方差分析
summary(fit)                                   # 报告结果
summary.aov(fit)                               # 对每个被解释变量报告结果
detach("UScereal")

3. 面板数据模型

3.1 模型描述

如果要分析多家上市公司历年的数据，就要用到面板数据模型：

$${y_{it}} = {\alpha _i} + {X_{it}}\beta + {\varepsilon_{it}}~~~~\left( {i = 1,2...N,~t = 1,2...T} \right)$$
假定：(1) $E\left(\varepsilon|X\right)=0$；(2) $\varepsilon \sim N\left( {0,{\sigma ^2}I_{N\times T}} \right)$；(3) $N>>T$

3.2 模型分类

• 混合面板（pooled）模型：
  $${\alpha _i} = {\alpha _j} = \alpha ~\text{或}~ Var\left( {{\alpha _i}} \right){\rm{ = }}0$$

• 固定效应（fixed effects）模型：
  

  $$Var\left( {{\alpha _i}} \right) \ne 0 ~\text{且}~ Cov\left( {{\alpha _i},{X_{it}}} \right) \ne 0$$
  比如在回归收入方差时，总有工作态度、工作能力等不可观测的个人因素在起作用，这将导致实际的残差项 ${u_{it}} = {\alpha _i} + {\varepsilon _{it}}$ 很可能与职位、受教育程度等其他解释变量相关，进而导致这些解释变量的系数估计有偏。

• 随机效应（random effects）模型：
  

  $$Var\left( {{\alpha _i}} \right) \ne 0 ~\text{且}~ Cov\left( {{\alpha _i},{X_{it}}} \right) = 0$$
  虽然这种情况不会造成估计结果有偏，但由于 $Cov\left( {{u_{it}},{u_{is}}} \right){\rm{ = }}\alpha _i^2 \ne 0$，将使参数估计不再有效，需要进行标准差调整。

3.3 程序实现

library(AER)
data(Grunfeld)
gr=subset(Grunfeld, firm %in% c("General Electric","General Motors", "IBM"))
## 选取三个公司的数据，设置为面板
library(plm)
pgr=pdata.frame(gr,c("firm","year"))                                 
## 混合数据回归(POOL)
gr_pool=plm(invest~value+capital, data=pgr, model="pooling")      
summary(gr_pool)
gr_ols=lm(invest~value+capital, data=gr)  # 等同于OLS回归
summary(gr_ols)
## 固定效应(FE)
gr_fe=plm(invest~value+capital, data=pgr, model="within")         
summary(gr_fe)
## 随机效应(RE)
gr_re=plm(invest~value+capital, data=pgr, model="random", random.method="walhus")   
summary(gr_re)
## Hausman检验
phtest(gr_re,gr_fe)  # 如果拒绝H0应选择固定效应模型

3.4 动态面板数据模型

当解释变量包含被解释变量的滞后项时，有动态面板数据模型：

$${y_{it}} = {\alpha _i} + \rho y_{i,t-1} + {X_{it}}\beta + {\varepsilon_{it}}~~~~\left( {i = 1,2...N,~t = 1,2...T} \right)$$
假定：(1) $E\left(\varepsilon|X\right)=0$；(2) $\varepsilon \sim N\left( {0,{\sigma ^2}I_{N\times T}} \right)$；(3) $N>>T$

• Arellano and Bond(1991)：差分广义矩（DIF-GMM）方法
  先对原模型做一阶差分：
  

  $${\Delta y_{it}} = y_{it}-y_{i,t-1} = \rho \Delta y_{i,t-1} + \Delta {X_{it}}\beta + {\Delta \varepsilon_{it}}$$
  由于 $E(y_{i,t-1},\varepsilon_{i,t-1}) \ne 0$，因此 $E(\Delta y_{i,t-1},\Delta \varepsilon_{it})\ne 0$，所以直接估计 $\rho$ 将有偏，但是：
  $$E(\Delta y_{i,t-1},y_{i,t-2})\ne 0,~~~E(y_{i,t-2},\Delta \varepsilon_{it})= 0$$
  因此可以用 $y_{i,t-2}$ 作为 $\Delta y_{i,t-1}$ 的工具变量估计差分方程。

• Blundell and Bond(1998)：系统广义矩（SYS-GMM）方法
  根据原方程有：
  

  $$\Delta{y_{i,t-1}} = {\alpha _i} + (\rho-1) y_{i,t-2} + {X_{i,t-1}}\beta + {\varepsilon_{i,t-1}}$$
  由于 $E(y_{i,t-2},\alpha_i) > 0$，因此如果回归上式将使 $\rho -1$ 被高估。这意味着即使 $\rho <0$，估计结果也会使 $\hat \rho -1 \to 0$， $y_{i,t-2}$ 好像和 $\Delta{y_{i,t-1}}$ 无关，即存在弱工具变量（weak IV）的问题。

对此Blundell and Bond(1998)提出可以把原方程和差分方程联立、同时回归。由于 $E(y_{i,t-1},\Delta {y_{i,t-1}}) \ne 0$，因此用 $\Delta {y_{i,t-1}}$ 作为原方程 $y_{i,t-1}$ 的工具变量。

• 程序实现：

data("EmplUK", package = "plm")
## Arellano and Bond (1991), Table 4b 
z1 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
           + log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99),
           data = EmplUK, effect = "twoways", model = "twosteps")
summary(z1)
## Blundell and Bond (1998) Table 4 (cf DPD for OX p.12 col.4)
z2 <- pgmm(log(emp) ~ lag(log(emp), 1)+ lag(log(wage), 0:1) +
             lag(log(capital), 0:1) | lag(log(emp), 2:99) +
             lag(log(wage), 2:99) + lag(log(capital), 2:99),        
           data = EmplUK, effect = "twoways", model = "onestep", 
           transformation = "ld")
summary(z2, robust = TRUE)

阅读讨论

1. 孔丹凤、孙宇辰、马驰骋、秦大忠，2015：信用风险转移工具真的转移了风险吗?——基于美国上市银行面板数据的实证研究，《金融研究》，2015:2
2. 张亦春、李晚春、彭江，2015：债权治理对企业投资效率的作用研究——来自中国上市公司的经验证据，《金融研究》，2015.7

参考资料

1. Kleiber, C. and A. Zeileis 2008: Applied Econometrics with R (Use R!), Springer
2. R.I. Kabacoff著：《R语言实战（第2版）》，王小宁、刘撷芯、黄俊文译，人民邮电出版社，2016
3. 薛毅、陈立萍编著，2012：《统计建模与R软件》，清华大学出版社