Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fixes mtest to pick up the correct 'B'; add many test cases/references; revert part of ab604aa4d75e9776dd3519149000d5c89291bfdc that claims to have fixed something (it did not fix) #64

Merged
merged 12 commits into from
Jan 5, 2025
1 change: 0 additions & 1 deletion NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@ subtitle: plm - Linear Models for Panel Data - A set of estimators and tests for

### Fixes:
* `pgmm`:
* for one-step GMM models: usual (non-robust) covariance matrix/standard errors fixed.
* argument `fsm` can now be set to influence the first step's weighing matrix.
* `vcovXX`: FD models with only one observation per group prior to
first-differencing errored ([#58](https://github.com/ycroissant/plm/issues/58)).
Expand Down
17 changes: 7 additions & 10 deletions R/est_gmm.R
Original file line number Diff line number Diff line change
Expand Up @@ -556,17 +556,11 @@ pgmm <- function(formula, data, subset, na.action,

residuals <- lapply(yX, function(x)
as.vector(x[ , 1L] - crossprod(t(x[ , -1L, drop = FALSE]), coefficients)))
outresid <- lapply(residuals, function(x) outer(x, x))

# non-robust variance-covariance matrix of one-step GMM:
# see Doornik/Arellano/Bond (2012), p. 31 (formula for V^hat1 with sig2 as in (4) on p. 30)
CPresid <- crossprod(unlist(residuals))
sig2 <- as.numeric(CPresid / (pdim$nT$N - NCOL(B1)))
vcov <- sig2 * B1


# A2 is also needed for "onestep" model in vcovHC.pgmm, hence calc. here and
# always include in model object
A2 <- mapply(function(x, y) crossprod(t(crossprod(x, y)), x), W, outresid, SIMPLIFY = FALSE)
A2 <- mapply(function(w, res) tcrossprod(crossprod(w, tcrossprod(res)), t(w)), W, residuals, SIMPLIFY = FALSE) # == mapply(function(w, res) t(w) %*% tcrossprod(res) %*% w, W, residuals, SIMPLIFY = FALSE)
A2 <- Reduce("+", A2)
minevA2 <- min(eigen(A2)$values)
A2 <- if (minevA2 < eps) {
Expand All @@ -578,7 +572,7 @@ pgmm <- function(formula, data, subset, na.action,
coef1s <- coefficients
t.CP.WX.A2 <- t(crossprod(WX, A2))
Y2 <- crossprod(t.CP.WX.A2, Wy)
vcov <- solve(crossprod(WX, t.CP.WX.A2))
vcov <- solve(crossprod(WX, t.CP.WX.A2)) # "B2"
coef2s <- as.numeric(crossprod(vcov, Y2))
names(coef2s) <- names.coef
coefficients <- list("step1" = coef1s, "step2" = coef2s)
Expand All @@ -589,6 +583,8 @@ pgmm <- function(formula, data, subset, na.action,
z <- as.vector(x[ , 1L] - crossprod(t(x[ , -1L, drop = FALSE]), coef2s))
names(z) <- nz
z})
} else {
vcov <- B1
}

rownames(vcov) <- colnames(vcov) <- names.coef
Expand Down Expand Up @@ -900,14 +896,15 @@ mtest.pgmm <- function(object, order = 1L, vcov = NULL, ...) {
X <- lapply(object$model, function(x) x[ , -1L, drop = FALSE])
W <- object$W
A <- if(model == "onestep") object$A1 else object$A2
B <- object$vcov # object$vcov is "B1" for one-step and "B2" for two-steps model
EX <- Reduce("+", mapply(crossprod, residl, X, SIMPLIFY = FALSE))
XZ <- Reduce("+", mapply(crossprod, W, X, SIMPLIFY = FALSE))
V <- mapply(tcrossprod, resid, SIMPLIFY = FALSE)
EVE <- Reduce("+", mapply(function(v, e) t(e) %*% v %*% e, V, residl, SIMPLIFY = FALSE))
ZVE <- Reduce("+", mapply(function(w, v, e) t(w) %*% v %*% e, W, V, residl, SIMPLIFY = FALSE))

num <- Reduce("+", mapply(crossprod, resid, residl, SIMPLIFY = FALSE))
denom <- EVE - 2 * EX %*% vcov(object) %*% t(XZ) %*% A %*% ZVE + EX %*% vv %*% t(EX)
denom <- EVE - 2 * EX %*% B %*% t(XZ) %*% A %*% ZVE + EX %*% vv %*% t(EX)
stat <- as.numeric(num / sqrt(denom))
names(stat) <- "normal"
if(!is.null(vcov)) vcov <- paste0(", vcov: ", deparse(substitute(vcov)))
Expand Down
7 changes: 4 additions & 3 deletions R/tool_vcovG.R
Original file line number Diff line number Diff line change
Expand Up @@ -1225,7 +1225,7 @@ vcovHC.pgmm <- function(x, ...) {
transformation <- describe(x, "transformation")
A1 <- x$A1
A2 <- x$A2
B1 <- x$B1
B1 <- x$B1 # needs to be B1 (from one-step model)

if(transformation == "ld") {
## yX <- lapply(x$model,function(x) rbind(diff(x),x))
Expand Down Expand Up @@ -1268,11 +1268,12 @@ vcovHC.pgmm <- function(x, ...) {
wexkw <- Reduce("+", mapply(function(x, y)
crossprod(x, crossprod(y, x)),
x$W, exk, SIMPLIFY = FALSE))
B2 <- vcov(x)
B2 <- x$vcov # is "B2" for two-step model
Dk <- -B2 %*% t(WX) %*% A2 %*% wexkw %*% A2 %*% We
D <- cbind(D, Dk)
}
# robust vcov for twosteps GMM model, Roodman (2019) form. (18)
# Windmeijer (2005) small sample bias correction for twosteps GMM model,
# see Windmeijer (2005), p. 33 formula (3.3); Roodman (2019) form. (18)
vcovr2s <- B2 + crossprod(t(D), B2) + t(crossprod(t(D), B2)) + D %*% vcovr1s %*% t(D)
}
if(model == "twosteps") vcovr2s else vcovr1s
Expand Down
12 changes: 6 additions & 6 deletions inst/tests/test_misc.Rout.save
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ alternative hypothesis: autocorrelation present
Arellano-Bond autocorrelation test of degree 2, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -0.36744, p-value = 0.7133
normal = -0.35166, p-value = 0.7251
alternative hypothesis: autocorrelation present

>
Expand Down Expand Up @@ -1493,8 +1493,8 @@ lag(log(capital), 0:1)1 -0.424393 0.058479 -7.2572 3.952e-13 ***
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sargan test: chisq(100) = 118.763 (p-value = 0.097096)
Autocorrelation test (1): normal = -5.519159 (p-value = 3.4063e-08)
Autocorrelation test (2): normal = -0.2524777 (p-value = 0.80067)
Autocorrelation test (1): normal = -4.808434 (p-value = 1.5212e-06)
Autocorrelation test (2): normal = -0.2800133 (p-value = 0.77947)
Wald test for coefficients: chisq(5) = 11174.82 (p-value = < 2.22e-16)
Wald test for time dummies: chisq(7) = 14.71138 (p-value = 0.039882)
>
Expand Down Expand Up @@ -1579,8 +1579,8 @@ lag(log(capital), 1) -0.424393 0.058479 -7.2572 3.952e-13 ***
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sargan test: chisq(100) = 118.763 (p-value = 0.097096)
Autocorrelation test (1): normal = -5.519159 (p-value = 3.4063e-08)
Autocorrelation test (2): normal = -0.2524777 (p-value = 0.80067)
Autocorrelation test (1): normal = -4.808434 (p-value = 1.5212e-06)
Autocorrelation test (2): normal = -0.2800133 (p-value = 0.77947)
Wald test for coefficients: chisq(5) = 11174.82 (p-value = < 2.22e-16)
Wald test for time dummies: chisq(7) = 14.71138 (p-value = 0.039882)
>
Expand Down Expand Up @@ -1891,4 +1891,4 @@ Coefficients:
>
> proc.time()
user system elapsed
9.45 0.90 10.32
9.40 1.07 10.51
36 changes: 20 additions & 16 deletions inst/tests/test_mtest.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,19 +9,24 @@ ab.a1 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
+ lag(log(capital), 0:2) + lag(log(output), 0:2) | lag(log(emp), 2:99),
data = EmplUK, effect = "twoways", model = "onestep")

mtest(ab.a1, 1, vcov = NULL)
mtest(ab.a1, 2, vcov = NULL)
mtest(ab.a1, 1, vcov = vcovHC)
mtest(ab.a1, 2, vcov = vcovHC) # -0.516 (reference)
mtest(ab.a1, 1, vcov = NULL) # -3.9394 (p = 0.0001) xtabond manual, example 1
mtest(ab.a1, 2, vcov = NULL) # -0.54239 (p = 0.5876) xtabond manual, example 1
mtest(ab.a1, 1, vcov = vcovHC) # -3.600 (p = 0.000) (reference DPD/Ox abest3.out) / -3.5996 (p = 0.0003) xtabond manual, example 2
mtest(ab.a1, 2, vcov = vcovHC) # -0.516 (p = 0.606) (reference A/B and DPD/Ox abest3.out) / -0.51603 (p = 0.6058) xtabond manual, example 2

# Windmeijer (2025), table 2, onestep with corrected std. err
wind.s1 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
# A/B (1991) model from col. b BUT as one-step model
ab.b.onestep <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
+ log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99),
data = EmplUK, effect = "twoways", model = "onestep")
mtest(wind.s1, 1, vcov = NULL)
mtest(wind.s1, 2, vcov = NULL)
mtest(wind.s1, 1, vcov = vcovHC) # -2.493 (reference)
mtest(wind.s1, 2, vcov = vcovHC) # -0.359 (reference)
# NB: coef match, non-robust std. errors do not match abest1.out

mtest(ab.b.onestep, 1, vcov = NULL) # -3.409 (p = 0.001) (reference in DPD/Ox abest1.out)
mtest(ab.b.onestep, 2, vcov = NULL) # -0.3695 (p = 0.712) (reference in DPD/Ox abest1.out)
mtest(ab.b.onestep, 1, vcov = vcovHC) # -2.493 (reference in Windmeijer) / -2.493 (p = 0.013) DPD/Ox abest1.out
mtest(ab.b.onestep, 2, vcov = vcovHC) # -0.359 (reference in Windmeijer) / -0.3594 (p = 0.719) DPD/Ox abest1.out



########################## two-steps ######################

Expand All @@ -30,18 +35,17 @@ ab.a2 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
+ lag(log(capital), 0:2) + lag(log(output), 0:2) | lag(log(emp), 2:99),
data = EmplUK, effect = "twoways", model = "twosteps")

mtest(ab.a2, 1, vcov = NULL)
mtest(ab.a2, 2, vcov = NULL) # -0.434 (reference)
mtest(ab.a2, 1, vcov = vcovHC)
mtest(ab.a2, 2, vcov = vcovHC)
mtest(ab.a2, 1, vcov = NULL) # -3.000 (p = 0.003) (reference DPD/Ox in abest3.out)
mtest(ab.a2, 2, vcov = NULL) # -0.434 (reference A/B) / DPD/Ox has -0.4158 (p = 0.678) in abest3.out, see comment in footnote 7, p. 32
mtest(ab.a2, 1, vcov = vcovHC) # -2.1255 (p = 0.0335) xtabond manual, example 4
mtest(ab.a2, 2, vcov = vcovHC) # -0.35166 (p = 0.7251) xtabond manual, example 4

## Arellano and Bond (1991), table 4 col. b / Windmeijer (2025), table 2
ab.b <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
+ log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99),
data = EmplUK, effect = "twoways", model = "twosteps")

