taller5

Librerias requeridas

Si no tiene alguna de estas librerías instaladas debe hacerlo antes de leerlas

r = getOption("repos")
r["CRAN"] = "http://cran.us.r-project.org"
options(repos = r)
library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(sandwich)
library(FSA)

## ## FSA v0.8.30. See citation('FSA') if used in publication.
## ## Run fishR() for related website and fishR('IFAR') for related book.

library(MASS)
library(nlme)

Ejercicio 1

Haga una rutina en R que implemente el método de GLS para este problema especíﬁco, simule unos datos y veriﬁque que su programa estima bien los parámetros para muestras grandes, T = 1000

En primer lugar simulamos los datos para el tamaño de muestra sugerido:

# Una variable Z_t estrictamente exógena y conocida (ponemos una distribución que no sea normal para que la varianza explote)
z_t<-rpois(1000,12)
# Calculamos unos errores normales preliminares
e_p<-rnorm(1000,0,0.6)
# Finalmente calculamos los errores del ejercicio
e_t<-e_p*z_t
# simulamos la matriz de X, para ello suponemos 3 variables
m_x<-c(14,15,16)
s_xx<-matrix(c(4,0,0,0,3,0,0,0,5), ncol=3)
X<-mvrnorm(1000,m_x,s_xx)
#suponesmos unos valores cualesquiera de Beta para hallar Y
betas1<-c(1.5,-0.3,0.6)
Y_t<-X%*%betas1+e_t
#guardamos los resultados en un data frame
dat1<-data.frame(Y_t,X,e_t)

Ahora obtenemos los betas estimados con el modelo simple

m1e1<-lm(Y_t~X1+X2+X3+0, data=dat1)
summary(m1e1)

## 
## Call:
## lm(formula = Y_t ~ X1 + X2 + X3 + 0, data = dat1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.9011  -4.6081   0.1107   4.6161  27.4846 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## X1   1.5510     0.1054  14.715  < 2e-16 ***
## X2  -0.3016     0.1060  -2.845  0.00453 ** 
## X3   0.5754     0.0920   6.255  5.9e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.647 on 997 degrees of freedom
## Multiple R-squared:  0.9235, Adjusted R-squared:  0.9232 
## F-statistic:  4010 on 3 and 997 DF,  p-value: < 2.2e-16

El modelo con errores robustos (empleando la matriz \(V^{-1}\) dentro del estimador)

coeftest(m1e1, vcov = vcovHC(m1e1, type = "HC0"))

## 
## t test of coefficients:
## 
##     Estimate Std. Error t value  Pr(>|t|)    
## X1  1.551010   0.103835  14.937 < 2.2e-16 ***
## X2 -0.301615   0.106279  -2.838  0.004632 ** 
## X3  0.575415   0.089894   6.401 2.372e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Empleando GLS alcanzables

m1fgls<-lm(Y_t~X1+X2+X3+0, data=dat1, weights = 1/m1e1$fitted.values^2)
summary(m1fgls)

## 
## Call:
## lm(formula = Y_t ~ X1 + X2 + X3 + 0, data = dat1, weights = 1/m1e1$fitted.values^2)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07317 -0.17327  0.00512  0.17815  1.21802 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## X1  1.52549    0.10564  14.440  < 2e-16 ***
## X2 -0.29431    0.10278  -2.863  0.00428 ** 
## X3  0.58868    0.09227   6.380  2.7e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3004 on 997 degrees of freedom
## Multiple R-squared:  0.9173, Adjusted R-squared:  0.917 
## F-statistic:  3686 on 3 and 997 DF,  p-value: < 2.2e-16

Empleando WLS

disegno_m1e1<-model.matrix(m1e1)
y_m<-as.matrix(dat1[,'Y_t'])
weight_m<-matrix(0,1000,1000)
diag(weight_m)<-z_t
betas1_wls<-solve(t(disegno_m1e1) %*% solve(weight_m) %*% disegno_m1e1) %*% (t(disegno_m1e1) %*% solve(weight_m) %*% y_m)
betas1_wls

