R语言计量：Newey-West调整

　　在计量经济学中，经常要对时间序列数据进行回归建模。时间序列数据通常具有异方差（Heteroscedasticity）和自相关（Autocorrelation）的性质，此时使用传统的最小二乘法（OLS）估计回归参数虽然仍可得到参数的无偏估计，但是传统方法计算出来的参数方差具有偏差，会导致参数的t检验不准确，常出现虚假显著的情况。为避免这种情况，计量经济学中常对上述参数的方差进行调整，最常用的是Newey-West调整（Newey and West，1987）。

　　在R语言中，对回归系数的t检验进行Newey-West调整可以使用AER包中的NeweyWest函数和coeftest函数（其实NeweyWest来自sandwich包，coeftest函数来自lmtest包，AER将他们合在了一起）。AER包是《Applied Econometrics with R》（Christian Kleiber和Achim Zeileis）一书的配套数据代码集，这本书介绍了常用计量方法的R语言实现，感兴趣的可以仔细读一读。

　　言归正传，下面介绍NeweyWest函数和coeftest函数的用法：

coeftest(x, vcov. = NULL, df = NULL, ...)

NeweyWest(x, lag = NULL, order.by = NULL, prewhite = TRUE, adjust = FALSE,

diagnostics = FALSE, sandwich = TRUE, ar.method = "ols", data = list(),

verbose = FALSE)

　　coeftest函数用于对回归系数进行检验（t检验或z检验），其参数x表示要进行检验的对象，一般需是一个回归模型（即lm类型数据）；vcov.表示参数的方差（矩阵），当取默认值NULL时，使用的就是传统回归的估计方法；df表示t检验的自由度，当取默认值NULL时，自由度为n-k（即样本量减参数量），而当df取负值时，检验方法将会从t检验（t分布假设）变为z检验（正态分布假设）。

　　NeweyWest函数用于产生经Newey-West法调整后的方差（矩阵），其参数x表示要进行检验的对象，一般需是一个回归模型（即lm类型数据）；lag表示带宽（详解见后文），取默认值NULL时程序会自动根据Newey and West (1994)计算出最优值；order.by表示排序，因为时间序列需按时间排序，默认值为NULL，即默认原始数据已经是按时间顺序排好了；prewhite表示是否进行prewhite处理（详解间后文），默认值是TRUE，即使用一阶自回归（AR(1)）进行prewhite处理；adjust表示是否对输出的方差（矩阵）进行调整，默认值为FALSE即不调整，如为TRUE，则结果将乘以n/(n - k)进行调整；其余参数解释见后文。

　　理论上说，对于一个回归模型x，计算其经过Newey-West法调整后的t检验可以直接使用下面的表达式：

coeftest(x, vcov. = NeweyWest)

　　但是在一般的计量经济学教材或金融工程教材中，如Turan G. Bali等人的《Empirical asset pricing: The cross section of stock returns》，通常不使用prewhite处理，且lag使用4*(n/100)^(2/9)向下取整，其中n是样本量。即：

