16 Respuestas a los ejercicios seleccionados
16.1 Capítulo 2
16.1.1 Función Cobb-Douglas
Partiendo de (2.10), y sin perder generalidad encontremos la derivada con respecto \(X_2\).
\[\begin{equation*} \frac{\partial Y_i }{\partial X_{2i}}=\alpha _2\left ( \alpha_0 X_{1i}^{\alpha _1}X_{2i}^{\alpha _2 - 1}X_{3i}^{\alpha _3}{\varepsilon _i} \right ) \end{equation*}\] Ahora, recociendo que \[\begin{equation*} Y_i =\alpha_0 X_{1i}^{\alpha _1}X_{2i}^{\alpha _2}X_{3i}^{\alpha _3}{\varepsilon _i} \end{equation*}\] podemos reexpresar la derivada de la siguiente manera: \[\begin{equation*} \frac{\partial Y_i }{\partial X_{2i}}=\alpha _2 \frac{Y_i}{X_{2i}}. \end{equation*}\] Ahora multiplicando a ambos lados por \(X_{2i}\) y dividiendo a ambos lados por \(Y_i\) obtenemos: \[\begin{equation*} \frac{\frac{\partial Y_i}{Y_i} }{\frac{\partial X_{2i}}{X_{2i}}}=\alpha _2 . \end{equation*}\] Ahora multiplicando por 100 el numerador y denominador de la mano derecha obtendremos: \[\begin{equation*} \frac{\frac{\partial Y_i}{Y_i} \dot 100}{\frac{\partial X_{2i}}{X_{2i}} \dot 100}=\alpha _2 . \end{equation*}\] Es fácil reconocer que el cambio porcentual de \(Y_i\) (\(\Delta\%Y_i\)) corresponde a \(\frac{\partial Y_i}{Y_i} \dot 100\). Así, tenemos que: \[\begin{equation*} \frac{\Delta\%Y_i}{\Delta\%X_{2i}}=\alpha _2. \end{equation*}\]
Es decir, \(\alpha _2\) representa el cambio porcentual en \(Y_i\) dado un cambio del 1% en \(X_{2i}\). Esto representa la elasticidad.
16.1.2 Ejercicio práctico
El primer paso es cargar los datos.
## 'data.frame': 34 obs. of 7 variables:
## $ Año : num 1988 1989 1990 1991 1992 ...
## $ CE (T) : num 1450 1477 1504 1531 1558 ...
## $ CD (T) : num 193 188 183 177 172 166 160 148 150 147 ...
## $ I (T) : num 71135 76495 81855 87215 92575 ...
## $ Ldies (T): num 931 905 879 853 827 802 775 741 733 691 ...
## $ LEl (T) : num 712 699 687 675 663 650 623 621 621 606 ...
## $ V (T) : num 36127 61829 74066 183746 220029 ...Los datos han sido leídos correctamente, pero los nombres de las variables no parecen los correctos.
## 'data.frame': 34 obs. of 7 variables:
## $ Año : num 1988 1989 1990 1991 1992 ...
## $ CE : num 1450 1477 1504 1531 1558 ...
## $ CD : num 193 188 183 177 172 166 160 148 150 147 ...
## $ I : num 71135 76495 81855 87215 92575 ...
## $ Ldies: num 931 905 879 853 827 802 775 741 733 691 ...
## $ LEI : num 712 699 687 675 663 650 623 621 621 606 ...
## $ V : num 36127 61829 74066 183746 220029 ...De esta manera, podemos proceder a estimar el modelo.
##
## Call:
## lm(formula = I ~ ., data = dejercicio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1925.0 -150.4 -20.9 110.5 3588.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.576e+04 2.244e+04 4.267 0.000204 ***
## CE 8.576e+01 8.391e+00 10.220 5.94e-11 ***
## CD 4.091e+01 7.722e+01 0.530 0.600490
## Ldies -5.518e+01 2.047e+01 -2.695 0.011754 *
## LEI -1.484e+02 2.393e+01 -6.201 1.07e-06 ***
## V 8.886e-04 1.417e-03 0.627 0.535758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 969.6 on 28 degrees of freedom
## Multiple R-squared: 0.9997, Adjusted R-squared: 0.9997
## F-statistic: 2.003e+04 on 5 and 28 DF, p-value: < 2.2e-16Ahora podemos interpretar los coeficientes. En este caso tendremos que:
- Intercepto: Si todas las otras variables del modelo fuesen cero, entonces el sector tendría egresos de 9.57552^{4} millones de dólares. Noten que este es un ejemplo en el que el intercepto no tiene sentido.
- Por cada millón adicional Kilovatios/hora de consumo de electricidad el ingreso del sector aumentará en 85.76 millones de dólares.
- Como se discutirá en el siguiente capítulo, el consumo de diésel no tiene efecto sobre el ingreso del sector (no es significativo).
- Por cada locomotora diésel adicional el ingreso del sector caerá en 55.18 millones de dólares.
- Por cada locomotora eléctrica adicional el ingreso del sector caerá en 148.37 millones de dólares.
- Como se discutirá en el siguiente capítulo, el número de viajeros no tiene efecto sobre el ingreso del sector (no es significativo).
16.2 Capítulo 4
16.2.1 Ejercicio práctico (continuación)
El primer paso es cargar los datos y hacer las modificaciones que ya habíamos efectuado.
dejercicio <- read.csv("./DataFinalCap/regmult.csv", sep=";")
str(dejercicio)
names(dejercicio) <- c("Año", "CE", "CD", "I", "Ldies", "LEI", "V")De esta manera, podemos proceder a estimar los modelos. Los resultados se reportan en los siguientes dos cuadros.
R1 <- lm(I ~ CE + CD + Ldies + LEI + V, dejercicio )
R2 <- lm(I ~ Ldies + LEI, dejercicio )
R3 <- lm(I ~ CE + LEI, dejercicio )
R4 <- lm(I ~ CE + CD + V , dejercicio )
R5 <- lm(I ~ CE + CD + Ldies + LEI , dejercicio )
if(knitr::is_latex_output()){
stargazer(R1, R2, R3, align=FALSE, t.auto=TRUE, p.auto = TRUE,
title="Modelos 1 a 3 estimados por MCO", label ="tab:ejercomp1", header= FALSE, table.placement="H", notes.align = "l", font.size = "footnotesize", type = "latex", column.labels= c("Modelo 1", "Modelo 2", "Modelo 3"))
}else{
# library("texreg")
htmlreg(list(R1, R2, R3, R4, R5), star.symbol = "\\*", doctype = FALSE, center = FALSE, caption = " (\\#tab:ejercomp1) Modelos estimados por MCO", bold = 0.05, caption.above = TRUE, custom.model.names = c("Modelo 1", "Modelo 2", "Modelo 3", "Modelo 4", "Modelo 5"))
}| Modelo 1 | Modelo 2 | Modelo 3 | Modelo 4 | Modelo 5 | |
|---|---|---|---|---|---|
| (Intercept) | 95755.20*** | 311892.45*** | 76169.70*** | -31642.69 | 92505.85*** |
| (22438.88) | (15639.65) | (20492.61) | (19588.98) | (21602.52) | |
| CE | 85.76*** | 101.62*** | 121.95*** | 87.63*** | |
| (8.39) | (6.80) | (8.46) | (7.76) | ||
| CD | 40.91 | -382.80*** | 32.51 | ||
| (77.22) | (39.67) | (75.25) | |||
| Ldies | -55.18* | -123.94*** | -52.43* | ||
| (20.47) | (27.36) | (19.78) | |||
| LEI | -148.37*** | -176.15** | -214.08*** | -148.84*** | |
| (23.93) | (57.66) | (14.99) | (23.66) | ||
| V | 0.00 | 0.00 | |||
| (0.00) | (0.00) | ||||
| R2 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 |
| Adj. R2 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 |
| Num. obs. | 34 | 34 | 34 | 34 | 34 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||||
Ahora calculemos las métricas de bondad de ajuste. Estas se resumen en el Cuadro 16.2.
| Modelo 1 | Modelo 2 | Modelo 3 | Modelo 4 | Modelo 5 | |
|---|---|---|---|---|---|
| R2.ajustado | 1.000 | 0.998 | 1.000 | 0.999 | 1.000 |
| AIC | 571.512 | 634.546 | 580.247 | 603.108 | 569.986 |
| BIC | 582.196 | 640.651 | 586.352 | 610.739 | 579.144 |
| Fuente: elaboración propia. |
Así el modelo 5 es el que tiene las mejores métricas de bondad de ajuste.
Dado que todos los modelos están anidados, procedamos a compararlos empleando pruebas de hipótesis. Puedes realizar todas las comparaciones posibles. En todos los casos podemos encontrar que el mejor modelo es el modelo (4.11) (Modelo 5).
## Analysis of Variance Table
##
## Model 1: I ~ CE + CD + Ldies + LEI
## Model 2: I ~ CE + CD + Ldies + LEI + V
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 29 26690966
## 2 28 26321440 1 369526 0.3931 0.535816.3 Capítulo 5
16.3.1 Continuación Ejercicio DBLANCOS
Tras estimar los modelos se obtienen los resultados reportados en el Cuadro 16.3.
res1a <- lm(q ~ D + inversión , datos.dummy)
datos.dummy$Dinversión <- datos.dummy$D*datos.dummy$inversión
res1b <- lm(q ~ inversión + Dinversión, datos.dummy)
if(knitr::is_latex_output()){
stargazer(res1, res1a, res1b, align=FALSE, t.auto=TRUE, p.auto = TRUE,
title="Modelos del ejercicio estimados por MCO", label ="tab:ejerdum1b", header= FALSE, table.placement="H",
notes.align = "l", font.size = "tiny", type = "latex")
}else{
# library("texreg")
htmlreg(list(res1, res1a, res1b), star.symbol = "\\*", doctype = FALSE, center = FALSE, caption = " (\\#tab:ejerdum1b) Modelos del ejercicio estimados por MCO", bold = 0.05, caption.above = TRUE)
}| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | -2299.81*** | 215805.52 | -42420.02 |
| (542.51) | (196090.46) | (191355.03) | |
| timedelta | 2.01*** | ||
| (0.29) | |||
| n_tokens_title | 113.10*** | ||
| (28.52) | |||
| num_hrefs | 31.98*** | ||
| (6.01) | |||
| num_self_hrefs | -62.93*** | ||
| (16.79) | |||
| num_imgs | 19.29* | ||
| (7.59) | |||
| average_token_length | -461.22*** | ||
| (91.19) | |||
| data_channel_is_lifestyle | -439.37 | ||
| (262.57) | |||
| data_channel_is_entertainment | -775.43*** | ||
| (156.01) | |||
| kw_min_max | -0.00* | ||
| (0.00) | |||
| kw_min_avg | -0.42*** | ||
| (0.07) | |||
| kw_max_avg | -0.21*** | ||
| (0.02) | |||
| kw_avg_avg | 1.84*** | ||
| (0.11) | |||
| self_reference_max_shares | 0.01*** | ||
| (0.00) | |||
| weekday_is_monday | 504.49** | ||
| (157.13) | |||
| global_subjectivity | 2420.73*** | ||
| (675.57) | |||
| avg_negative_polarity | -1777.66*** | ||
| (514.28) | |||
| is_weekend | 415.82* | ||
| (174.89) | |||
| D | -463504.94*** | ||
| (72900.21) | |||
| inversión | 4.84*** | 5.35*** | |
| (0.38) | (0.39) | ||
| Dinversión | -0.92*** | ||
| (0.14) | |||
| R2 | 0.02 | 0.48 | 0.48 |
| Adj. R2 | 0.02 | 0.47 | 0.47 |
| Num. obs. | 39644 | 228 | 228 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||
Es importante anotar que se requiere de un truco para estimar el modelo solo con interacción y sin intercepto. Para esto es necesario crear la variable de interacción aparte, pues la función lm() automáticamente se incluye el elemento de interacción, incluye el cambio en el intercepto. Por eso este modelo se estimó de la siguiente manera:
datos.dummy$Dinversión <- datos.dummy$D*datos.dummy$inversión
res1b <- lm(q ~ inversión + Dinversión, datos.dummy)Ahora comparemos estos modelos.
## Analysis of Variance Table
##
## Response: q
## Df Sum Sq Mean Sq F value Pr(>F)
## D 1 1.3319e+13 1.3319e+13 44.113 2.3e-10 ***
## inversión 1 4.8963e+13 4.8963e+13 162.163 < 2e-16 ***
## Residuals 225 6.7936e+13 3.0194e+11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: q
## Df Sum Sq Mean Sq F value Pr(>F)
## inversión 1 5.0076e+13 5.0076e+13 165.984 < 2.2e-16 ***
## Dinversión 1 1.2260e+13 1.2260e+13 40.638 1.025e-09 ***
## Residuals 225 6.7881e+13 3.0169e+11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: Lnih
## Df Sum Sq Mean Sq F value Pr(>F)
## yedu 1 306.49 306.486 911.905 < 2.2e-16 ***
## exp 1 66.29 66.286 197.226 < 2.2e-16 ***
## I(exp^2) 1 4.24 4.238 12.611 0.0003961 ***
## Residuals 1411 474.23 0.336
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: Lnih
## Df Sum Sq Mean Sq F value Pr(>F)
## yedu 1 306.49 306.486 911.905 < 2.2e-16 ***
## exp 1 66.29 66.286 197.226 < 2.2e-16 ***
## I(exp^2) 1 4.24 4.238 12.611 0.0003961 ***
## Residuals 1411 474.23 0.336
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: Lnih
## Df Sum Sq Mean Sq F value Pr(>F)
## yedu 1 306.49 306.486 911.905 < 2.2e-16 ***
## exp 1 66.29 66.286 197.226 < 2.2e-16 ***
## I(exp^2) 1 4.24 4.238 12.611 0.0003961 ***
## Residuals 1411 474.23 0.336
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Es decir, el modelo con solo cambio en el intercepto o solo cambio en pendiente es mejor que el modelo sin ningún cambio. Al igual que lo es el que tiene el cambio en ambos. No obstante, el modelo con cambio tanto en intercepto como en pendiente muestra coeficientes individualmente no significativos para los parámetros asociados con las variables dummy. Esto hace que esta decisión no sea fácil, pero claramente el modelo con cambio en pendiente e intercepto es mejor que los otros tres.
16.4 Capítulo 8
16.4.1 Ejercicio práctico: Fabricante de automóviles
El primer paso es cargar los datos.
## 'data.frame': 38 obs. of 9 variables:
## $ Edad : int 46 31 23 19 23 47 30 28 23 29 ...
## $ Peso : int 180 175 100 185 159 170 137 192 150 120 ...
## $ Altura_desde_zapato : num 187 168 154 190 178 ...
## $ Altura.desde.el.pie.descalzo: num 185 166 152 187 174 ...
## $ Altura.sentado : num 95.2 83.8 82.9 97.3 93.9 92.4 87.7 96.9 91.4 85.2 ...
## $ Longitud_brazo : num 36.1 32.9 26 37.4 29.5 36 32.5 35.8 29.4 26.6 ...
## $ Longitud_muslo : num 45.3 36.5 36.6 44.1 40.1 43.2 35.6 39.9 35.5 31 ...
## $ Longitud_piernainferior : num 41.3 35.9 31 41 36.9 37.4 36.2 43.1 33.4 32.8 ...
## $ Distancia_centro : num -206.3 -178.2 -71.7 -257.7 -173.2 ...Y procedamos a estimar el modelo y guardémolo en le objeto R1.
Ahora miremos si existe multicolinealidad
## Edad Peso
## 1.126843 3.458461
## Altura.desde.el.pie.descalzo Altura_desde_zapato
## 282.611282 278.854402# matriz X
XTX <- model.matrix(R1)
# se calculan los valores propios
e <- eigen(t(XTX) %*% XTX)
# se muestran los valores propios
e$val## [1] 3.206402e+06 2.027935e+04 8.958858e+03 8.322439e+00 7.090495e-02# se crea el valor propio mas grande
lambda.1 <- max(e$val)
# se crea el valor propio mas pequeño
lambda.k <- min(e$val)
# se calcula kappa
kappa <- sqrt(lambda.1/lambda.k)
kappa## [1] 6724.666Claramente las dos pruebas muestran un problema de multicolinealidad. En este caso la mejor opción para corregir el problema es eliminar variables con \(VIF\) grande. En este caso.
##
## Call:
## lm(formula = myForm, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -91.029 -25.249 0.294 25.200 54.717
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 532.877125 137.104715 3.887 0.000448 ***
## Edad 0.557597 0.406456 1.372 0.179094
## Peso -0.008688 0.310512 -0.028 0.977842
## Altura_desde_zapato -4.178042 0.996501 -4.193 0.000186 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.55 on 34 degrees of freedom
## Multiple R-squared: 0.6549, Adjusted R-squared: 0.6244
## F-statistic: 21.5 on 3 and 34 DF, p-value: 5.46e-08## Edad Peso Altura_desde_zapato
## 1.080473 3.418028 3.417264Nota que para un análisis de los resultados, sería necesario eliminar las variables que no son significativas.
16.5 Capítulo 9
16.5.1 Similación de Monte Carlo
El siguiente código permite realizar el experimento de Monte Carlo para una muestra de 20.
# se fija una semilla
set.seed(1234557)
# creamos un objeto nulo para guardar si se rechaza
# h0 o no para el coeficiente de X2
X2.rechaza <- matrix("NA", N, 3)
# h0 o no para el coeficiente de X3
X3.rechaza <- matrix("NA", N, 3)
# se fija el tamaño de cada muestra
n=20
# se fija el número de repeticiones en 10mil
N=10000
for (i in 1:N){
# creación de variable explicativa
X2 <- matrix(rnorm(n),n,1)
# creación de variable no significativa
X3 <- matrix(rnorm(n),n,1)
# error homoscedástico
err <- rnorm(n)
# error heteroscedástico
herr <- (X2^2)*err
# variable explicativa sin heteroscedasticidad
y1 <- 1 + X2*7 + err
# variable explicativa con heteroscedasticidad
y2 <- 1 + X2*7 + herr
# estimación del modelo con homoscedasticidad
res.homo <- summary(lm(y1 ~ X2 + X3))
# guardamos los resultados regresión homoscedasticidad
X2.rechaza[i,1] <- res.homo$coefficients[2,4] < 0.01
X3.rechaza[i,1] <- res.homo$coefficients[3,4] < 0.01
# estimación del modelo con heteroscedasticidad
res.hetero.1 <- lm(y2 ~ X2 + X3)
res.hetero <- summary(res.hetero.1)
# guardamos los resultados regresión heteroscedasticidad con MCO
X2.rechaza[i,2] <- res.hetero$coefficients[2,4]< 0.01
X3.rechaza[i,2] <- res.hetero$coefficients[3,4]< 0.01
# guardamos los resultados regresión heteroscedasticidad con HC4
coef.test.HC4<- coeftest(res.hetero.1, vcov = (vcovHC(res.hetero.1, "HC4")))
X2.rechaza[i,3] <- coef.test.HC4[2,4]< 0.01
X3.rechaza[i,3] <- coef.test.HC4[3,4]< 0.01
}
# poniendo el nombre de las columnas
X2.rechaza <- as.data.frame(X2.rechaza)
names(X2.rechaza) <- c("MCOsinHetero","MCOconHetero", "HC3")
X3.rechaza <- as.data.frame(X3.rechaza)
names(X3.rechaza) <- c("MCOsinHetero","MCOconHetero", "HC3")Ahora, miremos la proporción de rechazos de la hipótesis nula de no significancia para el caso de la pendiente que acompaña a \(x2\). Recuerden que en este caso la hipótesis nula es incorrecta.
X2.rechaza[,1] <- as.logical(X2.rechaza[,1])
X2.rechaza[,2] <- as.logical(X2.rechaza[,2])
X2.rechaza[,3] <- as.logical(X2.rechaza[,3])
apply(X2.rechaza,2,mean) ## MCOsinHetero MCOconHetero HC3
## 1.0000 0.9997 0.9706Ahora podemos replicar este resultado para los diferentes tamaños de muestra. ¿Qué concluyes?
16.6 Capítulo 11
16.6.1 Experimento de Monte Carlo 1
Partamos del siguiente modelo \[\begin{equation*} y_t = 1 + 7 x_{t} + \varepsilon_t \end{equation*}\] donde \(\varepsilon_t\) estará autocorrelacionado (AR1) o no. Adicionalmente, estimemos un modelo con dos variables explicativas, \(x_i\) y \(z_i\), para entender el comportamiento de las pruebas de significancia individual. Noten que la segunda variable no es significativa por construcción.
Consideremos el siguiente código que materializa este experimento.
set.seed(123)
# creamos un objeto nulo para guardar los coeficientes de x
x.estcoef <- NULL
# creamos un objeto nulo para guardar los coeficientes de z
z.estcoef <- NULL
n=50
N=10000
for (i in 1:N){
# Una variable explicativa
x <- matrix(rnorm(n),n,1)
#Una variable explicativa no significativa
z <- matrix(rnorm(n),n,1)
# error no.autocorrelacionado
err <- arima.sim(model = list(order = c(0, 0, 0)), n = n)
# error autocorrelacionado
Auto.err <- arima.sim(model = list(ar = 0.7), n = n)
y1 <- 1 + x*7 + err # variable explicativa sin auto
y2 <- 1 + x*7 + Auto.err # variable explicativa con auto
res.no.auto <- summary(lm(y1 ~ x + z))
x.cf.no.auto<-res.no.auto$coefficients[2,1]
z.cf.no.auto<-res.no.auto$coefficients[3,1]
res.auto <- summary(lm(y2 ~ x + z))
x.cf.auto <- res.auto$coefficients[2,1]
z.cf.auto <- res.auto$coefficients[3,1]
# guardando los resultados de la iteración
x.estcoef<-rbind(x.estcoef,cbind(x.cf.no.auto, x.cf.auto))
z.estcoef<-rbind(z.estcoef,cbind(z.cf.no.auto, z.cf.auto))
}Ahora, comparemos los resultados para el caso de un error no autocorrelacionado y uno si autocorrelacionado para la pendiente que acompaña a \(x_i\).
## x.cf.no.auto x.cf.auto
## 7 7## x.cf.no.auto x.cf.auto
## 0.14 0.20Los resultados muestran que para el estimador de la pendiente de \(x_i\):
- Los estimadores MCO siguen siendo insesgados en presencia de autocorrelación.
- El error estándar de los coeficientes es más grande en presencia de autocorrelación.
Ahora, hagamos lo mismo para la pendiente que acompaña a \(z_i\).
## z.cf.no.auto z.cf.auto
## 0 0## z.cf.no.auto z.cf.auto
## 0.15 0.20Los resultados son similares al caso de la variable que no debería estar en el modelo.
Por otro lado, miremos por un momento los histogramas de la distribución de los coeficientes estimados.
par(mfrow=c(2,2))
hist(x.estcoef[,1], 50, main = "Pendiente de x con no.auto", xlim = c(6, 8))
hist(x.estcoef[,2], 50, main = "Pendiente de x con auto", xlim = c(6, 8))
hist(z.estcoef[,1], 50, main = "Pendiente de z con no.auto", xlim = c(-1, 1))
hist(z.estcoef[,2], 50, main = "Pendiente de z con auto", xlim = c(-1, 1))Es interesante que la distribución sigue teniendo forma acampana. ¿Qué concluyes?
16.6.2 Experimento de Monte Carlo 2
Miremos lo que ocurre con la prueba de hipótesis individuales. Es decir, miremos que proporción de veces se rechaza la hipótesis nula de que cada una de las pendientes es igual a cero con y sin autocorrelación. Y finalmente el efecto de emplear una corrección H.A.C para el caso de errores con autocorrelación.
Usemos el mismo experimento del ejercicio anterior para ver esto, pero usemos una muestra de 500. (tu podrás reproducir este ejercicio para otros tamaños de muestra).
set.seed(123)
n=500
N=10000
# creamos un objeto nulo para guardar si se
# rechaza h0 o no para el coeficiente de z
z.rechaza <- matrix("NA", N, 3)
# creamos un objeto nulo para guardar si se
# rechaza h0 o no para el coeficiente de x.
# Las columnas serán los diferentes métodos
x.rechaza <- matrix("NA", N, 3)
library(sandwich)
library(AER)
library(lmtest)
for (i in 1:N){
# Una variable explicativa
x <- matrix(rnorm(n),n,1)
#Una variable explicativa no significativa
z <- matrix(rnorm(n),n,1)
# error no.autocorrelacionado
err <- arima.sim(model = list(order = c(0, 0, 0)), n = n)
# error autocorrelacionado
Auto.err <- arima.sim(model = list(ar = 0.7), n = n)
# variable explicativa sin auto
y1 <- 1 + x*7 + err
# variable explicativa con auto
y2 <- 1 + x*7 + Auto.err
modelo.no.auto <- lm(y1 ~ x + z)
res.no.auto <- summary(modelo.no.auto)
x.rechaza[i,1] <- res.no.auto$coefficients[2,4] < 0.01
z.rechaza[i,1]<- res.no.auto$coefficients[3,4] < 0.01
modelo.auto <- lm(y2 ~ x + z)
res.auto <- summary(lm(y2 ~ x + z))
x.rechaza[i,2] <- res.auto$coefficients[2,4] < 0.01
z.rechaza[i,2] <- res.auto$coefficients[3,3] < 0.01
# con corrección HCA NeweyWest
coef.test.NeweyWesy<- coeftest(modelo.auto,
vcov = NeweyWest(modelo.auto))
x.rechaza[i,3] <- coef.test.NeweyWesy[2,4] < 0.01
z.rechaza[i,3] <- coef.test.NeweyWesy[3,4] < 0.01
}
# poniendo el nombre de las columnas
x.rechaza <- as.data.frame(x.rechaza)
names(x.rechaza) <- c("OLS_sin_auto","OLS_con_auto", "HCA_NeweyWest")
z.rechaza <- as.data.frame(z.rechaza)
names(z.rechaza) <- c("OLS_sin_auto","OLS_con_auto", "HCA_NeweyWest")Ahora, miremos la proporción de rechazos de la hipótesis nula de no significancia para el caso de la pendiente que acompaña a \(x\). Recuerden que en este caso la hipótesis nula es incorrecta.
## 'data.frame': 10000 obs. of 3 variables:
## $ OLS_sin_auto : chr "TRUE" "TRUE" "TRUE" "TRUE" ...
## $ OLS_con_auto : chr "TRUE" "TRUE" "TRUE" "TRUE" ...
## $ HCA_NeweyWest: chr "TRUE" "TRUE" "TRUE" "TRUE" ...x.rechaza[,1] <- as.logical(x.rechaza[,1])
x.rechaza[,2] <- as.logical(x.rechaza[,2])
x.rechaza[,3] <- as.logical(x.rechaza[,3])
str(x.rechaza)## 'data.frame': 10000 obs. of 3 variables:
## $ OLS_sin_auto : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ OLS_con_auto : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
## $ HCA_NeweyWest: logi TRUE TRUE TRUE TRUE TRUE TRUE ...## OLS_sin_auto OLS_con_auto HCA_NeweyWest
## 1 1 1Los resultado muestran que independiente del método que se emplee en todos los casos se rechaza la hipótesis nula correctamente.
Finalmente, concentrémonos en el efecto sobre la prueba de significancia para el coeficiente asociado a \(z\). Recuerden que en este caso la hipótesis nula es correcta.
## 'data.frame': 10000 obs. of 3 variables:
## $ OLS_sin_auto : chr "FALSE" "FALSE" "FALSE" "FALSE" ...
## $ OLS_con_auto : chr "TRUE" "FALSE" "TRUE" "TRUE" ...
## $ HCA_NeweyWest: chr "FALSE" "FALSE" "FALSE" "FALSE" ...z.rechaza[,1] <- as.logical(z.rechaza[,1])
z.rechaza[,2] <- as.logical(z.rechaza[,2])
z.rechaza[,3] <- as.logical(z.rechaza[,3])
str(z.rechaza)## 'data.frame': 10000 obs. of 3 variables:
## $ OLS_sin_auto : logi FALSE FALSE FALSE FALSE FALSE TRUE ...
## $ OLS_con_auto : logi TRUE FALSE TRUE TRUE FALSE TRUE ...
## $ HCA_NeweyWest: logi FALSE FALSE FALSE FALSE FALSE FALSE ...## OLS_sin_auto OLS_con_auto HCA_NeweyWest
## 0.0107 0.4995 0.0116Ahora, noten que sin autocorrelación y empleando MCO, la proporción de veces que se rechaza incorrectamente la nula de no significancia es aproximadamente 1%. El nivel de significancia seleccionado.
Si se emplea erróneamente MCO en presencia de autocorrelación, se rechaza erróneamente el 50% de las veces la hipótesis nula. ¡Como tirar una moneda! Pero si se emplea la corrección H.A.C. de Newey West, la proporción de rechazos erróneos se aproxima al nivel de sigificancia empleado.
Noten que esto muestra que la solución de HCA mejora la inferencia en presencia de autocorrelación.
Ahora podemos replicar este resultado para los diferentes tamaños de muestra. ¿Qué concluyes?