##          [,1]
## X1  1.5744619
## X2 -0.3027422
## X3  0.5597089

Los resultados con WLS (el numeral C) son similares a los obtenidos con los errores robustos. no obstante, no son tan buenos como los GLS alcanzables.

Ejercicio 2

Haga una rutina en R que implemente el método de GLS para este problema especíﬁco, simule unos datos y veriﬁque que su programa estima bien los parámetros para muestras grandes, T = 1000

Para este caso cambiamos los valores de \(\epsilon_t\) suponiendo unos valores cualesquiera para \(\delta_0\) y \(\delta_1\)

e_t2<-e_p*(0.8+0.7*z_t)
# Reestimamos los valores de Y_t
Y_t2<-X%*%betas1+e_t2
#guardamos los resultados en un data frame
dat2<-data.frame(Y_t2,X,e_t2)

Calculamos el modelo lineal simple

m1e2<-lm(Y_t2~X1+X2+X3+0, data=dat2)
summary(m1e2)

## 
## Call:
## lm(formula = Y_t2 ~ X1 + X2 + X3 + 0, data = dat2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.4372  -3.5387   0.0818   3.5234  20.3960 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## X1  1.54089    0.08016  19.222  < 2e-16 ***
## X2 -0.30127    0.08063  -3.736 0.000197 ***
## X3  0.57999    0.06997   8.289 3.65e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.816 on 997 degrees of freedom
## Multiple R-squared:  0.954,  Adjusted R-squared:  0.9539 
## F-statistic:  6899 on 3 and 997 DF,  p-value: < 2.2e-16

El modelo con errores robustos (empleando la matriz \(V^{-1}\) dentro del estimador)

coeftest(m1e2, vcov = vcovHC(m1e2, type = "HC0"))

## 
## t test of coefficients:
## 
##     Estimate Std. Error t value  Pr(>|t|)    
## X1  1.540886   0.079037 19.4958 < 2.2e-16 ***
## X2 -0.301272   0.080802 -3.7285 0.0002034 ***
## X3  0.579991   0.068441  8.4744 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo WLS

disegno_m1e2<-model.matrix(m1e2)
y_m2<-as.matrix(dat2[,'Y_t2'])
weight_m<-matrix(0,1000,1000)
diag(weight_m)<-0.8+0.7*z_t
betas2_wls<-solve(t(disegno_m1e2) %*% solve(weight_m) %*% disegno_m1e2) %*% (t(disegno_m1e2) %*% solve(weight_m) %*% y_m2)
betas2_wls

##          [,1]
## X1  1.5577726
## X2 -0.3021528
## X3  0.5688189

Ejercicio 3

Haga una rutina en R que implemente el método de GLS para este problema especíﬁco, simule unos datos y veriﬁque que su programa estima bien los parámetros para muestras grandes T = 1000

Establecemos un valor para \(\rho\) cualquiera, y calculamos la matriz \(V\)

rho<-0.5
v_e3<-diag(1,nrow=1000,1000)
matriz_rho<-function(a,b){
  for(i in 2:1000){
    a[i, i-1]<-b
    a[i-1, i]<-b
  }
  return(a)
}
v_e3<-matriz_rho(v_e3,rho)

Con esta matriz calculamos los errores correlacionados en orden 1 y las demás variables

e_t3<-e_p%*%v_e3
e_t3<-t(e_t3)
# Reestimamos los valores de Y_t
Y_t3<-X%*%betas1+e_t3
#guardamos los resultados en un data frame
dat3<-data.frame(Y_t3,X,e_t3)

Calculamos el modelo lineal simple

m1e3<-lm(Y_t3~X1+X2+X3+0, data=dat3)
summary(m1e3)

