Buscar una ecuación OLS si un regresor está correlacionado con el error

Question

¿Cómo puedo expresar un estimador OLS si conozco la correlación, es decir, sé que $E(x_i u_i)=\rho$ (No estoy buscando IV o 2SLS).

Voy a explicar mi problema con un ejemplo:

En un problema sencillo $\ \ \ y_i = \beta_0 + \beta_1 x_i + u_i \ \ \ $ i asumir

$E(u_i)=0$

$E(x_i u_i) = \rho$ $ \ \ \ i=1,...,n$

Con la MM, puedo derivar las siguientes ecuaciones:

$\frac{1}{n}\sum_{i=1}^n(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = 0 \ \ \ \ \ \ \leftrightarrow \ \ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}$

$\frac{1}{n}\sum_{i=1}^n x_i(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i) = \rho \ \ \ \leftrightarrow \ \ \hat{\beta}_1 = \frac{\sum_{i=1}^n x_i(y_i - \bar{y})}{\sum_{i=1}^n x_i(x_i - \bar{x})} - \frac{\rho}{\sum_{i=1}^n x_i(x_i - \bar{x})}$

En términos del estimador ols y algún factor, ahora puedo escribir lo siguiente:

$\hat{\beta}_1 = {\hat{\beta}_1}_{ols} - \frac{\rho}{m_{xx}}$

$\hat{\beta}_0 = {\hat{\beta}_0}_{ols} + \bar{y}\frac{\rho}{m_{xx}}$

Ahora podemos ver por qué factor particular el $\boldsymbol{\beta}_{ols}$ el estimador está sesgado. Para este caso tan simple, pude derivar los estimadores muy fácilmente, pero para el caso más común con k regresores, no puedo encontrar un buen enfoque para expresar el nuevo estimador. Por ejemplo, pensé que podría utilizar

$\boldsymbol{X^T(y-X\beta) = \rho} \ \ $ en lugar de $ \ \ \boldsymbol{X^T(y-X\beta) = 0}$

con $ \ \boldsymbol{\rho} \ $ a $\ \ k\times 1$ -vector que contiene la información sobre la correlación entre los regresores y el término de error.

Pero la solución

$\boldsymbol{\hat{\beta} = (X^TX)^{-1}X^Ty - (X^TX)^{-1}\rho}$

no hace del todo el trabajo.

Answer 1

Aunque esta situación no se da en la práctica, está relacionada con el llamado enfoque de la función de control para tratar la endogeneidad.

Permítame reescribir su (simple) modelo $$ Y_i = \beta_0 + \beta_1X_i + U_i $$ junto con sus suposiciones $\mathbb{E}(U_i)=0$ y $\mathbb{E}(U_iX_i)=\rho$ .

Entonces $$ \mathbb{E}(X_i(U_i-\frac{\rho}{X_i}))=0 $$ de modo que si reescribo mi modelo $$ Y_i = \beta_0 + \beta_1X_i + \frac{\rho}{X_i} + \underbrace{U_i-\frac{\rho}{X_i}}_{\equiv V_i} $$ y estimar este modelo por OLS, restringiendo el coeficiente de la $\frac{1}{X_i}$ plazo para ser $\rho$ Debería obtener estimaciones consistentes de $\beta_1$ .

---

Así, consideremos las siguientes simulaciones de Stata

clear*
program simcont, rclass
    drop _all
    set obs 1000
    g x1 = rnormal()
    g x2 = rnormal()
    g u = x1 + rnormal()
    g x = x1 + x2

    g y = 2 + 3*x + u  // ols 
    reg y x
    mat mA = e(b)
    return scalar ols = el(mA, 1, colnumb(mA, "x"))

    g cont = 1/x
    constraint define 1 cont = 1  // true correlation between error and regressor
    cnsreg y x cont, constraints(1)  // constrained regression
    mat mA = e(b)
    return scalar cont = el(mA, 1, colnumb(mA, "x"))
end

simulate olsCoeff = r(ols) controlFuncCoeff=r(cont), reps(100): simcont 

