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基本分析流程).单细胞 ...
最新文章
- Spring源码分析【6】-ThreadLocal的使用和源码分析
- 材料成型计算机模拟第三版,材料成型计算机模拟考试复习资料.doc
- 赋值后页面不渲染_第七节：框架搭建之页面静态化的剖析
- [深入浅出WP8.1(Runtime)]Socket编程之UDP协议
- python调用sdk的文章_如何使用 python 接入虹软 ArcFace SDK
- python emoji 表情过滤
- oracle中的数据读取与查找
- win10安装oracle11g 服务端及配置详解
- 查询结果做缓存的例子
- 基础—数学—Exponential Family
- python和go哪个就业前景好_Python和Java就业前景对比
- 降维: 主成分分析(PCA） 局部线性嵌入（LLE）
- C++输出谢尔宾斯基三角形
- 计算机无法验证签名,计算机中win10/win7无法验证文件数字签名的解决方法
- Python如何爬取免费爬虫ip
- 新媒传信Java_新媒小课堂——多媒体、流媒体、富媒体
- Android手机录制屏幕及转GIF
- 要求公开华人程序员自杀真相，清华学霸被Facebook开除了
- 手机项目人力投入评估
- SQL注入及其危害、防御手段