Copula endogeneity correction of Breitung, Mayer and Wied (2024) (BMW)

Nonlinear endogeneity correction of Breitung, Mayer & Wied (2024). Each endogenous regressor z is regressed by OLS on the exogenous regressors x (always with an intercept, even when the structural model has none, because Assumption A4 requires E(e) = 0); the residual ehat = z - x'deltahat is transformed to a normal score etahat = Phi^-1(Fhat(ehat)) with Fhat(ehat_i) = rank(ehat_i) / (n + 1); and etahat then enters the structural regression of y on x and z. The order of operations is what separates BMW from 2sCOPE (CopReg2sCOPE()): BMW residualises the raw variables first and transforms them to normal scores second, while 2sCOPE transforms first and residualises second. With several endogenous regressors (Remark 2.1) each z is regressed on x on its own, never on the other endogenous regressors, and gets its own transform and its own rho. Without an exogenous part in formula the first stage has nothing to regress on, so the call is redirected to CopRegPG() with a warning.

Usage

CopRegBMW(
  formula,
  data,
  cdf = "rank.n1",
  ties = "max",
  nboots = 199,
  subset = NULL,
  contrasts = NULL,
  parallel = FALSE,
  ncores = NULL,
  verbose = interactive()
)

Arguments

formula: A two-part formula of the form y ~ endog_1 + endog_2 + ... | exog_1 + exog_2 + .... Position decides: terms written before the | are endogenous and receive a copula correction, terms written after it are exogenous and receive none. Everything lm() accepts is allowed: transformations, interactions, polynomials, factors, and -1 to drop the intercept. Endogenous regressors must be numeric. Missing values are handled by na.omit. Writing an interaction or a power of an endogenous regressor before the | gives it its own copula term and produces a warning, because such a term is strongly collinear with the regressor it is built from. BMW's first stage is linear in the exogenous regressors (Assumption A4); the authors' own remedy for a nonlinear relationship is to add nonlinear transformations of the exogenous regressors, e.g. w1:w2 or I(w^2), after the |, which then enter the first stage as written. Without an exogenous part the call is redirected to CopRegPG().
data: A data.frame containing the variables referenced in formula.
cdf: Character string naming the estimator of the marginal CDF of the first-stage residuals. Effectively fixed at "rank.n1", Rank / (n + 1), the transformation of the paper's Equation 2.3: unlike the other cross-sectional estimators in this package the choice is not free here, because the asymptotic theory of Proposition 3.1 is derived specifically for this normalisation. Supplying any other value still produces point estimates but triggers a warning that the reported standard errors are no longer covered by the published theory.
ties: Character string, "max" (default) for the counting function or "average" for midranks.
nboots: Number of bootstrap replicates used to compute the standard errors. Defaults to 199.
subset: An optional logical vector selecting the rows to use, as in lm(). Defaults to NULL (use every row).
contrasts: An optional list of contrasts for factor variables, passed on to model.matrix(). Defaults to NULL.
parallel: FALSE (default), TRUE, "multicore", or "snow" to spread the bootstrap replicates across worker processes.
ncores: Number of worker processes to use when parallel is not FALSE. NULL (default) uses one less than the number of detected cores.
verbose: Logical; show a progress message and progress bar while the bootstrap runs. Defaults to interactive().

Value

An object of class "copreg": a list including coefficients and std.error (point estimates and bootstrap standard errors for the structural and copula coefficients), vcov (the bootstrap covariance matrix of coefficients), rho and rho.se (the endogeneity measure corr(structural residual, copula data) for each endogenous regressor and its bootstrap standard error), residuals and fitted.values (structural model, copula terms excluded), residuals.augmented and fitted.augmented (including the copula terms), boot (raw bootstrap coefficient draws), ols.coefficients and ols.std.error (the uncorrected OLS benchmark from the same bootstrap resamples), icon (the resulting standard error inflation), first.stage.residuals (the ehat that the identification checks are run on), diagnostics, and call, method, cdf, ties. Methods exist for print(), summary(), coef(), vcov(), nobs(), residuals(), fitted(), predict(), confint(), formula(), and update().

Because Corollary 3.2 of Breitung, Mayer & Wied (2024) shows that, under the null rho = 0, the textbook t statistic built from classical OLS standard errors keeps a standard normal limit even though those standard errors are inconsistent for the structural coefficients, summary() on a fit from this function additionally reports that Durbin-Hausman-Wu test ($dhw on the summary object) next to the bootstrap-based one. Use validity() on the result for the identification checks; for this estimator the non-normality requirement that identification depends on applies to the first-stage residuals rather than to the endogenous regressor itself, and validity() tests them accordingly.

