diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 782da4d8dc..9447f51d87 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -29,7 +29,7 @@ sqrt, squeeze, ) -from numpy.linalg import eigvalsh, inv, matrix_rank, pinv +from numpy.linalg import eigvalsh, inv, matrix_rank, pinv, solve from pandas import DataFrame, Series, concat from scipy.optimize import minimize @@ -1044,10 +1044,44 @@ def _gmm_post_estimation( "weight_config": self._weight.config, "iterations": iters, "j_stat": self._j_statistic(params, weight_mat), + "centered_j_stat": self._centered_j_statistic(params), } return gmm_specific + def _centered_j_statistic( + self, params: linearmodels.typing.data.Float64Array + ) -> WaldTestStatistic | InvalidTestStatistic: + """ + J statistic that uses a centered estimate of the moment covariance + + The covariance of the moment conditions is estimated using the same + estimator as the model, i.e., the one selected by ``weight_type`` and + ``weight_config``, except that the mean of the moment conditions is + always subtracted, and it is evaluated at ``params``. + """ + y, x, z = self._wy, self._wx, self._wz + nobs, nvar, ninstr = y.shape[0], x.shape[1], z.shape[1] + name = "Centered J-test of overidentifying restrictions" + eps = y - x @ params + g_bar = (z * eps).mean(0) + + weight_matrix_estimator = WEIGHT_MATRICES[self._weight_type] + centered = weight_matrix_estimator(**{**self._weight_kwargs, "center": True}) + s_c = centered.weight_matrix(x, z, eps) + if matrix_rank(s_c, hermitian=True) < ninstr: + # The centered cluster sums add to zero, so a clustered estimator + # has rank at most one less than the number of clusters + reason = ( + "The centered covariance of the moment conditions is singular, " + "so the statistic is not defined. This occurs, for example, " + "when there are too few clusters." + ) + return InvalidTestStatistic(reason, name=name) + stat = float(nobs * g_bar @ solve(s_c, g_bar)) + null = "Expected moment conditions are equal to 0" + return WaldTestStatistic(stat, null, ninstr - nvar, name=name) + def _j_statistic( self, params: linearmodels.typing.data.Float64Array, @@ -1337,6 +1371,7 @@ def _gmm_post_estimation( "weight_config": self._weight.config, "iterations": iters, "j_stat": self._j_statistic(params, weight_mat), + "centered_j_stat": self._centered_j_statistic(params), } return gmm_specific diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index 78c7baa6fa..0f3c1593a5 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -1589,6 +1589,7 @@ def __init__( self._weight_config = results["weight_config"] self._iterations = results["iterations"] self._j_stat = results["j_stat"] + self._centered_j_stat = results["centered_j_stat"] @property def weight_matrix(self) -> linearmodels.typing.data.Float64Array: @@ -1648,6 +1649,75 @@ def j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: """ return self._j_stat + @property + def centered_j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: + r""" + J-test of overidentifying restrictions using a centered moment covariance + + Returns + ------- + WaldTestStatistic + J statistic test of overidentifying restrictions that uses a + covariance of the moment conditions that has been centered. + An InvalidTestStatistic is returned if that covariance is + singular, which can happen if there are too few clusters. + + Notes + ----- + The statistic is defined as + + .. math :: + + n \bar{g}'\hat{S}_c^{-1}\bar{g} \sim \chi^2_q + + where :math:`\bar{g} = n^{-1}\sum \hat{\epsilon}_i z_i`, + :math:`z_i` includes both the exogenous variables and instruments, + :math:`\hat{\epsilon}_i` are the model residuals and the degree of + freedom is :math:`q = n_{instr} - n_{endog}`. :math:`\hat{S}_c` is the + estimator of the variance of :math:`\sqrt{n}\bar{g}` selected by + ``weight_type`` and ``weight_config``, applied to the centered moment + conditions :math:`g_i - \bar{g}` and evaluated at the final parameter + estimates. For ``weight_type="robust"`` it is + + .. math :: + + \hat{S}_c = n^{-1}\sum_{i=1}^{n}(g_i - \bar{g})(g_i - \bar{g})' + + :math:`\hat{S}_c` and the :math:`W` of :attr:`j_stat` estimate the same + quantity when the model is correctly specified, and so the + statistics have the same asymptotic :math:`\chi^2_q` null + distribution. If the model is misspecified, :math:`E[g_i] \neq 0`, + both statistics are consistent tests and grow at rate :math:`n`, but + the uncentered estimator also contains :math:`\bar{g}\bar{g}'`. When + it is evaluated using weights from the final parameter estimates, as + in iterated GMM, this bounds the J statistic by :math:`n`, however + large the violation of the moment conditions, whereas the centered + statistic is not bounded. For ``weight_type="robust"`` and converged + iterated GMM, the two are related by + :math:`n \bar{g}'\hat{S}_c^{-1}\bar{g} = J / (1 - J / n)`, where + :math:`J` is the statistic in :attr:`j_stat`. This relationship is + :math:`J / (1 - s J / n)` if the estimator is scaled by + :math:`s = n / (n - k)` because ``debiased`` is ``True``. The + difference matters when :math:`J / n` is not small. + + Centering has no effect when ``weight_type`` is "unadjusted", + since that estimator uses the residual variance about its mean, and + the statistic then coincides with :attr:`j_stat` up to the step of + the estimation at which the weight matrix is evaluated. + It also has no effect in :class:`~linearmodels.iv.model.IVGMMCUE` + models unless ``center=False`` was set. + + References + ---------- + Hall, A. R. (2000). Covariance matrix estimation and the power of the + overidentifying restrictions test. *Econometrica*, 68(6), 1517-1528. + + Hansen, B. E. and Lee, S. (2021). Inference for iterated GMM under + misspecification. *Econometrica*, 89(3), 1419-1447. Theorem 1 and + Example 1 discuss the centered weight matrix. + """ + return self._centered_j_stat + def c_stat(self, variables: list[str] | str | None = None) -> WaldTestStatistic: r""" C-test of endogeneity diff --git a/linearmodels/tests/iv/results/centered-j-package-check.R b/linearmodels/tests/iv/results/centered-j-package-check.R new file mode 100644 index 0000000000..06af8ab945 --- /dev/null +++ b/linearmodels/tests/iv/results/centered-j-package-check.R @@ -0,0 +1,225 @@ +# Cross-check of centered-j-reference.csv against the R packages gmm and AER. +# +# Run from this directory: Rscript centered-j-package-check.R +# Requires: gmm (tested with 1.9.1) and AER (tested with 1.2.17) +# Output: a table of the agreement, and an error if the +# values in the file do not agree with the packages +# +# centered-j-reference.R, which produces the reference values, is written in +# base R so that nothing is shared with linearmodels. This script is not part +# of the tests and does not write anything. It checks that the reference values +# are what established R implementations compute, to the precision that they +# have, which is limited by their numerical optimizers and not by the +# statistics. +# +# * gmm::gmm with type="iterative" is the iterated GMM estimator. Its +# centeredVcov argument selects whether the covariance of the moment +# conditions is centered, and J is n times the objective. The J of a model +# whose estimator is not centered is gmm with centeredVcov=FALSE. The +# centered statistic is gmm::evalGmm with centeredVcov=TRUE evaluated at the +# estimates, with the weight matrix evaluated at the same estimates. +# * gmm::gmm with type="cue" is the continuously updating estimator. +# * vcov="iid" is the covariance of the moment conditions without any +# assumption on their variance, i.e., the heteroskedasticity robust +# estimator (evalGmm does not have vcov="MDS"), and vcov="HAC" with +# kernel="Bartlett" and prewhite=0 the kernel estimator. The +# bandwidth bw of gmm is the number of lags plus one, so that bw=5 is +# the 4 lags of linearmodels with weights 1 - j / 5. +# * AER::ivreg Sargan test is the J statistic with the unadjusted estimator +# of the covariance of the moments. It is the 2SLS J of the model, which +# is the iterated GMM estimator with that weight. +# +# What is not checked is what the packages do not offer: gmm does not have a +# one-step weight matrix that is the same as linearmodels' (so the two-step +# estimates are not compared, their first step is identity weighting), clustered +# covariances, or the small-sample scale of debiased estimators. Those are +# checked against theory by centered-j-reference.R. + +suppressMessages({ + library(foreign) + library(gmm) + library(AER) +}) + +ref <- read.csv("centered-j-reference.csv", stringsAsFactors = FALSE) +sim <- read.dta("simulated-data.dta") +sim$const <- 1 +hs <- read.csv("housing.csv") +hs$const <- 1 +reg <- model.matrix(~ region, hs)[, -1] +hs <- cbind(hs, reg) +mis <- read.csv("misspecified-data.csv") +mis$const <- 1 + +words <- function(s) strsplit(s, " ", fixed = TRUE)[[1]] + +# Data of a row of the reference table, weighted as linearmodels weights it +model_data <- function(r) { + d <- switch(r$dataset, sim = sim, housing = hs, mis = mis) + sw <- if (r$weighted) sqrt(d$weights / mean(d$weights)) else rep(1, nrow(d)) + exog <- words(r$exog) + endog <- words(r$endog) + instr <- words(r$instruments) + list(y = d[[r$dependent]] * sw, + X = as.matrix(d[, c(exog, endog), drop = FALSE]) * sw, + Z = as.matrix(d[, c(exog, instr), drop = FALSE]) * sw) +} + +# The moment function of gmm, with the data in a single matrix [y, X, Z] +moments <- function(k) { + function(theta, dat) { + e <- dat[, 1] - dat[, 2:(k + 1)] %*% theta + dat[, -(1:(k + 1)), drop = FALSE] * drop(e) + } +} + +gmm_args <- function(r) { + if (r$weight_type == "kernel") { + list(vcov = "HAC", kernel = "Bartlett", bw = as.integer(r$bandwidth) + 1, + prewhite = 0) + } else { + list(vcov = "iid") + } +} + +# Fit by gmm, returning J and the estimates +fit_gmm <- function(r, md, centered, type) { + k <- ncol(md$X) + dat <- cbind(md$y, md$X, md$Z) + n <- nrow(dat) + g <- moments(k) + # 2SLS starting values + start <- drop(solve(crossprod(md$X, md$Z) %*% solve(crossprod(md$Z)) %*% + crossprod(md$Z, md$X), + crossprod(md$X, md$Z) %*% solve(crossprod(md$Z)) %*% + crossprod(md$Z, md$y))) + run <- function(t0, type) { + args <- c(list(g = g, x = dat, t0 = t0, type = type, wmatrix = "optimal", + centeredVcov = centered, crit = 1e-12), + gmm_args(r)) + if (type == "iterative") args$itermax <- 1000 + do.call(gmm, args) + } + fit <- run(start, type) + if (type == "cue") { + # The CUE objective is hard to minimize. Use the best of two starts, the + # 2SLS estimates and the iterated GMM estimates + other <- run(coef(run(start, "iterative")), "cue") + if (other$objective < fit$objective) fit <- other + } + list(b = unname(coef(fit)), J = unname(fit$objective * n), n = n, dat = dat, g = g) +} + +# Objective of gmm times n at the parameters theta, with the weight matrix +# evaluated at theta_w, and with or without centering of the covariance of the +# moment conditions +eval_at <- function(r, md, theta, theta_w, centered) { + k <- ncol(md$X) + dat <- cbind(md$y, md$X, md$Z) + args <- c(list(g = moments(k), x = dat, t0 = theta, tetw = theta_w, + wmatrix = "optimal", centeredVcov = centered), gmm_args(r)) + unname(do.call(evalGmm, args)$objective * nrow(dat)) +} + +# Sargan statistic of the 2SLS fit, which is the J of the unadjusted estimator +sargan <- function(r, md) { + dat <- data.frame(y = md$y) + xcols <- colnames(md$X) + zcols <- colnames(md$Z) + d <- cbind(dat, as.data.frame(md$X), as.data.frame(md$Z[, setdiff(zcols, xcols), drop = FALSE])) + exog <- words(r$exog) + form <- as.formula(sprintf("y ~ 0 + %s | 0 + %s", + paste(xcols, collapse = " + "), + paste(zcols, collapse = " + "))) + fit <- ivreg(form, data = d) + unname(summary(fit, diagnostics = TRUE)$diagnostics["Sargan", "statistic"]) +} + +check <- function(label, got, want, rtol) { + err <- abs(got - want) / abs(want) + if (!is.finite(err) || err > rtol) { + stop(sprintf("%s: %.10g (package) vs %.10g (reference), relative difference %.2e", + label, got, want, err)) + } + err +} + +part_a <- list() +part_b <- list() +for (i in seq_len(nrow(ref))) { + r <- ref[i, ] + if (r$weight_type == "clustered" || r$debiased || is.na(r$centered_j)) next + md <- model_data(r) + label <- sprintf("%s [%s]", r$id, r$method) + + if (r$weight_type == "unadjusted") { + # 2SLS is the iterated GMM estimator with this weight, and its J is Sargan's + if (r$method == "two_step") next + s <- sargan(r, md) + part_a[[label]] <- c(j = check(paste(label, "J vs Sargan"), s, r$j, 1e-6), + centered_j = check(paste(label, "centered J vs Sargan"), s, + r$centered_j, 1e-6)) + next + } + + # A. The statistics at the estimates of the reference, which are the same + # parameters for both, so that only the statistics are compared. J uses the + # weight matrix that was used to estimate the parameters + theta <- as.numeric(words(r$params)) + theta_w <- as.numeric(words(r$weight_params)) + part_a[[label]] <- c( + j = check(paste(label, "J"), eval_at(r, md, theta, theta_w, r$center), r$j, 1e-8), + centered_j = check(paste(label, "centered J"), + eval_at(r, md, theta, theta, TRUE), r$centered_j, 1e-8) + ) + + # B. The estimators of the package. gmm does not have a two-step estimator + # that starts as linearmodels does, and its CUE is only compared for the + # robust estimator + if (r$method == "two_step" || (r$method == "cue" && r$weight_type != "robust")) next + fit <- fit_gmm(r, md, r$center, c(iterated = "iterative", cue = "cue")[[r$method]]) + if (r$method == "cue") { + # The reference value is a minimum, so that the package cannot be lower. + # The numerical optimizer of the package does not always find the minimum, + # in which case there is nothing to compare + if (fit$J < r$j * (1 - 1e-6)) { + stop(sprintf("%s: the package finds J=%.8g, below the reference minimum %.8g", + label, fit$J, r$j)) + } + if (fit$J > r$j * (1 + 5e-4)) { + cat(sprintf("note: %s: gmm stops at J=%.6g, above the minimum %.6g; not compared +", + label, fit$J, r$j)) + next + } + } + part_b[[label]] <- c( + j = check(paste(label, "J of the estimator"), fit$J, r$j, 5e-4), + params = max(abs(fit$b - theta) / pmax(abs(theta), 1e-2)) + ) + # The objective is flat in some directions, so that the estimates of the + # package are less precise than its J + if (part_b[[label]]["params"] >= 2e-2) { + stop(sprintf("%s: the estimates of the package differ by %.2e (relative) from the reference", + label, part_b[[label]]["params"])) + } +} + +cat(" +A. Statistics at the parameters of the reference (relative differences) +") +res_a <- do.call(rbind, part_a) +print(signif(res_a, 3)) +cat(" +B. Estimators of the package (relative differences in J and in the parameters) +") +res_b <- do.call(rbind, part_b) +print(signif(res_b, 3)) +cat(sprintf(paste0(" +A: %d models, largest relative difference %.2e in J and %.2e in ", + "the centered J +B: %d models, largest relative difference %.2e in J and ", + "%.2e in the parameters +"), + nrow(res_a), max(res_a[, "j"]), max(res_a[, "centered_j"]), + nrow(res_b), max(res_b[, "j"]), max(res_b[, "params"]))) diff --git a/linearmodels/tests/iv/results/centered-j-reference.R b/linearmodels/tests/iv/results/centered-j-reference.R new file mode 100644 index 0000000000..c739664263 --- /dev/null +++ b/linearmodels/tests/iv/results/centered-j-reference.R @@ -0,0 +1,423 @@ +# Reference values for IVGMMResults.centered_j_stat, the J statistic that uses +# a centered estimate of the covariance of the moment conditions. Written in +# base R so that nothing is shared with linearmodels. +# +# Run from this directory: Rscript centered-j-reference.R +# Output: centered-j-reference.csv (read by ../test_centered_j_stat.py) +# +# Definition (Hall 2000, Econometrica 68(6), 1517-1528; Hansen and Lee 2021, +# Econometrica 89(3), 1419-1447, Theorem 1 and Example 1): with moment +# conditions g_i = z_i e_i(b), where e_i(b) = y_i - x_i b, and gbar their mean +# +# J = n gbar' S^{-1} gbar S = covariance of g_i, not centered +# Jc = n gbar' S_c^{-1} gbar S_c = the same estimator applied to g_i - gbar +# +# Both are evaluated at the same parameters b. The standard J of linearmodels +# (and of Stata's ivregress gmm) uses the weight matrix that was used in the +# last estimation step, which is evaluated at the previous iterate, while Jc +# evaluates S_c at the final estimate. They are the same at convergence of the +# iterated estimator when the estimator of S is the same. +# +# GMM conventions follow linearmodels.iv.IVGMM and Stata's ivregress gmm. The +# first step uses W = (Z'Z/n)^{-1}, later steps use W = S(b_prev)^{-1} where S +# is an estimate of the covariance of the moment conditions based on the +# residuals of the previous step. "two_step" is the estimator after one update +# of W (iter_limit=2 in linearmodels), "iterated" repeats this until the +# estimates converge, and "cue" minimizes n gbar(b)' S(b)^{-1} gbar(b) over b. +# Four estimators of S are used, which are the same as Stata's wmatrix() +# options: +# robust S = (1/n) sum (z_i e_i)(z_i e_i)' +# unadjusted S = s2 Z'Z/n, with s2 the variance of e about its mean +# kernel Bartlett, S = (1/n) [G_0 + sum_j (1 - j/(m+1)) (G_j + G_j')] +# with G_j = sum_t (z_t e_t)(z_{t-j} e_{t-j})' and bandwidth m +# clustered S = (1/n) sum_g (sum_{i in g} z_i e_i)(sum_{i in g} z_i e_i)' +# "center" subtracts the mean of z_i e_i first, which is what the model's own +# estimator does when center=TRUE and has no effect on unadjusted, where the +# estimator is always centered. "debiased" uses the small-sample scale of +# linearmodels: n/(n-k) and (n-1)/(n-k) G/(G-1) for the clustered estimator, +# with k the number of regressors. Weighted models are estimated on data +# multiplied by the square root of the weights. +# +# Stata has no centered J, so there are no values from Stata for it. Instead +# the script +# 1. checks that its GMM and S estimators reproduce numbers produced by +# Stata that are already used in the tests, ivregress gmm on the housing +# and simulated data (stata-iv-housing-results.txt and +# stata-iv-simulated-results.txt), for all four estimators of S, and +# 2. checks properties of the centered statistic that follow from theory +# for every scenario that they apply to before writing it: +# * Hansen and Lee (2021) Theorem 1: the iterated GMM estimator is the +# same if S is centered, for the robust estimator, +# * Jc = J / (1 - s J / n) for the iterated estimator with the robust +# uncentered estimator, which follows from S = s (S_c + gbar gbar') and +# the Sherman-Morrison formula and implies that J < n / s. s = 1 +# unless the estimator is debiased, when s = n / (n - k), +# * Jc = J for unadjusted, whose estimator is always centered, +# * the CUE estimates are minima of the objective function. +# The values written were also compared with the R packages gmm and AER by +# centered-j-package-check.R, which needs those packages and is not needed to +# run this script or the tests. At the estimates that this script computes the +# packages' statistics agree with the values to about 1e-14. +# The scenarios include strongly misspecified models, where J is far from the +# chi2 distribution and close to its upper bound n. + +suppressMessages(library(foreign)) + +## ---------------------------------------------------------------- GMM ---- +wt_spec <- function(type = "robust", center = FALSE, debiased = FALSE, + bandwidth = 0, clusters = NULL) { + list(type = type, center = center, debiased = debiased, + bandwidth = bandwidth, clusters = clusters) +} + +# Estimate of the covariance of the moment conditions z_i e_i. nvar is the +# number of regressors +S_hat <- function(Z, e, wt, nvar) { + n <- nrow(Z) + if (wt$type == "unadjusted") { + S <- mean((e - mean(e))^2) * crossprod(Z) / n + return(if (wt$debiased) S * n / (n - nvar) else S) + } + ze <- Z * e + if (wt$center) ze <- sweep(ze, 2, colMeans(ze)) + if (wt$type == "robust") { + S <- crossprod(ze) / n + if (wt$debiased) S <- S * n / (n - nvar) + } else if (wt$type == "kernel") { + S <- crossprod(ze) + for (j in seq_len(wt$bandwidth)) { + op <- crossprod(ze[(j + 1):n, , drop = FALSE], ze[1:(n - j), , drop = FALSE]) + S <- S + (1 - j / (wt$bandwidth + 1)) * (op + t(op)) + } + S <- S / n + if (wt$debiased) S <- S * n / (n - nvar) + } else if (wt$type == "clustered") { + S <- matrix(0, ncol(ze), ncol(ze)) + for (g in unique(wt$clusters)) { + zb <- colSums(ze[wt$clusters == g, , drop = FALSE]) + S <- S + tcrossprod(zb) + } + S <- S / n + if (wt$debiased) { + ng <- length(unique(wt$clusters)) + S <- S * (n - 1) / (n - nvar) * ng / (ng - 1) + } + } else { + stop("unknown weight type") + } + S +} + +gmm_beta <- function(y, X, Z, W) { + A <- crossprod(X, Z) + drop(solve(A %*% W %*% t(A), A %*% W %*% crossprod(Z, y))) +} + +gmm_J <- function(y, X, Z, b, W) { + n <- nrow(Z) + g <- crossprod(Z, y - X %*% b) / n + drop(n * crossprod(g, W %*% g)) +} + +# Estimator after one update of the weight matrix +two_step <- function(y, X, Z, wt) { + n <- nrow(Z) + b1 <- gmm_beta(y, X, Z, solve(crossprod(Z) / n)) + W <- solve(S_hat(Z, drop(y - X %*% b1), wt, ncol(X))) + list(b = gmm_beta(y, X, Z, W), W = W, b1 = b1) +} + +# Iterate until the largest change in the parameters is tiny. The weight +# matrix returned was evaluated at the previous iterate, so that is is +# equal to the weight matrix of the final estimates up to the tolerance +iterated <- function(y, X, Z, wt, tol = 1e-13, itermax = 20000) { + n <- nrow(Z) + b <- gmm_beta(y, X, Z, solve(crossprod(Z) / n)) + for (it in seq_len(itermax)) { + W <- solve(S_hat(Z, drop(y - X %*% b), wt, ncol(X))) + b_new <- gmm_beta(y, X, Z, W) + converged <- max(abs(b_new - b)) <= tol * max(1, abs(b)) + b <- b_new + if (converged) break + } + stopifnot(converged) + list(b = b, W = W, iterations = it) +} + +# Objective function of the continuously updating estimator +cue_objective <- function(b, y, X, Z, wt) { + n <- nrow(Z) + e <- drop(y - X %*% b) + g <- colMeans(Z * e) + n * drop(crossprod(g, solve(S_hat(Z, e, wt, ncol(X)), g))) +} + +# The objective is flat near its minimum, which optim's BFGS with numerical +# gradients does not locate precisely enough, so use Newton's method with +# central difference derivatives. The Hessian is made positive definite and +# the step is shortened until the objective decreases +cue <- function(y, X, Z, wt, start, tol = 1e-12) { + f <- function(b) cue_objective(b, y, X, Z, wt) + k <- length(start) + unit <- diag(k) + grad <- function(b, h) { + sapply(seq_len(k), function(j) (f(b + h[j] * unit[, j]) - f(b - h[j] * unit[, j])) / (2 * h[j])) + } + b <- start + for (it in seq_len(200)) { + h <- 1e-5 * pmax(abs(b), 1e-2) + g <- grad(b, h) + H <- sapply(seq_len(k), function(j) { + (grad(b + h[j] * unit[, j], h) - grad(b - h[j] * unit[, j], h)) / (2 * h[j]) + }) + H <- (H + t(H)) / 2 + lambda <- 0 + while (inherits(try(chol(H + lambda * unit), silent = TRUE), "try-error")) { + lambda <- max(2 * lambda, 1e-6 * max(abs(diag(H)))) + } + step <- solve(H + lambda * unit, g) + size <- 1 + while (f(b - size * step) > f(b) && size > 1e-10) size <- size / 2 + b_new <- b - size * step + done <- max(abs(b_new - b)) <= tol * max(1, abs(b)) + b <- b_new + if (done) break + } + stopifnot(done) + list(b = b, J = f(b)) +} + +# J of the model, using S from the weight matrix that was used for estimation +# (W), and J that uses the centered estimator of S evaluated at b. The +# centered statistic is NA if its S is singular, which uses the rank +# tolerance of numpy.linalg.matrix_rank +centered_J <- function(y, X, Z, b, wt) { + n <- nrow(Z) + e <- drop(y - X %*% b) + g <- colMeans(Z * e) + wc <- wt + wc$center <- TRUE + S <- S_hat(Z, e, wc, ncol(X)) + d <- svd(S, nu = 0, nv = 0)$d + if (sum(d > max(d) * max(dim(S)) * .Machine$double.eps) < ncol(Z)) { + return(NA_real_) + } + n * drop(crossprod(g, solve(S, g))) +} + +near <- function(a, b, rtol, what) { + if (abs(a - b) > rtol * abs(b)) { + stop(sprintf("check failed for %s: %.10g vs %.10g", what, a, b)) + } + cat(sprintf(" ok %-62s %.10g (expected %.10g)\n", what, a, b)) +} + +## -------------------------------------------------------------- data ----- +sim <- read.dta("simulated-data.dta") +sim$const <- 1 +# Six blocks of 100 observations, which is as many clusters as moments in the +# scenario with the clusters "block6" +sim$block6 <- (seq_len(nrow(sim)) - 1) %/% 100 +hs <- read.csv("housing.csv") +hs$const <- 1 +reg <- model.matrix(~ region, hs)[, -1] +hs <- cbind(hs, reg) +reg_cols <- colnames(reg) +mis <- read.csv("misspecified-data.csv") +mis$const <- 1 + +mat <- function(d, cols) as.matrix(d[, cols, drop = FALSE]) + +## ------------------------------------------------------- calibration ----- +cat("Calibration against Stata results already used in the tests\n") +# ivregress gmm rent pcturban (hsngval = faminc i.region), wmatrix(...) from +# stata-iv-housing-results.txt: gmm-{robust,unadjusted,cluster}-asymptotic +Xh <- mat(hs, c("hsngval", "pcturban", "const")) +Zh <- mat(hs, c("pcturban", "const", "faminc", reg_cols)) +fit <- two_step(hs$rent, Xh, Zh, wt_spec()) +near(fit$b[1], 0.00146432787, 1e-6, "housing GMM b[hsngval]") +near(fit$b[2], 0.76154815601, 1e-6, "housing GMM b[pcturban]") +near(fit$b[3], 112.12271295, 1e-6, "housing GMM b[_cons]") +near(gmm_J(hs$rent, Xh, Zh, fit$b, fit$W), 6.8364006463, 1e-6, + "housing Hansen J, robust, chi2(3)") +hj <- function(wt) { + f <- two_step(hs$rent, Xh, Zh, wt) + gmm_J(hs$rent, Xh, Zh, f$b, f$W) +} +near(hj(wt_spec("unadjusted")), 11.287665072, 1e-6, "housing Hansen J, unadjusted") +near(hj(wt_spec("clustered", clusters = as.integer(factor(hs$division)))), + 3.6677640329, 1e-6, "housing Hansen J, clustered by division") + +# ivregress gmm y x3 x4 x5 (x1 = z1 z2), wmatrix(...) on the simulated data. The +# dependent variable differs by weight matrix, as in the Stata run +Xs <- mat(sim, c("x1", "const", "x3", "x4", "x5")) +Zs <- mat(sim, c("const", "x3", "x4", "x5", "z1", "z2")) +sj <- function(dep, wt) { + f <- two_step(sim[[dep]], Xs, Zs, wt) + gmm_J(sim[[dep]], Xs, Zs, f$b, f$W) +} +cl <- sim$cluster_id +near(sj("y_unadjusted", wt_spec("unadjusted")), 0.38076944209, 1e-6, + "simulated Hansen J, unadjusted") +near(sj("y_robust", wt_spec("robust")), 0.22164820274, 1e-6, + "simulated Hansen J, robust") +near(sj("y_robust", wt_spec("robust", center = TRUE)), 0.22173011287, 1e-6, + "simulated Hansen J, robust, centered") +near(sj("y_clustered", wt_spec("clustered", clusters = cl)), 0.40598636292, 1e-6, + "simulated Hansen J, clustered") +near(sj("y_clustered", wt_spec("clustered", center = TRUE, clusters = cl)), + 0.40736456674, 1e-6, "simulated Hansen J, clustered, centered") +near(sj("y_kernel", wt_spec("kernel", bandwidth = 12)), 0.42317006975, 1e-6, + "simulated Hansen J, Bartlett kernel with 12 lags") +near(sj("y_kernel", wt_spec("kernel", center = TRUE, bandwidth = 12)), + 0.42592978028, 1e-6, "simulated Hansen J, Bartlett kernel, centered") + +## ----------------------------------------------------------- scenarios --- +# dataset, dependent, exog, endog and instruments are space separated name +# lists. weighted uses the column "weights" of the dataset. bandwidth is the +# number of lags of a kernel weight matrix, and clusters the name of the +# column that defines the clusters of a clustered weight matrix. iterated and +# cue select the estimators that are computed in addition to two_step. +scen <- function(id, dataset, dependent, exog, endog, instruments, + weight_type = "robust", center = FALSE, debiased = FALSE, + weighted = FALSE, bandwidth = "", clusters = "", + iterated = TRUE, cue = FALSE) { + data.frame(id = id, dataset = dataset, dependent = dependent, exog = exog, + endog = endog, instruments = instruments, weighted = weighted, + weight_type = weight_type, center = center, debiased = debiased, + bandwidth = bandwidth, clusters = clusters, iterated = iterated, + cue = cue, stringsAsFactors = FALSE) +} +hinstr <- paste(c("faminc", reg_cols), collapse = " ") +e2 <- list(exog = "const x3", endog = "x1 x2", instruments = "z1 z2 x4 x5") +e1 <- list(exog = "const x3 x5", endog = "x1", instruments = "z1 z2 x4") +sc <- rbind( + # Heteroskedasticity robust + scen("sim_e1_i3", "sim", "y_robust", e1$exog, e1$endog, e1$instruments), + scen("sim_e2_i4", "sim", "y_robust", e2$exog, e2$endog, e2$instruments, cue = TRUE), + scen("sim_e2_i4_weighted", "sim", "y_robust", e2$exog, e2$endog, e2$instruments, + weighted = TRUE), + scen("sim_e2_i4_center", "sim", "y_robust", e2$exog, e2$endog, e2$instruments, + center = TRUE, cue = TRUE), + scen("sim_e2_i4_debiased", "sim", "y_robust", e2$exog, e2$endog, e2$instruments, + debiased = TRUE), + scen("housing", "housing", "rent", "const pcturban", "hsngval", hinstr), + scen("mis_robust", "mis", "y", "const", "x", "z1 z2 zb", cue = TRUE), + scen("mis_robust_weighted", "mis", "y", "const", "x", "z1 z2 zb", weighted = TRUE), + scen("mis_robust_center", "mis", "y", "const", "x", "z1 z2 zb", center = TRUE, + cue = TRUE), + scen("mis_robust_debiased", "mis", "y", "const", "x", "z1 z2 zb", debiased = TRUE), + # Unadjusted + scen("sim_e2_i4_unadjusted", "sim", "y_unadjusted", e2$exog, e2$endog, + e2$instruments, weight_type = "unadjusted"), + scen("housing_unadjusted", "housing", "rent", "const pcturban", "hsngval", hinstr, + weight_type = "unadjusted"), + scen("mis_unadjusted", "mis", "y", "const", "x", "z1 z2 zb", + weight_type = "unadjusted"), + # Kernel + scen("sim_e2_i4_kernel6", "sim", "y_kernel", e2$exog, e2$endog, e2$instruments, + weight_type = "kernel", bandwidth = 6), + scen("sim_e1_i3_kernel3_center", "sim", "y_kernel", e1$exog, e1$endog, + e1$instruments, weight_type = "kernel", center = TRUE, bandwidth = 3), + scen("sim_e2_i4_kernel12_weighted", "sim", "y_kernel", e2$exog, e2$endog, + e2$instruments, weight_type = "kernel", weighted = TRUE, bandwidth = 12), + scen("mis_kernel4", "mis", "y", "const", "x", "z1 z2 zb", weight_type = "kernel", + bandwidth = 4), + # Clustered + scen("sim_e2_i4_clustered", "sim", "y_clustered", e2$exog, e2$endog, + e2$instruments, weight_type = "clustered", clusters = "cluster_id"), + scen("sim_e1_i3_clustered_center", "sim", "y_clustered", e1$exog, e1$endog, + e1$instruments, weight_type = "clustered", center = TRUE, + clusters = "cluster_id"), + scen("housing_clustered", "housing", "rent", "const pcturban", "hsngval", hinstr, + weight_type = "clustered", clusters = "division"), + scen("mis_clustered", "mis", "y", "const", "x", "z1 z2 zb", + weight_type = "clustered", clusters = "group", cue = TRUE), + scen("mis_clustered_debiased", "mis", "y", "const", "x", "z1 z2 zb", + weight_type = "clustered", debiased = TRUE, clusters = "group"), + # As many clusters as moment conditions, so that the centered covariance has + # rank 5 and not 6 and the centered statistic is not defined + scen("sim_e1_i3_six_clusters", "sim", "y_clustered", e1$exog, e1$endog, + e1$instruments, weight_type = "clustered", clusters = "block6", + iterated = FALSE) +) + +words <- function(s) strsplit(s, " ", fixed = TRUE)[[1]] + +run_scenario <- function(r) { + d <- switch(r$dataset, sim = sim, housing = hs, mis = mis) + sw <- if (r$weighted) sqrt(d$weights / mean(d$weights)) else rep(1, nrow(d)) + y <- d[[r$dependent]] * sw + X <- mat(d, c(words(r$exog), words(r$endog))) * sw + Z <- mat(d, c(words(r$exog), words(r$instruments))) * sw + n <- nrow(Z) + wt <- wt_spec(r$weight_type, r$center, r$debiased, + if (r$bandwidth == "") 0 else as.integer(r$bandwidth), + if (r$clusters == "") NULL else as.integer(factor(d[[r$clusters]]))) + df <- ncol(Z) - ncol(X) + # b is the estimate, and wb the parameters that the weight matrix of J was + # evaluated at, which are those of the previous step of a two step estimator + fmt <- function(b) paste(sprintf("%.17g", b), collapse = " ") + row <- function(method, b, j, wb = b) { + data.frame(r[, c("id", "dataset", "dependent", "exog", "endog", "instruments", + "weighted", "weight_type", "center", "debiased", "bandwidth", + "clusters")], + method = method, nobs = n, df = df, j = j, + centered_j = centered_J(y, X, Z, b, wt), params = fmt(b), + weight_params = fmt(wb), stringsAsFactors = FALSE) + } + out <- list() + + ts <- two_step(y, X, Z, wt) + out$two_step <- row("two_step", ts$b, gmm_J(y, X, Z, ts$b, ts$W), wb = ts$b1) + + if (r$iterated) { + it <- iterated(y, X, Z, wt) + j <- gmm_J(y, X, Z, it$b, it$W) + res <- row("iterated", it$b, j) + out$iterated <- res + jc <- res$centered_j + if (r$weight_type == "robust") { + # Hansen and Lee (2021) Theorem 1: the iterated estimator does not + # change if S is centered + other <- wt + other$center <- !wt$center + b_other <- iterated(y, X, Z, other)$b + stopifnot(max(abs(b_other - it$b)) <= 1e-8 * max(1, abs(it$b))) + if (!r$center) { + # S = s (S_c + gbar gbar') with the small-sample scale s of the + # estimator, and Sherman-Morrison, which implies J < n / s + s <- if (r$debiased) n / (n - ncol(X)) else 1 + stopifnot(j < n / s) + near(jc, j / (1 - s * j / n), 1e-8, + paste(r$id, "Jc = J / (1 - s J / n)")) + } + } + if (r$weight_type == "unadjusted") { + near(jc, j, 1e-8, paste(r$id, "Jc = J for unadjusted")) + } + } + + if (r$cue) { + fit <- cue(y, X, Z, wt, ts$b) + # Another start finds the same minimum + again <- cue(y, X, Z, wt, ts$b * (1 + 0.05 * c(1, -1, 1, -1)[seq_along(ts$b)])) + near(again$J, fit$J, 1e-7, paste(r$id, "CUE minimum is the same from other start")) + stopifnot(fit$J <= gmm_J(y, X, Z, ts$b, ts$W) + 1e-8) + out$cue <- row("cue", fit$b, fit$J) + } + do.call(rbind, out) +} + +cat("Scenarios and checks of the properties of the statistic\n") +res <- do.call(rbind, lapply(seq_len(nrow(sc)), function(i) run_scenario(sc[i, ]))) +rownames(res) <- NULL +for (i in seq_len(nrow(res))) { + cat(sprintf("%-30s %-9s J=%-12.8g Jc=%-12.8g df=%d\n", res$id[i], res$method[i], + res$j[i], res$centered_j[i], res$df[i])) +} +# The scenarios are useful only if the statistics differ in some of them +stopifnot(any(res$centered_j > 1.1 * res$j, na.rm = TRUE)) +stopifnot(any(is.na(res$centered_j))) +write.csv(res, "centered-j-reference.csv", row.names = FALSE, quote = TRUE) +cat("wrote centered-j-reference.csv\n") diff --git a/linearmodels/tests/iv/results/centered-j-reference.csv b/linearmodels/tests/iv/results/centered-j-reference.csv new file mode 100644 index 0000000000..b8e6841e2a --- /dev/null +++ b/linearmodels/tests/iv/results/centered-j-reference.csv @@ -0,0 +1,51 @@ +"id","dataset","dependent","exog","endog","instruments","weighted","weight_type","center","debiased","bandwidth","clusters","method","nobs","df","j","centered_j","params","weight_params" +"sim_e1_i3","sim","y_robust","const x3 x5","x1","z1 z2 x4",FALSE,"robust",FALSE,FALSE,"","","two_step",600,2,34.447915388289,34.6464711949478,"1.0387363975923256 0.19783480893277858 0.72729673669628914 2.7408293475914358","1.021816507727221 0.19914956452741089 0.75655796213391535 2.6458175928313792" +"sim_e1_i3","sim","y_robust","const x3 x5","x1","z1 z2 x4",FALSE,"robust",FALSE,FALSE,"","","iterated",600,2,32.7908483734886,34.6865154902307,"1.0455094471113475 0.20262777324545789 0.72701332478259306 2.7366609258144501","1.0455094471113475 0.20262777324545789 0.72701332478259306 2.7366609258144501" +"sim_e2_i4","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",FALSE,FALSE,"","","two_step",600,2,76.8900475843116,72.7530359165364,"0.9408575627669552 -0.020539359597322622 2.446918401298948 1.3079653233086512","0.87760165534261225 0.048978217611413744 1.7717561888575446 1.8632094501989847" +"sim_e2_i4","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",FALSE,FALSE,"","","iterated",600,2,67.3949199186077,75.9229557949153,"0.95142714972837494 -0.0043767829885885028 2.333398448788143 1.4205869914513498","0.95142714972837494 -0.0043767829885885028 2.333398448788143 1.4205869914513498" +"sim_e2_i4","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",FALSE,FALSE,"","","cue",600,2,1.28649068948938,1.28925504718034,"31.72372594879841 2.3355679604859185 555.12059706633113 -533.20016370919905","31.72372594879841 2.3355679604859185 555.12059706633113 -533.20016370919905" +"sim_e2_i4_weighted","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",TRUE,"robust",FALSE,FALSE,"","","two_step",600,2,31.267260384654,46.8042960576167,"0.87335629751578259 -0.028979848924509905 0.87921377054415206 2.7957275405373188","0.83678744559873286 0.010070709871668278 0.25809104395926868 3.3317816178426503" +"sim_e2_i4_weighted","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",TRUE,"robust",FALSE,FALSE,"","","iterated",600,2,46.1521354805488,49.9979923410109,"0.90691273923123794 -0.064548062387158767 1.086967975674006 2.6512700656698764","0.90691273923123794 -0.064548062387158767 1.086967975674006 2.6512700656698764" +"sim_e2_i4_center","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",TRUE,FALSE,"","","two_step",600,2,88.1918387091334,69.4795148047788,"0.95015532074190212 -0.030757497858393448 2.5461580628691833 1.2263519891704622","0.87760165534261225 0.048978217611413744 1.7717561888575446 1.8632094501989847" +"sim_e2_i4_center","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",TRUE,FALSE,"","","iterated",600,2,75.9229557949236,75.922955794921,"0.95142714972836206 -0.004376782988609316 2.3333984487879187 1.4205869914515856","0.95142714972836206 -0.004376782988609316 2.3333984487879187 1.4205869914515856" +"sim_e2_i4_center","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",TRUE,FALSE,"","","cue",600,2,1.28925504718064,1.28925504718064,"31.724011799382083 2.3355806077397281 555.12576773913645 -533.20515269592761","31.724011799382083 2.3355806077397281 555.12576773913645 -533.20515269592761" +"sim_e2_i4_debiased","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",FALSE,TRUE,"","","two_step",600,2,76.3774472670827,72.2680156770978,"0.94085756276694621 -0.020539359597323049 2.4469184012987815 1.3079653233088127","0.87760165534261225 0.048978217611413744 1.7717561888575446 1.8632094501989847" +"sim_e2_i4_debiased","sim","y_robust","const x3","x1 x2","z1 z2 x4 x5",FALSE,"robust",FALSE,TRUE,"","","iterated",600,2,66.9456204524933,75.4168027562837,"0.95142714972837017 -0.0043767829886057694 2.3333984487880852 1.4205869914514218","0.95142714972837017 -0.0043767829886057694 2.3333984487880852 1.4205869914514218" +"housing","housing","rent","const pcturban","hsngval","faminc regionNE regionSouth regionWest",FALSE,"robust",FALSE,FALSE,"","","two_step",50,3,6.83640064627001,8.10619997080815,"112.12271295048038 0.76154815600611991 0.001464327867106802","120.70651453650972 0.081515968185227675 0.0022398329844168334" +"housing","housing","rent","const pcturban","hsngval","faminc regionNE regionSouth regionWest",FALSE,"robust",FALSE,FALSE,"","","iterated",50,3,3.85790639827124,4.18046310552144,"113.77151906725848 1.133841782934728 0.00084401328302554073","113.77151906725848 1.133841782934728 0.00084401328302554073" +"mis_robust","mis","y","const","x","z1 z2 zb",FALSE,"robust",FALSE,FALSE,"","","two_step",300,2,80.750068052939,109.038231451155,"0.80419228898649464 1.543770895293592","0.74780817604971728 1.6534283755456953" +"mis_robust","mis","y","const","x","z1 z2 zb",FALSE,"robust",FALSE,FALSE,"","","iterated",300,2,78.6336127431786,106.565789482687,"0.87780510620591945 1.4412073515961232","0.87780510620591945 1.4412073515961232" +"mis_robust","mis","y","const","x","z1 z2 zb",FALSE,"robust",FALSE,FALSE,"","","cue",300,2,74.3965010577646,98.9299830099001,"0.81526183679751796 0.50272373802602077","0.81526183679751796 0.50272373802602077" +"mis_robust_weighted","mis","y","const","x","z1 z2 zb",TRUE,"robust",FALSE,FALSE,"","","two_step",300,2,69.9274730795802,90.212761493308,"0.68049955551076058 1.4414358920434238","0.75822457446680103 1.5761259557894858" +"mis_robust_weighted","mis","y","const","x","z1 z2 zb",TRUE,"robust",FALSE,FALSE,"","","iterated",300,2,67.7949647633942,87.5884944023458,"0.60100578081470113 1.2890041135384014","0.60100578081470113 1.2890041135384014" +"mis_robust_center","mis","y","const","x","z1 z2 zb",FALSE,"robust",TRUE,FALSE,"","","two_step",300,2,110.490435279729,108.090710931535,"0.82495863959064697 1.5033838862111475","0.74780817604971728 1.6534283755456953" +"mis_robust_center","mis","y","const","x","z1 z2 zb",FALSE,"robust",TRUE,FALSE,"","","iterated",300,2,106.565789482688,106.565789482687,"0.8778051062059441 1.4412073515961228","0.8778051062059441 1.4412073515961228" +"mis_robust_center","mis","y","const","x","z1 z2 zb",FALSE,"robust",TRUE,FALSE,"","","cue",300,2,98.9299830099001,98.9299830099001,"0.81526178847968311 0.50272367657004624","0.81526178847968311 0.50272367657004624" +"mis_robust_debiased","mis","y","const","x","z1 z2 zb",FALSE,"robust",FALSE,TRUE,"","","two_step",300,2,80.2117342659194,108.311309908147,"0.80419228898649453 1.5437708952935916","0.74780817604971728 1.6534283755456953" +"mis_robust_debiased","mis","y","const","x","z1 z2 zb",FALSE,"robust",FALSE,TRUE,"","","iterated",300,2,78.1093886582241,105.855350886136,"0.87780510620591923 1.441207351596123","0.87780510620591923 1.441207351596123" +"sim_e2_i4_unadjusted","sim","y_unadjusted","const x3","x1 x2","z1 z2 x4 x5",FALSE,"unadjusted",FALSE,FALSE,"","","two_step",600,2,223.330987412239,223.330987412239,"0.91426856040012738 0.11478559796790451 1.5186069645655693 1.7709213590354511","0.9142685604001286 0.11478559796790429 1.5186069645655877 1.7709213590354336" +"sim_e2_i4_unadjusted","sim","y_unadjusted","const x3","x1 x2","z1 z2 x4 x5",FALSE,"unadjusted",FALSE,FALSE,"","","iterated",600,2,223.330987412239,223.330987412239,"0.91426856040012738 0.11478559796790451 1.5186069645655693 1.7709213590354511","0.91426856040012738 0.11478559796790451 1.5186069645655693 1.7709213590354511" +"housing_unadjusted","housing","rent","const pcturban","hsngval","faminc regionNE regionSouth regionWest",FALSE,"unadjusted",FALSE,FALSE,"","","two_step",50,3,11.2876650716273,11.2876650716275,"120.70651453650946 0.081515968185257623 0.0022398329844167974","120.70651453650972 0.081515968185227675 0.0022398329844168334" +"housing_unadjusted","housing","rent","const pcturban","hsngval","faminc regionNE regionSouth regionWest",FALSE,"unadjusted",FALSE,FALSE,"","","iterated",50,3,11.2876650716273,11.2876650716275,"120.70651453650946 0.081515968185257623 0.0022398329844167974","120.70651453650946 0.081515968185257623 0.0022398329844167974" +"mis_unadjusted","mis","y","const","x","z1 z2 zb",FALSE,"unadjusted",FALSE,FALSE,"","","two_step",300,2,211.372404390446,211.372404390446,"0.74780817604971739 1.6534283755456956","0.74780817604971728 1.6534283755456953" +"mis_unadjusted","mis","y","const","x","z1 z2 zb",FALSE,"unadjusted",FALSE,FALSE,"","","iterated",300,2,211.372404390446,211.372404390446,"0.74780817604971739 1.6534283755456956","0.74780817604971739 1.6534283755456956" +"sim_e2_i4_kernel6","sim","y_kernel","const x3","x1 x2","z1 z2 x4 x5",FALSE,"kernel",FALSE,FALSE,"6","","two_step",600,2,41.3984143198759,96.1216891975662,"0.9806921937056059 0.11741820312078874 1.4689711544327031 1.7402995197293119","0.93145558950296836 0.160755674825863 0.74007895622013231 2.4338825837448415" +"sim_e2_i4_kernel6","sim","y_kernel","const x3","x1 x2","z1 z2 x4 x5",FALSE,"kernel",FALSE,FALSE,"6","","iterated",600,2,39.4114686426462,72.4907094432892,"1.0498953805989797 0.20016077086503781 2.3124330691735411 0.82479472509501139","1.0498953805989797 0.20016077086503781 2.3124330691735411 0.82479472509501139" +"sim_e1_i3_kernel3_center","sim","y_kernel","const x3 x5","x1","z1 z2 x4",FALSE,"kernel",TRUE,FALSE,"3","","two_step",600,2,41.5508448637388,41.0743822965028,"1.1414903877247398 0.31488481349462794 0.74149574484729364 2.1847103385480025","1.1094335279001255 0.32193463168905173 0.78053420201140611 2.1649598362582192" +"sim_e1_i3_kernel3_center","sim","y_kernel","const x3 x5","x1","z1 z2 x4",FALSE,"kernel",TRUE,FALSE,"3","","iterated",600,2,41.0690669745437,41.0690669745453,"1.1428515389669041 0.3125079784234468 0.73991748597279927 2.1843709473650539","1.1428515389669041 0.3125079784234468 0.73991748597279927 2.1843709473650539" +"sim_e2_i4_kernel12_weighted","sim","y_kernel","const x3","x1 x2","z1 z2 x4 x5",TRUE,"kernel",FALSE,FALSE,"12","","two_step",600,2,18.59404242854,34.8188249206781,"0.94676116696388724 0.053168988231793016 0.0044164983153378327 3.2638401885633641","0.92232960734112135 0.10463072159863709 -0.21836367049057767 3.4281497164123005" +"sim_e2_i4_kernel12_weighted","sim","y_kernel","const x3","x1 x2","z1 z2 x4 x5",TRUE,"kernel",FALSE,FALSE,"12","","iterated",600,2,22.1131918074663,42.191744571811,"0.94444881453433283 -0.0036517315982342504 0.34890844484208844 2.9678444126460026","0.94444881453433283 -0.0036517315982342504 0.34890844484208844 2.9678444126460026" +"mis_kernel4","mis","y","const","x","z1 z2 zb",FALSE,"kernel",FALSE,FALSE,"4","","two_step",300,2,39.770174244194,118.543999270778,"0.7392863881812014 1.6016660444778315","0.74780817604971728 1.6534283755456953" +"mis_kernel4","mis","y","const","x","z1 z2 zb",FALSE,"kernel",FALSE,FALSE,"4","","iterated",300,2,39.3835763550875,115.250838336216,"0.65548306871914863 1.4110162295991544","0.65548306871914863 1.4110162295991544" +"sim_e2_i4_clustered","sim","y_clustered","const x3","x1 x2","z1 z2 x4 x5",FALSE,"clustered",FALSE,FALSE,"","cluster_id","two_step",600,2,62.2483755732058,94.1394589720762,"0.99806333351181897 0.13270714998462066 2.4578648920759503 0.81581702220814822","0.9163634823363811 0.11201216540664237 1.5199951179972002 1.7819111916644939" +"sim_e2_i4_clustered","sim","y_clustered","const x3","x1 x2","z1 z2 x4 x5",FALSE,"clustered",FALSE,FALSE,"","cluster_id","iterated",600,2,51.3114922869231,89.642056283282,"1.0066838407525809 0.16884019545567766 2.5268541279760379 0.69441663333333525","1.0066838407525809 0.16884019545567766 2.5268541279760379 0.69441663333333525" +"sim_e1_i3_clustered_center","sim","y_clustered","const x3 x5","x1","z1 z2 x4",FALSE,"clustered",TRUE,FALSE,"","cluster_id","two_step",600,2,76.8502495855277,71.4977190631225,"1.0530214773280551 0.26530909072874681 0.64426664330335692 2.4043632312450396","1.0562112223258302 0.26210387413359709 0.76198062545416534 2.303542128137714" +"sim_e1_i3_clustered_center","sim","y_clustered","const x3 x5","x1","z1 z2 x4",FALSE,"clustered",TRUE,FALSE,"","cluster_id","iterated",600,2,71.3055086921529,71.3055086921562,"1.0537108122009007 0.26578783651094978 0.63628665259002115 2.4073749718913571","1.0537108122009007 0.26578783651094978 0.63628665259002115 2.4073749718913571" +"housing_clustered","housing","rent","const pcturban","hsngval","faminc regionNE regionSouth regionWest",FALSE,"clustered",FALSE,FALSE,"","division","two_step",50,3,3.66776403286596,74.2870360851392,"100.87359600312622 0.63495496019525044 0.00192718161301284","120.70651453650972 0.081515968185227675 0.0022398329844168334" +"housing_clustered","housing","rent","const pcturban","hsngval","faminc regionNE regionSouth regionWest",FALSE,"clustered",FALSE,FALSE,"","division","iterated",50,3,3.29010053586092,4.60884847815156,"121.416047124421 1.6823708566717663 -0.00017292415191441058","121.416047124421 1.6823708566717663 -0.00017292415191441058" +"mis_clustered","mis","y","const","x","z1 z2 zb",FALSE,"clustered",FALSE,FALSE,"","group","two_step",300,2,40.5643574093102,124.996028398726,"0.78022708697624199 1.6341622302702141","0.74780817604971728 1.6534283755456953" +"mis_clustered","mis","y","const","x","z1 z2 zb",FALSE,"clustered",FALSE,FALSE,"","group","iterated",300,2,40.3240274709715,122.964272525224,"0.89531875044019082 1.5822768530447593","0.89531875044019082 1.5822768530447593" +"mis_clustered","mis","y","const","x","z1 z2 zb",FALSE,"clustered",FALSE,FALSE,"","group","cue",300,2,38.7233743057053,109.199761829026,"-0.7619316201980153 -0.7897820714187106","-0.7619316201980153 -0.7897820714187106" +"mis_clustered_debiased","mis","y","const","x","z1 z2 zb",FALSE,"clustered",FALSE,TRUE,"","group","two_step",300,2,39.7548791510865,122.501681789654,"0.78022708697624188 1.6341622302702143","0.74780817604971728 1.6534283755456953" +"mis_clustered_debiased","mis","y","const","x","z1 z2 zb",FALSE,"clustered",FALSE,TRUE,"","group","iterated",300,2,39.5193450944605,120.510470431354,"0.89531875044019138 1.5822768530447595","0.89531875044019138 1.5822768530447595" +"sim_e1_i3_six_clusters","sim","y_clustered","const x3 x5","x1","z1 z2 x4",FALSE,"clustered",FALSE,FALSE,"","block6","two_step",600,2,5.71934733454157,NA,"1.0325290142041232 0.26003681823006841 0.72986933824949174 2.3646741217275","1.0562112223258302 0.26210387413359709 0.76198062545416534 2.303542128137714" diff --git a/linearmodels/tests/iv/results/misspecified-data.csv b/linearmodels/tests/iv/results/misspecified-data.csv new file mode 100644 index 0000000000..22efd0e8fc --- /dev/null +++ b/linearmodels/tests/iv/results/misspecified-data.csv @@ -0,0 +1,301 @@ +y,x,z1,z2,zb,group,weights +2.958269403663809,1.915710415483861,0.7303700527552083,0.9465944302291408,-0.5418142555300646,0,0.7437036062309645 +0.8329547650674989,-0.10852328420915325,-0.21381523798065996,-0.516388702287081,-0.6045941477566886,0,0.5935889614343164 +1.9616616399723164,-0.4940112818018646,-1.5126036083820427,0.250701629659096,-0.06965958416659856,0,1.2593447206523778 +2.913993190600483,-0.08054418158805787,-0.1732102597741807,-1.284425339598826,0.46034786976896325,0,1.7674705329008589 +0.7367793493619592,1.7880750393593732,1.3043997025709384,-0.13850895919233727,-1.3379672766634154,0,1.4092977559921365 +5.239544514178503,0.5607702438669324,1.3077163143119432,1.0469153516308232,1.4668146690107366,1,1.7008814166843462 +-2.4285438708052314,-3.741403635122796,-2.0119317872413487,-0.6237351330081786,0.877252377982353,1,1.6136206949560232 +4.959145205773169,-1.243064377092407,-1.1166313035532758,-0.8846953680095427,1.5867133024579132,1,0.3938012581003548 +4.955808843553651,1.051997524729542,1.5354413871651544,0.5120736671972522,-1.0719991346777977,1,1.538620691836892 +0.6937691764664002,0.5359029732896249,0.30483441760118724,0.4480195052525989,-1.4747397447608415,1,0.9838427565543996 +1.1023623770462374,1.167054439170946,0.0716586127680766,0.5373648879764925,-0.6200519416313377,2,1.764963111969013 +-7.659732837595938,-2.0312718363509945,-0.5370594248937794,0.2007034711156063,-2.4706426700488757,2,0.6746205146456019 +2.6725054622816344,-1.9049727787293989,-0.7416388918945219,-1.5740144536456735,1.3222504130063062,2,1.5353026686613247 +-2.586692789779767,-0.7248976911391105,2.003525565912089,-0.9098756793925423,-0.5728181235222058,2,0.28218675141947813 +5.705155200695566,1.6118038930810732,1.1828478926834436,-2.1733541281405504,0.9403741289009797,2,0.5951186554987755 +-1.1635892267615724,-0.7169302684479273,-1.7561734737161978,0.6491478357360618,-0.20193860604557629,3,1.7798618927716443 +-3.976245510045868,-0.30985725549483156,-0.6527831555636565,0.038984833174710694,-1.5216290087771944,3,1.6246845672261294 +-1.2100297095792296,-0.7862646779675071,0.31648398743622336,0.7126815271314701,-0.1689157616600646,3,0.38914000018319717 +-0.05642278845435428,0.6329295545443289,0.1742137791075811,0.33197420008089346,-1.2574491404468284,3,1.7165017300018686 +1.756026939855871,0.3375443671406677,1.1849178834779082,-0.41331111002887166,1.2500920403429188,3,0.6707155664480049 +2.683035891110621,2.105245296690994,0.0848817431495182,1.2579523719638201,0.05593581299426142,4,1.8817574113506719 +-1.4808665602512066,-1.061358886944535,0.8149433496438384,-0.4945268234975099,-0.19483068026711564,4,1.9771141283133224 +1.4980208519255054,-0.5272235202226212,-0.005866515227590213,-1.1050231439891915,-0.13583436730303952,4,1.6480921195004727 +7.717316737363743,-0.06747841899700702,-1.5816742214515522,1.387283402401902,3.3572660017180516,4,1.9892754793616658 +1.33825308240635,-0.48035144516491507,-1.004295315432542,0.22997880166584145,-0.3031174180585578,4,0.8498443138038592 +5.53648452825758,1.5533498174307816,0.47219152414978993,-0.2460139393388056,0.7045998471315106,5,1.6950540620442294 +-0.48297777628399946,-1.021851359646012,-0.7088896260968297,0.1345497911283033,-0.6360809868624419,5,0.9618333595234898 +10.73248090675222,3.2770655916076974,0.9560576653693094,0.8330588836329301,2.750169828075492,5,1.7382006179010143 +-3.2878970581663975,-0.6949885285765827,-0.40330583408782217,1.5138175302864794,-1.3686741220384573,5,1.0509409715587321 +-2.1808189674615437,0.13587790314066636,0.38043512024827153,0.8471863905883188,-1.0994883599122016,5,0.6878755428021809 +-0.02199649250204344,0.8823839687976198,0.3215169895554431,0.20718511068405546,-0.1693528266990382,6,1.4864456659897987 +-3.0215136611274978,-2.14656987004925,-0.5914929661201958,-0.9561268161020038,-0.6157935455225435,6,0.8524013505876884 +-1.4905678283425792,-1.216114719196303,-0.08404079929191595,-1.0158082426055868,-0.21760304190360832,6,1.0586328902006714 +-6.6186798221202015,-3.9757030961753888,-1.9480721198925557,-1.976740133051463,-1.5323183119501242,6,1.0080548815553616 +3.3214505549870643,0.6902908488639572,0.24333945803320584,-0.7259884607908632,0.953248352287193,6,1.5318082754726023 +2.1541604770069895,-0.1591035202567608,-1.1399497020374056,0.9778329882212438,0.20965530842587354,7,0.7059333011830795 +3.0753930146066475,-0.48971476736232944,-0.3299050716835647,1.0379571290736098,0.7332764132971796,7,0.350251881162834 +2.9418358801189903,0.7769107705932368,-0.5265320423242844,1.2294979556970282,0.0988549137898959,7,1.681956682515771 +9.134317141800878,0.7723168808107467,-0.8398428936338217,0.24186956820386282,1.9774833327878385,7,1.0663462000322175 +-7.739845773267658,-4.170534948135549,-1.7229073032352038,-0.6959361266428196,-0.4879188019060541,7,1.1734276483271096 +2.187788467532581,-0.532931861319037,-0.27396447370556976,-0.962552523524673,0.5021304599301235,8,1.2368964382855943 +-1.4499814248306877,-2.0499994499674337,0.02134121222320837,-1.4940344406028248,-1.0585570219353782,8,1.5672475633955543 +-1.123754717532035,-1.2726469672144445,-0.6183329599934674,-0.636556352896624,-0.43252652566383953,8,1.6122801519124685 +-3.206862713316379,-2.723376501891806,-0.9761346372981895,-0.8501941582492976,-0.7622385486674632,8,0.4567809208766068 +1.7685998238990601,0.3073913883313429,0.07484155412904495,0.6769061536096508,-0.30899310689199455,8,0.5431156336920409 +-3.841122562676031,-1.819916986239762,0.9549759094267388,0.44216967889788067,0.020288762796722152,9,1.433053932008368 +-1.2147927194215207,0.7604071447567975,-0.9219743227017724,0.08228873018718033,-1.5720144774512335,9,1.3840494295984498 +1.1879603272809653,0.14382545087185655,1.0788184241554817,-0.4952071353931553,-0.32604297244801356,9,0.5218509470390867 +1.1949373820426152,1.902343662520812,1.0048278013804426,0.682194273562486,0.08874122692860963,9,1.4942147260718288 +2.583040111892981,-0.05660745979360281,0.27350550643716637,-1.0779105205442923,0.6482693850101949,9,0.8477185899700719 +-0.8639054204579744,0.010779063727825361,0.2863231520281442,-0.5994808684255928,-0.3516373718573524,10,0.7024262032331474 +2.83908887098378,-0.3183861799117526,0.11834262998395179,-0.5650257546365541,0.7792733646645184,10,1.9770979197974023 +5.02237426881751,0.9303883972261859,0.39528432105309713,-1.2360719222929948,0.7967106011817288,10,1.1725945675249636 +-5.709185111917011,-0.6597928632875054,0.3792309508133216,-0.8367503080253135,-2.5000913624289116,10,1.0409096430700928 +-7.255066041521814,-0.5273646714342332,-0.1515545930857107,0.696563416420119,-2.819202202174404,10,0.91617398988349 +-0.8759032830750588,-0.8718754363661335,-0.9244558503784825,-0.18194343257987478,-0.028694224780999984,11,1.6461307618381988 +3.8629111903480466,-0.10906964791206337,-0.4271452597930619,-0.27740886364487927,1.1894760079260656,11,0.827687026916424 +3.156185192765769,2.5239187562970784,1.0056704456584105,1.0263363323021022,-0.4539774754754613,11,0.9946172928214545 +-0.7474792081447086,-1.4025146589001523,0.008716743682122545,-1.2775585646341714,-0.27357163793555495,11,0.97128518409894 +2.676804963471636,1.3727444532934205,-0.3505146920201209,1.689204306088344,-0.45126720460485537,11,1.7753957198158805 +2.1682677897872007,0.9224368949200572,1.2710771168619477,-0.5573536641302073,-0.8446960112791175,12,0.7463829180248113 +-1.689740600936338,-1.4154352610007512,-1.495859565838478,0.30577814714910717,-0.8565515886147606,12,1.6810971327900186 +1.8897395701979889,-0.7815514965182655,0.30532602350490995,-0.4806756689898086,0.7165893211134127,12,1.4083979855829998 +0.14584908270646246,0.2615633756188306,-0.46270851924038126,1.8106978733799526,-0.46247015522978485,12,0.5944847964884951 +-1.4042738182832393,-0.11684134661774605,1.2622581505307584,-1.0662898095898612,-1.2824699572188512,12,0.487513441404545 +-5.244465195066562,-1.2129964464149066,1.330095582143609,-0.4949953447154995,-1.5827750894654862,13,1.942262481645788 +-2.2572150290956814,0.4922761620252976,0.6428671170277993,-0.4209845564040565,-0.3106606120888969,13,0.8407516706148545 +-0.09239812725498853,1.1649413517643268,0.2961586289520059,1.5880861167131346,-1.0646861061282682,13,1.048160343912389 +2.9828763019201925,1.0389886288419419,-0.11333166454283633,0.6958405719937386,-0.023440520069429514,13,0.981178711109005 +3.657696655586323,-0.25418718012685615,-0.5010672611843601,-0.47957647202646503,1.1184020284819,13,0.3580739612066969 +-6.052974149573677,-1.8800425694026348,-1.5878012793647969,-1.3610997309816681,-0.6135701845640176,14,1.5925189288728094 +0.9676535653726783,-1.7537987789317173,-1.6672607263015693,-0.5227885584181439,-0.3005743537754992,14,1.3798986503925916 +2.309310240823281,-0.4748370479463045,-0.47168795768132643,0.06343257154615001,0.23231131868181207,14,0.9516133692768602 +1.758989983941735,0.530641524876688,0.7055413065845478,0.9636198477380411,-0.44631687380469837,14,1.1576861396416787 +6.022772590616309,1.6899536642515656,1.4643125946813216,-0.8008080322121723,0.6329695360218472,14,1.3051146898073007 +1.9489072800597869,-1.799911949597624,-1.9853575045657152,0.36850166022893677,0.9570215817658836,15,1.8239625154052592 +1.3139166471404167,0.31037633051456837,0.44156345874939956,0.5495304372025912,-0.26118922959721136,15,1.1118088397047363 +-0.11494341139248743,0.15432944763951867,0.9011706305068803,-0.8237318783985146,-1.5248893018448968,15,1.9406311224774 +-1.906205204923348,-1.7240529634638302,-1.7688672701651957,0.2556628775094522,-0.34305942259562244,15,1.119771292560712 +3.661217747135879,0.9009145729867991,0.5256422755371501,-0.2930897437797102,-0.34160077154416796,15,1.9679090789149645 +1.9251994298097133,-1.2094535103382391,-1.352735454779728,0.14160542287193054,1.4254542138197857,16,0.28913585110561174 +0.7947983269462622,0.1445255106499168,-0.3312380549736514,-0.33088661277934184,0.7956350031812963,16,1.7949733535189614 +0.4953380831991976,1.357228666206045,-0.77990126116703,1.0488944256378545,-0.4395738166572779,16,1.5923661992322884 +0.41255254606866765,1.3425253022436958,0.3758439975960465,-0.6069348116022925,-0.35423400733857546,16,0.9138917257371248 +-3.8820456054901236,1.5611128965984793,1.195163267996701,0.6939163884025127,-2.04953467129589,16,0.7902376907117855 +4.525277669673589,1.2009964414608818,-0.05271684643720397,1.6252802241690487,-0.23981233332207563,17,1.1510875993659266 +1.287402539347271,1.1421751131035762,-0.6205753660675222,0.456409408012839,-0.6047010331305173,17,1.3574950534745927 +-0.7909720987421349,-0.12675999854389963,0.2614475741462509,1.128673729781687,-0.5009589125948162,17,1.247851350242334 +-2.6237007305134106,-0.12513704146310567,0.39714804146598065,1.2900698528345285,-0.5117330146385948,17,1.6362885243363252 +0.6226063711267851,-1.04371455485327,-0.5322356250673066,-0.09631039443601601,-0.3497494331958686,17,0.9868012580730767 +4.187731686239691,-0.06643728624499073,-0.6973748166366792,-0.8121955904201338,1.0545087480939057,18,1.460576555515775 +1.2466267887192095,-0.5372306445604031,-1.125903364855824,-0.6332231681549062,0.555875446957435,18,1.0872254890712318 +1.3369975249234094,1.1841164756404934,1.0389355598812964,1.1258434542731033,0.7941464780897706,18,1.7945395764108034 +1.069899535227543,-1.4722331392838754,0.3951341328034816,-1.034078106667284,0.00533672311560644,18,0.9521764386587247 +1.8073664890606782,0.18123830332133883,-0.3175603483187787,-0.7992566344922983,-0.33418366063779,18,1.4677424307159956 +2.2213717419692944,1.1867422479329708,0.9993499480010414,0.2992377625117017,-0.4797284354491165,19,0.5632680218826159 +0.6916760251566612,-0.4506030810090258,-1.6856843896466558,0.03179453542448021,0.24317412967889346,19,1.6551613306380388 +1.6618660908537064,-0.7049298230977961,-0.7505465379519102,1.261883564341513,0.6526349339696002,19,0.5640554337314281 +11.625148898371076,2.4101512786529167,-1.042432873746889,0.749270650937372,1.612038258512965,19,0.3343640057709053 +6.5062331962611415,1.8471486445633323,0.18486772169001758,0.048129142657954045,0.9583687333203612,19,1.6153500758028403 +6.884861493784417,0.28103058871064956,0.18244535826153857,0.1299130184483145,1.7348029464757158,20,0.9098731528202394 +-6.385101853243487,-1.7871764485475605,-0.3416284336920741,0.584215918386603,-2.242242577505365,20,0.6527884598069624 +-3.0678900106232465,-0.8811076108756211,-0.4389705682377842,0.6323810738437057,-0.8394125847203557,20,0.35263762749935634 +-3.0258573435465177,-1.817865275337763,0.4403396205105493,-0.20383269406818813,-0.6260620848645042,20,1.1329993543735852 +1.690629741616742,-0.6636908581391319,-1.218174751295915,0.304274606404671,-0.563360395147635,20,1.5174159915168839 +-2.6478942461214365,-2.5809675297690653,-1.7783180385990425,-0.3517011613372201,0.37162736053903866,21,0.924966123600714 +2.55445536565185,-0.612281991205873,-2.2000123411064822,-1.3904303087597665,0.06375094369859083,21,0.48364815280772067 +0.5564887789622652,-0.8192433893688775,-0.3296992672298233,-0.7183270235067828,1.0575461503430779,21,0.2783251661730618 +-0.49550282004135227,1.117843950586921,0.9226522729910517,0.7185126499504312,0.47080153674852426,21,0.6471631092187962 +-9.566769148251744,-2.528722684186441,0.9297486940142597,-0.2572341183913677,-2.1377871966138424,21,0.7460689107531396 +2.4750760941376684,1.7426249608330302,1.518766334807177,2.7382393569028616,-0.23368896847456896,22,1.7884577970031423 +2.8709816307735636,-0.2673740826207912,-0.18881883720680706,-0.7632264460806376,1.1708629395301626,22,0.2697620089148394 +7.866754799646828,2.530301071080695,1.4737151731515785,0.5781067434083211,0.9607310440960617,22,0.3031464442276401 +5.27411582208624,0.43640380976669446,-1.3128189757604696,0.4393421619402547,0.6038642570832362,22,1.2599590395922105 +-1.3321533599274982,0.8831108063957337,0.940699268020265,-0.3060359975420481,-0.8467242662479798,22,1.444824385424271 +8.464883885405554,2.483458270892269,1.7018911510333001,0.48088549736014397,1.2645912450095473,23,0.8819629828328642 +2.5426397901068754,1.1605857357716696,-0.012004826133123703,0.18248709864787196,-0.03268161613343976,23,0.7978062746589233 +1.5234231434489593,0.1240996107160599,-1.3783682448329218,-0.976929364129759,-0.016570418593109336,23,1.124723792631013 +-3.652566838503973,-1.820133991426799,-0.29037442970703536,0.041493062219172094,-1.0690076327998437,23,1.447642530265149 +3.3511536548100973,1.4147146189285094,1.1670566340784954,1.819239190415962,1.3972919619989541,23,0.9537783612295092 +-2.7840355898987545,-1.7672135254249501,-1.9893447387667267,-0.04821728765763476,-0.6993672947188914,24,0.6690384228554732 +-3.474757737738534,0.5867553059532915,1.8584181084776343,1.3468158602131868,-1.3360435250396931,24,1.0918189189970273 +-3.5569441686636467,-1.3321756781735246,0.9316195328404271,0.2168038544377314,-1.0199266672386937,24,1.550534303141136 +0.38708381022151417,0.6749593309867109,0.5025533526031402,-1.426929750905541,-0.6753778625402143,24,1.1871484332829412 +4.428334043738305,-0.710889271527996,-1.2952045214859662,-1.2028216972489947,1.1361220891719663,24,1.8939380914001902 +7.3355415436010185,3.8868616538283165,1.0652270795783119,2.38700854658317,0.10833072902051977,25,1.0545535303540419 +5.522128636010948,0.04409106488706987,-0.5663948583519773,-0.27129918776321954,1.721749621622646,25,0.4813479022226952 +-4.506272877310426,-1.6807724230056396,-0.8672385298506878,0.4767250415702574,-0.5433392678863881,25,0.3954131164942386 +-0.9043057503536266,0.06664337734086745,-0.6900138631157284,0.1938036059532553,-0.17598326360529884,25,1.524444912210235 +2.31403609246989,0.21678670394456268,-0.05036354701728859,1.5162389420969513,0.4966525424815544,25,0.4991587124645881 +-3.0197648958121035,-0.7512092050389922,-0.5006774959221059,0.05850239546374439,-1.1221883400494146,26,0.6012824231671631 +-6.827596419751689,-2.6526833246771813,-1.3180606667768042,-0.36135576653645984,-0.564322150884393,26,1.3313376166946127 +1.4886618326297765,-1.4050121738432848,-0.42648501758997626,-0.5127135441166234,0.941667821176125,26,1.2628941372411397 +-8.30929728442482,-1.8278745465922674,0.6296098171555479,-1.1414719878373067,-1.6935628597265313,26,1.084654025170713 +1.795611775321138,0.4511669908130097,0.35014246442008523,1.8271435938713774,0.49713146611618836,26,1.6619827985247322 +-4.600787390576942,-0.3980823341101012,0.02699353776594693,-0.053317480376366404,-1.9124151634071218,27,1.6699635104796668 +0.5858450699144563,-0.08637809983912093,-0.03854299544313712,-0.8349172278922764,-1.2411565066201022,27,0.9084360400606828 +4.018486031087461,-0.7070757081557745,-0.6488863890952399,0.8243527892812863,1.2313969033246235,27,0.32640213985690847 +6.283684321566685,0.7046266064823534,-0.15928562545469804,0.8557438387271517,1.8101324356255462,27,1.3313161908207465 +-1.2269115412309366,-1.2743370233291336,-1.7708253747379483,-0.07205951156669704,-0.5354335516108314,27,0.28635156281060914 +1.1479279951225774,1.1067004333465158,-0.37424052700261723,-0.34021726025251775,-0.9713808282838574,28,1.1428221437198398 +6.177243418433327,1.4742090404237516,1.0245796557515758,0.3805505295757566,1.5418400537790602,28,1.4020920196537487 +0.9995103498505842,-0.38753623441610485,-0.0663615649157623,-0.7084310659122245,0.167684402310061,28,1.2264028216552687 +2.397248777550938,1.386148351693962,1.5168891305301522,0.565713387598629,-0.37826733010230523,28,1.2882736693821808 +-1.7161686160606555,-3.2918877165694886,0.45085702822802937,-1.798510531940951,1.0823608994227083,28,1.7529124402104532 +-2.4882661365981424,-1.8306266200073158,-2.1995915658267315,-0.13278874601279006,-1.0505984085054574,29,1.4681944903573252 +-0.8709149917351522,-1.8029031230604078,-0.19100269756225322,0.5193995398106401,0.0951202881572464,29,1.7537220166942034 +2.4752479729145227,0.05832924664294831,0.1451068391466729,-0.5788078571319755,-0.17469831263216104,29,1.0111687130173972 +3.1886737508943677,0.2820733275466238,0.5790165409714535,0.11007055110046843,0.4239714797766376,29,0.6218401407088382 +3.7924118936326017,-0.24696188394776586,-1.2413909424214062,-0.2191224895361873,-0.5425361012427711,29,1.6413354665814888 +0.4967587598545684,-0.49351717125956895,-0.5482081102201689,-1.2283083323017474,0.4111211721693767,30,1.3031846029627803 +1.3162858444970598,0.43054480281801233,-0.33187876419206924,0.10477131129042651,-0.46421175080593335,30,0.8004143960948834 +-5.971504216337894,-1.6048497937030153,0.7392226604181471,-1.9866281481594055,-1.000944575429672,30,1.007827229663853 +1.0426510049365623,0.8916791278234543,-0.3868251559173378,1.6030746207032154,0.05613826718516512,30,0.8503898198198158 +1.2857515142380158,0.7379502486936043,-0.2465680892147891,0.04205747138095975,-0.6266088808863725,30,0.682461615314174 +6.414756404291507,3.1286009289453416,1.73185040273552,-0.39611978484626553,0.7528719023355489,31,1.5127067755715287 +-6.510366635055064,0.0021317422622247895,0.8795090246717956,-0.5693955599130082,-3.306687283368152,31,1.4141229442246974 +-1.3957247502793773,-0.5634244229959043,-0.10518196626865614,-0.23892336036657807,-1.4164759717967066,31,0.543762409762581 +-1.2636221118280058,-1.526185831610099,-0.955380245656014,0.04991691526086433,-0.6867283679382477,31,1.704084287359247 +5.052506668183763,0.7500005708884747,0.5142086819448389,0.8350771490817094,1.3489011128149766,31,0.25304041939201893 +-2.2209167978103603,-1.509775258474062,-0.6352018740297543,-1.0907764913795124,-1.2418652246393675,32,1.15657743498669 +0.6556348351617173,0.8179329828016897,1.3788080999410215,0.7566500208244716,-0.6224504820616833,32,0.5299683938506706 +-1.9922237424618385,-1.4248733530859017,-1.6004590305006017,-0.22079581486553834,0.8236612499082416,32,1.98941935061742 +-4.816802220039019,-3.8548650974199545,-2.0679709950965006,0.035641668839981454,-0.026326573019543695,32,0.38025264435123196 +2.4801034604546435,0.6265530589052375,-1.3347239109348175,0.5723701698037224,0.8433487697741786,32,1.3238988683112924 +-2.6339559275092173,-1.371129926337003,-1.8796940335376988,0.5299033029279521,-0.22100746045576064,33,1.5484666855368132 +2.738981721786967,0.697056903793119,1.127833848790633,0.4417128686489708,0.40541527041848635,33,0.5073654009257306 +-3.714735578482114,0.11107106032039954,0.5138011391583975,0.5481632683028156,-0.7724186915625151,33,1.9803562028033475 +0.5932688277610381,0.5503471908990123,-1.0989788703764964,-0.4623270833031166,-0.9050035334922378,33,0.37979148934012175 +-5.783259896479529,-1.630996746811028,1.0206751699925425,-0.9107529215336012,-0.509138953625345,33,0.5903371320967248 +-1.8998139151667504,-0.6311359844025606,0.16385852847889537,-0.2817286024535976,-0.5051267014367744,34,0.9724240488652776 +4.8413604343062,-0.19367476819257434,-1.5907246387282596,-0.5197441797921192,0.59548759322631,34,1.521735097811945 +5.0930260274245915,4.051061559224024,0.856706129305663,1.9636107361349704,0.2814023788155969,34,1.51991157525788 +-4.091521614147816,0.6370605004656008,1.7174926333836207,-0.1528978179445471,-1.5583152829338953,34,0.8956500037123452 +-0.9386772350711514,0.21365734199611908,0.45386152571859484,1.941204808781373,-0.6734464269743918,34,1.2137794083732447 +3.144420456334498,1.3871673556894535,1.9415297020990891,-0.33951562455713924,0.4005843891133172,35,0.8576543215109097 +2.7873644880934596,0.6835145989621223,0.13216941342952604,1.2113246121832957,0.5449926337905644,35,1.6015304975293838 +2.3055860729083033,1.4185466525950063,0.7675122670316653,0.5150486734231124,-0.021551331727591844,35,1.1225985610014884 +-0.4286133268455466,-1.894143697697265,-0.5640111869298329,-1.6701212898400648,-0.5115430662716092,35,1.5561261336910372 +2.186230058489132,1.4398852582574735,0.26746232322929175,1.2827893172151272,-0.026888638399345438,35,1.7135602647476749 +6.526438151894126,1.839980520629568,0.37750157814005514,0.09779456944361518,1.7404404294273872,36,1.8930216984017247 +2.3134488557373354,0.2245536629418532,0.506784108206313,-1.1284816625130412,0.8811802376474023,36,0.4336480758504833 +-2.393897307718211,-1.1772833673522736,-1.0062625009788675,0.23468701347821067,0.04141640635625313,36,0.28038464927109064 +0.8713351141697798,-0.6070327201783865,1.4120769516062535,-1.4432888712873622,0.4830672014188797,36,0.6706354285471448 +3.5201254943330085,0.2683032803133399,-0.8776747472953214,0.3612058197618119,0.46823051100270113,36,0.8136045746943357 +5.478743436416232,0.4721190205150302,-0.7373691658865668,-0.42442735045662555,0.615801798911929,37,1.9771155748277094 +-1.3884263787174194,-2.2510686530789834,-1.3740919816046713,-0.209719812972025,0.2568975870621886,37,0.6618159273422525 +2.998903471621987,-0.8459474199539849,-1.36055811180906,-0.1890946836236865,-0.03743779057737937,37,0.5470958892318193 +-2.339749068976543,-1.3130682226941655,-0.401032900287155,-1.3605208921856324,-0.3901701476580143,37,1.4163913981840488 +-1.1466713876647237,-1.6319463497738709,0.24896193188211158,-1.575530797371647,-0.4495190420994745,37,0.9816908304013633 +-1.243176815867205,-1.8775536159385582,-0.1660608246933138,-1.5804780123717725,-0.08428585412733673,38,0.9689244877054256 +0.7701272216054957,-0.36727139311283136,-0.795584021235714,-0.6268811334383697,0.3746113036009825,38,0.9204083089209835 +4.8718325524639345,-0.6271057092041404,-1.1545884704020546,0.7650979474795238,1.7036248605175843,38,0.47940216449348166 +1.6698034331624976,2.3383373993623944,1.2084986129216344,1.3473391082018202,-0.5925593079547427,38,1.0800009970137099 +6.897765535022611,1.8457556134294173,-0.16248141271170635,0.5810705542771605,0.7765814024319553,38,1.8029506852560735 +-0.07990844203979242,-0.9706729711061535,-1.6904021874219677,-0.032903551929972526,-0.31010320561934485,39,1.4757989280712147 +-1.4247057157601288,-0.6471924841286812,1.6089227507071475,-1.80307684108756,-0.03962932459113808,39,0.3724012760881715 +-0.5381963970403275,-0.8133330141960637,-0.5085224369397449,-1.3876873138315784,-0.6162027491181034,39,1.2226578845655742 +-1.3477884403431024,-1.0861695744073199,-0.8404034583450445,-0.4289185026337298,-0.09348043928709296,39,1.242779325697338 +-2.690111923304227,-0.9856874194290158,0.5603349365589029,-0.9122731035737668,-1.6047032428367873,39,0.9142495170157954 +-1.2321725514307205,1.1668158742835728,0.9530620582535467,0.4210271650712709,0.2423748435227483,40,0.27165096885414153 +6.239019846677261,2.800705402953729,0.7971644770829803,0.40563099925028634,1.235709125670592,40,1.714473105913164 +-0.007618532600873751,-1.9475011103081394,-1.6483128356071628,-1.1514774477481735,0.6720369777100385,40,1.7245921161143571 +-8.785890871653441,-1.5944713486889355,0.0933897808834858,0.48599525852730424,-2.909702072318526,40,0.9443249185838204 +1.7205179142772975,0.8865623090850334,1.6358003907155216,1.6354853068711444,0.7548852995702088,40,1.8053032389682304 +0.4929541285406479,-1.303686765024481,-0.7462896729235644,-2.2808933920354377,0.19677746509725427,41,0.8162358508998595 +-4.373364226618549,-0.22723676565139905,1.0993528905896428,-1.8965573017294486,-1.9885764765224399,41,0.38705729299467645 +3.8381077752425536,0.602130102710569,1.0731616596321611,0.16107546689425103,1.0279494673127199,41,0.9580431950837303 +4.121983289872635,1.1248945620345008,1.1916963585158757,-0.742235890703221,0.514886599555834,41,0.6243504868178903 +5.611760307616764,1.2182466779513936,1.1350255009509471,-1.161282285007108,1.142806163802987,41,0.7161605392392139 +3.9100496337974544,1.5009564806769262,-0.16612256878798212,0.23562923382244247,0.7735706643725259,42,1.7044046359379605 +1.0145438145198928,-0.6420162634926119,-0.4740532862069214,-1.1136120212640725,1.3883634534043112,42,1.9029395274568117 +0.6022437560036675,-2.1789096720088237,-1.152483943218276,-2.475734282950941,1.240093549848944,42,1.2398592953550325 +0.46724807628394305,0.2922317511924343,0.545520518597859,-1.9891628085026174,-0.035023882360477744,42,1.0953814082661157 +-1.911789189852102,-0.9529220327923733,-0.16475101442038903,-1.7569827930663346,-0.3253091982708507,42,0.319496554237089 +-4.892445772098854,-1.2648558020846563,0.7966895829687445,0.7030103314497222,-0.13445817987796957,43,0.6970162277049119 +2.8622778522561125,-0.13630690131890394,-0.26834250411441557,-0.6516695999838324,1.175738344019974,43,1.2256356941525848 +-3.858941193899293,-0.9879604285306847,-0.9636977566479143,-0.5289022573291949,0.014403073942865167,43,1.4578588097310923 +2.8201048637810517,2.214453289441427,0.8674999738652763,0.2529462017920496,0.11580024475141427,43,0.5068985002493539 +-4.992177721240758,-1.609996761381251,-1.4273254448535762,0.9405709843164768,-1.2027454321250315,43,0.5605331082436859 +4.355102471457346,1.9534514507957603,0.21811073470446082,-1.1094527156601515,0.09717301218741618,44,1.343639522020031 +-1.904788412849844,-1.1143519646099525,0.17898134155823422,0.41263082030327475,-0.3981201406329871,44,0.708042822788882 +-5.16670349868338,-1.2998847643970155,-0.5965280642818794,-0.04688845093593669,-1.0243860961068005,44,1.97868743822958 +-3.442563503296976,-2.6365256532496026,-2.1496238055727135,-1.060136908630993,0.2541686667742621,44,1.4361854660188185 +0.9349781897792396,0.7917397382313568,-0.6142900351519119,0.43737817719107164,-0.6045893773027468,44,0.5052021494552403 +-2.0115012699677406,0.3391327925869765,-0.7083619573983271,1.1356301753878941,-2.2331434714781486,45,0.8495348782409666 +0.8629880145159039,-0.7662060748548022,-0.21681274674335713,-0.9481573697140789,-0.37848397075203954,45,0.6989067071930125 +-1.693754640676681,-1.4590496322788564,-0.164591586581964,-0.17852144838257958,-0.49843732639733923,45,1.4999683368539636 +3.215494234265117,0.8123677893510237,0.3013606719982764,0.06052379632760912,-0.8419738019126738,45,1.8011778644213896 +-0.8700307216006806,0.50748128681761,0.18079012547111326,-0.0022209374475656777,-1.3065954371118635,45,1.8040844428233653 +-5.392010501454382,-2.149969006170661,-1.112957089322431,2.052489445178177,-1.413120491628515,46,1.1337626838292592 +7.377227684876863,2.328563140153082,1.8249185132629642,1.1982705388086639,0.38747195654796035,46,0.756565200164873 +0.25604429283083696,0.677275336688257,1.0554202565486248,-1.3599791538692199,-0.1445379441873209,46,1.7972195973463552 +1.4496440416039134,-0.980503302395497,-0.6212806319783595,0.6710024416983313,0.18254017372313475,46,1.0410766657164356 +2.70446810180254,-0.2699289675593234,0.26796795186316186,-0.30440079457928276,0.6824460049778844,46,0.7457896496529525 +3.543810781935966,-0.40557424609619996,-1.0247746579750339,0.41942763449445275,1.209611978451004,47,0.9942052692952679 +-0.7410569545711146,-1.9242243218168005,-1.4909379722704446,0.5137714890206001,-0.05472902365413854,47,1.6359390918152672 +1.8579640770930337,0.5261977610590112,0.5777261758905913,0.5818473652704421,0.3981383826147005,47,1.8180554013907764 +-2.6722268152321162,-1.998393887029826,1.076704382246489,-1.4000688814844495,0.9598448881782858,47,1.5354255620282695 +-3.5029230105403846,-0.6240577498316174,-1.454729465942609,0.852958580911974,-1.6685614508250604,47,1.1544781609085053 +0.9091157161103938,-0.7152034779416973,-0.8382870106873932,-0.6023820062163783,0.2679063820276917,48,1.7874805232090698 +3.1237510786746565,1.3338028520931848,0.3513397786476463,1.0900003514552132,0.24743912688502245,48,0.6594130043136364 +-6.814471648408313,0.38115586118524203,2.09600802737824,0.20981393182726657,-1.4363145198051428,48,1.4708651848830643 +-8.701235019274554,-3.635510572598299,-1.240613090774368,-0.9118940426839416,-0.6843845062861176,48,1.3581496701089057 +-3.9724256980137866,-0.9052403969612929,-0.30227154085215935,-0.814817746901905,-1.0033478052194953,48,1.442464017795692 +-1.4051051236291698,0.04045459304102056,-0.3257031841209866,-0.04852060877483198,-0.11905735550405369,49,1.448912579910415 +0.5542895363576127,-1.1695454658068962,-0.3566530453951007,-1.6185440464511336,0.6659538866622805,49,1.006272840130744 +-0.17108758631112608,-0.21315684402329652,-0.6510200708176037,0.265136152543015,0.28713864729869004,49,1.119173618880999 +8.977966162114514,2.1230971985507967,0.6078293958390719,0.07482264482233784,2.8031190764962095,49,0.5420577790784393 +2.3238125237424945,2.8778473961701807,2.2448361970325283,0.519295015085277,0.09901736240529331,49,1.4692055327526687 +5.68908894767185,1.2937672897369257,0.6809618806551178,-0.2387060075224627,-0.016949495324513664,50,1.1056317341675626 +2.689321758666659,-0.9698331418289579,-2.3811710352449134,-0.6872347134320678,-0.6153297846888606,50,1.9826038789183105 +-3.748068848725099,-3.2268474938375973,-0.5888255427820462,-1.6545292209778635,-0.6031740459896611,50,1.5814550393484077 +8.689883763393988,1.6406725849666035,0.5352070807387453,0.46497183665822994,1.0006992581569747,50,1.3978101421469433 +9.139753744898728,0.5821385551484934,-2.403264575370268,-0.26246219661028347,1.1603184874754868,50,1.8564370821878002 +2.5958699159027567,1.4961147795689906,1.0410018290273593,-0.16009524497964253,-0.029328440048627172,51,0.34711629903813634 +0.7020367680843776,0.08164737117060628,-0.5510628828912335,0.6307225222712023,-1.1528412589405401,51,1.8399033101509825 +6.391998357119437,-0.058388286799849975,-1.8352187512769618,-0.8091856145446071,0.7339656794017347,51,1.583745346391374 +4.769128870672328,0.21748780855122896,0.4925827732839163,-0.2766154367264763,1.042646842991677,51,1.5402194168233803 +2.1674026849801984,-2.849047429733865,-1.6882506954190402,-1.5791184146165067,1.570839875840887,51,1.3767365260501674 +-1.4805016365618595,1.0997086296795273,1.6562685135870694,-0.1119581410696217,-2.4156120933892806,52,1.6867340433265066 +2.9047405083723272,0.6376117300207831,1.3384696796582944,-0.8058188736831878,0.5708016642234405,52,1.945477912558558 +0.15308588215028976,-0.0003385032538973931,-0.6209502519426773,1.2897746441442113,-0.8905582697814965,52,0.907596937968671 +-5.034376550991284,-2.775798245432158,-0.3319692744744976,-0.6926072122964692,-0.858849376606259,52,0.598211615017449 +-0.08484310914864543,0.5489833867642129,0.8070923073198584,0.9509039319545639,-0.6806687977792959,52,0.8604240246854583 +2.64030833122583,-1.4221125239661518,-1.4979742374525198,-0.7458892410951394,0.6640481338870257,53,0.4201052708612847 +4.756978352726562,2.6968726249176083,0.5487107618463026,1.8220101955200685,0.2742427677458051,53,1.893908485632363 +2.477220122051279,1.3020262206693674,-1.2351943921140671,1.6318201960251957,-0.6792085175743192,53,0.2874028765011786 +4.527176130644375,0.36582049921731585,-0.13558152945965646,0.5779599621846817,1.6034210603315222,53,1.4259490022414283 +1.4009658671628111,-0.1565960525547,-0.5040728244391428,1.6427131579308971,0.2814352481578166,53,0.9770263158129711 +-1.9600935861494024,-0.5320432533250777,0.10291834413897047,-0.39250555240758744,-1.0667717930845262,54,1.2057043427463459 +-2.4318096457111245,1.0368018825473246,0.5267387479957301,0.7543508334473852,-1.8453420107589362,54,0.5059228018431419 +1.8168124545173154,1.629999658834049,0.5863944186577773,1.097032012159575,-1.200153691056126,54,0.6152081698960475 +0.26315407315775174,1.045642316400072,-1.544827439801125,0.6543595519088503,-0.6053233269655924,54,0.9807400600012495 +-5.175497019241831,-1.7389084576099387,-0.3442070300377855,1.437332051712372,-0.9768082254110931,54,1.2167241137217055 +1.3276406792249429,1.104040395526299,1.7513379000157343,-0.5731081347314152,-1.0015948009063924,55,1.2361902600935517 +-2.2969868776577567,-0.09931464263766132,0.2880178356014904,-0.5214322226033576,-1.8625246738148322,55,1.5612391899009912 +-1.3206015423644626,-0.10216513632020768,-0.35006545540343287,0.7324209656196093,-0.643951650068028,55,1.3555706222486867 +2.6657899885902223,-0.1882057763615128,-1.119089037677491,0.7640505665777328,-0.10884935419944351,55,0.39032413543678246 +-9.460497272390601,-2.15607929067495,-1.6563266245739119,0.061175879001195255,-2.337466657182652,55,1.0021552124015694 +-2.402936448157052,-2.2270751964211852,-1.6646556237776047,0.40909938244913446,0.6855773685894414,56,1.0978358950512785 +-1.7323779066780358,-1.0771709817991515,-0.4839346787481199,-0.7426637799215418,-0.07258400905730082,56,0.28203930323815835 +2.1139740682335977,1.4137971144535568,1.6892310636777486,-0.026511306507893587,-0.29540580553344375,56,0.5058646645549704 +-7.753336445432895,-4.444799462253192,-1.896677951874554,-0.5345880800715095,-0.6370659063354339,56,1.5312343076508415 +2.685002208221723,0.6463579198672376,1.1026696008135948,-0.1659575304275721,0.796257071143577,56,0.6980101477029222 +1.036119214211098,1.0882589392145414,0.6876438284191397,-1.0901355664745063,-0.5290251597061963,57,1.8625943942163103 +-0.5809148514038966,-1.5844824771164352,0.0963785992729772,-0.5777512918851638,-0.15844935658592782,57,1.3063186377342444 +1.9209526373346275,0.6120715216253134,-0.4335191046359218,0.16471960849365547,-0.8019584498673994,57,1.1847098269435539 +3.9309328328856226,2.1806855412024384,1.0498470167482439,-0.5350534531512635,0.3229259543953575,57,0.3549939468936314 +0.10284566675372858,0.5487181860470249,0.7728837223526349,0.6999400432050135,-0.8376891228981093,57,1.4410947772885068 +0.8352868499284816,-0.9151323558003258,-0.4008436843272913,-0.510346353910523,1.3872420589607772,58,1.595286323156781 +-1.7182639852429005,-0.004637739893034162,0.9460220607225864,-1.4349020474613015,-0.30422261464977385,58,1.0901122389554647 +-7.847957421729675,-2.587495129810671,-0.6088972848215375,-1.3404169287016776,-1.854602834051177,58,1.7585420822896238 +1.876526579373426,-0.22379256570973938,-1.583195525130429,0.41142334984097023,1.3073950325973713,58,0.36828649046530154 +4.334019865086071,1.1710003043447947,0.4450350453488521,0.8096621595519858,0.7869907837176355,58,0.9007818821905731 +-2.24929829145176,-0.5434602566644398,1.1574219739196652,-0.04922082469498444,-0.22680636688093317,59,1.5624645229352452 +2.0272590920044875,1.576251449678759,-0.09540638546458373,-0.040271698733606844,0.3448300447491122,59,1.9414578472323252 +-1.824354878951385,-0.5553514626755143,-0.3036222658225003,-1.0208577120014233,-0.2964973053874775,59,1.8963932348032277 +-1.543423779203957,-0.08261261065342684,0.3223263863411442,0.6957585157288493,0.21477601026103457,59,0.8765948698668644 +-4.6321345858895215,-1.6304986777681814,-0.4495616616593726,-0.5951422085516106,-0.9562845948066991,59,1.1584432491714636 diff --git a/linearmodels/tests/iv/test_centered_j_stat.py b/linearmodels/tests/iv/test_centered_j_stat.py new file mode 100644 index 0000000000..1181d2c58c --- /dev/null +++ b/linearmodels/tests/iv/test_centered_j_stat.py @@ -0,0 +1,514 @@ +""" +Tests of the centered J statistic, IVGMMResults.centered_j_stat + +The reference values in results/centered-j-reference.csv were produced by +results/centered-j-reference.R, a base R implementation of the definition in +Hall (2000) and Hansen and Lee (2021, Theorem 1 and Example 1). Before it +writes the reference values the script checks that it reproduces results from +Stata that are used elsewhere in the tests: ivregress gmm on the housing and +simulated data, including the Hansen J statistic of overidentified models for +the robust, unadjusted, clustered and kernel weight matrices with and without +centering. It also checks the properties of the statistic that follow from +theory in every scenario where they apply. Stata does not have a centered J +statistic, so the values of the statistic are not from Stata. + +The values were also checked against the R packages gmm and AER by +results/centered-j-package-check.R, which is not part of the tests. At the +estimates in the reference the statistics are identical to numerical precision +for all of the models that the packages can estimate. + +The scenarios include misspecified models, where the J statistic that uses an +uncentered covariance is bounded by the sample size and the centered statistic +is not, and every estimator of the covariance of the moment conditions with and +without centering and small-sample adjustments, weights, two-step, iterated and +continuously updating GMM. + +References +---------- +Hall, A. R. (2000). Covariance matrix estimation and the power of the +overidentifying restrictions test. Econometrica, 68(6), 1517-1528. + +Hansen, B. E. and Lee, S. (2021). Inference for iterated GMM under +misspecification. Econometrica, 89(3), 1419-1447. +""" + +import os + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest +from scipy import stats + +from linearmodels.iv import IVGMM, IVGMMCUE +from linearmodels.iv.gmm import KernelWeightMatrix +from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic + +CWD = os.path.split(os.path.abspath(__file__))[0] +RESULTS = os.path.join(CWD, "results") + +SIMULATED = pd.read_stata(os.path.join(RESULTS, "simulated-data.dta")) +SIMULATED["const"] = 1.0 +# Same clusters as the R script, as many as there are moments in some models +SIMULATED["block6"] = np.arange(SIMULATED.shape[0]) // 100 +HOUSING = pd.read_csv(os.path.join(RESULTS, "housing.csv"), index_col=0) +HOUSING["const"] = 1.0 +# Same names as the columns created by model.matrix in the R script +REGIONS = pd.get_dummies(HOUSING.region, drop_first=True).add_prefix("region") +HOUSING = pd.concat([HOUSING, REGIONS.astype(float)], axis=1) +# The relationship between y and x is misspecified, see the generating script +MISSPECIFIED = pd.read_csv(os.path.join(RESULTS, "misspecified-data.csv")) +MISSPECIFIED["const"] = 1.0 +DATASETS = {"sim": SIMULATED, "housing": HOUSING, "mis": MISSPECIFIED} + +REFERENCE = pd.read_csv( + os.path.join(RESULTS, "centered-j-reference.csv"), + keep_default_na=False, + na_values={"centered_j": ["NA"]}, +) +RTOL = {"two_step": 1e-6, "iterated": 1e-6, "cue": 1e-6} +# The CUE objective is flat near its minimum in some scenarios, where the +# estimates are not as precise as the statistics +RTOL_PARAMS = {"two_step": 1e-8, "iterated": 1e-6, "cue": 1e-4} +# The iterated estimator converges to the same estimates as in R, which are +# found with a tolerance of 1e-13. A tolerance of 1e-14 here is reached within +# 50 iterations in all scenarios, whereas 1e-20 requires 1500 in one +ITERATED_TOL = 1e-14 +# The default optimizer of IVGMMCUE, BFGS with numerical derivatives, does not +# find the minimum in some of the scenarios, e.g., it stops at J=1.66 in +# sim_e2_i4 where the minimum is 1.29 whatever the tolerance. These settings +# find the minima that R finds +CUE_OPTIONS = { + "method": "Nelder-Mead", + "options": {"xatol": 1e-10, "fatol": 1e-14, "maxiter": 20000, "maxfev": 40000}, +} + + +def build_model(row, model=IVGMM): + """Create the model described by a row of the reference table""" + data = DATASETS[row.dataset] + weights = data.weights if row.weighted else None + # The way that the covariance of the moment conditions is estimated + config = { + "weight_type": row.weight_type, + "center": bool(row.center), + "debiased": bool(row.debiased), + } + if row.weight_type == "kernel": + config.update({"kernel": "bartlett", "bandwidth": int(row.bandwidth)}) + elif row.weight_type == "clustered": + config["clusters"] = pd.factorize(data[row.clusters])[0] + return model( + data[row.dependent], + data[row.exog.split()], + data[row.endog.split()], + data[row.instruments.split()], + weights=weights, + **config, + ) + + +def fit_model(row): + """Estimate the model of a row of the reference table""" + if row.method == "cue": + return build_model(row, IVGMMCUE).fit(opt_options=CUE_OPTIONS) + if row.method == "iterated": + res = build_model(row).fit(iter_limit=1000, tol=ITERATED_TOL) + assert res.iterations < 1000 + return res + return build_model(row).fit() + + +def reference_rows(): + return [ + pytest.param(row, id=f"{row.id}-{row.method}") + for row in REFERENCE.itertuples(index=False) + ] + + +@pytest.mark.parametrize("row", reference_rows()) +def test_centered_j_against_r(row): + res = fit_model(row) + rtol = RTOL[row.method] + # The estimates, so that the statistics are those of the same estimator + params = np.array(row.params.split(), dtype=float) + assert_allclose(np.asarray(res.params), params, rtol=RTOL_PARAMS[row.method]) + # The J statistic is the one of the same model and estimator + assert_allclose(res.j_stat.stat, row.j, rtol=rtol) + assert res.j_stat.df == row.df + centered = res.centered_j_stat + if np.isnan(row.centered_j): + assert isinstance(centered, InvalidTestStatistic) + assert np.isnan(centered.stat) + assert np.isnan(centered.pval) + return + assert type(centered) is WaldTestStatistic + assert centered.df == row.df + assert_allclose(centered.stat, row.centered_j, rtol=rtol) + assert_allclose( + centered.pval, stats.chi2.sf(row.centered_j, row.df), rtol=max(rtol, 1e-5) + ) + assert centered.dist_name == f"chi2({row.df})" + + +def test_reference_covers_the_cases_that_matter(): + # Guard against the table losing the cases that distinguish the statistics + assert set(REFERENCE.method) == {"two_step", "iterated", "cue"} + assert set(REFERENCE.weight_type) == {"robust", "unadjusted", "kernel", "clustered"} + assert set(REFERENCE.dataset) == {"sim", "housing", "mis"} + for column in ("weighted", "center", "debiased"): + assert REFERENCE[column].any() + assert not REFERENCE[column].all() + defined = REFERENCE.dropna(subset=["centered_j"]) + # Strongly different, in both directions, and strongly rejecting + assert (defined.centered_j > 1.5 * defined.j).sum() >= 10 + assert (defined.centered_j < 0.95 * defined.j).sum() >= 3 + pvals = stats.chi2.sf(defined.centered_j, defined.df) + assert (pvals < 1e-6).sum() >= 10 + assert (pvals > 0.1).sum() >= 3 + # J is bounded by the number of observations, and the new statistic is not + # Unadjusted models are excluded since their statistic is not + # affected by the mean of the moment conditions + iterated = defined[ + (defined.method == "iterated") & (defined.weight_type != "unadjusted") + ] + assert (iterated.j < iterated.nobs).all() + assert REFERENCE.centered_j.isna().sum() == 1 + + +# ---------------------------------------------------------------------------- +# Independent calculations that do not use any of the weight matrix classes +# ---------------------------------------------------------------------------- + + +def moments(res, data, dep, exog, endog, instruments, weighted=False): + """Moment conditions at the estimates, computed from the data""" + sw = np.ones(data.shape[0]) + if weighted: + sw = np.sqrt(data.weights / data.weights.mean()).to_numpy() + y = data[dep].to_numpy() * sw + x = data[exog + endog].to_numpy() * sw[:, None] + z = data[exog + instruments].to_numpy() * sw[:, None] + e = y - x @ np.asarray(res.params) + return z * e[:, None] + + +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("iter_limit", [1, 2, 3, 100]) +def test_centered_j_matches_covariance_of_the_moments(iter_limit, weighted): + data = MISSPECIFIED + spec = ("y", ["const"], ["x"], ["z1", "z2", "zb"]) + weights = data.weights if weighted else None + res = IVGMM(data.y, data[spec[1]], data[spec[2]], data[spec[3]], weights=weights) + res = res.fit(iter_limit=iter_limit) + g = moments(res, data, *spec, weighted=weighted) + gbar = g.mean(0) + # np.cov centers the moments and does not adjust for the number of terms + s_c = np.cov(g, rowvar=False, bias=True) + expected = g.shape[0] * gbar @ np.linalg.solve(s_c, gbar) + assert_allclose(res.centered_j_stat.stat, expected, rtol=1e-10) + + +@pytest.mark.parametrize("weighted", [False, True]) +def test_centered_j_cluster_formula(weighted): + data = MISSPECIFIED + spec = ("y", ["const"], ["x"], ["z1", "z2", "zb"]) + weights = data.weights if weighted else None + clusters = data.group.to_numpy() + mod = IVGMM( + data.y, + data[spec[1]], + data[spec[2]], + data[spec[3]], + weights=weights, + weight_type="clustered", + clusters=clusters, + ) + res = mod.fit() + g = moments(res, data, *spec, weighted=weighted) + gbar = g.mean(0) + sums = pd.DataFrame(g - gbar).groupby(clusters).sum().to_numpy() + s_c = sums.T @ sums / g.shape[0] + expected = g.shape[0] * gbar @ np.linalg.solve(s_c, gbar) + assert_allclose(res.centered_j_stat.stat, expected, rtol=1e-10) + # and this differs from using the heteroskedasticity robust estimator + robust = IVGMM(data.y, data[spec[1]], data[spec[2]], data[spec[3]]).fit() + assert abs(res.centered_j_stat.stat - robust.centered_j_stat.stat) > 1 + + +@pytest.mark.parametrize("weight_type", ["robust", "unadjusted", "kernel", "clustered"]) +@pytest.mark.parametrize( + ("cov_type", "cov_config"), + [ + ("robust", {}), + ("unadjusted", {}), + ("kernel", {"kernel": "parzen", "bandwidth": 2}), + ("clustered", {"clusters": "group"}), + ], +) +def test_centered_j_uses_the_weight_estimator_not_the_covariance_estimator( + weight_type, cov_type, cov_config +): + # The statistic tests the moment conditions with the covariance estimator + # of the weight matrix, as j_stat does, and not with the estimator used for + # the covariance of the parameters + data = MISSPECIFIED + config = {"weight_type": weight_type} + if weight_type == "kernel": + config.update({"kernel": "bartlett", "bandwidth": 3}) + elif weight_type == "clustered": + config["clusters"] = data.group.to_numpy() + cov_config = { + key: data[value].to_numpy() if key == "clusters" else value + for key, value in cov_config.items() + } + mod = IVGMM(data.y, data.const, data.x, data[["z1", "z2", "zb"]], **config) + base = mod.fit() + res = mod.fit(cov_type=cov_type, **cov_config) + assert_allclose(res.centered_j_stat.stat, base.centered_j_stat.stat, rtol=1e-12) + assert_allclose(res.j_stat.stat, base.j_stat.stat, rtol=1e-12) + + +def test_centered_j_does_not_change_the_model(): + data = MISSPECIFIED + mod = IVGMM( + data.y, + data.const, + data.x, + data[["z1", "z2", "zb"]], + weight_type="kernel", + bandwidth=3, + ) + first = mod.fit() + config = first.weight_config + assert config["center"] is False + second = mod.fit() + assert second.weight_config == config + assert second.weight_config["bandwidth"] == 3 + assert_allclose(second.centered_j_stat.stat, first.centered_j_stat.stat) + assert_allclose(second.j_stat.stat, first.j_stat.stat) + + +@pytest.mark.parametrize("kernel", ["bartlett", "parzen", "qs"]) +def test_centered_j_keeps_the_options_of_the_weight_estimator(kernel): + # A bandwidth that is chosen using the data, optimal_bw, is chosen for the + # centered estimator too. This is not reported by the configuration of + # the estimator, which has the bandwidth that it selected + data = SIMULATED + dep, exog = "y_kernel", ["const", "x3"] + endog, instr = ["x1", "x2"], ["z1", "z2", "x4", "x5"] + results = {} + for optimal_bw in (True, False): + mod = IVGMM( + data[dep], + data[exog], + data[endog], + data[instr], + weight_type="kernel", + kernel=kernel, + bandwidth=None, + optimal_bw=optimal_bw, + ) + results[optimal_bw] = res = mod.fit() + g = moments(res, data, dep, exog, endog, instr) + estimator = KernelWeightMatrix( + kernel=kernel, bandwidth=None, center=True, optimal_bw=optimal_bw + ) + x = data[exog + endog].to_numpy() + z = data[exog + instr].to_numpy() + eps = (data[dep].to_numpy() - x @ np.asarray(res.params))[:, None] + s_c = estimator.weight_matrix(x, z, eps) + gbar = g.mean(0) + expected = g.shape[0] * gbar @ np.linalg.solve(s_c, gbar) + assert_allclose(res.centered_j_stat.stat, expected, rtol=1e-10) + # The option is relevant for these data + chosen, default = results[True], results[False] + assert abs(chosen.centered_j_stat.stat - default.centered_j_stat.stat) > 1e-3 + + +# ---------------------------------------------------------------------------- +# Properties +# ---------------------------------------------------------------------------- + + +@pytest.mark.parametrize("weighted", [False, True]) +@pytest.mark.parametrize("data", [MISSPECIFIED, SIMULATED], ids=["mis", "sim"]) +def test_iterated_gmm_identity(data, weighted): + # With iterated GMM and the robust estimator, S = S_c + gbar gbar', and + # so, by the Sherman-Morrison formula, J_c = J / (1 - J / n). J can then + # not exceed n, and J_c is not bounded + if data is MISSPECIFIED: + spec = ("y", ["const"], ["x"], ["z1", "z2", "zb"]) + else: + spec = ("y_robust", ["const", "x3"], ["x1", "x2"], ["z1", "z2", "x4", "x5"]) + weights = data.weights if weighted else None + dep, exog, endog, instr = spec + mod = IVGMM(data[dep], data[exog], data[endog], data[instr], weights=weights) + res = mod.fit(iter_limit=1000, tol=ITERATED_TOL) + nobs = data.shape[0] + j = res.j_stat.stat + assert j < nobs + assert res.centered_j_stat.stat > j + assert_allclose(res.centered_j_stat.stat, j / (1 - j / nobs), rtol=1e-8) + + +def test_both_statistics_are_proportional_to_the_sample_size(): + # With a fixed violation of a moment condition both J / n and J_c / n have + # limits, which are a / (1 + a) and a, and so both statistics diverge at + # rate n. J < n always, but this is a bound on the statistic and not on how + # it grows with the sample size + rng = np.random.default_rng(20260628) + nobs_values = (1000, 4000, 16000) + ratios = [] + for nobs in nobs_values: + z = rng.standard_normal((nobs, 3)) + v = rng.standard_normal(nobs) + u = 0.5 * v + rng.standard_normal(nobs) + x = z @ np.array([0.8, 0.5, 0.4]) + v + y = 1 + x + 0.5 * z[:, 2] + u + res = IVGMM(y, np.ones((nobs, 1)), x, z).fit(iter_limit=1000, tol=ITERATED_TOL) + j, centered = res.j_stat.stat, res.centered_j_stat.stat + assert j < nobs + ratios.append((j / nobs, centered / nobs)) + ratios = np.array(ratios) + # Not changing with the sample size, within sampling error + assert np.all(ratios.max(0) < 1.3 * ratios.min(0)) + # a / (1 + a) and a + assert_allclose(ratios[:, 0], ratios[:, 1] / (1 + ratios[:, 1]), rtol=1e-6) + # which is a large statistic in the largest sample + assert ratios[-1, 1] * nobs_values[-1] > 100 + + +@pytest.mark.parametrize("debiased", [False, True]) +def test_unadjusted_is_not_changed_by_centering(debiased): + # The estimator of the covariance is already centered + data = SIMULATED + mod = IVGMM( + data.y_unadjusted, + data[["const", "x3"]], + data[["x1", "x2"]], + data[["z1", "z2", "x4", "x5"]], + weight_type="unadjusted", + debiased=debiased, + ) + res = mod.fit(iter_limit=100) + assert_allclose(res.centered_j_stat.stat, res.j_stat.stat, rtol=1e-8) + + +@pytest.mark.parametrize("center", [True, False]) +def test_cue(center): + data = MISSPECIFIED + mod = IVGMMCUE(data.y, data.const, data.x, data[["z1", "z2", "zb"]], center=center) + res = mod.fit() + j = res.j_stat.stat + nobs = data.shape[0] + if center: + # The weight matrix of the estimator is the one of the statistic + assert_allclose(res.centered_j_stat.stat, j, rtol=1e-10) + else: + # The weight and the moments are evaluated at the same parameters, + # and so this holds exactly and not only at a minimum + assert_allclose(res.centered_j_stat.stat, j / (1 - j / nobs), rtol=1e-10) + # The default of the CUE estimator is to center + default = IVGMMCUE(data.y, data.const, data.x, data[["z1", "z2", "zb"]]).fit() + assert_allclose(default.centered_j_stat.stat, default.j_stat.stat, rtol=1e-10) + + +def test_result_types(): + data = MISSPECIFIED + exog, endog, instr = data.const, data.x, data[["z1", "z2", "zb"]] + for res in ( + IVGMM(data.y, exog, endog, instr), + IVGMMCUE(data.y, exog, endog, instr), + ): + centered = res.fit().centered_j_stat + assert type(centered) is WaldTestStatistic + assert centered.df == 2 + assert 0 <= centered.pval <= 1 + assert centered.dist_name == "chi2(2)" + assert "Expected moment conditions are equal to 0" in str(centered) + assert "Centered J-test" in str(centered) + assert set(centered.critical_values) == {"10%", "5%", "1%"} + + +def test_exactly_identified(): + # As for j_stat, there is nothing to test. The statistic is 0 up to + # numerical error and has no degrees of freedom + data = MISSPECIFIED + res = IVGMM(data.y, data.const, data.x, data[["z1"]]).fit() + assert res.centered_j_stat.df == res.j_stat.df == 0 + assert_allclose(res.centered_j_stat.stat, 0, atol=1e-8) + assert_allclose(res.j_stat.stat, 0, atol=1e-8) + assert "chi2(0)" in str(res.centered_j_stat) + assert "R-squared" in str(res.summary) + + +@pytest.mark.parametrize("nclusters", [4, 5, 20]) +def test_too_few_clusters(nclusters): + # The centered cluster sums add to zero, so that the covariance has at + # most nclusters - 1 non-zero eigenvalues and the model has 4 moments + data = MISSPECIFIED + clusters = np.arange(data.shape[0]) % nclusters + mod = IVGMM( + data.y, + data.const, + data.x, + data[["z1", "z2", "zb"]], + weight_type="clustered", + clusters=clusters, + ) + res = mod.fit() + centered = res.centered_j_stat + assert np.isfinite(res.j_stat.stat) + if nclusters > 4: + assert type(centered) is WaldTestStatistic + assert np.isfinite(centered.stat) + return + assert isinstance(centered, InvalidTestStatistic) + assert np.isnan(centered.stat) + assert np.isnan(centered.pval) + assert "singular" in str(centered) + assert "Centered J-test" in str(centered) + + +def test_summary_is_not_changed(): + # The statistic is an attribute and does not replace anything in the + # header of the summary, which has the R-squared and F-statistic of the + # model as for other IV models + data = MISSPECIFIED + res = IVGMM(data.y, data.const, data.x, data[["z1", "z2", "zb"]]).fit() + summary = str(res.summary) + for label in ("R-squared:", "Adj. R-squared:", "F-statistic:", "P-value (F-stat)"): + assert label in summary + assert "Centered" not in summary + assert "Endogenous: x" in summary + assert "Instruments: z1, z2, zb" in summary + assert max(len(line) for line in summary.split("\n")) <= 78 + + +@pytest.mark.slow +def test_size_and_power_of_the_test(): + # The statistic has the right size, and is not less powerful than J + rng = np.random.default_rng(20260628) + nobs, reps = 250, 1000 + reject = {violation: np.zeros((reps, 2), dtype=bool) for violation in (0.0, 0.3)} + critical = stats.chi2.ppf(0.95, 2) + for rep in range(reps): + z = rng.standard_normal((nobs, 3)) + v = rng.standard_normal(nobs) + u = (1 + 0.5 * np.abs(z[:, 0])) * (0.5 * v + rng.standard_normal(nobs)) + x = z @ np.array([0.6, 0.5, 0.4]) + v + for violation, rejected in reject.items(): + y = 1 + x + violation * z[:, 2] + u + res = IVGMM(y, np.ones((nobs, 1)), x, z).fit(iter_limit=100) + rejected[rep] = [ + res.j_stat.stat > critical, + res.centered_j_stat.stat > critical, + ] + size = reject[0.0].mean(0) + power = reject[0.3].mean(0) + # The standard error of a rejection frequency is about 0.007 + assert 0.03 < size[1] < 0.075 + assert 0.03 < size[0] < 0.075 + assert power[1] >= power[0] - 0.01 + assert power[1] > 0.5