## 
## Call:
## lm(formula = Y_t3 ~ X1 + X2 + X3 + 0, data = dat3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.34250 -0.49006  0.00357  0.49497  2.43256 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## X1  1.50548    0.00999  150.70   <2e-16 ***
## X2 -0.29766    0.01005  -29.62   <2e-16 ***
## X3  0.59690    0.00872   68.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7248 on 997 degrees of freedom
## Multiple R-squared:  0.9992, Adjusted R-squared:  0.9992 
## F-statistic: 4.384e+05 on 3 and 997 DF,  p-value: < 2.2e-16

El modelo con errores robustos (empleando la matriz \(V^{-1}\) dentro del estimador)

coeftest(m1e3, vcov = vcovHC(m1e3, type = "HC0"))

## 
## t test of coefficients:
## 
##      Estimate Std. Error t value  Pr(>|t|)    
## X1  1.5054764  0.0099546 151.234 < 2.2e-16 ***
## X2 -0.2976557  0.0102102 -29.153 < 2.2e-16 ***
## X3  0.5969010  0.0087957  67.863 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo WLS

disegno_m1e3<-model.matrix(m1e3)
y_m3<-as.matrix(dat3[,'Y_t3'])
betas3_wls<-solve(t(disegno_m1e3) %*% solve(v_e3) %*% disegno_m1e3) %*% (t(disegno_m1e3) %*% solve(v_e3) %*% y_m3)
betas3_wls

##          [,1]
## X1  1.5005059
## X2 -0.2993831
## X3  0.6004992

Ejercicio 4

Haga una rutina en R que implemente el método de GLS para este problema especíﬁco, simule unos datos y veriﬁque que su programa estima bien los parámetros para muestras grandes T = 1000

En este caso, partiendo de la matriz \(V\) del punto 3, se multiplica por \(\delta_0+\delta_1*Z_t\)

v_e4<-(0.8+0.7*z_t)*v_e3

con esta matriz se recalculan las variables y se estiman los modelos:

e_t4<-e_p%*%v_e4
e_t4<-t(e_t4)
# Reestimamos los valores de Y_t
Y_t4<-X%*%betas1+e_t4
#guardamos los resultados en un data frame
dat4<-data.frame(Y_t4,X,e_t4)

Calculamos el modelo lineal simple

m1e4<-lm(Y_t4~X1+X2+X3+0, data=dat4)
summary(m1e4)

## 
## Call:
## lm(formula = Y_t4 ~ X1 + X2 + X3 + 0, data = dat4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.9078  -4.6931   0.1413   4.6436  20.0703 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## X1  1.55567    0.09704  16.032  < 2e-16 ***
## X2 -0.29121    0.09761  -2.984  0.00292 ** 
## X3  0.57146    0.08470   6.747 2.55e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.041 on 997 degrees of freedom
## Multiple R-squared:  0.9351, Adjusted R-squared:  0.9349 
## F-statistic:  4787 on 3 and 997 DF,  p-value: < 2.2e-16

El modelo con errores robustos (empleando la matriz \(V^{-1}\) dentro del estimador)

coeftest(m1e4, vcov = vcovHC(m1e4, type = "HC0"))

## 
## t test of coefficients:
## 
##     Estimate Std. Error t value  Pr(>|t|)    
## X1  1.555670   0.096499 16.1212 < 2.2e-16 ***
## X2 -0.291211   0.099876 -2.9157  0.003628 ** 
## X3  0.571464   0.085857  6.6560 4.638e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo WLS

disegno_m1e4<-model.matrix(m1e4)
y_m4<-as.matrix(dat4[,'Y_t4'])
betas4_wls<-solve(t(disegno_m1e4) %*% solve(v_e4) %*% disegno_m1e4) %*% (t(disegno_m1e4) %*% solve(v_e4) %*% y_m4)
betas4_wls

##         [,1]
## X1  1.217003
## X2 -1.640373
## X3  2.049504