mtest(ab.b, 1, vcov = NULL) # -2.826 (reference about A/B's m1 in Windmeijer)
mtest(ab.b, 2, vcov = NULL) # -0.327 (reference in A/B)
mtest(ab.b, 1, vcov = NULL) # -2.826 (reference about A/B's m1 in Windmeijer) / DPD/Ox has -2.428 (p = 0.015) in abest1.out
mtest(ab.b, 2, vcov = NULL) # -0.327 (reference in A/B and in Windmeijer) / DPD/Ox has -0.3325 (p = 0.739) in abest1.out, explanation in footnote 7 on p. 32
mtest(ab.b, 1, vcov = vcovHC) # -1.999 (reference in Windmeijer)
mtest(ab.b, 2, vcov = vcovHC) # -0.316 (reference in Windmeijer)

64 changes: 30 additions & 34 deletions inst/tests/test_mtest.Rout.save
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ### test results with plm 2.6-4
> ### test results with plm 2.6-4 (none of the literature's results were reproduced)
> library(plm)
> data("EmplUK", package = "plm")
>
Expand All @@ -26,83 +26,80 @@ Type 'q()' to quit R.
+ + lag(log(capital), 0:2) + lag(log(output), 0:2) | lag(log(emp), 2:99),
+ data = EmplUK, effect = "twoways", model = "onestep")
>
> mtest(ab.a1, 1, vcov = NULL)
> mtest(ab.a1, 1, vcov = NULL) # -3.9394 (p = 0.0001) xtabond manual, example 1

