在计量经济学中，经常要对时间序列数据进行回归建模。时间序列数据通常具有异方差（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 ‘ ’ 1
Residual 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 ‘ ’ 1Newey-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) %*% x
for (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 <- 0
for (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_b
a 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.