En este caso el moelo WLS no se ajusta tan bien a los betas supuestos inicialmente. Los ejercicios en cada uno de estos puntos, como se puede notar, tienen los mismos pasos, por lo que se pueden meter en una rutina function pero en este caso se dejé el paso a paso para que quedara claro el proceso de los diferentes casos de GLS.

Ejercicio 5

Para los ejercicios de Wooldridge debemos en primer lugar instalar la librería, y así acceder a la base de datos:

library(wooldridge)

## Warning: package 'wooldridge' was built under R version 3.6.3

## 
## Attaching package: 'wooldridge'

## The following object is masked from 'package:MASS':
## 
##     cement

data('gpa3')

Posteriormente, el ejercicio nos dice que se estima una ecuación para el promedio académico a término (trmgpa) para los estudiantes de segundo semestre de otoño. Por tanto, en primer lugar se deben obtener los subconjuntos de datos, de lo contrario, se obtendría una estimación completamente diferente a la del libro.

Los estudiantes de otoño se obtienen cuando la variable “fall” es igual a 1, mientras que los estudiantes de segundo semestre se obtienen cuando la variable “frstsem” es igual a cero, por tanto:

w4<-subset(subset(gpa3, term==1), frstsem==0)

Con estos datos podemos estimar el modelo del libro con un simple lm para los errores estándar simples, y con coeftest obtener los errores robustos.

m1<-lm(trmgpa~crsgpa+cumgpa+tothrs+sat+hsperc+female+season, data=w4)
print("Modelo con errores estándar simples")

## [1] "Modelo con errores estándar simples"

summary(m1)

## 
## Call:
## lm(formula = trmgpa ~ crsgpa + cumgpa + tothrs + sat + hsperc + 
##     female + season, data = w4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56399 -0.32599  0.05126  0.38069  1.42707 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.120089   0.551935  -3.841 0.000154 ***
## crsgpa       0.899780   0.174806   5.147 5.20e-07 ***
## cumgpa       0.192810   0.064036   3.011 0.002860 ** 
## tothrs       0.001378   0.001246   1.106 0.269898    
## sat          0.001793   0.000212   8.457 1.98e-15 ***
## hsperc      -0.003944   0.001814  -2.174 0.030600 *  
## female       0.351496   0.085411   4.115 5.19e-05 ***
## season      -0.157044   0.097678  -1.608 0.109095    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5408 on 261 degrees of freedom
## Multiple R-squared:  0.4648, Adjusted R-squared:  0.4505 
## F-statistic: 32.39 on 7 and 261 DF,  p-value: < 2.2e-16

m2<-coeftest(m1, vcov = vcovHC(m1, type = "HC0"))
print("Errores robustos")

## [1] "Errores robustos"

m2

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -2.12008930  0.54637485 -3.8803 0.0001322 ***
## crsgpa       0.89978033  0.16552415  5.4359 1.251e-07 ***
## cumgpa       0.19280964  0.07427735  2.5958 0.0099714 ** 
## tothrs       0.00137801  0.00120345  1.1450 0.2532389    
## sat          0.00179270  0.00019394  9.2434 < 2.2e-16 ***
## hsperc      -0.00394381  0.00190132 -2.0742 0.0390368 *  
## female       0.35149557  0.07867221  4.4678 1.178e-05 ***
## season      -0.15704384  0.07950372 -1.9753 0.0492874 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En el libro, los errores simples se encuentran entre paréntesis mientras que los errores robustos se encuentran entre cajas. Teniendo en cuenta lo anterior, los resultados concuerdan perfectamente.

5.1

¿Las variables crsgpa, cumgpa y tothrs tienen los efectos estimados esperados? ¿Cuáles de estas variables son estadísticamente significativas al nivel del 5%? ¿Importa qué errores estándar se utilizan?

Respuesta

Estas variables tienen los signos que se esperarían con respecto a su relación lógica con el promédio académico a término.

Si un estudiante toma cursos en los que las calificaciones son, en promedio, más altas (crsgpa más altas) entonces su promedio a término será más alto.