Arellano-Bond autocorrelation test of degree 1

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -4.1906, p-value = 2.782e-05
alternative hypothesis: autocorrelation present

> mtest(ab.a1, 2, vcov = NULL)
> mtest(ab.a1, 2, vcov = NULL) # -0.54239 (p = 0.5876) xtabond manual, example 1

Arellano-Bond autocorrelation test of degree 2

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -0.64241, p-value = 0.5206
alternative hypothesis: autocorrelation present

> mtest(ab.a1, 1, vcov = vcovHC)
> mtest(ab.a1, 1, vcov = vcovHC) # -3.600 (p = 0.000) (reference DPD/Ox abest3.out) / -3.5996 (p = 0.0003) xtabond manual, example 2

Arellano-Bond autocorrelation test of degree 1, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = NaN, p-value = NA
normal = -3.5996, p-value = 0.0003187
alternative hypothesis: autocorrelation present

Warning message:
In sqrt(denom) : NaNs produced
> mtest(ab.a1, 2, vcov = vcovHC) # -0.516 (reference)
> mtest(ab.a1, 2, vcov = vcovHC) # -0.516 (p = 0.606) (reference A/B and DPD/Ox abest3.out) / -0.51603 (p = 0.6058) xtabond manual, example 2

Arellano-Bond autocorrelation test of degree 2, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = NaN, p-value = NA
normal = -0.51603, p-value = 0.6058
alternative hypothesis: autocorrelation present

