Gaussian Copula for Endogeneity Correction

An interactive tour of Park and Gupta’s instrument-free control function, and the assumptions that decide whether it works

What a Gaussian copula actually is, how it turns an endogeneity problem into an extra regressor, and why the whole thing collapses when your endogenous variable is normally distributed.
econometrics
endogeneity
copula
interactive
Author

Ashwin Malshe

Published

July 31, 2026

Modified

July 31, 2026

Note

This post is entirely written by Anthropic’s Opus 5 model. My key objective for using AI was to create the visualizations, which are a little beyond my expertise. However, I also asked it to write an explanation for the Gaussian copula and provide more details to complete this post.

Every empirical paper that regresses sales on price eventually runs into the same referee comment: price is endogenous. Firms do not set prices at random. They set them knowing things about demand that you, the analyst, never observe. Whatever those things are, they sit in your error term and they are correlated with the regressor you care most about.

The textbook answer is an instrumental variable. The practical problem with the textbook answer is that good instruments are scarce, and a weak instrument leaves you worse off than doing nothing. So there is real appeal in a method that needs no instrument at all.

Park and Gupta (2012) proposed exactly that. If the endogenous regressor is non-normally distributed, you can construct the correction out of the regressor itself. No exclusion restriction, no outside variable, one extra column in your design matrix.

That sounds too good to be true, and the honest answer is that it is conditionally true. The method rests on assumptions that are easy to state and easy to violate, and the literature since 2012 has been mostly about mapping out where it breaks. This post builds the method up from scratch, with widgets you can push around, and then makes every assumption explicit and shows you what each one looks like when it fails.

Note

The R code uses endogCopula, my package implementing eight copula-based estimators behind a single copreg() entry point.

remotes::install_github("ashgreat/endogCopula")

What a copula actually is

Start with the definition, because “copula” gets used loosely and the loose version will not help you here.

Take any two random variables \(P\) and \(\xi\) with a joint distribution \(F(p, x)\) and marginals \(F_P\) and \(F_\xi\). Sklar’s theorem says that joint distribution always factors into two separate pieces:

\[F(p, x) = C\big(F_P(p),\; F_\xi(x)\big)\]

The marginals \(F_P\) and \(F_\xi\) describe each variable on its own. The function \(C\) describes only how they move together. That function is the copula.

The mechanism is the probability integral transform. If \(P\) is continuous, then \(U = F_P(P)\) is Uniform\((0,1)\) no matter what \(P\) looked like. Same for \(V = F_\xi(\xi)\). The copula is simply the joint distribution of that pair \((U, V)\), both marginals flattened out to uniform, so nothing is left but the dependence.

The consequence worth internalizing: the copula is invariant to strictly increasing transformations of the marginals. Take logs of \(P\), square it, rescale it in dollars instead of euros. The copula does not move. Push the marginal selector below and watch the left panel change shape completely while the middle and right panels stay frozen.

Widget 1 — the copula is the dependence, and nothing else

400 draws. Every marginal here is a strictly increasing transform of the same underlying normal draw, which is why only the first panel can move.

The middle panel is the copula. The right panel is the same thing viewed through \(\Phi^{-1}\), which turns the unit square back into a plane and, when the copula is Gaussian, turns the cloud into a tidy ellipse.

Why the Gaussian copula gives you a regressor

A Gaussian copula is the copula you get by taking a bivariate normal and flattening its marginals:

\[C(u, v; \rho) = \Phi_\rho\big(\Phi^{-1}(u),\, \Phi^{-1}(v)\big)\]

It has exactly one parameter, \(\rho\). Assuming it is a real assumption about the shape of dependence, not a free lunch, and we will come back to that.

Now put it to work. The structural model is

\[y = \beta_0 + \beta_1 P + \beta_2 W + \xi\]

where \(P\) is endogenous and \(W\) is exogenous. Define the normal scores

\[P^* = \Phi^{-1}\big(H(P)\big), \qquad \xi^* = \xi / \sigma_\xi\]

where \(H\) is the marginal CDF of \(P\). Assume (i) \(\xi\) is normal and (ii) the dependence between \(P\) and \(\xi\) is Gaussian-copular. Together these say precisely that \((P^*, \xi^*)\) is bivariate standard normal with correlation \(\rho\).

For a bivariate normal, the conditional mean is linear. So