kdensity olsCoeff, xline(3, lcolor(green)) addplot(kdensity controlFuncCoeff) ///
    legend(label(1 "KDE of OLS coeff. estimates") label(2 "KDE of control function coeff. estimates")) ///
    xtitle("estimates") ytitle("density") title("Comparison of OLS and control function approaches") 

que produce la siguiente imagen (basada en 100 réplicas)

Sin embargo, yo me lo pensaría muy bien y experimentaría con más regresores y, en general, con más configuraciones de datos antes de poner en práctica esta estrategia de estimación con datos reales.

---

Preguntas de seguimiento:

El OP ha pedido algunas aclaraciones en los comentarios para las que proporciono una respuesta actualizada.

Déjeme reescribir su modelo $$ Y_i = \beta_0 + \beta_1 Z_{1i} + \beta_2 Z_{2i} + \beta_3 X_i + U_i $$ donde $Z_{1i}$ y $Z_{2i}$ son exógenos, y $X_{i}$ es endógena, es decir $\mathbb{E}(X_iU_i) = \rho$ . Además, quiere que la variable $Z_{2i}$ que se construirá como $$ Z_{2i} = \mathbf{1}_{[i\text{ is odd.}]} $$ Está simulando las estimaciones OLS del coeficiente de $X_i$ Es decir $\beta_3$ . Aquí hay un pequeño Stata script para hacerlo.

