R:使用 dynlm 计算 waldtest

Question

我运行了以下两个回归:

library("dynlm")
library("lmtest")
zoop <- (test1[, -1])
f <- d(y) ~ d(x)+ d(z) + d(m) + d(log(p))
m1 <- dynlm(f, data = zoop, start = 1,end = 15)
coeftest(m1, vcov=NeweyWest)
m2 <- dynlm(f, data = zoop, start = 16,end = 31)
coeftest(m2, vcov=NeweyWest)

这给了我 m1 的输出:

               Estimate  Std. Error t value Pr(>|t|)  
(Intercept) -0.00068055  0.00021691 -3.1374  0.01056 *
d(x)         0.27475798  0.10605395  2.5907  0.02692 *
d(z)         0.00046720  0.00129363  0.3612  0.72550  
d(m)         0.00047590  0.00024276  1.9604  0.07838 .
d(log(p))    0.01876845  0.00829852  2.2617  0.04723 *

和 m2 的输出:

               Estimate  Std. Error t value Pr(>|t|)  
(Intercept) -0.00037592  0.00023431 -1.6044  0.14309  
d(x)         0.29475934  0.12946162  2.2768  0.04882 *
d(z)        -0.00514108  0.00219475 -2.3424  0.04384 *
d(m)        -0.00011535  0.00065369 -0.1765  0.86383  
d(log(p))   -0.00501189  0.03535847 -0.1417  0.89040  