Cuanto mejor haya sido el desempeño del estudiante en el pasado, según lo medido por cumgpa, mejor será el estudiante (en promedio) en el semestre actual.

Finalmente, tothrs es una medida de la experiencia, y su coeficiente indica un retorno creciente a la experiencia.

Las variables crsgpa y cumgpa son significativas al 5% con cualquiera de los dos tipos de errores estándar, mientras que tothrs no es significativo con ninguno de los dos. En este caso no importa el tipo de errores que se utilizan.

5.2

¿Por qué tiene sentido la hipótesis H0: \(\beta_{crgspa}=1\) ? Pruebe esta hipótesis a dos colas con un nivel de significancia del 5%, utilizando ambos errores estándar. Describa sus conclusiones.

Respuesta

Este test se puede hacer gracias a la librería FSA. Aplicamos el test al modelo con errores simples:

hoCoef(m1, 2, 1, alt = "two.sided")

##  term Ho Value  Estimate Std. Error          T  df   p value
##     2        1 0.8997803  0.1748062 -0.5733187 261 0.5669226

Y posteriormente al modelo con errores robustos. En este caso, como m2 no es un objeto lm se puede calcular a mano simplemente como: \(t_k=(b_k-1)/SE\)

t_k<-(m2[2,1]-1)/m2[2,2]
t_k

## [1] -0.6054685

# para el p-value
2*pt(-abs(t_k),df=nrow(w4)-8)

## [1] 0.5453938

La hipótesis tiene sentido porque, sin información adicional del estudiante, el mejor predictor del promedio académico al término del semestre es el promedio académico del semestre cursado. En ambos casos se acepta la hipótesis nula.

5.3

Pruebe si el hecho de que el estudiante participe en una actividad deportiva de la temporada tiene efecto sobre el promedio de calificaciones, utilizando ambos errores estándar. ¿El nivel de significancia al que se puede rechazar el valor nulo depende del error estándar utilizado?

Respuesta

En este caso no es necesario realizar la hipótesis de si el estimador es estádisticamente cero pues, el modelo lm y el Coeftest nos arrojan ese resultado en el p-value. En este caso sí importa el tipo de errores seleccionados pues, con errores robustos el pvalue es 0.049 (pasa raspando) mientras que con errores estándar simples es de 0.10 (no es significativo al 5%)

Ejercicio 6

Cargamos los datos del punto y estimamos el primer modelo fijando el subset male==1

data('econmath')
m1_6<-lm(score~colgpa+act, data=subset(econmath, male==1))
summary(m1_6)

## 
## Call:
## lm(formula = score ~ colgpa + act, data = subset(econmath, male == 
##     1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.576  -5.668   0.988   6.578  27.425 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.5185     3.7245   5.509 6.45e-08 ***
## colgpa       13.6047     0.9424  14.437  < 2e-16 ***
## act           0.6703     0.1503   4.459 1.07e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.81 on 403 degrees of freedom
##   (22 observations deleted due to missingness)
## Multiple R-squared:  0.4205, Adjusted R-squared:  0.4176 
## F-statistic: 146.2 on 2 and 403 DF,  p-value: < 2.2e-16

El modelo para las mujeres sería con el subset male==0

m2_6<-lm(score~colgpa+act, data=subset(econmath, male==0))
summary(m2_6)

## 
## Call:
## lm(formula = score ~ colgpa + act, data = subset(econmath, male == 
##     0))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.804  -6.571   0.391   7.706  32.396 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   13.786      4.111   3.354 0.000872 ***
## colgpa        11.885      1.093  10.878  < 2e-16 ***
## act            1.034      0.177   5.841 1.06e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.89 on 405 degrees of freedom
##   (20 observations deleted due to missingness)
## Multiple R-squared:  0.3666, Adjusted R-squared:  0.3635 
## F-statistic: 117.2 on 2 and 405 DF,  p-value: < 2.2e-16