Warning message:
In sqrt(denom) : NaNs produced
>
> # Windmeijer (2025), table 2, onestep with corrected std. err
> wind.s1 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
> # A/B (1991) model from col. b BUT as one-step model
> ab.b.onestep <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
+ + log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99),
+ data = EmplUK, effect = "twoways", model = "onestep")
> mtest(wind.s1, 1, vcov = NULL)
> # NB: coef match, non-robust std. errors do not match abest1.out
>
> mtest(ab.b.onestep, 1, vcov = NULL) # -3.409 (p = 0.001) (reference in DPD/Ox abest1.out)

Arellano-Bond autocorrelation test of degree 1

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = -3.9009, p-value = 9.583e-05
alternative hypothesis: autocorrelation present

> mtest(wind.s1, 2, vcov = NULL)
> mtest(ab.b.onestep, 2, vcov = NULL) # -0.3695 (p = 0.712) (reference in DPD/Ox abest1.out)

Arellano-Bond autocorrelation test of degree 2

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = -0.46764, p-value = 0.64
alternative hypothesis: autocorrelation present

> mtest(wind.s1, 1, vcov = vcovHC) # -2.493 (reference)
> mtest(ab.b.onestep, 1, vcov = vcovHC) # -2.493 (reference in Windmeijer) / -2.493 (p = 0.013) DPD/Ox abest1.out

Arellano-Bond autocorrelation test of degree 1, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = NaN, p-value = NA
normal = -2.4934, p-value = 0.01265
alternative hypothesis: autocorrelation present

Warning message:
In sqrt(denom) : NaNs produced
> mtest(wind.s1, 2, vcov = vcovHC) # -0.359 (reference)
> mtest(ab.b.onestep, 2, vcov = vcovHC) # -0.359 (reference in Windmeijer) / -0.3594 (p = 0.719) DPD/Ox abest1.out

Arellano-Bond autocorrelation test of degree 2, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = NaN, p-value = NA
normal = -0.35945, p-value = 0.7193
alternative hypothesis: autocorrelation present

Warning message:
In sqrt(denom) : NaNs produced
>
>
>
> ########################## two-steps ######################
>
Expand All @@ -111,36 +108,36 @@ In sqrt(denom) : NaNs produced
+ + lag(log(capital), 0:2) + lag(log(output), 0:2) | lag(log(emp), 2:99),
+ data = EmplUK, effect = "twoways", model = "twosteps")
>
> mtest(ab.a2, 1, vcov = NULL)
> mtest(ab.a2, 1, vcov = NULL) # -3.000 (p = 0.003) (reference DPD/Ox in abest3.out)

Arellano-Bond autocorrelation test of degree 1

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -2.9998, p-value = 0.002702
alternative hypothesis: autocorrelation present

> mtest(ab.a2, 2, vcov = NULL) # -0.434 (reference)
> mtest(ab.a2, 2, vcov = NULL) # -0.434 (reference A/B) / DPD/Ox has -0.4158 (p = 0.678) in abest3.out, see comment in footnote 7, p. 32