\[E[\xi^* \mid P^*] = \rho P^* \quad\Longrightarrow\quad E[\xi \mid P] = \rho\,\sigma_\xi\,\Phi^{-1}\big(H(P)\big)\]

Conditioning on \(P\) and on \(P^*\) is the same thing, since \(H\) is monotone. And that is the whole method. The part of the error that depends on \(P\) has a closed form, it is proportional to \(P^*\), so put \(P^*\) in the regression:

\[y = \beta_0 + \beta_1 P + \beta_2 W + \gamma P^* + u, \qquad \gamma = \rho\,\sigma_\xi\]

The endogeneity has been converted into an omitted variable, and then un-omitted. This is a control function, in the same family as Heckman’s inverse Mills ratio, except the control function is built from the regressor’s own empirical distribution instead of from an outside instrument.

Two things fall out for free. \(\hat\gamma\) estimates \(\rho\sigma_\xi\), so testing \(H_0: \gamma = 0\) is a test for endogeneity. And \(\hat\gamma\)’s sign tells you which way naive OLS was wrong.

Widget 2 — the bias, and the correction

True slope is 2. Drag ρ and watch naive OLS walk away from it while the corrected estimate stays put.

Note that \(\hat\gamma\) tracks \(\rho\sigma_\xi\) as you drag, which is the theory above showing up numerically. Note also that at \(\rho = 0\) the two estimates coincide: the correction costs you nothing in bias when there was no endogeneity, only a little precision.

The assumptions

This is the part that matters. The method is not hard to run; it is hard to justify. endogCopula ships a validity() method that walks these one at a time, so I have organised this section to match its four steps.

[1] The endogenous regressor must be non-normal

This is not a regularity condition. It is the entire source of identification, and it is worth seeing why.

Suppose \(P \sim N(\mu, \sigma_P^2)\). Then \(H(p) = \Phi\!\big((p - \mu)/\sigma_P\big)\), so

\[P^* = \Phi^{-1}\big(H(P)\big) = \frac{P - \mu}{\sigma_P}\]

The copula term is an exact linear function of \(P\). Adding it to a regression that already contains \(P\) and an intercept adds a perfectly collinear column. There is nothing to estimate. The design matrix is singular, and no amount of data will fix it.

So the method runs on departure from normality. The further \(P\) is from normal, the more \(\Phi^{-1}(H(\cdot))\) bends away from a straight line, and the more independent variation the copula term carries. Drag the slider below from right to left and watch the curve straighten.

Widget 3 — the identification cliff

P is lognormal with log-scale σ. As σ → 0 it becomes normal, the copula term becomes a straight line in P, and the estimator dies.

That \(\omega\) readout is the same quantity endogCopula prints under “Collinearity of the copula terms”. It is one minus the \(R^2\) from regressing the copula term on everything else in the model, so it measures the independent variation the correction has left to work with. Near zero means weakly identified, and you will see it approach zero long before the design matrix is literally singular.

How non-normal is non-normal enough? Becker, Proksch and Ringle (2022) ran the simulations and give boundary conditions: at \(n = 300\) and 80% power you want \(|\text{skewness}| \geq 1.932\), or an Anderson–Darling statistic above roughly 19. A significant Shapiro–Wilk test is not sufficient. With enough observations, that test rejects normality for departures far too small to identify anything.

One subtlety that trips people up: the non-normality has to survive conditioning on your exogenous regressors. If \(P\) looks skewed only because it inherits skewness from a covariate already in the model, the copula term has nothing left to contribute once \(W\) is partialled out. validity() checks the regressor itself; the residual version is what actually matters, and the \(\omega\) diagnostic is the honest summary.

[2] The copula term must be uncorrelated with the exogenous regressors

Park and Gupta’s derivation assumes the copula terms are uncorrelated with the exogenous regressors \(W\). When that fails, the PG estimator is inconsistent.

This is not exotic. If \(W\) and \(P\) are correlated, which is the ordinary case, then \(P^*\) and \(W\) generally will be too. validity() reports both the individual correlations with a Holm correction and a joint \(F\) test, and tells you which way to go:

  • Assumption satisfied → PG is consistent, and more efficient than the alternatives.
  • Assumption violated → use 2sCOPE (Yang, Qian and Xie 2025), which is built for exactly this case. In endogCopula, method = "2scope".

[3] The structural error should be normal

The derivation needs \(\xi\) normal to get \(E[\xi \mid P] = \rho\sigma_\xi P^*\). Drop normality and the conditional mean is no longer proportional to \(P^*\), so the added regressor is the wrong control function.

