Variables Macroeconomicas

A lo largo del trabajo se usará el software de programación R, el cual es considerado como un lenguaje estadístico, el proyecto fue iniciado por Robert Gentleman y Ross Ihaka en el departamento de estadística de la Universidad de Auckland en 1995, actualmente este lenguaje es de software libre, por lo que no se requiere una licencia y es desarrollado por voluntarios. En esta sección se trabajará con la base precargada, “state.x77”, la cual contiene 50 renglones correspondientes a información económica de cada uno estados que conforman la Unión Americana, con 8 columnas que hacen referencia a: - Population: Es la población estimada del estado correspondiente al primero de julio de 1975. - Income: Indice per capita en 1974. - Illiteracy: Porcentaje de alfabetización del estado correspondiente en 1974. - Life Exp: Esperanza de vida en años de los habitantes del estado correspondiente en 1969-71. - Murder: Tasa de homicidios y asesinatos del estado correspondiente por cada 100,000 habitantes. - HS Grad: Porcentaje de universitarios del estado correspondiente en 1970. - Frost: El promedio de días con bajas temperaturas de 1931 a 1960 del estado correspondiente. - Area: Área del estado en millas cuadradas.

Estos datos pueden ser visualizados mediante el siguiente código, en el presente trabajo sólo se muestran los primeros 5 resultados.

datos = as.data.frame(state.x77) ##Accede a los datos y se cambia el formato a una tabla

Se desea analizar la relación lineal que existe entre la tasa de asesinatos y la esperanza de vida, es decir, se ajustará el comportamiento de la tasa de homicidios, la cual se denotará con la variable \(y\), mediante la variable relacionada a la esperanza de vida en un determinado estado, la cual se denominará \(x\). Para ello se conformará un conjunto de vectores \((x_i, y_i)\) para cada pareja de observaciones, los cuales se grafican en el eje cartesiano, observando que la muestra sigue una tendencia inversamente proporcional, es decir, entre mayor sea la esperanza de vida en un estado entonces la tasa de homicidios tiende a disminuir de forma lineal.

library(ggplot2)
X=datos[,4]   ##renglon 4 correspondiente a la esperanza de vida
Y=datos[,5]  ##renglon 5 correspondiente a la tasa de homicidios

##Se agregan la información en una tabla
Observaciones=data.frame(X,Y)  

##Se procede a graficar los datos
##De la matriz observaciones toma a X y Y
plot = ggplot(Observaciones, mapping=aes(X, Y))+  
geom_point(color="royalblue", alpha=0.7, size=5)+ ##realiza la grafica de puntos
labs( title= "Observaciones de esperanza de vida vs. tasa de homicidios", 
x="Esperanza de Vida",
y="Tasa de Homicidios")+
geom_smooth(method="lm", color="darkred")+  ##Marca la linea de tendencia
theme_minimal()+
theme(legend.title.align=0.5)
plot
## `geom_smooth()` using formula = 'y ~ x'

La línea marcada de color rojo es realizada con el valor estimado de \(y\) utilizando un primer modelo de regresión lineal simple automatizado de la forma: \[\hat{ y }=\hat{ \beta }_0+\hat{ \beta }_1 x.\]

Se sabe que los estimadores \(\hat{\beta}_0\), \(\hat{\beta}_1\) están dados de la siguiente manera:

  1. \(\hat{\beta}_0=\bar{ y }-\hat{\beta}_1\bar{ x }\).

  2. \(\hat{\beta}_1=\frac{S_{xy}}{S_{xx}}=\frac{\sum_{ i=1 }^n (x_i-\bar{ x })(y_i-\bar{ y })}{\sum_{i=1}^n (x_i-\bar{ x })^2}\).

Por lo que realizando las operaciones mediante el siguiente código se sabe que: \(\hat{\beta}_1=-2.147\) y $_0= 159.576 $:

xbarra=mean(X)  ##calcular la media muestral de X
ybarra=mean(Y)  ##calcular la media muestral de Y

Sxy=sum((X-xbarra)*(Y-ybarra)) #Calculo de Sxy
Sxx=sum((X-xbarra)^2)     #Calculo de Sxx


B1=Sxy/Sxx           #Estimador de B1
B0=ybarra-xbarra*B1  #Estimador de B0

##Imprimir resultados
B0
## [1] 159.5757
B1
## [1] -2.147301