Arellano-Bond autocorrelation test of degree 2

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -0.41575, p-value = 0.6776
alternative hypothesis: autocorrelation present

> mtest(ab.a2, 1, vcov = vcovHC)
> mtest(ab.a2, 1, vcov = vcovHC) # -2.1255 (p = 0.0335) xtabond manual, example 4

Arellano-Bond autocorrelation test of degree 1, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -2.3456, p-value = 0.01899
normal = -2.1255, p-value = 0.03355
alternative hypothesis: autocorrelation present

> mtest(ab.a2, 2, vcov = vcovHC)
> mtest(ab.a2, 2, vcov = vcovHC) # -0.35166 (p = 0.7251) xtabond manual, example 4

Arellano-Bond autocorrelation test of degree 2, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + lag(log(capital), ...
normal = -0.36744, p-value = 0.7133
normal = -0.35166, p-value = 0.7251
alternative hypothesis: autocorrelation present

>
Expand All @@ -149,15 +146,15 @@ alternative hypothesis: autocorrelation present
+ + log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99),
+ data = EmplUK, effect = "twoways", model = "twosteps")
>
> mtest(ab.b, 1, vcov = NULL) # -2.826 (reference about A/B's m1 in Windmeijer)
> mtest(ab.b, 1, vcov = NULL) # -2.826 (reference about A/B's m1 in Windmeijer) / DPD/Ox has -2.428 (p = 0.015) in abest1.out

Arellano-Bond autocorrelation test of degree 1

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = -2.4278, p-value = 0.01519
alternative hypothesis: autocorrelation present

> mtest(ab.b, 2, vcov = NULL) # -0.327 (reference in A/B)
> mtest(ab.b, 2, vcov = NULL) # -0.327 (reference in A/B and in Windmeijer) / DPD/Ox has -0.3325 (p = 0.739) in abest1.out, explanation in footnote 7 on p. 32

Arellano-Bond autocorrelation test of degree 2

Expand All @@ -170,19 +167,18 @@ alternative hypothesis: autocorrelation present
Arellano-Bond autocorrelation test of degree 1, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = -1.5357, p-value = 0.1246
normal = -1.5385, p-value = 0.1239
alternative hypothesis: autocorrelation present

> mtest(ab.b, 2, vcov = vcovHC) # -0.316 (reference in Windmeijer)

Arellano-Bond autocorrelation test of degree 2, vcov: vcovHC

data: log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + log(capital) + ...
normal = -0.30388, p-value = 0.7612
normal = -0.27968, p-value = 0.7797
alternative hypothesis: autocorrelation present

>
>
> proc.time()
user system elapsed
3.89 0.48 4.35
2.39 0.43 2.76
7 changes: 4 additions & 3 deletions inst/tests/test_pgmm.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ ab.a1 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1)
+ lag(log(capital), 0:2) + lag(log(output), 0:2) | lag(log(emp), 2:99),
data = EmplUK, effect = "twoways", model = "onestep")

(s.ab.a1 <- summary(ab.a1, robust = FALSE))
(s.ab.a1 <- summary(ab.a1, robust = FALSE)) # xtabond manual, example 1, has slightly different values for non-robust standard errors
(s.ab.a1r <- summary(ab.a1, robust = TRUE)) # as tabulated by Arellano/Bond

### Arellano-Bond test, Sargan test and Wald test yield different results for onestep
Expand Down Expand Up @@ -291,8 +291,9 @@ if(pdynmc.avail) {
# check coefficients
stopifnot(isTRUE(all.equal(round( s.ab.a1[[1L]][ , 1], 5), round( ab.a1.pdynmc$coefficients[1:10], 5), check.attributes = FALSE)))

# check standard errors (non-robust and robust)
stopifnot(isTRUE(all.equal(round( s.ab.a1[[1L]][ , 2], 5), round( ab.a1.pdynmc$stderr$step1[1:10], 5), check.attributes = FALSE)))
# check standard errors (non-robust and robust)
# pdynmc 2025-01-05 does not seem to produce correct non-robust SEs in one-step case
# stopifnot(isTRUE(all.equal(round( s.ab.a1[[1L]][ , 2], 5), round( ab.a1.pdynmc$stderr$step1[1:10], 5), check.attributes = FALSE)))
stopifnot(isTRUE(all.equal(round(s.ab.a1r[[1L]][ , 2], 5), round(ab.a1r.pdynmc$stderr$step1[1:10], 5), check.attributes = FALSE)))

## two-steps GMM
Expand Down
Loading
Loading