Becker et al. find that with a non-normal error, the approach is no longer consistent in models with an intercept. Yang et al. and Qian, Koschmann and Xie do permit a non-normal error, but only under a decomposition \(\xi = U + V\) where \(U\) is normal and carries all the dependence with the regressors, and \(V\) is independent and can be anything.

Here is the honest difficulty, and endogCopula says so in its own output rather than papering over it: your residuals estimate \(\xi\), not \(U\). So residual skewness neither establishes nor rules out a violation. Whether the decomposition holds is an argument you have to make from the suspected source of the endogeneity, not from a diagnostic plot. That is uncomfortable, and it is correct.

[4] Continuous regressor, independent draws, and the price of the correction

Three smaller conditions and one tradeoff.

\(P\) must be continuous. The transform \(\Phi^{-1}(H(P))\) is built on the probability integral transform, which needs a continuous CDF. With a discrete or coarsely rounded regressor, \(H\) steps rather than climbs, ties pile up, and the copula term degrades. Discrete extensions exist and are not well behaved.

Observations must be independent, and \(n\) must be reasonably large. The correction is estimated, not known, so it inherits sampling error, and the bootstrap that endogCopula uses for standard errors needs room to work. Panel data needs the panel formulation, not this one.

Finally, the correction is not free. Adding a near-collinear regressor inflates standard errors. Qian, Koschmann and Xie call this ICON, the ratio of the corrected standard error to the uncorrected one, and suggest 6 as a rough ceiling. Above that you are buying unbiasedness at a price in precision that may not be worth paying.

The Monte Carlo below shows both halves of that trade at once. Naive OLS is tight and wrong; the correction is centred and wider. Push \(\sigma\) down toward normality and watch the corrected histogram spread out until it is useless.

Widget 4 — sampling distributions over 400 replications

One draw can flatter any estimator. This is what actually matters: where the estimates land in repeated samples.

Doing it in R

Simulate the same design: a lognormal endogenous regressor, an exogenous covariate, and a structural error correlated with the regressor at \(\rho = 0.6\).

library(endogCopula)

set.seed(20260731)
n   <- 500
rho <- 0.6

z1 <- rnorm(n)
z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n)   # Gaussian copula, by construction

p  <- as.numeric(scale(exp(z1)))               # lognormal: strongly skewed
w  <- rnorm(n)                                 # exogenous
xi <- z2                                       # structural error
y  <- 1 + 2 * p + 0.5 * w + xi                 # true beta_1 is 2

dat <- data.frame(y, p, w)

Naive OLS first, so we know what we are correcting:

coef(summary(lm(y ~ p + w)))
             Estimate Std. Error  t value      Pr(>|t|)
(Intercept) 0.9672314 0.03609615 26.79597  1.557385e-98
p           2.4795912 0.03612938 68.63087 1.105879e-255
w           0.5088887 0.03723142 13.66826  2.479346e-36

Badly biased, and confident about it. The point estimate is off by about 24%, with a standard error small enough that the true value of 2 is many standard errors away. This is the characteristic damage of endogeneity: not noise, but precision pointed in the wrong direction.

Now the copula correction. The formula is two-part: endogenous regressors before the |, exogenous after.

fit <- copreg(y ~ p | w, data = dat, method = "pg", nboots = 199)
summary(fit)

Copula endogeneity correction: PG (Park & Gupta 2012) 

Call:
copreg(formula = y ~ p | w, data = dat, method = "pg", nboots = 199)

Residuals of the augmented regression (u = xi - C gamma):
     Min       1Q   Median       3Q      Max 
-2.28588 -0.49861 -0.02098  0.47458  2.76658 

Residuals of the structural model (xi = y - mu - P alpha - W beta):
      Min        1Q    Median        3Q       Max 
-3.004899 -0.631978  0.002421  0.562421  3.464059 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.95902    0.04329  22.154  < 2e-16 ***
p            2.03431    0.07721  26.348  < 2e-16 ***
w            0.50285    0.03777  13.313  < 2e-16 ***
p_cop        0.54230    0.07226   7.505 6.13e-14 ***
---
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.56614    0.05282   10.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fit, on 496 residual degrees of freedom:
                        augmented structural
