R: compute waldtest using dynlm
I have performed the following two regressions:
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)
which gives me the output for 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 *
and the output for m2 is:
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
I would like to compute waldtest for a parameter, for example d (z) variables from both models, i.e. d (z) -d (z) = 0 (where the first d (z) is from module m1 and the second d (z) is m2. How can this be done with R? And a second similar question: I would also like to compute waldtest eg for a model like d (z) -d (x) = 0 "Thanks a lot in advance!
Sample data:
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
source to share
The easiest way is to accommodate a nested model with interactions rather than two separate models. So you can first generate a coefficient that encodes two segments:
fac <- factor(as.numeric(time(zoop) > as.Date("2005-01-24")))
fac <- zoo(fac, time(zoop))
And then you can put a model where all the coefficients must be equal ( M0
) and where they differ ( M1
). The latter is equivalent to the separate models M1
and m2
that you put above:
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"))
Then you can easily check all coefficients separately for differences by looking
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
The line fac1:d(z)
has a slope difference test d(z)
. And if you want to test all coefficients together, you can do:
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
For complete reproducibility, I include the code I used to read the data:
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)
source to share