Econometria: regressão por matrizes

Abstract

This is an undergrad student level exercise for class use. We estimate an equation using matrices.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Econometria: regressão por matrizes. Campo Grande-MS,Brasil: RStudio/Rpubs, 2020. Disponível em http://rpubs.com/amrofi/estimation_matrix.

1 Chamando os dados:

# os dados são chamados fazendo a rotina abaixo, saída do `dput(dados)`:
dados <- structure(list(y = c(70, 65, 90, 95, 110, 115, 120, 140, 155, 150), x = c(80, 
    100, 120, 140, 160, 180, 200, 220, 240, 260), obs = c(1, 2, 3, 4, 5, 6, 7, 8, 
    9, 10)), row.names = c(NA, -10L), class = c("tbl_df", "tbl", "data.frame"))
print(dados)  # funcao para visualizar os dados na tela do RStudio

     y   x obs
1   70  80   1
2   65 100   2
3   90 120   3
4   95 140   4
5  110 160   5
6  115 180   6
7  120 200   7
8  140 220   8
9  155 240   9
10 150 260  10

library(knitr)  # pacote para a funcao kable
attach(dados)  # funcao para que o RStudio entenda os rótulos de 'dados'

# Estatisticas Descritivas

summary(dados)

       y                x            obs       
 Min.   : 65.00   Min.   : 80   Min.   : 1.00  
 1st Qu.: 91.25   1st Qu.:125   1st Qu.: 3.25  
 Median :112.50   Median :170   Median : 5.50  
 Mean   :111.00   Mean   :170   Mean   : 5.50  
 3rd Qu.:135.00   3rd Qu.:215   3rd Qu.: 7.75  
 Max.   :155.00   Max.   :260   Max.   :10.00  

2 Estimação por matrizes

2.1 Criar matrizes X e y1

# Para que o R entenda os rótulos, usarei attach(dados)
attach(dados)
# Para criar a matriz X coluna 1 com valores unitários e x1
X <- as.matrix.data.frame(cbind(1, dados[, 2]))
print(X)

      [,1] [,2]
 [1,]    1   80
 [2,]    1  100
 [3,]    1  120
 [4,]    1  140
 [5,]    1  160
 [6,]    1  180
 [7,]    1  200
 [8,]    1  220
 [9,]    1  240
[10,]    1  260

# Para transformar a variavel (y) em vetor:
y1 <- as.matrix(dados[, 1])
print(y1)

      [,1]
 [1,]   70
 [2,]   65
 [3,]   90
 [4,]   95
 [5,]  110
 [6,]  115
 [7,]  120
 [8,]  140
 [9,]  155
[10,]  150

2.2 Estimar beta

# Para estimar o vetor (X'X) da equacao, primeiro se obtem o parametro pelo
# seguintes passos: 1) transposta de X: (X')
trX <- (t(X))
# 2) Produto da transposta de X por X, com o codigo %*%: (X'X)
X_X <- trX %*% X
X_X

     [,1]   [,2]
[1,]   10   1700
[2,] 1700 322000

# Para obter a inversa (X'X)-1 se deve primeiro ativar o pacote library(matlib),
# e usar inv() para inverter a matriz
det(X_X)

[1] 330000

library(matlib)
invX_X <- (inv(X_X))
invX_X  # (X'X)-1

            [,1]        [,2]
[1,]  0.97575758 -0.00515152
[2,] -0.00515152  0.00003030

# Uma vez que se tem a inversa (X'X)-1, se procede o produto X'y:
Xy <- trX %*% y1
Xy  # X'y

       [,1]
[1,]   1110
[2,] 205500

# Para calcular o vetor beta:
beta <- invX_X %*% Xy  # beta = (X'X)-1 X'y
beta

           [,1]
[1,] 24.4535538
[2,]  0.5084628

2.2.1 Pela função nativa do R lm()

attach(dados)
reg1 <- lm(y ~ x, data = dados)  # comparar beta com coefficients do summary   
# sumario de resultados do objeto reg1: saída da regressao
summary(reg1)