Residual standard error 0.7609    0.9231    
R-squared               0.9179    0.8792    
Adjusted R-squared      0.9174    0.8784    
  sigma above is the standard error of the structural model, the one 
  entering xi* = xi / sigma.
Standard errors from 199 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 38.12 3.7e-24 1.258e-16

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  4.177e-05    w   0.9993 1.745e-09  0.9993
 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.8545 0.2696

Read that output against the derivation:

  • p is back near 2, where it belongs.
  • p_cop is \(\hat\gamma\), the coefficient on the copula term. It is significant, which is your endogeneity test.
  • rho(p*, xi*) recovers \(\rho\) directly, close to the 0.6 we simulated. Multiply it by the structural residual standard error and you get p_cop back, since \(\gamma = \rho\sigma_\xi\).
  • The standard error on p is about twice the OLS one. That is the price, and it is a reasonable one here.

Then check the assumptions rather than assuming them:

validity(fit)

Validity check for PG (Park & Gupta 2012)
n = 500, intercept: yes, target power 80%
Sources: Becker, Proksch & Ringle (2022); Yang, Qian & Xie (2025);
         Qian, Koschmann & Xie (2025)

[1] Nonnormality of the endogenous regressors
  skewness ex.kurtosis    AD  CvM     KS p Yang ok Becker ok
p    3.072       13.95 38.12 6.97 1.26e-16    TRUE      TRUE
    Becker et al. at n = 500: |skewness| >= 1.932, or AD > 18.964, or CvM > 3.488
    Yang et al.: KS p < .05

[2] Assumption: correlation of the copula transformation term
    with the exogenous regressors
  corr(W, CTT) p (Holm)
w            0    0.999
    Joint test: R2 = 1.745e-09, F = 8.69e-07, p = 0.9993
    => satisfied. Park & Gupta is consistent here and more efficient.

[3] Error term: structural residuals xi
    skewness = 0.1015, excess kurtosis = 0.4097, AD = 0.3778 (p = 0.407)
    Becker et al. (2022) find that with a nonnormal error the approach
    is no longer consistent in models with an intercept. Yang et al.
    (2025) and Qian et al. (2025) do permit a nonnormal error, but only
    under the decomposition xi = U + V into a normally distributed
    endogenous part U, which carries the entire dependence with the
    regressors, and an independent nonnormal V. The residuals show xi,
    not U, so their skewness neither establishes nor rules out a
    violation; whether that decomposition holds has to be argued from
    the suspected sources of endogeneity.

[4] ICON: standard error inflation relative to uncorrected OLS
            SE (corrected) SE (uncorrected)   ICON
(Intercept)        0.04329          0.03658 1.1835
p                  0.07721          0.05799 1.3313
w                  0.03777          0.03831 0.9859
    => largest ICON = 1.331, below the threshold of 6.

The blocks map onto the sections above, though which ones appear depends on the estimator’s own identification condition. The non-normality block carries Becker et al.’s thresholds. The next tests whether the copula terms are uncorrelated with w, and tells you to switch estimators if they are not. The error-term block reports diagnostics with the caveat about \(\xi\) versus \(U\). The last is ICON.

The negative control is worth running once so you recognise the failure when you meet it. Make the regressor nearly normal and everything falls apart:

set.seed(20260731)
p_norm <- as.numeric(scale(exp(0.02 * z1)))    # sigma = 0.02: essentially normal
dat2   <- data.frame(y = 1 + 2 * p_norm + 0.5 * w + xi, p = p_norm, w = w)

fit2 <- copreg(y ~ p | w, data = dat2, method = "pg", nboots = 199)
summary(fit2)$coefficients
             Estimate Std. Error   z value     Pr(>|z|)
(Intercept) 0.9673344 0.05330599 18.146824 1.360425e-73
p           1.6542152 0.76655840  2.157977 3.092964e-02
w           0.5029865 0.03533345 14.235419 5.523112e-46
p_cop       0.8996036 0.77691952  1.157911 2.469004e-01

The point estimate wanders off and the standard error on p inflates by roughly a factor of ten against the skewed case, exactly as widget 3 predicted. The copula term is no longer distinguishable from zero, because it is no longer distinguishable from p.

Nothing errored. Nothing warned you at the console. The model fit, returned numbers, and one of them even cleared the 5% threshold. It just did not identify anything. validity() is what tells you so:

validity(fit2)

Validity check for PG (Park & Gupta 2012)
n = 500, intercept: yes, target power 80%
Sources: Becker, Proksch & Ringle (2022); Yang, Qian & Xie (2025);
         Qian, Koschmann & Xie (2025)