En el tercer y cuarto modelo se estima para el conjunto de datos completo

m3_6<-lm(score~colgpa+act, data=econmath)
summary(m3_6)

## 
## Call:
## lm(formula = score ~ colgpa + act, data = econmath)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.492  -6.228   0.302   7.052  31.120 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  17.1189     2.8063   6.100 1.64e-09 ***
## colgpa       12.5096     0.7230  17.302  < 2e-16 ***
## act           0.8791     0.1164   7.552 1.16e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.48 on 811 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.3805, Adjusted R-squared:  0.379 
## F-statistic: 249.1 on 2 and 811 DF,  p-value: < 2.2e-16

m4_6<-lm(score~male+colgpa+act+male*colgpa+male*act, data=econmath)
summary(m4_6)

## 
## Call:
## lm(formula = score ~ male + colgpa + act + male * colgpa + male * 
##     act, data = econmath)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.576  -6.076   0.634   7.194  32.396 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.7862     3.9124   3.524 0.000449 ***
## male          6.7323     5.5493   1.213 0.225419    
## colgpa       11.8853     1.0400  11.428  < 2e-16 ***
## act           1.0338     0.1685   6.137 1.32e-09 ***
## male:colgpa   1.7194     1.4398   1.194 0.232755    
## male:act     -0.3635     0.2315  -1.570 0.116761    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.37 on 808 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.3968, Adjusted R-squared:  0.393 
## F-statistic: 106.3 on 5 and 808 DF,  p-value: < 2.2e-16

6.1

Calcule el estadístico de Chow habitual para probar la hipótesis nula de que las ecuaciones de regresión son las mismas para hombres y mujeres. Encuentra el valor p de la prueba.

Respuesta

En primer lugar creamos un objeto con la suma de los residuos al cuadrado de los modelos separados y el modelos general (hombres y mujeres)

SSR=NULL
SSR$total<-m3_6$residuals^2
SSR$H<-m1_6$residuals^2
SSR$M<-m2_6$residuals^2

Calculamos el k

K<-m3_6$rank

Estimamos la prueba de Chow

numerador<-(sum(SSR$total)-(sum(SSR$H)+sum(SSR$M)))/K
denominador<-(sum(SSR$H)+sum(SSR$M))/(nrow(econmath)-2*K)
chow<-numerador/denominador
chow

## [1] 7.629867

El valor p. En este caso recordemos que el estadístico de Chow es una F

1-pf(chow, K, (nrow(econmath)-2*K))

## [1] 4.892474e-05

El resultado muestra que las ecuaciones son diferentes.

6.2

Calcule el estadístico Chow habitual para probar la hipótesis nula de que los coeficientes de pendiente son los mismos para hombres y mujeres, e informe el valor p.

Respuesta

En este casos e realiza la misma prueba pero empleando el modelo que incluye diferencias en las pendientes dependiendo el genero (interacciones entre las variables)

SSR2=NULL
SSR2$total<-m4_6$residuals^2
SSR2$H<-m1_6$residuals^2
SSR2$M<-m2_6$residuals^2

Calculamos el k

K2<-m4_6$rank

Estimamos la prueba de Chow

numerador2<-(sum(SSR2$total)-(sum(SSR2$H)+sum(SSR2$M)))/K2
denominador2<-(sum(SSR2$H)+sum(SSR2$M))/(nrow(econmath)-2*K2)
chow2<-numerador2/denominador2
chow2

## [1] -4.715911e-14

El valor p. En este caso recordemos que el estadístico de Chow es una F

1-pf(chow2, K2, (nrow(econmath)-2*K))

## [1] 1

En este caso se acepta la hipótesis nula de que las pendientes son las mismas, por lo que el cambio está en el intercepto.

6.3

¿Tiene suficiente información para calcular versiones robustas de heteroscedasticidad de las pruebas en (ii) y (iii)? Explique.

Respuesta