coeftest(x, vcov. = NeweyWest(x, lag = floor(4*(n/100)^(2/9)), prewhite = F)

　　而如果要对Fama-Macbeth回归的回归系数进行检验，我们只需构建一个因变量为所检验系数，自变量为1的回归方程，然后套入上述表达式中即可：

coeftest(lm(x~1), vcov. = NeweyWest)
coeftest(lm(x~1), vcov. = NeweyWest(lm(x~1), lag = floor(4*(n/100)^(2/9)), prewhite = F)

　　下面我们举一个例子：

library(AER)
# 随机生成样本量n为30的自变量x和因变量y
set.seed(20190331)
n <- 30
y <- rnorm(n) * 100
x <- runif(n) * 100
# 构建回归模型
lm_temp <- lm(y~x)
# 经典回归的参数检验
summary(lm_temp)
# coeftest函数默认的参数检验
coeftest(lm_temp)
# 默认的经Newey-West调整后的参数检验
coeftest(lm_temp, vcov. = NeweyWest)
# 一般计量经济学或金融工程教材中使用的Newey-West调整
coeftest(lm_temp, vcov. = NeweyWest(lm_temp, lag = floor(4*(n/100)^(2/9))))

　　其结果如下。

> # 经典回归的参数检验
> summary(lm_temp)Call:
lm(formula = y ~ x)Residuals:Min      1Q  Median      3Q     Max
-172.22  -36.29  -10.27   44.38  174.31 Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -93.7898    28.7032  -3.268  0.00287 **
x             0.6483     0.4612   1.406  0.17075
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 80.07 on 28 degrees of freedom
Multiple R-squared:  0.06594,   Adjusted R-squared:  0.03258
F-statistic: 1.977 on 1 and 28 DF,  p-value: 0.1708> # coeftest函数默认的参数检验
> coeftest(lm_temp)t test of coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -93.78980   28.70322 -3.2676 0.002868 **
x             0.64834    0.46115  1.4059 0.170752
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1> # 默认的经Newey-West调整后的参数检验
> coeftest(lm_temp, vcov. = NeweyWest)t test of coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -93.78980   37.33587 -2.5121  0.01804 *
x             0.64834    0.53002  1.2232  0.23144
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1> # 一般计量经济学或金融工程教材中使用的Newey-West调整
> coeftest(lm_temp, vcov. = NeweyWest(lm_temp, lag = floor(4*(n/100)^(2/9))))t test of coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -93.78980   37.34376 -2.5115  0.01807 *
x             0.64834    0.54410  1.1916  0.24343
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

　　Newey-West调整的R语言实现就介绍完毕了，下面我们简单介绍一下Newey-West调整的原理，做到不仅知其然，而且知其所以然。部分内容参考了Zeileis (2004)，感兴趣的可以去阅读原文。

　　对于一个时间序列的线性多元回归方程：

　　　　　　　　　　　　　　，（1）

　　可以写成矩阵的形式：

　　　　　　　　　　　　　　　　　　　　，（2）

　　也就是：

　　　　　　　　　　　　，（3）

　　使用OLS进行参数估计，得：

　　　　　　　　　　　　　　　　　　　，（4）

　　将式（2）代入式（4）得：

　　　　　　　　　　　　　　　　　　，（5）

　　根据回归模型的基本假设：

　　　　　　　　　　　　　　　　　　　　　，（6）

　　可得：

　　　　　　　　　　　　　　　　　　　　　，（7）

　　即回归参数beta的估计是无偏估计。进而，其协方差矩阵可以写成：

　　　　　　　　　　　　，（8）

　　其中：

　　　　　　　　　　　　　　　　　，（9）

　　表示和的协方差。

　　在经典的回归模型中，被认为是独立同分布，因此式（9）可以写成：

　　　　　　　　　　　　　　　　　　，（10）

　　可以使用残差估计出来（注意自由度是n – k – 1）：

　　　　　　　　　　　　　　　　　　　　　，（11）

　　但是对于时间序列数据来说，独立同分布的残差假设通常不符合实际，此时协方差矩阵如何估计？在仅考虑异方差的情况下，White (1980)指出，式（12）的协方差矩阵估计方法是渐进无偏的：

　　　　　　　　　　　　　　　　　　，（12）

　　既然如此，我们自然就会想到，在同时考虑异方差和自相关的情况下，协方差矩阵式（9）是否可以使用式（13）的形式估计？

　　　　　　　　　　　　　　　　，（13）

　　答案是不可以，通常情况下，按这种方式估计出来的矩阵不正定，不符合协方差矩阵的半正定性质。为了满足半正定性质，我们通常只考虑的前几阶自相关，即式（14）的中间红色部分。这种考虑当然也是合理的，因为实际情况中时间序列的自相关随阶数增加而衰减，也就是说只要时间间隔m足够大，和可以看成是独立的。

　　　　　　　　　　　，（14）

　　式（14）是一个带状矩阵，只有中间带状部分是非0的，其宽度由L决定，称为带宽，也就是前面提到的NeweyWest函数的lag参数。在对进行估计时，也不能直接使用，Newey and West (1987)给出一种估计方法：

　　　　　　　　　　　　　　　　，（15）

　　即权重系数随着m的增加线性衰减。当然，除了Newey and West (1987)提出的这种线性衰减法（即Bartlett法），还有许多其他权重赋值方法，详见Zeileis (2004)。

　　至此，同时考虑异方差和自相关情况的协方差就可以估计了，这就是著名的Newey-West调整。还有一个问题，L取多大合适？Newey and West (1994)给出以下建议。

　　首先构建：

　　　　　　　　　　　　　　　　　　　　　　　，（16）

其中，是标量，代表残差向量的第t个元素。

　　令：

　　　　　　　　　　　　　　　　　　　　　　，（17）

　　　　　　　　　　　　　　　　　　　　　，（18）

　　　　　　　　　　　　　　　　，（19）

　　　　　　　　　　　　　　　　　　　　　，（20）

　　　　　　　　　　　　　　　　　　　　　，（21）

　　　　　　　　　　　　　　　　　　　　，（22）

　　　　　　　　　　　　　　　　　　　　　　，（23）

　　由此可见不是带宽L，而是计算L所用的一个参数，计量经济学和金融工程课本中的表述可能有误。

　　带宽L的计算，可以使用bwNeweyWest函数。

bwNeweyWest(x, order.by = NULL, kernel = c("Bartlett", "Parzen",

"Quadratic Spectral", "Truncated", "Tukey-Hanning"), weights = NULL,

prewhite = 1, ar.method = "ols", data = list(), ...)

　　参数x表示要进行检验的对象，一般需是一个回归模型（即lm类型数据），kernel是不同的核函数，决定了式（18）和式（22）的计算，weights就是式（17）的权重。

　　可以看到，在NeweyWest函数和bwNeweyWest中，都存在参数prewhite，所谓的prewhite处理是怎么回事？

　　我们首先回到式（8），将其改造后得：

　　　　　　　　　　　　　　，（24）

　　方程右端的中间部分可以看作的协方差矩阵。因此之前我们直接对进行的估计，也可以转化为对进行估计，即估计的协方差矩阵。而所谓prewhite处理就是，不直接使用，而是先对做向量自回归处理（通常使用VAR(1)），然后取其残差估计。不过，经过prewhite处理得到的矩阵不能直接使用，要经过如下处理：

　　　　　　　　　　　　　　　　　　　，（25）

其中，I是单位矩阵，A是向量自回归的系数矩阵。

　　根据以上原理，我们也可以自己计算Newey-West调整后的协方差矩阵，顺便对AER包中的函数进行检验。

　　带宽L：

# 计算带宽L
n_temp <- floor(4*(30/100)^(2/9))
ww <- u * x
ss_0 <- function (x, n) {ss_sum <- t(x) %*% xfor (i in 1:n) {ss_temp <- t(x[1:(length(x) - i)]) %*% x[(1 + i):length(x)]ss_sum <- ss_sum + 2 * ss_temp}return(ss_sum)
}
ss_1 <- function (x, n) {ss_sum <- 0for (i in 1:n) {ss_temp <- t(x[1:(length(x) - i)]) %*% x[(1 + i):length(x)]ss_sum <- ss_sum + 2 * i * ss_temp}return(ss_sum)
}
gamma <- 1.1447 * ((ss_1(ww, n_temp)/ss_0(ww, n_temp))^2)^(1/3)
L <- gamma * n^(1/3)
# 我们自己计算的L
L
# bwNeweyWest函数计算的L
bwNeweyWest(lm_temp, prewhite = F)

> L[,1]
[1,] 11.24111
> # bwNeweyWest函数计算的L
> bwNeweyWest(lm_temp, prewhite = F)
[1] 11.24111

　　协方差矩阵：

## Newey-West协方差矩阵估计
ee <- u %*% t(u)
s_0 <- diag(diag(ee), ncol = length(u), nrow = length(u))
temp <- matrix(0, ncol = length(u), nrow = length(u))
L <- 11
for (i in 1:nrow(temp)) {for (j in 1:L) {if (i + j > ncol(temp)) {next()}temp[i, i + j] <- (1 - j/(L + 1)) * ee[i, i + j]}
}
s_L <- temp + t(temp)
s <- s_0 + s_L
var_b <- solve(t(xx) %*% xx) %*% t(xx) %*% s %*% xx %*% (solve(t(xx) %*% xx))
# 我们计算的协方差矩阵
var_b
# NeweyWest计算的协方差矩阵
NeweyWest(lm_temp, lag = 11, prewhite = F)

> var_ba           x
a 868.83744 -12.1655102
x -12.16551   0.1800566
> # NeweyWest计算的协方差矩阵
> NeweyWest(lm_temp, lag = 11, prewhite = F)(Intercept)           x
(Intercept)   868.83744 -12.1655102
x             -12.16551   0.1800566

参考文献

Bali T G, Engle R F, Murray S. Empirical asset pricing: The cross section of stock returns[M]. John Wiley & Sons, 2016.

Kleiber C, Zeileis A. Applied econometrics with R[M]. Springer Science & Business Media, 2008.

Newey W K, West K D. A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance 　　Matrix[J]. Econometrica, 1987, 55(3): 703-708.

Newey W K, West K D. Automatic lag selection in covariance matrix estimation[J]. The Review of Economic Studies, 1994, 61(4): 631-653.

White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity[J]. econometrica, 1980, 48(4): 817-838.

Zeileis A. Econometric Computing with HC and HAC Covariance Matrix Estimators[J]. Journal of Statistical Software, 2004, 11(10): 1-17.

R语言计量：Newey-West调整相关推荐

R语言计量（一）：一元线性回归与多元线性回归分析
文章目录一.数据调用与预处理二.一元线性回归分析三.多元线性回归分析 (一)解释变量的多重共线性检测 (二)多元回归 1. 多元最小二乘回归 2. 逐步回归 (三)回归诊断四.模型评价-常用的 ...
关于R语言计量检验之三平稳性（单位根）
单位根检验是指检验序列中是否存在单位根,因为存在单位根就是非平稳时间序列了.单位根就是指单位根过程,可以证明,序列中存在单位根过程就不平稳,会使回归分析中存在伪回归. ----- 2 单位根检验方法- ...
生信——R语言：1.windows软件安装与配置
跨专业搞生信一.安装软件 1.安装R语言直接在下面网址下载安装R语言,windows直接下一步无脑安装下载适用于 Windows 的 R-4.2.1.用于统计计算的 R 项目. (r-projec ...
r语言面板数据回归_R语言——伍德里奇计量经济导论案例实践第十三章横截面与面板数据（一）...
哈喽,停更了大概有三周的计量笔记又要重新开始啦!虽然美国的疫情没有停歇的迹象,可是依旧阻挡不了大学开学的热情.从8月3号开始上课到现在,也经历了很多事情,每天都是抱着死猪不怕开水烫的心情,暗地里安慰自 ...
r语言找不到cochrane函数_R语言——伍德里奇计量经济导论案例实践第十二章时间序列的序列相关和异方差问题...
在上一章节的复习笔记中,我们介绍了时间序列比较常见的AR模型和随机游走序列.在对时间序列进行回归时,我们和横截面数据一样做了很多假设,但是上一章内容没有回答如何解决误差项之间的序列相关性 (seria ...
R语言编写自定义函数使用Wilcoxon符号秩检验（Wilcoxon signed rank）实现多分组非参数成对检验（pairwise）、并使用p.adjust函数调整概率值
R语言编写自定义函数使用Wilcoxon符号秩检验(Wilcoxon signed rank)实现多分组非参数成对检验(Nonparametric pairwise multiple comparis ...
R语言使用ggplot2包geom_jitter()函数绘制分组（strip plot，一维散点图）带状图（双分类变量分组：色彩配置、添加箱图、位置参数调整）实战
R语言使用ggplot2包geom_jitter()函数绘制分组(strip plot,一维散点图)带状图(双分类变量分组:色彩配置.添加箱图.位置参数调整)实战目录 R语言使用ggplot2包ge ...
R语言可视化包ggplot2包调整线条粗细实战（Adjust Line Thickness）
R语言可视化包ggplot2包调整线条粗细实战(Adjust Line Thickness) 目录 R语言可视化包ggplot2包调整线条粗细实战(Adjust Line Thickness)
R语言学习 - 热图美化（数值标准化和调整坐标轴顺序）
生物信息学习的正确姿势 NGS系列文章包括NGS基础.在线绘图.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞 ...

R语言计量：Newey-West调整

R语言计量：Newey-West调整相关推荐

最新文章

热门文章