R：使用 dynlm 计算 waldtest
R: compute waldtest using dynlm
我运行以下两个回归：
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
我想计算 parameterequality 的 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
最简单的方法是拟合具有交互作用的嵌套模型，而不是两个单独的模型。所以可以先生成一个对这两个段进行编码的因子：
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)
问题 2：
m3<- update(m1, formula = . ~ . - d(x))
waldtest(m1,m3)
根据 Achim Zeileis 的评论，这应该是：
linearHypothesis(m1, "d(z) = d(x)")
我运行以下两个回归：
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
我想计算 parameterequality 的 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
最简单的方法是拟合具有交互作用的嵌套模型，而不是两个单独的模型。所以可以先生成一个对这两个段进行编码的因子：
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)
问题 2：
m3<- update(m1, formula = . ~ . - d(x))
waldtest(m1,m3)
根据 Achim Zeileis 的评论，这应该是：
linearHypothesis(m1, "d(z) = d(x)")