What endogeneity is
A regressor is endogenous when it is correlated with the error term of the regression it appears in. That correlation arises for several familiar reasons: a variable that drives both the regressor and the outcome is left out of the model, the regressor and the outcome are determined together rather than one causing the other, or the regressor is measured with error. Whatever the source, ordinary least squares then estimates the wrong thing. The coefficient on the endogenous regressor absorbs part of the correlation between the regressor and the error, so it no longer estimates the causal effect, and the bias does not go away as the sample grows.
The standard fix is an instrumental variable, something that shifts
the endogenous regressor but has no direct effect on the outcome. Good
instruments are hard to find in practice. Many candidates fail the
exclusion restriction, and even a valid instrument can be so weakly
correlated with the regressor that the resulting estimates are imprecise
to the point of being useless. Park and Gupta (2012) proposed a way to
correct for endogeneity without an instrument at all, and
endogCopula implements that approach and several later
refinements of it.
The copula idea in plain terms
The instrument-free approach works by modeling the joint distribution of the endogenous regressor and the structural error directly, using a Gaussian copula to describe how the two move together. A copula separates the dependence between two variables from their individual shapes, so the correlation between the regressor and the error can be estimated from the data without ever observing the error itself. The method turns the regressor’s own estimated cumulative distribution function into a normal score, adds that score as an extra “copula term” in the regression, and lets its coefficient absorb the part of the error that the regressor was correlated with. What is left in the regressor’s own coefficient is the part of its variation that has nothing to do with the error, which is the part we want.
The method only works under one identifying assumption. The
endogenous regressor must not be normally distributed, while the
structural error is normal. If both were normal, the copula term and the
regressor itself would carry the same information and there would be
nothing to separate them by. A skewed or heavy tailed regressor, such as
a price, a duration, or a count, gives the method something to work
with. validity(), covered in the third article, checks this
assumption directly on a fitted model.
The two-part formula
Every estimator in this package takes the same formula shape:
y ~ endogenous_1 + endogenous_2 | exogenous_1 + exogenous_2Position decides which side of the model a regressor is on. A term
written before the | is treated as endogenous and gets a
copula correction; the same term written after the | is
treated as exogenous and gets none. Everything lm() accepts
works on either side: transformations such as log(x),
interactions, polynomials, factors, and - 1 to drop the
intercept. Endogenous regressors have to be numeric, because each one
needs its own estimated marginal distribution.
A worked example
We simulate data with a known true effect so that we can check
whether each method recovers it. The endogenous regressor p
is drawn from a chi-squared distribution, which is skewed and therefore
satisfies the nonnormality assumption, and it shares a Gaussian noise
component with the error u, which is what makes it
endogenous. The exogenous regressor w is unrelated to the
error.
set.seed(1)
n <- 300
w <- rnorm(n)
u <- rnorm(n)
p <- rchisq(n, df = 3) + 0.3 * u
y <- 1 + 2 * p + 0.5 * w + u
dat <- data.frame(y = y, p = p, w = w)The coefficient on p in the data generating process is
exactly 2. A naive OLS regression that ignores the endogeneity gets it
wrong:
ols <- lm(y ~ p + w, data = dat)
coef(ols)
#> (Intercept) p w
#> 0.7433303 2.0842374 0.5330108
confint(ols)["p", ]
#> 2.5 % 97.5 %
#> 2.035910 2.132564The 95% confidence interval for p does not cover 2, even
though 2 is the true value, and that gap is what endogeneity bias looks
like in practice: not a point estimate that merely happens to be off,
but an interval that is wrong about how sure it should be.
CopRegPG() fits the Park and Gupta (2012) correction. It
adds a copula term for p, built from p’s
estimated marginal CDF, and reports bootstrap standard errors. We set a
seed first because the bootstrap is random and no seed is set
internally.
set.seed(1)
fit <- CopRegPG(y ~ p | w, data = dat, nboots = 99)
coef(fit)
#> (Intercept) p w p_cop
#> 1.1144394 1.9555238 0.5283453 0.3381301
confint(fit, parm = "p")
#> 2.5 % 97.5 %
#> p 1.797888 2.113159The corrected estimate of p sits close to 2, and its
confidence interval covers the true value. The correction did what it is
meant to do. It moved the estimate away from the OLS bias and gave an
interval that is honest about where the true effect lies.
summary()
summary() gives the full picture: the coefficient table
for the augmented regression, the endogeneity measure rho(P,
xi) which is the estimated correlation between the regressor’s
normal score and the structural error’s normal score, and fit statistics
for both the augmented model (with the copula term) and the structural
model (without it).
summary(fit)
#>
#> Copula endogeneity correction: PG (Park & Gupta 2012)
#>
#> Call:
#> CopRegPG(formula = y ~ p | w, data = dat, nboots = 99)
#>
#> Residuals of the augmented regression (u = xi - C gamma):
#> Min 1Q Median 3Q Max
#> -3.230572 -0.650312 -0.008736 0.698965 3.425732
#>
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#> Min 1Q Median 3Q Max
#> -2.85058 -0.69539 -0.09088 0.73995 4.16745
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.11444 0.24097 4.625 3.75e-06 ***
#> p 1.95552 0.08043 24.314 < 2e-16 ***
#> w 0.52835 0.06338 8.336 < 2e-16 ***
#> p_cop 0.33813 0.19418 1.741 0.0816 .
#> ---
#> 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(p*, xi*) 0.3077 0.1549 1.986 0.047 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Fit, on 296 residual degrees of freedom:
#> augmented structural
#> Residual standard error 1.0230 1.0752
#> R-squared 0.9610 0.9569
#> Adjusted R-squared 0.9606 0.9565
#> sigma above is the standard error of the structural model, the one
#> entering xi* = xi / sigma.
#> Standard errors from 99 bootstrap replicates; cdf = "kde.silverman", 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
#> p 9.866 6.729e-24 2.419e-05
#>
#> Correlation of the copula terms with the exogenous regressors
#> (Park & Gupta assume this is zero; 'joint' tests all of them at once,
#> the Holm p value refers to the single largest correlation):
#> max |corr| with p (Holm) joint R2 joint p
#> p -0.01978 w 0.7332 0.0003912 0.733
#> full matrix in summary(object)$diagnostics$exog.correlation.matrix
#>
#> Collinearity of the copula terms (omega near 0 = weakly identified):
#> corr(P, C) omega
#> p_cop 0.9457 0.1054A rho far from zero says the regressor was indeed correlated with the
error, which is why the correction moved the estimate. The
identification diagnostics printed at the bottom of
summary() are covered in depth in the third article,
“Checking the identifying assumptions”.
Working with the fitted model
A fitted copreg object supports the same extractors as
an lm object, plus a few of its own. coef()
returns the augmented regression’s coefficients, structural regressors
first and copula terms last.
coef(fit)
#> (Intercept) p w p_cop
#> 1.1144394 1.9555238 0.5283453 0.3381301confint() gives a confidence interval either from the
normal approximation using the bootstrap standard error
(type = "normal", the default) or from the empirical
quantiles of the bootstrap draws themselves
(type = "percentile"):
confint(fit, parm = "p", level = 0.95, type = "normal")
#> 2.5 % 97.5 %
#> p 1.797888 2.113159
confint(fit, parm = "p", level = 0.95, type = "percentile")
#> 2.5 % 97.5 %
#> p 1.783432 2.07983residuals() has two flavors.
type = "structural" (the default) gives xi = y - mu - P
alpha - W beta, the residual of the causal model before the copula term
is subtracted off, which is what the identification diagnostics in
validity() are built from. type = "augmented"
gives the residual of the regression that was actually fitted, copula
term included.
head(residuals(fit, type = "structural"))
#> 1 2 3 4 5 6
#> 0.8701159 -1.0475642 2.1105549 -0.3381369 1.6868280 1.5952111
head(residuals(fit, type = "augmented"))
#> 1 2 3 4 5 6
#> 0.9809027 -1.1163678 1.7670547 -0.6325607 1.5237370 1.3819524predict() always uses the structural model, y = mu + P
alpha + W beta. Copula terms are controls for endogeneity, not part of
the causal model whose predictions we usually want, so they never enter
a prediction, in or out of sample.
Where to go next
This article covered one estimator, CopRegPG(), on data
built to satisfy its assumptions. Real data are less obliging in two
ways: the endogenous and exogenous regressors are often correlated with
each other, which breaks an assumption PG relies on, and the identifying
assumptions themselves need checking rather than taking on faith. The
article “Choosing an estimator” walks through the estimators that relax
PG’s assumption and explains when to use each one, and the article
“Checking the identifying assumptions” goes through
validity() check by check, with an example where the
assumptions hold and one where they do not.