[1] Nonnormality of the endogenous regressors
  skewness ex.kurtosis    AD   CvM  KS p Yang ok Becker ok
p   -0.099      -0.068 0.398 0.066 0.799   FALSE     FALSE
    Becker et al. at n = 500: |skewness| >= 1.932, or AD > 18.964, or CvM > 3.488
    Yang et al.: KS p < .05

[2] Assumption: correlation of the copula transformation term
    with the exogenous regressors
  corr(W, CTT) p (Holm)
w      -0.0011    0.981
    Joint test: R2 = 1.199e-06, F = 0.0005972, p = 0.9805
    => satisfied. Park & Gupta is consistent here and more efficient.

[3] Exogenous regressors as identifying variation
    (continuous, KS p < .001, first-stage F > 10)
    Not available: this estimator has no first stage.

[4] Error term: structural residuals xi
    skewness = 0.1051, excess kurtosis = 0.281, AD = 0.4072 (p = 0.347)
    Becker et al. (2022) find that with a nonnormal error the approach
    is no longer consistent in models with an intercept. Yang et al.
    (2025) and Qian et al. (2025) do permit a nonnormal error, but only
    under the decomposition xi = U + V into a normally distributed
    endogenous part U, which carries the entire dependence with the
    regressors, and an independent nonnormal V. The residuals show xi,
    not U, so their skewness neither establishes nor rules out a
    violation; whether that decomposition holds has to be argued from
    the suspected sources of endogeneity.

[5] ICON: standard error inflation relative to uncorrected OLS
            SE (corrected) SE (uncorrected)   ICON
(Intercept)        0.05331          0.03311  1.610
p                  0.76656          0.03996 19.181
w                  0.03533          0.03516  1.005
    => largest ICON = 19.18, above the threshold of 6: weak
       identification or a misspecified dependence model.

When not to use it

Reach for something else when the endogenous regressor is close to normal, discrete, or coarsely measured; when your sample is small; when the copula terms correlate with your exogenous regressors (use 2sCOPE); when you have reason to think the structural error is non-normal in a way that carries the dependence; or when you have a genuinely good instrument, in which case use it.

And treat the Gaussian copula itself as an assumption. It is one specific dependence structure among many. It has no tail dependence, and it is symmetric. If you believe the endogeneity is concentrated in the upper tail, say firms cut prices hardest exactly when they see a demand shock, a Gaussian copula is the wrong shape and no diagnostic in this post will catch it.

endogCopula implements seven other estimators through the same copreg() interface for when PG’s conditions do not hold: 2scope and np for the correlated-exogenous case, bmw for a first-stage-residual identification condition, plus ima and jams. The validity() decision tree points at the right one.

References

Becker, J.-M., D. Proksch, and C. M. Ringle (2022). Revisiting Gaussian copulas to handle endogenous regressors. Journal of the Academy of Marketing Science 50, 46–66. https://doi.org/10.1007/s11747-021-00805-y

Breitung, J., A. Mayer, and D. Wied (2024). Asymptotic properties of endogeneity corrections using nonlinear transformations. The Econometrics Journal 27(3), 362–383. https://doi.org/10.1093/ectj/utae002

Haschka, R. E. (2025). Endogeneity corrections using between-regressor correlations. IMA Journal of Management Mathematics. https://doi.org/10.1093/imaman/dpae018

Hu, X., Y. Qian, and H. Xie (2025). Correcting endogeneity via nonparametric copula control functions. NBER Working Paper 33607. https://doi.org/10.3386/w33607

Park, S. and S. Gupta (2012). Handling endogenous regressors by joint estimation using copulas. Marketing Science 31(4), 567–586. https://doi.org/10.1287/mksc.1120.0718

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. https://doi.org/10.1177/00222437241296453

They run entirely in your browser using Quarto’s built-in Observable JS, with the numerics in a small stats.js module beside this post. No R is involved once the page loads.

One deliberate difference from the R code: the widgets estimate \(H\) with the rescaled empirical CDF, \(\text{rank}/(n+1)\), which is the cdf = "rank.n1" option in endogCopula. The R code uses the package default, a kernel-smoothed CDF with Silverman’s bandwidth, which is what Park and Gupta use in the original article. The distinction does not matter for anything argued here, but it will produce small numerical differences if you compare the two directly.