No se tiene suficiente información porque, si bien el estimador es el mismo, lo que cambia es el error estándar, afectando la varianza del error de los residuales, pero no se obtiene por defecto un valor de los residuales en cada observación.

Ejercicio 7

Considere el siguiente modelo para explicar el comportamiento del sueño:

\(sleep=\beta_0+\beta_1totwrk+\beta_2educ+\beta_3age+\beta_4age^2+\beta_5yngkid+\beta_6male+u\)

Escriba un modelo que permita a la varianza de \(u\) diferir entre hombres y mujeres. La varianza no depende de otros factores.

Respuesta

la varianza estaría dada por: \(Var(u|totwrk, educ, age, yngkid,male)=Var(u|male)=\delta_0+\delta_1male\) Dónde la varianza para las mujeres sería \(\delta_0\) y para los hombres sería \(\delta_0+\delta_1\)

Use los datos en SLEEP75 para estimar los parámetros del modelo para heterocedasticidad. (Primero tiene que estimar la ecuación de sleep por OLS para obtener los residuos de OLS). ¿Es la varianza estimada de u mayor para hombres o mujeres?

Respuesta

Cargamos los datos y estimamos el modelo:

data('sleep75')
m1_7<-lm(sleep~totwrk+educ+age+age^2+yngkid+male,data=sleep75)
summary(m1_7)

## 
## Call:
## lm(formula = sleep ~ totwrk + educ + age + age^2 + yngkid + male, 
##     data = sleep75)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2349.05  -241.44     5.82   265.69  1345.44 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3640.23376  114.33203  31.839   <2e-16 ***
## totwrk        -0.16569    0.01801  -9.202   <2e-16 ***
## educ         -11.76532    5.87132  -2.004   0.0455 *  
## age            2.00994    1.52083   1.322   0.1867    
## yngkid         4.78424   50.01991   0.096   0.9238    
## male          87.54557   34.66501   2.525   0.0118 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 418 on 700 degrees of freedom
## Multiple R-squared:  0.1216, Adjusted R-squared:  0.1153 
## F-statistic: 19.38 on 5 and 700 DF,  p-value: < 2.2e-16

Estimamos un modelo para los errores al cuadrado

sleep2<-cbind(sleep75,m1_7$residuals)
m2_7<-lm((m1_7$residuals^2)~male,data=sleep2)
summary(m2_7)

## 
## Call:
## lm(formula = (m1_7$residuals^2) ~ male, data = sleep2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -190365 -155648 -113078   26824 5327646 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   190369      20394   9.335   <2e-16 ***
## male          -30235      27094  -1.116    0.265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 356700 on 704 degrees of freedom
## Multiple R-squared:  0.001766,   Adjusted R-squared:  0.0003478 
## F-statistic: 1.245 on 1 and 704 DF,  p-value: 0.2648

Como el coeficiente de “male” es negativo, la varianza es mayor para las mujeres

¿Es la varianza de u estadísticamente diferente para hombres y mujeres?

Respuesta

No, dado que en el segundo modelo la variable “male” no es estadísticamente significativa, no existe diferencia en la varianza entre hombres y mujeres.

Ejercicio 8

Use VOTE1 para este ejercicio.

Estime un modelo con \(voteA\) como variable dependiente y \(prtystrA\), \(democA\), \(log(expendA)\) y \(log(expendB)\) como variables independientes. Obtenga los residuos de OLS, \(\hat{u}_i\), y regréselos en todas las variables independientes. Explique por qué obtiene \(R^2=0\).

data('vote1')
m1_8<-lm(voteA~prtystrA+democA+log(expendA)+log(expendB), data=vote1)
summary(m1_8)