References

Breitung, J., A. Mayer, and D. Wied (2024). Asymptotic properties of endogeneity corrections using nonlinear transformations. The Econometrics Journal 27(3), 362-383.

Park, S. and S. Gupta (2012). Handling endogenous regressors by joint estimation using copulas. Marketing Science 31(4), 567-586.

Qian, Y., A. Koschmann, and H. Xie (2025). A practical guide to endogeneity correction using copulas. Journal of Marketing.

Yang, F., Y. Qian, and H. Xie (2025). Addressing endogeneity using a two-stage copula generated regressor approach. Journal of Marketing Research 62(4), 601-623.

Zhao, Y., I. Gijbels, and I. Van Keilegom (2020). Inference for semiparametric Gaussian copula model adjusted for linear regression using residual ranks. Bernoulli 26(4), 2815-2846.

Examples

set.seed(2)
n <- 200
x <- rnorm(n)
w <- rnorm(n)
e <- rexp(n, rate = 1) - 1               # nonnormal first-stage error
z <- 1 + 0.5 * x + 0.5 * w + e           # endogenous regressor
y <- 1 + 2 * z + 0.5 * x + 0.3 * e + rnorm(n, sd = 0.5)
dat <- data.frame(y = y, z = z, x = x, w = w)

fit <- CopRegBMW(y ~ z | x + w, data = dat, nboots = 25)
summary(fit)
#> 
#> Copula endogeneity correction: BMW (Breitung, Mayer & Wied 2024) 
#> 
#> Call:
#> CopRegBMW(formula = y ~ z | x + w, data = dat, nboots = 25)
#> 
#> Residuals of the augmented regression (u = xi - C gamma):
#>      Min       1Q   Median       3Q      Max 
#> -1.47991 -0.32525 -0.00841  0.37096  1.33434 
#> 
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#>      Min       1Q   Median       3Q      Max 
#> -1.38606 -0.30332 -0.01385  0.39400  1.26246 
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.65213    0.12916   5.049 4.45e-07 ***
#> z            2.32327    0.10476  22.177  < 2e-16 ***
#> x            0.31253    0.05731   5.453 4.95e-08 ***
#> w           -0.07146    0.07458  -0.958    0.338    
#> z_cop        0.04984    0.11839   0.421    0.674    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Endogeneity: rho(P*, xi*) is the correlation between the normal score 
#>   of an endogenous regressor and that of the structural error, xi* = xi / sigma.
#>              Estimate Std. Error z value Pr(>|z|)
#> rho(z*, xi*)  0.09064    0.20507   0.442    0.658
#> 
#> Durbin-Hausman-Wu test of no endogeneity, rho = 0. Under that null the
#>   textbook t statistic stays valid even though the copula terms are
#>   generated regressors (Breitung, Mayer & Wied 2024, Corollary 3.2), so
#>   these standard errors are classical rather than bootstrapped:
#>       Estimate Std. Error t value Pr(>|t|)
#> z_cop  0.04984    0.10733   0.464    0.642
#> 
#> Fit, on 195 residual degrees of freedom:
#>                         augmented structural
#> Residual standard error 0.5420    0.5442    
#> R-squared               0.9680    0.9677    
#> Adjusted R-squared      0.9673    0.9671    
#>   sigma above is the standard error of the structural model, the one 
#>   entering xi* = xi / sigma.
#> Standard errors from 25 bootstrap replicates; cdf = "rank.n1", ties = "max".
#> Pr(>|z|) in both tables: Wald test using the normal approximation,
#>   z = Estimate / Std. Error, with the bootstrap standard error.
#>   See confint(object, type = "percentile") for bootstrap percentile intervals.
#> 
#> --- Identification diagnostics ------------------------------------
#> 
#> Non-normality of the endogenous regressors (small p = non-normal, good):
#>      AD     AD p  KS p
#> z 1.316 0.001994 0.371
#> 
#> Correlation of the copula terms with the exogenous regressors
#> (this estimator does not assume it is zero; it projects the correlation
#>  out in the first stage):
#>   max |corr| with p (Holm) joint R2 joint p
#> z    0.08317    x    0.484   0.0114  0.3233
#>  full matrix in summary(object)$diagnostics$exog.correlation.matrix
#> 
#> Collinearity of the copula terms (omega near 0 = weakly identified):
#>       corr(P, C)  omega
#> z_cop     0.7665 0.1335
#>