Sustituyendo los valores de \(\hat{\beta}_0\) y \(\hat{\beta}_1\), se tiene que el modelo de regresión lineal o de valores ajustados está dado por la siguiente ecuación:

\[\hat{ y } = 159.576 - 2.147x.\]

El anterior resultado puede realizarse de una manera más directa usando la función lm cuya implementación y resultado sería:

Modelo = lm(Y ~ X, data = Observaciones) 
summary(Modelo)
## 
## Call:
## lm(formula = Y ~ X, data = Observaciones)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7272 -1.6733 -0.1734  1.4909  4.8680 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  159.576     17.579   9.078 5.45e-12 ***
## X             -2.147      0.248  -8.660 2.26e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.33 on 48 degrees of freedom
## Multiple R-squared:  0.6097, Adjusted R-squared:  0.6016 
## F-statistic: 74.99 on 1 and 48 DF,  p-value: 2.26e-11

Este proceso automatizado de R es de gran importancia, ya que resume y muestra toda la información relacionada al modelo de regresión lineal. Para poder visualizar estos valores, se debe ingresar el comando summary. De esta manera se mostrará un análisis cuantitativo de los residuales al mostrar los quantiles, así como el mínimo y máximo de estos, el valor numérico de los coeficientes asociados (\(\hat{\beta}_0\), \(\hat{\beta}_1\)) así como la desviación estándar de los estimadores y el p-value de la hipótesis \(H_0: \hat{\beta}_j = 0\), para \(j=0,1\). Para finalizar, la función summary también muestra el error estándar del modelo \(\hat{\sigma}\), el coeficiente de determinación \(R^2\), así como el p-value de la significancia del modelo.

summary(Modelo)
## 
## Call:
## lm(formula = Y ~ X, data = Observaciones)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7272 -1.6733 -0.1734  1.4909  4.8680 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  159.576     17.579   9.078 5.45e-12 ***
## X             -2.147      0.248  -8.660 2.26e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.33 on 48 degrees of freedom
## Multiple R-squared:  0.6097, Adjusted R-squared:  0.6016 
## F-statistic: 74.99 on 1 and 48 DF,  p-value: 2.26e-11

Aunque esta técnica es muy efectiva, se presentará la forma en la que se construye cada una de estas estadísticas para que se verifiquen los resultados mostrados a lo largo del capítulo. Para ello, se comprobará el teorema \(\ref{teo:3}\), el cual establece que \(\sum_{i=1}^n e_i = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = 0\), esto con la ayuda del siguiente código:

Suma_ei = sum(Y - B0 - B1 * X)
Suma_ei
## [1] -3.979039e-13
###El anterior código puede ser realizado de forma alternativa como:
Suma_eiM = sum(Modelo$residuals) ## Modelo$residuals forma de acceder a los residuales
Suma_eiM
## [1] -4.218847e-15

Se observa que en ambos casos, se obtiene que la suma de los residuales es cercana a 0, \(Suma_{ei} = -2.842171e-14\) y \(Suma_{eiM} = 3.441691e-15\).

De forma teórica, debería converger a 0, sin embargo, por cuestiones de decimales, muchas veces se obtienen aproximaciones cercanas a 0.

Es natural cuestionarse por la varianza del modelo. Para ello, se sabe que el estimador de \(\sigma^2\) es:

\[\hat{\sigma}^2 = \frac{\sum_{i=1}^n (y_i - \hat{y})^2}{n-2}.\]

Además, se sabe por el teorema \(\ref{teo:5}\) que las varianzas de la estimación en los parámetros son:

  1. \(Var(\hat{\beta}_0) = \left(\frac{1}{n} + \frac{\bar{X}^2}{S_{xx}}\right) \sigma^2\).

  2. \(Var(\hat{\beta}_1) = \frac{1}{S_{xx}} \sigma^2\).

De esta manera, realizando el código que se muestra a continuación, se conoce el valor numérico de las siguientes varianzas:

n = length(X) ##número de observaciones
sigma2 = sum(Modelo$residuals^2) / (n-2) ##cálculo de la varianza del modelo
sigma = sqrt(sigma2)    ##Desviación estándar del modelo

Var_B1 = sigma2 / Sxx      ##Varianza del estimador B1
Var_B0 = (1/n + xbarra^2 / Sxx) * sigma2  ##Varianza del estimador B0