## 
## Call:
## lm(formula = voteA ~ prtystrA + democA + log(expendA) + log(expendB), 
##     data = vote1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.576  -4.864  -1.146   4.903  24.566 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.66141    4.73604   7.952 2.56e-13 ***
## prtystrA      0.25192    0.07129   3.534  0.00053 ***
## democA        3.79294    1.40652   2.697  0.00772 ** 
## log(expendA)  5.77929    0.39182  14.750  < 2e-16 ***
## log(expendB) -6.23784    0.39746 -15.694  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.573 on 168 degrees of freedom
## Multiple R-squared:  0.8012, Adjusted R-squared:  0.7964 
## F-statistic: 169.2 on 4 and 168 DF,  p-value: < 2.2e-16

Estimamos el modelo de los residuales

vote1_2<-cbind(vote1,m1_8$residuals)
m2_8<-lm(m1_8$residuals~prtystrA+democA+log(expendA)+log(expendB), data=vote1_2)
summary(m2_8)

## 
## Call:
## lm(formula = m1_8$residuals ~ prtystrA + democA + log(expendA) + 
##     log(expendB), data = vote1_2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.576  -4.864  -1.146   4.903  24.566 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -2.705e-15  4.736e+00       0        1
## prtystrA      1.563e-16  7.129e-02       0        1
## democA       -4.577e-16  1.407e+00       0        1
## log(expendA) -3.685e-16  3.918e-01       0        1
## log(expendB) -5.594e-16  3.975e-01       0        1
## 
## Residual standard error: 7.573 on 168 degrees of freedom
## Multiple R-squared:  5.211e-32,  Adjusted R-squared:  -0.02381 
## F-statistic: 2.189e-30 on 4 and 168 DF,  p-value: 1

Aunque el \(R^2\) no es exactamente cero (puede parecerlo por redondeo de la computadora), este es estimado en los OLS para que no se correlacione con las variables regresoras, hecho que explica por qué son cero incluso los \(\beta\) estimados.

Ahora, calcule la prueba de Breusch-Pagan para detectar heterocedasticidad. Use la versión estadística F e informe el valor p.

bptest(m1_8)

## 
##  studentized Breusch-Pagan test
## 
## data:  m1_8
## BP = 9.0934, df = 4, p-value = 0.05881

a un nivel de significancia del 5%, existe presencia de heterocedasticidad.

Calcule el caso especial de la prueba de White para detectar heterocedasticidad, nuevamente utilizando la forma estadística F. ¿Qué tan fuerte es la evidencia de heterocedasticidad ahora?

Esta prueba se puede obener para un objeto lm creando una base de datos con los residuales y los valores estimados del modelo 1, y regresando los residuos al cuadrado en los valores del \(y\) estimado e \(y\) estimado al cuadrado.

vote1_3<-cbind(m1_8$residuals,data=m1_8$fitted.values)
vote1_3<-data.frame(vote1_3)
colnames(vote1_3)<-c("u","voteAg")
vote1_3$voteAg2<-vote1_3$voteAg^2
m3_8<-lm((u^2)~voteAg+voteAg2, data=vote1_3)
summary(m3_8)

## 
## Call:
## lm(formula = (u^2) ~ voteAg + voteAg2, data = vote1_3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -107.65  -44.80  -29.88   23.81  539.52 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 171.85840   53.14214   3.234  0.00147 **
## voteAg       -4.26368    2.16653  -1.968  0.05070 . 
## voteAg2       0.03574    0.02124   1.682  0.09434 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 81.66 on 170 degrees of freedom
## Multiple R-squared:  0.03173,    Adjusted R-squared:  0.02034 
## F-statistic: 2.786 on 2 and 170 DF,  p-value: 0.0645

El valor del estadístico F confirma un poco más de heterocedasticidad, aunque no difiere mucho del B-P test.

taller5

Mateo Linares

31/5/2020

Librerias requeridas

Ejercicio 1

Ejercicio 2

Ejercicio 3

Ejercicio 4

Ejercicio 5

5.1

Respuesta

5.2

Respuesta

5.3

Respuesta

Ejercicio 6

6.1

Respuesta

6.2

Respuesta

6.3

Respuesta

Ejercicio 7

Respuesta

Respuesta

Respuesta

Ejercicio 8