Actividad 3
Ricardo Cortés
2023-09-29
library(pacman)
p_load(data.table, fixest, magrittr, sandwich, mvtnorm, MASS)
# La función p_load de pacman nos permite cargar múltiples paquetes a la vez. Aún así, los tenemos que tener instalados.
set.seed(123)1. Sesgo por variables omitidas
# Load data
library(wooldridge)
data("bwght")
bwght$cigs2 <- bwght$cigs*bwght$cigs
# adds a column with squares
# Run regressions
est_0 = feols(bwght ~ cigs, data = bwght)
est_1 = feols(bwght ~ cigs + motheduc, data = bwght)
est_2 = feols(bwght ~ cigs + motheduc + faminc, data = bwght)
est_3 = feols(bwght ~ cigs + cigs2 + motheduc + faminc, data = bwght)
#Nicer tables
## Standard errors (using fixest package)
etable(est_0, est_1, est_2, est_3, se = "standard")## est_0 est_1 est_2
## Dependent Var.: bwght bwght bwght
##
## Constant 119.8*** (0.5723) 115.4*** (3.107) 116.8*** (3.138)
## cigs -0.5138*** (0.0905) -0.4862*** (0.0926) -0.4633*** (0.0927)
## motheduc 0.3308 (0.2328) 0.0143 (0.2580)
## faminc 0.0915** (0.0325)
## cigs2
## _______________ ___________________ ___________________ ___________________
## S.E. type IID IID IID
## Observations 1,388 1,387 1,387
## R2 0.02273 0.02420 0.02977
## Adj. R2 0.02202 0.02279 0.02767
##
## est_3
## Dependent Var.: bwght
##
## Constant 117.3*** (3.150)
## cigs -0.7787*** (0.2099)
## motheduc -0.0074 (0.2581)
## faminc 0.0898** (0.0324)
## cigs2 0.0123. (0.0073)
## _______________ ___________________
## S.E. type IID
## Observations 1,387
## R2 0.03174
## Adj. R2 0.02894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1a. Usando est_2 prediga el valor en onzas del peso de un niño cuya madre fumó 2 cigarros por día durante el embarazo, con todo lo demás en su promedio.
# Limpiar base de datos
# Hay NA?
sum(is.na(bwght))## [1] 197# Por tanto, hay 197 datos no registrados en el data set. Podemos eliminar los NA
# Checamos si tenemos NA en nuestras variables de interés
c(sum(is.na(bwght$motheduc)), sum(is.na(bwght$faminc)))## [1] 1 0Para no perder más observaciones de las necesarias, quitaremos sólo la fila de educación de la madre con NA.
Creando un nuevo data set para manipular:
# Quitamos el NA de motheduc y definimos lo demás como indica la pregunta:
new_data <- data.frame(cigs = 2, motheduc = mean(na.omit(bwght$motheduc)), faminc = mean(bwght$faminc))
predicted_weight <- predict(est_2, new_data)
cat("El peso estimado en onzas del bebé es:", predicted_weight, "onzas\n")## El peso estimado en onzas del bebé es: 118.7477 onzasb. ¿Cuál es el signo de la correlación entre educación de la madre y número de cigarros fumados al día durante el embarazo? ¿Es lo que esperaba? Explique. (Vea Tabla 3.2. en Wooldridge ed. 7)
Probaremos otro método de limpieza de NA:
bwght <- na.omit(bwght)
cor(bwght$motheduc, bwght$cigs)## [1] -0.2167217Como tenemos una correlación negativa, implica que entre mayor la educación de la madre, menor su consumo de cigarros.
c. ¿Cuál es el efecto de fumar un cigarro extra en el embarazo en el peso en onzas considerando una relación no lineal? ¿Por qué utilizaría una relación no lineal? Explique.
Usaría una relación no lineal porque el efecto que tiene fumar un cigarro no es el mismo al que tiene fumar una cajetilla diaria. Me imagino que el efecto puede incrementar mucho conforme la cantidad consumida es grande, pero después el efecto se vuelve marginal.
\[bwght = \beta_{0} + \beta_{1}cigs + \beta_{2} cigs^2 + \beta_{3}motheduc + \beta_{4}faminc + \epsilon \] \[\therefore \frac{\partial bwght}{\partial cigs} = \beta_{1} + 2 \beta_{2}cigs = -0-7787 + (2*0.0123)cigs\]
Esto implica que el cambio de fumar un cigarro adicional sobre el peso en onzas depende de un nivel especifico de cigarros. Por ejemplo, el impacto sobre el peso en onzas de fumar con dos, cinco y treinta cigarros es:
dos_cig <- -0.7787 + (2*0.0123)*2
cinco_cig <- -0.7787 + (2*0.0123)*5
treinta_cig <- -0.7787 + (2*0.0123)*30
c(dos_cig, cinco_cig, treinta_cig)## [1] -0.7295 -0.6557 -0.0407d. Interprete \(R^2\) y \(R^2\) adj en est_3. Formalmente ¿por qué \(R^2\) adj es menor?
r2 <- etable(est_3)[11,2]
adj_r2 <- etable(est_3)[12,2]
c(r2, adj_r2)## [1] "0.03174" "0.02894"\(R^2\) adj es menor porque est_3 tiene 4 regresores y estos castigan el valor en su fórmula.
2. Varianza de los estimadores
# Standard errors (using fixest package)
# Table 1:
# Under MLR5
etable(est_0,est_1,est_2,est_3, se = "standard")## est_0 est_1 est_2
## Dependent Var.: bwght bwght bwght
##
## Constant 119.8*** (0.5723) 115.4*** (3.107) 116.8*** (3.138)
## cigs -0.5138*** (0.0905) -0.4862*** (0.0926) -0.4633*** (0.0927)
## motheduc 0.3308 (0.2328) 0.0143 (0.2580)
## faminc 0.0915** (0.0325)
## cigs2
## _______________ ___________________ ___________________ ___________________
## S.E. type IID IID IID
## Observations 1,388 1,387 1,387
## R2 0.02273 0.02420 0.02977
## Adj. R2 0.02202 0.02279 0.02767
##
## est_3
## Dependent Var.: bwght
##
## Constant 117.3*** (3.150)
## cigs -0.7787*** (0.2099)
## motheduc -0.0074 (0.2581)
## faminc 0.0898** (0.0324)
## cigs2 0.0123. (0.0073)
## _______________ ___________________
## S.E. type IID
## Observations 1,387
## R2 0.03174
## Adj. R2 0.02894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Table 2:
# Fixes Heteroskedasticity (HC1 type)
etable(est_0,est_1,est_2,est_3, se = "hetero")## est_0 est_1 est_2
## Dependent Var.: bwght bwght bwght
##
## Constant 119.8*** (0.5745) 115.4*** (2.928) 116.8*** (2.945)
## cigs -0.5138*** (0.0877) -0.4862*** (0.0890) -0.4633*** (0.0893)
## motheduc 0.3308 (0.2202) 0.0143 (0.2408)
## faminc 0.0915** (0.0312)
## cigs2
## _______________ ___________________ ___________________ ___________________
## S.E. type Heteroskedast.-rob. Heteroskedast.-rob. Heteroskedast.-rob.
## Observations 1,388 1,387 1,387
## R2 0.02273 0.02420 0.02977
## Adj. R2 0.02202 0.02279 0.02767
##
## est_3
## Dependent Var.: bwght
##
## Constant 117.3*** (2.962)
## cigs -0.7787*** (0.2018)
## motheduc -0.0074 (0.2410)
## faminc 0.0898** (0.0312)
## cigs2 0.0123. (0.0070)
## _______________ ___________________
## S.E. type Heteroskedast.-rob.
## Observations 1,387
## R2 0.03174
## Adj. R2 0.02894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1a. Demuestre matemáticamente por qué el SE en la segunda tabla es distinto al SE la Tabla 1, al menos para el estimador asociado a cigs. (i.e derive la varianza del error homoscedastica y heteroscedástica).
\[\hat\sigma^2 = \frac{\sum^{n}_{i=1} \hat u^2_{i}}{(n-k-1)}\] Esta fórmula es un término del error estándar, el término de arriba es la suma cuadrada de los residuos. La varianza homocedastica es menor porque su varianza es constante, pues reduce los residuos al hacer la varianza constante y el error estandar también es constante. En el caso heterocedastico, el error estándar es más grande porque los datos presentan más dispersión.
b. ¿Qué supuesto asumimos como vulnerado en la estimación de Tabla 2?
Asumimos como vulnerado el MLR5 que dice que los datos son homocedasticos. Es decir, que la varianza de los errores debe ser constante a lo largo de todas las observaciones, lo que significa que la dispersión de los errores debe ser la misma en todos los niveles de las variables independientes.
c. Dado el siguiente código y considerando est_0 (en ejercicio 1) ¿aumenta la varianza del estimador cigs al incluirse mothereduc en est_1? Discuta lo anterior utilizando las fórmulas de la varianza de beta para un modelo simple y para uno múltiple.
summary(lm(bwght$cigs ~ bwght$motheduc)) ##
## Call:
## lm(formula = bwght$cigs ~ bwght$motheduc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.099 -2.308 -1.350 -0.392 37.692
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.05689 0.83521 9.647 < 2e-16 ***
## bwght$motheduc -0.47907 0.06258 -7.655 3.99e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.219 on 1189 degrees of freedom
## Multiple R-squared: 0.04697, Adjusted R-squared: 0.04617
## F-statistic: 58.6 on 1 and 1189 DF, p-value: 3.986e-14La varianza de un estimador es: \(\hat var(\hat\beta_{i}) = \frac{\sigma^2}{\sum(x_{i}- \bar{x})^2}\)
d. ¿Cuál es el VIF Tras agregar motheduc? Discuta, considerando el trade-off sesgo varianza si debemos incluir motheduc.
VIF <- 1/(1- 0.02420)
cat(VIF)## 1.0248e. ¿A qué concepto nos referimos en el inciso c?
El concepto al cuál nos referimos es el error estándar
3. The data BEAUTY (Wooldridge Library in R) contains data on wages and an index of people’s beauty.
a. Explore las variables con ?beauty y agregue una columna a los datos que incluya una aproximación de edad para cada i:
data("beauty")
beauty$age <- beauty$educ+beauty$exper+6b. Estime una regresión simple: de log(wage) con looks como explicativa con errores estándar robustos a heteroscedasticidad (HC1): pegue el resultado abajo e interprete
regresion_1 <- feols(lwage ~ looks, data = beauty)
etable(regresion_1, se = "white")## regresion_1
## Dependent Var.: lwage
##
## Constant 1.496*** (0.0818)
## looks 0.0510* (0.0252)
## _______________ _________________
## S.E. type Heteroskeda.-rob.
## Observations 1,260
## R2 0.00345
## Adj. R2 0.00265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Esto significa que cuando la belleza aumente en una unidad [1:5], el salario esperado debería de crecer un 5.10%
c. Agregue como controles a la regresión en b. la experiencia, la experiencia al cuadrado, la educación en años y la edad ¿Qué pasa con su regresión y por qué? ¿Qué supuesto se puede estar violando?
beauty$exper2 <- (beauty$exper)^2
regresion_2 <- feols(lwage ~ looks + exper + exper2 + educ + age, data = beauty)## The variable 'age' has been removed because of collinearity (see $collin.var).etable(regresion_2, se = "white")## regresion_2
## Dependent Var.: lwage
##
## Constant 0.0463 (0.1084)
## looks 0.0621** (0.0222)
## exper 0.0487*** (0.0046)
## exper2 -0.0007*** (0.0001)
## educ 0.0684*** (0.0059)
## _______________ ___________________
## S.E. type Heteroskedast.-rob.
## Observations 1,260
## R2 0.23225
## Adj. R2 0.22981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como podemos observar, el coeficiente de edad (age) no se incluye, lo cual tiene sentido dado que (como diseñamos anterioremente) la edad está compuesta (combineación lineal) por la educación y la experiencia, lo cual implica que hay colinearidad entre los regresores de educacion y experiencia. R nos indica que la borra.
d. Solucione el problema en c. y vuelva a estimar su regresión. Pegue el resultado abajo e interprete
El software elimina automáticament la varible age por la colinearidad, de tal modo que si los ‘looks’ aumentan en una unidad, el salario aumentará un 6.21%
4. Which of the following can cause OLS estimators to be biased? Why?
**a. Heteroskedasticity.* b. Omitting an important variable. *c. A sample correlation coefficient of .95 between two independent variables both included in the model.**
La heteroscedasticidad se refiere a la situación en la que la varianza del término de error (ε) en un modelo de regresión no es constante en diferentes valores de las variables independientes. En presencia de heteroscedasticidad, los estimadores OLS siguen siendo insesgados pero ya no son eficiente.
Omitir una variable importante de un modelo de regresión puede efectivamente causar que los estimadores OLS estén sesgados. Esto se conoce como sesgo por variable omitida. Cuando se excluye una variable importante del modelo, su efecto en la variable dependiente se atribuye al término de error. Como resultado, los estimadores OLS de las variables incluidas pueden estar sesgados porque no tienen en cuenta la influencia de la variable omitida.
Una alta correlación entre variables independientes (multicolinealidad) puede generar problemas en la estimación OLS, pero no necesariamente resulta en sesgo en los estimadores. En cambio, puede hacer que las estimaciones sean menos precisas (aumentar sus errores estándar), lo que dificulta distinguir los efectos individuales de las variables correlacionadas. Sin embargo, siempre y cuando la correlación no cause multicolinealidad perfecta (donde una variable puede expresarse como una combinación lineal de las otras), los estimadores OLS aún pueden ser insesgados.
Por tanto, la respuesta es b.
6. Confirm the partialling out interpretation of the OLS estimates by explicitly doing the partialling out for Example 3.2 (see Wooldridge 7th Edition Section 3.2f).
This first requires regressing educ on exper and tenure and saving the residuals, r1. Then, regress log(wage) on r1. Compare the coefficient on r1 with the coefficient on educ in the regression of log(wage) on educ, exper, and tenure.
a. Pegue abajo el resultado de las dos regresiones de interés y explique intuitivamente por qué se observa que el coeficiente de r1 es igual al de educ cuando controla por exper y tenure.
data("wage1")
first_regression <- lm(educ ~ exper + tenure, data = wage1)
r1 <- first_regression$residuals
second_regression <- lm(lwage ~ r1, data = wage1)
third_regression <- lm(lwage ~ educ + exper + tenure, data = wage1)
#Usaremos stargazer para una table más bonita:
#Si queremos que en Markdown se imprima dentro de nuestro mismo documento:
#{r, results = 'asis'}
# stargazer(second_regression, third_regression, type = "html")La intuición detrás de esto es que cuando se hace una regresión de educ sobre exper y tenure, esencialmente se está “particiendo” la parte de educ que puede explicarse mediante exper y tenure. Los residuos r1 representan la parte inexplicada de educ después de tomar en cuenta exper y tenure.
Entonces, cuando haces una regresión de log(wage) en r1, estás observando la relación entre la parte inexplicada de educ y log(wage), controlando por exper y tenure Si el coeficiente de r1 es similar al coeficiente de educ en el modelo completo [log(wage) con educ, exper y tenure], sugiere que el efecto de educ sobre log(wage) es capturado en gran medida por la parte de educ que no tiene relación con exper y tenure.
Este es un enfoque común para comprender la contribución única de una regresor después de tener en cuenta otras variables relevantes en el modelo.
7. Follow the next script step by step and answer the questions below.
# set seed
set.seed(123)
#Example before computation:
#we define a matrix X with two vectors (regressors x1 and x2) from random normal values with mean [5,10] and a matrix "sigma" of variances [3,6], and covariance 2
mu <- c(5,10)
sigma <- matrix(c(3,2,2,6), 2, 2) # Se define una matriz sigma que representa la matriz de covarianza entre las dos variables aleatorias. Los valores en la matriz determinan las varianzas y la covarianza entre las variables.
X <- rmvnorm(1000, mean = mu, sigma = sigma)
# Se genera una matriz X de dimensiones 1000x2 utilizando la función rmvnorm, que genera números aleatorios multivariados a partir de una distribución normal multivariada.
head(X)## [,1] [,2]
## [1,] 3.955929 9.171603
## [2,] 7.623002 10.937024
## [3,] 6.059552 14.178891
## [4,] 5.142173 7.191611
## [5,] 3.639927 8.592305
## [6,] 7.209863 11.466354summary(X)## V1 V2
## Min. : 0.06447 Min. : 3.454
## 1st Qu.: 3.92522 1st Qu.: 8.403
## Median : 5.16410 Median :10.059
## Mean : 5.07445 Mean :10.066
## 3rd Qu.: 6.11615 3rd Qu.:11.795
## Max. :11.04032 Max. :19.072var(X)## [,1] [,2]
## [1,] 2.768283 1.727895
## [2,] 1.727895 6.114899#plot vectors 1 and 2 in our matrix X (with covariance 2!)
plot(X[,1], X[,2])
##end of example.
# Set number of observations
n <- 500
#Another example of what is going to happen... Here we create a matrix X
#with two vectors (regressors x1 and x2) with median 50 and 100 each, and cov(X_1,X_2) = 0.25
X <-rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
# Aquí, se está creando una nueva matriz X. Esta matriz tiene dos columnas que representan dos variables aleatorias (regresores), x1 y x2. Se están generando datos aleatorios para estas dos variables de manera que la media de x1 sea 50, la media de x2 sea 100 y la covarianza entre x1 y x2 sea 0.25. Se están generando n observaciones.
summary (X)## V1 V2
## Min. :41.58 Min. : 92.43
## 1st Qu.:47.90 1st Qu.: 97.96
## Median :50.00 Median : 99.71
## Mean :49.98 Mean : 99.88
## 3rd Qu.:51.98 3rd Qu.:102.04
## Max. :59.31 Max. :110.06plot(X[,1], X[,2]) 
#Now the loop to do this 10000 times, fasten your seat belts:
# Empty vectors of coefficients to be filled
coefs1 <- cbind("beta_hat_1" = numeric(10000), "beta_hat_2" = numeric(10000))
coefs2 <- coefs1
# Se están creando dos matrices coefs1 y coefs2 que se utilizarán para almacenar los coeficientes de regresión estimados en un bucle posterior. Estas matrices tienen 10,000 filas y dos columnas llamadas "beta_hat_1" y "beta_hat_2".
###From here###
# Loop sampling and estimation (be patient it takes a while)
for (i in 1:10000) { # Se inicia un bucle for que se ejecutará 10,000 veces.
# for median 50 and 100 and cov(X_1,X_2) = 0.25
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
# En cada iteración del bucle, se generan nuevos datos aleatorios para X con las mismas características que se describieron anteriormente: medias de 50 y 100 para x1 y x2, y covarianza de 0.25.
u <- rnorm(n, sd = 5)
# Se genera un vector u de errores aleatorios con distribución normal con una desviación estándar de 5. Luego, se define una variable Y como una función lineal de x1 y x2 más los errores u.
#We define PRF
Y <- 5 + 2.5 * X[, 1] + 3.5 * X[, 2] + u
#we compute SRF 10000 times using our n values in matrix X
coefs1[i, ] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
# Se realiza una regresión lineal de Y sobre x1 y x2 utilizando la función lm. Se almacenan los coeficientes estimados (sin incluir la intersección) en la matriz coefs1 en la fila i.
# for cov(X_1,X_2) = 0.85
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))
Y <- 5 + 2.5 * X[, 1] + 3.5 * X[, 2] + u
coefs2[i, ] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
# Se repiten los mismos pasos, pero esta vez se generan datos con una covarianza diferente (0.85) y se almacenan los coeficientes en la matriz coefs2.
}
#Histogram of coefficients x1 and x2 for the two 10000s loops of regressions.
hist(coefs1)
hist(coefs2)
# Se crean histogramas para visualizar la distribución de los coeficientes estimados para x1 y x2 en ambos conjuntos de datos generados.
# Obtain variance estimates
diag(var(coefs1))## beta_hat_1 beta_hat_2
## 0.005266405 0.005323088diag(var(coefs2))## beta_hat_1 beta_hat_2
## 0.01817545 0.01825663 # Se calculan estimaciones de la varianza para los coeficientes estimados en ambas matrices coefs1 y coefs2. Las diagonales de estas matrices representan las varianzas de los coeficientes.a. A partir de la marca ###From here###. Explique cada línea pertinente del código Pongo la explicación de cada linea en el codigo arriba con #comentarios
b. Escriba la matriz de varianza – covarianza de los estimadores 1 y 2.
sigma <- matrix(c(3,2,2,6), 2, 2) # Se define una matriz sigma que representa la matriz de covarianza entre las dos variables aleatorias. Los valores en la matriz determinan las varianzas y la covarianza entre las variables.c. ¿Por qué los histogramas de coefs1 y coefs2 son
distintos? Los histogramas de coefs1 y
coefs2 son distintos porque los datos utilizados para
estimar los coeficientes de regresión en ambos casos son generados con
diferentes matrices de covarianza. Vamos a analizarlo en detalle:
Matriz de Covarianza Distinta: En el código, se generan dos conjuntos de datos diferentes,
X, con dos matrices de covarianza distintas. La primera matriz de covarianza tiene una covarianza de 0.25 entrex1yx2, mientras que la segunda matriz de covarianza tiene una covarianza de 0.85 entrex1yx2. Estas diferencias en las covarianzas implican que las variablesx1yx2en los dos conjuntos de datos tienen relaciones diferentes entre sí.Efecto en la Regresión: Cuando se ajusta un modelo de regresión lineal a estos datos, la estructura de correlación entre las variables predictoras (
x1yx2) y la variable de respuesta (Y) afectará los coeficientes estimados. En otras palabras, la relación entrex1,x2yYserá diferente en los dos conjuntos de datos debido a las diferencias en las covarianzas.Estimación de Coeficientes: Los coeficientes de regresión estimados en el bucle
forse calculan utilizando diferentes datos en cada iteración. En cada iteración, se generan nuevos valores deXy, por lo tanto, se obtienen diferentes estimaciones de los coeficientes de regresión.Histogramas Diferentes: Los histogramas que se crean al final del código representan la distribución de los coeficientes de regresión estimados a lo largo de las 10,000 iteraciones para cada conjunto de datos. Dado que los datos son diferentes debido a las matrices de covarianza diferentes, las distribuciones de los coeficientes estimados también serán diferentes. Los valores y la dispersión de los coeficientes variarán en función de la estructura de correlación en los datos.
d. ¿Por qué las varianzas de coefs1 y coefs2 son
distintas? Las varianzas de coefs1 y
coefs2 son distintas debido a las diferencias en la
estructura de los datos utilizados para estimar los coeficientes de
regresión en cada caso. Aquí está la explicación detallada:
Diferentes Matrices de Covarianza: En el código, se generan dos conjuntos de datos diferentes,
X, con dos matrices de covarianza distintas. La primera matriz de covarianza tiene una covarianza de 0.25 entrex1yx2, mientras que la segunda matriz de covarianza tiene una covarianza de 0.85 entrex1yx2. Estas diferencias en las covarianzas indican que las variablesx1yx2en los dos conjuntos de datos tienen estructuras de correlación diferentes entre sí.Efecto en la Regresión: Cuando se ajusta un modelo de regresión lineal a estos datos, la estructura de correlación entre las variables predictoras (
x1yx2) y la variable de respuesta (Y) afectará los coeficientes estimados. Las diferencias en las covarianzas pueden influir en la variabilidad de los coeficientes estimados, especialmente si existe una correlación más fuerte o débil entrex1yx2.Estimación de Coeficientes: Los coeficientes de regresión estimados en el bucle
forse calculan utilizando diferentes datos en cada iteración. En cada iteración, se generan nuevos valores deXy, por lo tanto, se obtienen diferentes estimaciones de los coeficientes de regresión. La estructura de correlación diferente en los datos tiene un impacto directo en la variabilidad de las estimaciones de coeficientes.Varianza de Coeficientes: La varianza de los coeficientes estimados se calcula al final del código para ambos conjuntos de datos. Debido a las diferencias en la estructura de correlación y las covarianzas en los datos, las estimaciones de coeficientes tendrán diferentes varianzas. La varianza refleja la dispersión de los valores estimados alrededor de su valor esperado, y esta dispersión será diferente en función de cómo están relacionadas
x1yx2en cada conjunto de datos.
e. ¿qué parte del código genera estas diferencias? Las diferencias en las varianzas de coefs1 y coefs2 se generan principalmente debido a las diferencias en la estructura de correlación de los datos utilizados en las dos partes del bucle for que estiman los coeficientes de regresión. Concretamente, las diferencias se producen en la generación de los datos aleatorios (X) y, en última instancia, en cómo se relacionan las variables predictoras (x1 y x2) en cada conjunto de datos. Aquí está la parte del código que genera estas diferencias:
# for median 50 and 100 and cov(X_1,X_2) = 0.25
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
#...
# for cov(X_1,X_2) = 0.85
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))En las líneas anteriores, se generan los datos aleatorios (X) para cada conjunto de datos dentro del bucle for. Las diferencias en la estructura de correlación y covarianza entre x1 y x2 se especifican en las matrices sigma que se pasan como argumento a la función rmvnorm. En particular:
En la primera parte del bucle (for median 50 and 100 and cov(X_1,X_2) = 0.25), se genera X con una covarianza de 0.25 entre x1 y x2, lo que indica una relación más débil entre las dos variables.
En la segunda parte del bucle (for cov(X_1,X_2) = 0.85), se genera X con una covarianza de 0.85 entre x1 y x2, lo que indica una relación más fuerte entre las dos variables.
Estas diferencias en la estructura de correlación y covarianza entre x1 y x2 en los datos aleatorios generados son la causa principal de las diferencias en las varianzas de los coeficientes estimados (coefs1 y coefs2) en el resultado final del código.
f. ¿Qué concepto/supuesto central estamos simulando en este ejercicio? El ejercicio que se presenta en el código tiene como objetivo simular el concepto central de la econometría y la estadística conocido como “Homocedasticidad” o “Homogeneidad de la varianza”. Este concepto es fundamental en el contexto de la regresión lineal y tiene importantes implicaciones para la validez de las inferencias estadísticas.
La homocedasticidad se refiere a la suposición de que la varianza de los errores en un modelo de regresión es constante en todos los niveles de las variables predictoras. En otras palabras, se asume que la dispersión de los errores no cambia a medida que las variables independientes (en este caso, x1 y x2) varían. Esto es un supuesto importante para que las estimaciones de los coeficientes de regresión sean eficientes y no sesgadas.
En el ejercicio, se simulan dos escenarios diferentes para ilustrar cómo la covarianza entre las variables predictoras x1 y x2 (representada en las matrices de covarianza sigma) puede afectar la homocedasticidad y, por lo tanto, la varianza de los coeficientes estimados en un modelo de regresión lineal.
En el primer escenario, se genera un conjunto de datos con una covarianza baja (0.25) entre x1 y x2. Esto representa una relación más débil entre las variables predictoras y, en general, un escenario de mayor homocedasticidad.
En el segundo escenario, se genera otro conjunto de datos con una covarianza más alta (0.85) entre x1 y x2. Esto representa una relación más fuerte entre las variables predictoras y, en general, un escenario de menor homocedasticidad.
Al comparar los resultados de estos dos escenarios, el ejercicio ilustra cómo la estructura de correlación y covarianza entre las variables predictoras puede influir en la variabilidad de los coeficientes estimados y, por lo tanto, en la validez de las inferencias realizadas en un modelo de regresión. En última instancia, esto destaca la importancia de considerar la homocedasticidad al realizar análisis de regresión y evaluar la validez de los resultados.
g. Ahora cambie n a 1000, ¿qué sucede con los histogramas y las varianzas? Explique intuitivamente, pero con el uso de alguna fórmula ¿por qué se da este cambio?
# set seed
set.seed(123)
# Set number of observations (changed to 1000)
n <- 1000
#Another example of what is going to happen... Here we create a matrix X
#with two vectors (regressors x1 and x2) with median 50 and 100 each, and cov(X_1,X_2) = 0.25
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
summary (X)## V1 V2
## Min. :41.18 Min. : 91.54
## 1st Qu.:48.04 1st Qu.: 97.92
## Median :50.20 Median :100.06
## Mean :50.13 Mean :100.08
## 3rd Qu.:52.02 3rd Qu.:102.37
## Max. :59.90 Max. :111.39plot(X[,1], X[,2]) 
#Now the loop to do this 10000 times, fasten your seat belts:
# Empty vectors of coefficients to be filled
coefs1 <- cbind("beta_hat_1" = numeric(10000), "beta_hat_2" = numeric(10000))
coefs2 <- coefs1
###From here###
# Loop sampling and estimation (be patient it takes a while)
for (i in 1:10000) {
# for median 50 and 100 and cov(X_1,X_2) = 0.25
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
u <- rnorm(n, sd = 5)
#We define PRF
Y <- 5 + 2.5 * X[, 1] + 3.5 * X[, 2] + u
#we compute SRF 10000 times using our n values in matrix X
coefs1[i, ] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
# for cov(X_1,X_2) = 0.85
X <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))
Y <- 5 + 2.5 * X[, 1] + 3.5 * X[, 2] + u
coefs2[i, ] <- lm(Y ~ X[, 1] + X[, 2])$coefficients[-1]
}
#Histogram of coefficients x1 and x2 for the two 10000s loops of regressions.
hist(coefs1)
hist(coefs2)
# Obtain variance estimates
diag(var(coefs1))## beta_hat_1 beta_hat_2
## 0.002618979 0.002669144diag(var(coefs2))## beta_hat_1 beta_hat_2
## 0.008932250 0.008922246Cuando cambias el valor de n a 1000 (anteriormente era
500), esto significa que estás generando un conjunto de datos más grande
para cada escenario. Esto tiene varias implicaciones en cuanto a lo que
sucede con los histogramas y las varianzas de los coeficientes
estimados:
Histogramas:
Los histogramas se volverán más suaves y tenderán a seguir una distribución más cercana a la distribución normal, especialmente cuando
nes lo suficientemente grande. Esto se debe al Teorema del Límite Central, que establece que cuando se toman muestras suficientemente grandes de una población, la distribución de la media de esas muestras se aproximará a una distribución normal, independientemente de la forma de la distribución poblacional original.En el caso de \(n = 1000\), los histogramas de \(coefs1\) y \(coefs2\) serán más simétricos y se asemejarán más a campanas de Gauss (distribuciones normales) en comparación con el caso anterior con \(n = 500\).
Varianzas:
Con un tamaño de muestra mayor (\(n = 1000\)), las estimaciones de los coeficientes de regresión tendrán varianzas más pequeñas. Esto se debe a la fórmula para la varianza de la media muestral, que indica que a medida que aumenta el tamaño de la muestra, la varianza de la media disminuye. La fórmula es:
\[ \text{Var}(\bar{X}) = \frac{\sigma^2}{n} \]
Donde:
- \(\text{Var}(\bar{X})\) es la varianza de la media muestral.
- \(\sigma^2\) es la varianza de la población.
- \(n\) es el tamaño de la muestra.
Por lo tanto, cuando \(n\) se duplica de 500 a 1000, la varianza de las estimaciones de los coeficientes (como las de \(coefs1\) y \(coefs2\)) se reduce a la mitad aproximadamente, lo que significa que las estimaciones son más precisas y tienen menos variabilidad.
En resumen, aumentar el tamaño de la muestra (\(n\)) a 1000 tiene el efecto de suavizar los histogramas y reducir las varianzas de las estimaciones de los coeficientes. Esto se debe a la ley de los grandes números y al Teorema del Límite Central, que indican que a medida que aumenta el tamaño de la muestra, las estimaciones se vuelven más precisas y se aproximan más a una distribución normal.
h. Aún bajo el supuesto MRL4 ¿Se puede obtener un estimador inconsistente en una muestra aleatoria de individuos? De un ejemplo con el uso de los histogramas.
El supuesto MRL4 (Media de los errores igual a cero) establece que la esperanza matemática de los errores de un modelo de regresión lineal simple es igual a cero. Este supuesto es fundamental en la regresión lineal, y si no se cumple, los estimadores de los coeficientes de regresión pueden ser sesgados. Sin embargo, incluso bajo este supuesto, es posible obtener estimadores inconsistentes en una muestra aleatoria de individuos si se violan otros supuestos clave de la regresión lineal.
Uno de los supuestos críticos en la regresión lineal es el supuesto de homocedasticidad, que establece que la varianza de los errores es constante en todos los niveles de las variables predictoras. Si este supuesto se viola, es decir, si la varianza de los errores no es constante en todos los niveles de las variables predictoras, entonces los estimadores de los coeficientes de regresión pueden ser inconsistentes. Esto significa que, a medida que aumenta el tamaño de la muestra, los estimadores no convergen al valor verdadero del parámetro, lo que lleva a estimaciones sesgadas y poco confiables.
Para ilustrar esto con un ejemplo utilizando histogramas, podemos simular datos en los que se viole el supuesto de homocedasticidad. A continuación, se muestra un ejemplo de código en R que genera datos con heterocedasticidad y luego ajusta un modelo de regresión lineal:
# Establecer una semilla para reproducibilidad
set.seed(123)
# Crear datos simulados con heterocedasticidad
n <- 100
x <- seq(1, n)
error <- rnorm(n, mean = 0, sd = x) # Varianza creciente con x
y <- 2 * x + 5 + error
# Ajustar un modelo de regresión lineal
model <- lm(y ~ x)
# Ver los coeficientes estimados
summary(model)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -179.394 -27.067 1.914 20.882 191.507
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.5531 10.7540 -0.33 0.742
## x 2.3012 0.1849 12.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.37 on 98 degrees of freedom
## Multiple R-squared: 0.6125, Adjusted R-squared: 0.6086
## F-statistic: 154.9 on 1 and 98 DF, p-value: < 2.2e-16# Crear histograma de los residuos
hist(model$residuals, main = "Histograma de Residuos", xlab = "Residuos")
En este ejemplo, se generan datos en los que la varianza de los errores aumenta a medida que la variable \(x\) aumenta. Esto viola el supuesto de homocedasticidad. Cuando ajustamos un modelo de regresión lineal a estos datos y examinamos el histograma de los residuos, podemos ver que los residuos no tienen una distribución normal y la dispersión de los residuos cambia con x. Esto indica que el supuesto de homocedasticidad se ha violado y, como resultado, los estimadores de los coeficientes de regresión pueden ser inconsistentes, incluso si el supuesto MRL4 se cumple.