clear*
program simcont, rclass
    syntax [, errorVariance(real 1.0)]
    drop _all
    set obs 1000
    scalar beta0 = 5
    scalar beta1 = 1
    scalar beta2 = 2
    scalar beta3 = 3
    scalar rho = 0.1

    g z1 = rnormal()
    g z2 = mod(_n, 2)
    g u = sqrt(`errorVariance')*rnormal() 
    g x = rho*u/`errorVariance' + rnormal()

    g y = beta0 + beta1*z1 + beta2*z2 + beta3*x + u
    reg y z1 z2 x
    mat mA = e(b)
    return scalar ols = el(mA, 1, colnumb(mA, "x"))

    // return the results of the heteroskedasticity test
    estat hettest, rhs iid
    return scalar hettestPValues = r(p)
end

// simulate with error variance = 1
simulate olsCoeff = r(ols) hettestPValues = r(hettestPValues), reps(100): simcont 
su olsCoeff hettestPValues // p-values have the correct mean; no heteroskedasticity

cap mat drop biasBeta
forvalues errorVariance = 2(1)6 {
    simulate olsCoeff = r(ols), reps(100): simcont, errorVariance(`errorVariance') 
    qui su olsCoeff
    mat biasBeta = (nullmat(biasBeta), r(mean) - 3)
    local colNames "`colNames' errVar:`errorVariance' "
}
mat colnames biasBeta = `colNames'
mat list biasBeta

Ahora bien, hay que tener en cuenta que aquí hay al menos dos discrepancias.

• Lo primero es que no hay heteroscedasticidad en el modelo tal y como lo has escrito. La media de los valores p de una prueba LM de homoscedasticidad de las simulaciones indica que se están dibujando bajo la nula.

. su olsCoeff hettestPValues // p-values have the correct mean; no heteroskedasticity

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
    olsCoeff |       100    3.094182    .0326048   3.007202   3.179945
hettestPVa~s |       100    .4762952    .2885692   .0039824   .9844939

• Lo siguiente que hay que tener en cuenta es que el sesgo asintótico es el mismo para todos los valores de la varianza del error, siempre que el grado de endogeneidad sea el mismo. Aquí están los resultados de mi simulación

. mat list biasBeta

biasBeta[1,5]
       errVar:    errVar:    errVar:    errVar:    errVar:
            2          3          4          5          6
r1  .09913641  .10274499   .0912232  .11309942  .09122764

Si obtienes resultados diferentes a estos, entonces debes mostrarnos tu código y podremos comparar los dos.

Actualización del seguimiento

Déjeme reescribir su modelo una vez más $$ Y_i = \beta_0 + \beta_1Z_{1i} + \beta_2 Z_{2i} + \beta_3 X_i + U_{1i} $$ donde como antes, $Z_{1i}$ y $Z_{2i}$ son exógenos y $X_i$ es endógena. Se supone que el proceso de generación de datos $$ \begin{align} Z_{2i} &= \mathbf{1}_{[i\text{ is odd}]}\\ X_i &= \alpha_0 + \alpha_1Z_{1i} + U_{2i} \end{align} $$ Se introduce la endogeneidad asumiendo que $U_{1i}$ y $U_{2i}$ están correlacionados. $$ \mathbb{C}(U_{1i}, U_{2i}) = \rho $$

Además, asume que los errores en la ecuación de forma reducida son heterocedásticos: $$ \begin{align} \mathbb{V}(U_{2i}\mid Z_{2i} = 0) &= 1 \\ \mathbb{V}(U_{2i}\mid Z_{2i} = 1) &= \tfrac{1}{q} \\ \end{align} $$

Aquí está el código de Stata para simular este DGP.

clear*
program simcont, rclass
    syntax [, Q(real 1.0)]

    drop _all
    set obs 1000
    scalar beta0 = 5
    scalar beta1 = 1
    scalar beta2 = 2
    scalar beta3 = 3

    scalar alpha0 = 1
    scalar alpha1 = 4

    scalar rho = 0.1

    g z1 = rnormal()
    g z2 = mod(_n, 2)
    scalar a11 = 1
    scalar a12 = rho*sqrt(1)*sqrt(1)
    scalar a13 = rho*sqrt(1)*sqrt(1/`q')
    scalar a21 = rho*sqrt(1)*sqrt(1)
    scalar a22 = 1
    scalar a23 = 0
    scalar a31 = rho*sqrt(1)*sqrt(1/`q')
    scalar a32 = 0
    scalar a33 = 1/`q'

    mat corrMatrix= (a11, a12, a13 \ a21, a22, a23 \a31, a32, a33)
    drawnorm u1 u3 u4, cov(corrMatrix)
    g u2 = cond(z2, u3, u4)

    g x = alpha0 + alpha1*z2 + u2
    g y = beta0 + beta1*z1 + beta2*z2 + beta3*x + u1
    reg y z1 z2 x
    mat mA = e(b)
    return scalar ols = el(mA, 1, colnumb(mA, "x"))

    // return the results of the heteroskedasticity test
    qui reg x z2
    estat hettest, rhs iid
    return scalar hettestPValues = r(p)
end

// simulate to check for heteroskedasticity
simulate olsCoeff = r(ols) hettestPValues = r(hettestPValues), reps(100): simcont, q(5) 
su hettestPValues // strong evidence of heteroskedasticity in the reduced form

cap mat drop biasBeta
forvalues q = 2(1)5 {
    simulate olsCoeff = r(ols), reps(1000): simcont, q(`q') 
    qui su olsCoeff
    mat biasBeta = (nullmat(biasBeta), r(mean) - 3)
    local colNames "`colNames' q:`q' "
}
mat colnames biasBeta = `colNames'
mat list biasBeta

• Tenga en cuenta que la heteroscedasticidad está donde usted la puso en el modelo, y los resultados de la simulación son ahora capaces de encontrarla.

. su hettestPValues // strong evidence of heteroskedasticity in the reduced form

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
hettestPVa~s |       100    1.53e-26    8.12e-26   3.09e-38   6.22e-25

• Además, puedo confirmar que lo que afirmas es realmente cierto, que a medida que aumenta la heteroscedasticidad en la ecuación de forma reducida, el sesgo en el estimador OLS también aumenta.

biasBeta[1,4]
            q:         q:         q:         q:
            2          3          4          5
r1  .11242463  .11808594  .12122701  .12210455

Una explicación sencilla de esto es que como $q$ aumenta, el condicional (sobre $Z_{2i}=1$ ) y, por tanto, la varianza incondicional del error de la forma reducida $U_{2i}$ disminuye (compruébelo mediante una descomposición de la varianza), lo que significa que la varianza del regresor endógeno disminuye, lo que significa que el $(\mathbf{X}'\mathbf{X})^{-1}$ aumenta en magnitud, y el sesgo general aumenta (esta es una descripción muy aproximada, estoy seguro de que se puede formalizar).