Copula endogeneity correction of Park and Gupta (2012) (PG)

Fits the linear control-function estimator of Park and Gupta (2012): each endogenous regressor is transformed to a normal score using an estimate of its marginal CDF, and the transformed variable (the "copula term") is added to the structural regression alongside the untransformed regressors. Exogenous regressors are not transformed and do not enter the copula term. Standard errors are obtained by a nonparametric bootstrap.

Usage

CopRegPG(
  formula,
  data,
  cdf = "kde.silverman",
  ties = "max",
  nboots = 199,
  subset = NULL,
  contrasts = NULL,
  parallel = FALSE,
  ncores = NULL,
  verbose = interactive()
)

Arguments

formula: A two-part formula y ~ endogenous_1 + endogenous_2 + ... | exogenous_1 + exogenous_2 + .... Position decides: a term written before the | is endogenous and receives a copula correction, the same term written after it is exogenous and receives none. Everything lm() accepts is allowed: transformations (log(p)), interactions (p*m), polynomials (I(p^2)), factors (as.factor(g)), and -1 to drop the intercept. The exogenous part may be omitted entirely. Endogenous regressors must be numeric. Missing values are handled by na.omit.
data: A data.frame containing the variables referenced in formula.
cdf: CDF estimator used for the copula transformation of the endogenous regressors. One of "kde.silverman" (default; integral of a Gaussian kernel density with Silverman's bandwidth, the estimator Park and Gupta use in the original article), "kde.cv" (same, with a cross-validated bandwidth), "kde.plugin" (same, with a plug-in bandwidth), "ecdf.fixed" (empirical CDF with a replaced boundary), "ecdf.adj" (adjusted empirical CDF), "rank.n" (rescaled empirical CDF, Rank/n with a correction), or "rank.n1" (Rank/(n+1)).
ties: How tied values are ranked when computing the empirical CDF: "max" (default, the counting function) or "average" (midranks).
nboots: Number of bootstrap replicates used for the standard errors.
subset: An optional logical vector specifying a subset of observations to use, as in lm(). NULL (default) uses all rows.
contrasts: An optional list of contrasts passed to model.matrix(). NULL (default) uses the defaults in options().
parallel: FALSE (default), TRUE, "multicore" or "snow" to spread the bootstrap over ncores workers. No seed is set internally either way; call set.seed() beforehand for reproducible standard errors.
ncores: Number of workers used when parallel is not FALSE. NULL (default) uses one less than the number of detected cores.
verbose: Show a progress bar during the bootstrap. Defaults to interactive().

Value

An object of class "copreg", a list including the structural coefficients and their bootstrap standard errors, the matched call, the uncorrected OLS benchmark from the same bootstrap resamples (ols.coefficients, ols.std.error), the raw bootstrap draws (boot), identification diagnostics (diagnostics), and the standard error inflation relative to OLS (icon). Methods exist for coef(), vcov(), confint(), residuals(), fitted(), predict(), print() and summary().

References

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

Examples

set.seed(42)
n <- 300
w <- rnorm(n)
u <- rnorm(n)
p <- rchisq(n, df = 3) + 0.3 * u   # non-normal endogenous regressor
y <- 1 + 2 * p + 0.5 * w + u
dat <- data.frame(y = y, p = p, w = w)

fit <- CopRegPG(y ~ p | w, data = dat, nboots = 25)
summary(fit)
#> 
#> Copula endogeneity correction: PG (Park & Gupta 2012) 
#> 
#> Call:
#> CopRegPG(formula = y ~ p | w, data = dat, nboots = 25)
#> 
#> Residuals of the augmented regression (u = xi - C gamma):
#>     Min      1Q  Median      3Q     Max 
#> -2.6892 -0.6745 -0.0929  0.6247  3.0787 
#> 
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#>      Min       1Q   Median       3Q      Max 
#> -2.53982 -0.62128 -0.06731  0.61248  3.31175 
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.35084    0.43959   3.073  0.00212 ** 
#> p            1.86613    0.16466  11.333  < 2e-16 ***
#> w            0.53260    0.05763   9.242  < 2e-16 ***
#> p_cop        0.49819    0.33814   1.473  0.14066    
#> ---
#> 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.4512     0.1810   2.493   0.0127 *
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Fit, on 296 residual degrees of freedom:
#>                         augmented structural
#> Residual standard error 0.9668    1.0834    
#> R-squared               0.9645    0.9555    
#> Adjusted R-squared      0.9642    0.9550    
#>   sigma above is the standard error of the structural model, the one 
#>   entering xi* = xi / sigma.
#> Standard errors from 25 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 8.608 5.845e-21 0.001517
#> 
#> 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.02702    w   0.6414 0.0007299  0.6412
#>  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.9214  0.15
#>