Taller I - Análisis de regresión
Angélica Méndez - Tomás Ávila - Álvaro García
20 de diciembre de 2024
\(~\)
Con los datos del archivo “ExpImpDatos.txt”:
Primer punto
Estimar el modelo de inversión extranjera (Invext) como variable dependiente y, como independientes, el tiempo (mes), importaciones (Import), exportaciones (Export) e Ipc. Hacerlo con el modelo completo y con el modelo particionado; para este último, la matriz \(\mathbf{X}_2=[\mbox{Export}\quad\mbox{Ipc}]\). El modelo es
\[ \mbox{Invext}=\beta_0+\beta_1\mbox{mes}+\beta_2\mbox{Import}+\beta_3\mbox{Export}+\beta_4\mbox{Ipc}+\varepsilon. \]
Solución
En primer lugar, creamos la matriz \(\mathbf{X}\) del modelo
\[ \vec{Y}=\mathbf{X}\vec{\beta}+\vec{\varepsilon}, \]
donde \(\vec{Y}=[\mbox{Invext}]\), \(\vec{\beta}=(\beta_0, \beta_1, \beta_2, \beta_3, \beta_4)\), \(\vec{\varepsilon}=(\varepsilon_0, \varepsilon_1, \varepsilon_2, \varepsilon_3, \varepsilon_4)\) y \(\mathbf{X}=[\mathbf{1}\quad\mbox{mes}\quad\mbox{Import}\quad\mbox{Export}\quad\mbox{Ipc}]\). De esta manera, \(\mathbf{X}\in\mathcal{M}_{n\times p}(\mathbb{R})\) con \(n\) el número de observaciones, \(p:=k+1\) el número de parámetros involucrados en el modelo y \(k\) el número de regresores. En este caso, \(\mathbf{X}\in\mathcal{M}_{24\times5}(\mathbb{R})\):
ExpImp <- read.table("C:/Users/USER/Downloads/ExpImpDatos.txt", header = T)
attach(ExpImp)
#Asignación de variables:
Y = Invext
n = length(Y)
uno <- rep(1,n) #Vector de unos
mes = mes
import = Import
export = Export
ipc = Ipc
#Construcción de la matriz X:
X = cbind(uno, mes, import, export, ipc); X## uno mes import export ipc
## [1,] 1 1 2193.2 1087.5 1.56
## [2,] 1 2 2242.6 921.5 1.22
## [3,] 1 3 1725.2 941.6 0.47
## [4,] 1 4 2159.2 826.7 0.03
## [5,] 1 5 1899.0 924.1 0.29
## [6,] 1 6 2015.4 884.8 0.35
## [7,] 1 7 1646.5 855.8 0.17
## [8,] 1 8 1736.5 952.6 0.91
## [9,] 1 9 1412.6 711.9 2.21
## [10,] 1 10 1297.5 763.3 1.70
## [11,] 1 11 2330.9 946.1 0.94
## [12,] 1 12 1343.3 831.8 0.78
## [13,] 1 13 1398.5 926.1 0.48
## [14,] 1 14 1552.1 1038.8 0.28
## [15,] 1 15 1260.5 1018.8 0.31
## [16,] 1 16 1273.4 1000.7 0.50
## [17,] 1 17 1177.1 1061.1 0.33
## [18,] 1 18 1074.0 964.5 0.35
## [19,] 1 19 1173.4 1068.9 0.48
## [20,] 1 20 1189.3 1243.4 0.53
## [21,] 1 21 1184.9 1012.8 1.29
## [22,] 1 22 1374.5 1059.4 2.30
## [23,] 1 23 1273.0 1068.5 1.71
## [24,] 1 24 2664.6 892.8 1.00Luego, estimamos los coeficientes que conforman al vector \(\vec{\beta}\) del modelo usando que \(\widehat{\vec{\beta}}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\vec{Y}\):
#Cálculo de beta-estimado:
Be = solve(t(X)%*%X)%*%t(X)%*%Y #Estimación mediante operaciones matriciales
round(t(Be), 3)## uno mes import export ipc
## [1,] 15598.56 225.139 -0.564 1.064 -138.261Para visualizar los valores de \(\mbox{Invext}\) que arroja el modelo con estos coeficientes (las entradas del vector \(\widehat{\vec{Y}}\)) respecto de los observados en la base de datos (las entradas del vector \(\vec{Y}\)), calculamos \(\widehat{\vec{Y}}\) como \(\mathbf{X}\widehat{\vec{\beta}}\):
Ye = X%*%Be #Y estimado
plot(Y, col = "darkseagreen", xlim = c(0,25), main = "Modelo ajustado", xlab = "Mes", ylab = "Invext", pch = 19, bty = "n", yaxt="n", xaxt="n")
lines(Ye, lwd = 2, col = "darkslateblue")
axis(2, at = c(15000, 17000, 19000, 21000))
axis(1, at = c(0, 4, 8, 12, 16, 20, 24))
legend(1, 20700, legend = c("Invext observado", "Invext estimado"), pch = c(19, NA), lty = c(NA, 1), lwd = c(NA, 2), col = c("darkseagreen", "darkslateblue"), bty = "n")
Ahora, calculamos los coeficientes de \(\widehat{\vec{\beta}}\) mediante el comando \(\verb|lm|\) de R:
##
## Call:
## lm(formula = Invext ~ mes + Import + Export + Ipc, data = ExpImp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1081.98 -277.83 -26.47 231.19 913.43
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15598.5580 1016.9621 15.338 3.72e-12 ***
## mes 225.1392 17.1544 13.124 5.63e-11 ***
## Import -0.5640 0.2445 -2.307 0.0325 *
## Export 1.0639 0.9433 1.128 0.2734
## Ipc -138.2609 153.1738 -0.903 0.3780
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 461.6 on 19 degrees of freedom
## Multiple R-squared: 0.9471, Adjusted R-squared: 0.936
## F-statistic: 85.1 on 4 and 19 DF, p-value: 7.417e-12## (Intercept) mes Import Export Ipc
## 15598.558 225.139 -0.564 1.064 -138.261Viendo ambas estimaciones juntas:
## uno mes import export ipc
## [1,] 15598.56 225.139 -0.564 1.064 -138.261
## [2,] 15598.56 225.139 -0.564 1.064 -138.261Finalmente, estimamos los coeficientes de \(\vec{\beta}\) mediante el modelo particionado, donde
\[ \mathbf{X}=[\mathbf{X_1}\quad|\quad\mathbf{X}_2]:=[\mathbf{1}\quad\mbox{mes}\quad\mbox{Import}\quad|\quad\mbox{Export}\quad\mbox{Ipc}] \]
y
\[ \vec{\beta}=(\vec{\beta}_1,\vec{\beta}_2):=\begin{bmatrix} \beta_0\\ \beta_1\\ \beta_2\\ \rule[.5ex]{2em}{.4pt}\\ \beta_3\\ \beta_4\end{bmatrix}. \]
Para hacer las estimaciones matricialmente, utilizamos que
\[ \widehat{\vec{\beta}}_1=\widetilde{\vec{\beta}}_1-(\mathbf{X}_1^T\mathbf{X}_1)^{-1}\mathbf{X}_1^T\mathbf{X}_2\widehat{\vec{\beta}}_2\quad\mbox{y}\quad\widehat{\vec{\beta}}_2=(\mathbf{X}_2^T\mathbf{M}_1\mathbf{X}_2)^{-1}\mathbf{X}_2^T\mathbf{M}_1\vec{Y}, \] donde \(\mathbf{M}_1\) y \(\mathbf{M}_2\) son las matrices creadoras residuales asociadas a \(\mathbf{X}_1\) y \(\mathbf{X}_2\) respectivamente, y \(\widetilde{\vec{\beta}}_1:=(\mathbf{X}_1^T\mathbf{X}_1)^{-1}\mathbf{X}_1^T\vec{Y}\).
#Construcción de algunas matrices:
X1 = cbind(uno, mes, Import) #Matriz X1
X2 = cbind(Export, Ipc) #Matriz X2
I = diag(n) #Matriz identidad de orden 24
H1 = X1%*%solve(t(X1)%*%X1)%*%t(X1) #Matriz hat asociada a X1
M1 = I-H1 #Matriz creadora residual asociada a X1
#Estimación del vector beta-2 mediante operaciones matriciales:
B2e = solve(t(X2)%*%M1%*%X2)%*%t(X2)%*%M1%*%Y
#Construcción del vector beta-1 moño:
B1tilde = solve(t(X1)%*%X1)%*%t(X1)%*%Y
#Estimación del vector beta-1 mediante operaciones matriciales:
B1e = B1tilde-solve(t(X1)%*%X1)%*%t(X1)%*%X2%*%B2eLas estimaciones son:
## 15598.56 225.139 -0.564 | 1.064 -138.261De esta manera, \(\widehat{\vec{\beta}}=\left(\widehat{\vec{\beta}}_1,\widehat{\vec{\beta}}_2\right)\) es (escrito de forma transpuesta):
## uno mes Import Export Ipc
## [1,] 15598.56 225.139 -0.564 1.064 -138.261Cuyos coeficientes coinciden con los calculados en el modelo completo:
## uno mes import export ipc
## [1,] 15598.56 225.139 -0.564 1.064 -138.261
## [2,] 15598.56 225.139 -0.564 1.064 -138.261
## [3,] 15598.56 225.139 -0.564 1.064 -138.261Segundo punto
Para el modelo completo, probar la hipótesis
\[ H_0\colon\mathbf{K}^T\vec\beta=\vec{m}\quad\mbox{vs.}\quad H_1\colon\mathbf{K}^T\vec\beta\neq\vec{m} \]
con \(\vec{\beta}^T=[\beta_0\quad\beta_1\quad\beta_2\quad\beta_3\quad\beta_4]\) y \(\vec{m}^T=[0\quad0]\). En particular, el sistema de hipótesis
\[ \begin{matrix}H_0\colon\beta_3=0,\\\phantom{H_0\colon}\beta_4=0\end{matrix}\quad\mbox{vs.}\quad\begin{matrix}H_1\colon\beta_3\neq0,\\\phantom{H_1\colon}\beta_4\neq0\end{matrix}. \]
Definir \(\mathbf{K}\), \(\vec{\beta}\) y \(\vec{m}\) para formular los sistemas de hipótesis.
Solución
La ecuación matricial que nos permite obtener el sistema de hipótesis planteado es:
\[ \underbrace{\begin{bmatrix} 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix}}_{\mathbf{K}^T} \underbrace{\begin{bmatrix} \beta_0\\ \beta_1\\ \beta_2\\ \beta_3\\ \beta_4 \end{bmatrix}}_{\vec{\beta}}=\underbrace{\begin{bmatrix} 0 \\ 0 \end{bmatrix}}_{\vec{m}}, \]
pues, multiplicando el lado izquierdo de la igualdad:
\[ \underbrace{\begin{bmatrix} \beta_3\\ \beta_4 \end{bmatrix}}_{\mathbf{K}^T\vec{\beta}}=\underbrace{\begin{bmatrix} 0 \\ 0 \end{bmatrix}}_{\vec{m}}. \]
Esto último es \(H_0\).
La estadística de prueba que empleamos para juzgar \(H_0\) frente a \(H_1\) es
\[ F:=\frac{(\mathbf{K}^T\widehat{\vec{\beta}}-\vec{m})^T[\mathbf{K}^T(\mathbf{X}^T\mathbf{X})^{-1}\textbf{K}]^{-1}(\mathbf{K}^T\widehat{\vec{\beta}}-\vec{m})}{s\widehat{\sigma}^2}\overset{H_0}{\sim}\mathrm{F}_{n-k-1}^s, \]
donde \(s:=\mathrm{rg}(\mathbf{K})\), \(\mathrm{F}_{n-k-1}^s\) es la distribución F con \(s\) grados de libertad en el numerador y \(n-k-1\) grados de libertad en el denominador, y
\[ \widehat{\sigma}^2=\mathrm{CME}=\frac{\mathrm{SCE}}{df\mathrm{E}}=\frac{\mathrm{SCE}}{n-k-1}=\frac{\vec{Y}^T\mathbf{M}\vec{Y}}{n-k-1}. \]
Para el caso particular de este ejercicio, la estadística de prueba queda de la siguiente manera:
\[ \begin{align*} F&\overset{H_0}{=}\frac{\left(\begin{bmatrix}\widehat{\beta}_3\\\widehat{\beta}_4\end{bmatrix}-\begin{bmatrix}0\\0\end{bmatrix}\right)^T[\mathbf{K}^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{K}]^{-1}\left(\begin{bmatrix}\widehat{\beta}_3\\\widehat{\beta}_4\end{bmatrix}-\begin{bmatrix}0\\0\end{bmatrix}\right)}{2\frac{\vec{Y}^T\mathbf{M}\vec{Y}}{24-4-1}}\sim\mathrm{F_{24-4-1}^2}\\ &=\frac{\begin{bmatrix}\widehat{\beta}_3&\widehat{\beta}_4\end{bmatrix}[\mathbf{K}^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{K}]^{-1}\begin{bmatrix}\widehat{\beta}_3\\\widehat{\beta}_4\end{bmatrix}}{\frac{2}{19}{\vec{Y}^T\mathbf{M}\vec{Y}}}\sim\mathrm{F_{19}^2}. \end{align*} \]
La decisión será rechazar \(H_0\) si un valor observado de \(F\), digamos \(F_c\), cumple que \(F_c>\mathrm{F}_{19, 1-\alpha}^2\), donde \(\mathrm{F}_{19, 1-\alpha}^2\) es el percentil \(1-\alpha\) de la distribución \(\mathrm{F}_{19}^2\). Calculamos \(F_c\) con los coeficientes \(\widehat{\beta}_3\) y \(\widehat{\beta}_4\) obtenidos en el primer punto:
#Cálculo de Fc:
H = X%*%solve(t(X)%*%X)%*%t(X) #Matriz hat asociada a X
K <- matrix(c(0, 0, 0, 1, 0, 0, 0, 0, 0, 1), nrow = 5, ncol = 2)
s = qr(K)$rank #Rango de K
bet34 = Be[c(4,5)] #Vector de beta_3 y beta_4 estimados
W = solve(t(K)%*%solve(t(X)%*%X)%*%K)
M = I-H #Matriz creadora residual asociada a X
SCE = t(Y)%*%M%*%Y
Fc = as.numeric(t(bet34)%*%W%*%bet34/(s/19)%*%SCE); Fc #F calculado## [1] 1.300997Si establecemos un nivel de significancia \(\alpha=0.05\), entonces \(\mathrm{F}_{19, 1-\alpha}^2\) vale:
## [1] 3.521893Esta situación se ilustra en la siguiente gráfica:
lb = 0; ub = qf(1-0.05, 2, 19)
x <- seq(0, 4, length=10000)
hx <- df(x, 2, 19)
plot(x, hx, type="n", xlab="valores de X", ylab="", xaxt = "n", axes=FALSE, main = TeX("$P(0<X<F_{19, 0.95}^2)=0.95$"))
i <- x >= lb & x <= ub
lines(x, hx)
polygon(c(lb,x[i],ub), c(0,hx[i],0), col="darkseagreen")
text(0.75, 0.6, TeX("$F_{19}^2$"))
text(0.7, 0.2, "0.95")
axis(1, at = c(0, 1, Fc, 2, 3, ub, 4), labels = c(0, 1, TeX("$F_c$"), 2, 3, TeX("$F_{19, 0.95}^2$"), 4))
abline(v = Fc, lty = 2, lwd = 2, col = "darkslateblue")
Por lo tanto, ya que \(F_c<\mathrm{F_{19, 0.95}^2}\), con una confiabilidad del 95%, no se rechaza \(H_0\). Es decir, existe suficiente evidencia empírica para no rechazar que \(\beta_3=0\) y \(\beta_4=0\).
Tercer punto
Para el modelo particionado, probar
\[ H_0\colon\vec{\beta}_2=\mathbf{0}\quad\mbox{vs.}\quad H_1\colon\vec{\beta}_2\neq\mathbf{0}. \]
Solución
Como se mencionó en la segunda parte del primer punto, \(\vec{\beta}_2=(\beta_3, \beta_4)\).
La estadística de prueba que empleamos para juzgar \(H_0\) frente a \(H_1\) es
\[ F:=\frac{\mathrm{SCR}(\vec{\beta}_2\,|\,\vec{\beta}_1)/r}{\widehat{\sigma}^2}\overset{H_0}{\sim}\mathrm{F}_{n-k-1}^r, \]
donde \(r\) es el número de variables en \(\mathbf{X}_2\) y la suma extra de cuadrados debida a \(\vec{\beta}_2\) es \(\mathrm{SCR}(\vec{\beta}_2\,|\,\vec{\beta}_1)=\mathrm{SCR}(\vec{\beta})-\mathrm{SCR}(\widetilde{\vec{\beta}}_1)\). Estas sumas de cuadrados de la regresión se calculan como
\[ \mathrm{SCR}(\vec{\beta})=\vec{Y}^T\left(\mathbf{H}-\bar{\mathbf{J}}_n\right)\vec{Y}\quad\mbox{y}\quad\mathrm{SCR}(\widetilde{\vec{\beta}}_1)=\vec{Y}^T\left(\mathbf{H}_1-\bar{\mathbf{J}}_n\right)\vec{Y}, \] donde \(\bar{\mathbf{J}}_n:=1/n\,\mathbf{J}_n:=1/n\,\mathbf{1}_n\mathbf{1}_n^T\) con \(\mathbf{1}_n\) vector columna formado por \(n\) unos.
Es fácil ver que \(\mathrm{SCR}(\vec{\beta}_2\,|\,\vec{\beta}_1)=\vec{Y}^T(\mathbf{H-\mathbf{H}_1)}\vec{Y}\). De esta manera, calculamos \(\mathrm{SCR}(\vec{\beta}_2\,|\,\vec{\beta}_1)\):
Para el caso particular de este ejercicio, la estadística de prueba queda como
\[ F=\frac{\mathrm{SCR}(\vec{\beta}_2\,|\,\vec{\beta}_1)}{2\widehat{\sigma}^2}=\frac{\vec{Y}^T(\mathbf{H}-\mathbf{H}_1)\vec{Y}}{2\widehat{\sigma}^2}\overset{H_0}{\sim}\mathrm{F}_{19}^2. \]
Como en el punto anterior, la decisión será rechazar \(H_0\) si un valor observado de \(F\) es mayor a \(F_{19,1-\alpha}^2\). Teniendo en cuenta el último cálculo efectuado, hallamos el valor observado de \(F\):
## [1] 1.300997Que resulta siendo el mismo del modelo sin particionar. Si establecemos de nuevo un nivel de significancia de 0.05, entonces el percentil en esta prueba nos queda igual al de la anterior; i. e., 3.522. Como el valor observado es menor que esta cantidad, con una confiabilidad del 95%, no se rechaza \(H_0\) .
En síntesis, el percentil \(\mathrm{F}_{19,0.95}^2\), el valor calculado para la estadística de prueba en el modelo completo, y el valor calculado para la estadística de prueba en el modelo particionado son, respectivamente:
## [,1]
## 3.521893
## Fc 1.300997
## Fc2 1.300997\(~\)