Call:
lm(formula = y ~ x, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.364  -4.977   1.409   4.364   8.364 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 24.45455    6.41382   3.813  0.00514 ** 
x            0.50909    0.03574  14.243 5.75e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.493 on 8 degrees of freedom
Multiple R-squared:  0.9621,    Adjusted R-squared:  0.9573 
F-statistic: 202.9 on 1 and 8 DF,  p-value: 5.753e-07

3 Cálculos para testes dos parâmetros

3.1 Desvios-padrões e testes t

# agora calculamos os desvios-padrões obtencao de sigma quadrado estimado
yly <- t(y) %*% y
blXy <- t(beta) %*% Xy
ulu <- yly - blXy
ulu  # Soma dos Quadrados dos Residuos

         [,1]
[1,] 337.2727

uhat <- y - X %*% (beta)
uhat  # residuos estimados

             [,1]
 [1,]   4.8181818
 [2,] -10.3636364
 [3,]   4.4545455
 [4,]  -0.7272727
 [5,]   4.0909091
 [6,]  -1.0909091
 [7,]  -6.2727273
 [8,]   3.5454545
 [9,]   8.3636364
[10,]  -6.8181818

# obter sigma quadrado para gl=n-p=8
n <- nrow(dados)  # numero de observacoes
k <- ncol(X)  # numero de parametros a estimar
sigsqhat <- as.numeric(ulu/(n - k))
sigsqhat  # s2 - variancia estimada do residuo

[1] 42.15909

# variancia-covariancia de beta
varcovbeta <- sigsqhat * invX_X
varcovbeta

           [,1]         [,2]
[1,] 41.1370523 -0.217183196
[2,] -0.2171832  0.001277548

# obtendo a raiz da variancia da diagonal de varcovbeta desvio padrao de beta
sebeta <- sqrt(diag(varcovbeta))
sebeta

[1] 6.41381730 0.03574281

# parametros beta
beta

           [,1]
[1,] 24.4545455
[2,]  0.5090909

# estatistica t dos parametros
tbeta <- beta/sebeta
tbeta

          [,1]
[1,]  3.812791
[2,] 14.243171

# obtendo a probabilidade de tbeta para 5% e df=n-k
pvalue <- 2 * pt(-abs(tbeta), n - k)
pvalue  # probabilidade de t dos parametros

             [,1]
[1,] 5.142172e-03
[2,] 5.752746e-07

4 Obtendo \(R^2\) e \(R^2 ajustado\)

# soma dos quadrados totais
sqtot <- yly - n * mean(y1)^2
sqtot

     [,1]
[1,] 8890

# soma dos quadrados dos residuos
sqres <- ulu
sqres

         [,1]
[1,] 337.2727

# soma dos quadrados da regressao
sqreg <- sqtot - sqres
sqreg

         [,1]
[1,] 8552.727

# coeficiente de determinação R2
r2 <- sqreg/sqtot
r2

          [,1]
[1,] 0.9620616

# coeficiente de determinação ajustado R2-ajustado
r2aj <- 1 - (sqres/(n - k))/(sqtot/(n - 1))
r2aj

          [,1]
[1,] 0.9573193

5 Coeficiente F de significação

F <- (sqreg/(k - 1))/(sqres/(n - k))
F

         [,1]
[1,] 202.8679

# probabilidade de F para H0: betas iguais a zero
probF <- 1 - pf(F, k - 1, n - k)
probF

             [,1]
[1,] 5.752746e-07

6 Comparar com saída da função lm

summary(reg1)

Call:
lm(formula = y ~ x, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.364  -4.977   1.409   4.364   8.364 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 24.45455    6.41382   3.813  0.00514 ** 
x            0.50909    0.03574  14.243 5.75e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.493 on 8 degrees of freedom
Multiple R-squared:  0.9621,    Adjusted R-squared:  0.9573 
F-statistic: 202.9 on 1 and 8 DF,  p-value: 5.753e-07