####imprimir resultados
sigma2
## [1] 5.429329
Var_B0
## [1] 309.0105
Var_B1
## [1] 0.06148799

Con estos valores se construyen los intervalos de confianza para \(\beta_0\) y \(\beta_1\) con un nivel de significancia \(\alpha\). En este ejercicio se trabajará con un nivel de significancia \(\alpha = 0.05\). Por lo que observando las secciones 1.5.1 y 1.5.2 se tiene los siguientes intervalos de confianza:

  1. \(\beta_{0} \in \left(\hat{\beta}_{0} - t_{n-2}^{\alpha/2}\sqrt{\hat{Var}(\hat{\beta}_{0})}, \hat{\beta}_{0} + t_{n-2}^{\alpha/2}\sqrt{\hat{Var}(\hat{\beta}_{0})}\right)\).

  2. \(\beta_{1} \in \left(\hat{\beta}_{1} - t_{n-2}^{\alpha/2}\sqrt{\hat{Var}(\hat{\beta}_{1})}, \hat{\beta}_{1} + t_{n-2}^{\alpha/2}\sqrt{\hat{Var}(\hat{\beta}_{1})}\right)\).

Los cuales son calculados en el siguiente código:

alpha = 0.05   ##nivel de significancia
t = qt(1 - alpha/2, df = n - 2)  ##Obtención del cuantil
intervalo_B0 = B0 + c(-1, 1) * t * sqrt(Var_B0)  #cálculo de intervalo de B0
intervalo_B1 = B1 + c(-1, 1) * t * sqrt(Var_B1)  #cálculo de intervalo de B1

##Se agrega la información a una tabla del tipo data frame
Intervalos = as.data.frame(rbind(intervalo_B0, intervalo_B1))
colnames(Intervalos) = c("Inferior", "Superior") #cambio de nombres columna
rownames(Intervalos) = c("B0", "B1")  #cambio de nombre del renglón

##Imprimir resultados
Intervalos
##      Inferior   Superior
## B0 124.231359 194.920026
## B1  -2.645874  -1.648729

Sin embargo, existe una manera directa de obtener los intervalos de confianza y es mediante la función confint(Modelo, level = 0.95), la cual recibe como parámetros el modelo obtenido con la función lm y el nivel de confianza del intervalo.

confint(Modelo, level = 0.95)
##                  2.5 %     97.5 %
## (Intercept) 124.231359 194.920026
## X            -2.645874  -1.648729

Como se observa, los resultados son similares al implementar ambos códigos.

Finalmente, se calcula la ANOVA para medir la asociación lineal entre la variable respuesta \(Y\) y la variable explicativa \(x\).

GL = c(1, n - 2, n - 1) ##Grados de Libertad

###Calculo de la suma de cuadrados
SC_reg = B1 * Sxy     
SC_total = sum((Y - ybarra)^2)
SC_error = SC_total - SC_reg
SC = c(SC_reg, SC_error, SC_total)

##Cuadrados Medios
CM = c(SC_reg, SC_error / (n - 2), 0)

##Estadístico F
F = c(SC_reg / (SC_error / (n - 2)), 0, 0)

###Acomodar información en tabla y cambiar nombres
ANOVA = data.frame(GL, SC, CM, F)
rownames(ANOVA) = c("Regresión", "Error", "Total")
colnames(ANOVA) = c("GL", "SC", "CM", "F")

##imprimir resultado
ANOVA
##           GL       SC         CM        F
## Regresión  1 407.1380 407.138028 74.98865
## Error     48 260.6078   5.429329  0.00000
## Total     49 667.7458   0.000000  0.00000
#cálculo del p-value
Pivalue = pf(F[1], df1 = 1, df2 = n - 2, lower.tail = FALSE)
Pivalue
## [1] 2.26007e-11

Por lo que de esta forma se observa que se rechaza la hipótesis nula de hipótesis, con un nivel de confianza el 99%, por lo que se tiene evidencia de que \(\beta_1 \neq 0\), por lo que existe una asociación lineal entre la variable respuesta \(y\) y la variable explicativa \(x\).

El anterior código puede ser simplificado con el siguiente resultado, en el cual ambos métodos obtienen resultados similares:

anova(Modelo)
## Analysis of Variance Table
## 
## Response: Y
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## X          1 407.14  407.14  74.989 2.26e-11 ***
## Residuals 48 260.61    5.43                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1