我想计算参数相等的 waldtest 例如对于来自两个模型的变量 d(z)，即:d(z) -d(z) = 0(其中第一个 d(z) 来自模型 m1，第二个d(z) 来自 m2。如何使用 R 来做到这一点？还有第二个类似的问题:我也想计算模型一的 waldtest，例如 d(z)-d(x) = 0？非常感谢提前!

数据样本:

Date         y       x      m       z       p
03.01.2005  2.154   2.089   14.47   17.938  344999
04.01.2005  2.151   2.084   14.51   17.886  344999
05.01.2005  2.151   2.087   14.42   17.95   333998
06.01.2005  2.15    2.085   13.8    17.95   333998
07.01.2005  2.146   2.086   13.57   17.913  333998
10.01.2005  2.146   2.087   12.92   17.958  333998
11.01.2005  2.146   2.089   13.68   17.962  333998
12.01.2005  2.145   2.085   14.05   17.886  339999
13.01.2005  2.144   2.084   13.64   17.568  339999
14.01.2005  2.144   2.085   13.57   17.471  339999
17.01.2005  2.143   2.085   13.2    17.365  339999
18.01.2005  2.144   2.085   13.17   17.214  347999
19.01.2005  2.143   2.086   13.63   17.143  354499
20.01.2005  2.144   2.087   14.17   17.125  354499
21.01.2005  2.143   2.087   13.96   17.193  354499
24.01.2005  2.143   2.086   14.11   17.283  354499
25.01.2005  2.144   2.086   13.63   17.083  354499
26.01.2005  2.143   2.086   13.32   17.348  347999
27.01.2005  2.144   2.085   12.46   17.295  352998
28.01.2005  2.144   2.084   12.81   17.219  352998
31.01.2005  2.142   2.084   12.72   17.143  352998
01.02.2005  2.142   2.083   12.36   17.125  352998
02.02.2005  2.141   2.083   12.25   17  357499
03.02.2005  2.144   2.088   12.38   16.808  357499
04.02.2005  2.142   2.084   11.6    16.817  357499
07.02.2005  2.142   2.084   11.99   16.798  359999
08.02.2005  2.141   2.083   11.92   16.804  355500
09.02.2005  2.142   2.08    12.19   16.589  355500
10.02.2005  2.14    2.08    12.04   16.5    355500
11.02.2005  2.14    2.078   11.99   16.429  355500

Answer 1

最简单的方法是拟合具有交互作用的嵌套模型，而不是两个单独的模型。因此，您可以首先生成对两个段进行编码的因子:

fac <- factor(as.numeric(time(zoop) > as.Date("2005-01-24")))
fac <- zoo(fac, time(zoop))

然后你可以拟合一个模型，其中所有系数都被约束为相等( M0 )和一个它们不同的模型( M1 )。后者相当于单独的模型 m1和 m2你适合上面:

M0 <- dynlm(d(y) ~ d(x)+ d(z) + d(m) + d(log(p)),
  data = zoop, start = as.Date("2005-01-04"))
M1 <- dynlm(d(y) ~ fac * (d(x)+ d(z) + d(m) + d(log(p))),
  data = zoop, start = as.Date("2005-01-04"))

然后，您可以通过查看来轻松地单独测试所有系数的差异

coeftest(M1, vcov = NeweyWest)
## t test of coefficients:
## 
##                   Estimate  Std. Error t value Pr(>|t|)   
## (Intercept)    -0.00068055  0.00020683 -3.2904 0.003847 **
## fac1            0.00030463  0.00027961  1.0895 0.289561   
## d(x)            0.27475798  0.09503821  2.8910 0.009361 **
## d(z)            0.00046720  0.00129029  0.3621 0.721280   
## d(m)            0.00047590  0.00028483  1.6708 0.111147   
## d(log(p))       0.01876845  0.01134940  1.6537 0.114617   
## fac1:d(x)       0.02000136  0.16417829  0.1218 0.904315   
## fac1:d(z)      -0.00560828  0.00237869 -2.3577 0.029261 * 
## fac1:d(m)      -0.00059126  0.00073696 -0.8023 0.432305   
## fac1:d(log(p)) -0.02378034  0.03454593 -0.6884 0.499540   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fac1:d(z)的线路对 d(z) 中的差异进行了测试连续下坡。如果你想联合测试所有系数，你可以这样做:

waldtest(M0, M1, vcov = NeweyWest)
## Wald test
## 
## Model 1: d(y) ~ d(x) + d(z) + d(m) + d(log(p))
## Model 2: d(y) ~ fac * (d(x) + d(z) + d(m) + d(log(p)))
##   Res.Df Df      F  Pr(>F)  
## 1     24                    
## 2     19  5 3.7079 0.01644 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

为了完全重现，我包含了我用来读取数据的代码:

zoop <- read.zoo(textConnection("Date         y       x      m       z       p
03.01.2005  2.154   2.089   14.47   17.938  344999
04.01.2005  2.151   2.084   14.51   17.886  344999
05.01.2005  2.151   2.087   14.42   17.95   333998
06.01.2005  2.15    2.085   13.8    17.95   333998
07.01.2005  2.146   2.086   13.57   17.913  333998
10.01.2005  2.146   2.087   12.92   17.958  333998
11.01.2005  2.146   2.089   13.68   17.962  333998
12.01.2005  2.145   2.085   14.05   17.886  339999
13.01.2005  2.144   2.084   13.64   17.568  339999
14.01.2005  2.144   2.085   13.57   17.471  339999
17.01.2005  2.143   2.085   13.2    17.365  339999
18.01.2005  2.144   2.085   13.17   17.214  347999
19.01.2005  2.143   2.086   13.63   17.143  354499
20.01.2005  2.144   2.087   14.17   17.125  354499
21.01.2005  2.143   2.087   13.96   17.193  354499
24.01.2005  2.143   2.086   14.11   17.283  354499
25.01.2005  2.144   2.086   13.63   17.083  354499
26.01.2005  2.143   2.086   13.32   17.348  347999
27.01.2005  2.144   2.085   12.46   17.295  352998
28.01.2005  2.144   2.084   12.81   17.219  352998
31.01.2005  2.142   2.084   12.72   17.143  352998
01.02.2005  2.142   2.083   12.36   17.125  352998
02.02.2005  2.141   2.083   12.25   17  357499
03.02.2005  2.144   2.088   12.38   16.808  357499
04.02.2005  2.142   2.084   11.6    16.817  357499
07.02.2005  2.142   2.084   11.99   16.798  359999
08.02.2005  2.141   2.083   11.92   16.804  355500
09.02.2005  2.142   2.08    12.19   16.589  355500
10.02.2005  2.14    2.08    12.04   16.5    355500
11.02.2005  2.14    2.078   11.99   16.429  355500"),
format = "%d.%m.%Y", header = TRUE)