Econometria: regressão linear por matrizes

Abstract

This is an undergrad student level instruction for class use.

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 linear por matrizes. Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em http://rpubs.com/amrofi/matrix_reg.

1 Introdução

Neste arquivo utilizo uma série fictícia da Tabela 2.4 do livro de Gujarati e Porter (2011: p.67) das Despesas familiares de consumo semanal Y e renda familiar semanal X. É um exercício rápido de regressão com uso de matrizes e comparativamente à saída da função lm() nativa do R.

2 Dados

O leitor pode carregar direto na Tabela 2.4 do livro de Gujarati e Porter (2011: p.67) como abaixo ou usar a saída do dput(dados) conforme colocado abaixo.

library(readxl)    # pacote para chamar os dados do Excel
dados <- read_excel("exercicio_estimacao_excel_anexo.xlsx", 
                  sheet = "dados")

# os dados também podem ser chamados fazendo a rotina abaixo `dput(dados)`:
# A execução da linha abaixo fornece exatamente o mesmo que a linha 2
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"))
View(dados)  # funcao para visualizar os dados na tela do RStudio
library(knitr)  # pacote para a funcao kable
kable(dados)  # funcao para visualizar os dados em uma forma de tabela

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

attach(dados)  # funcao para que o RStudio entenda os rótulos de 'dados'

# Estatisticas Descritivas

kable(summary(dados), caption = "Estatísticas descritivas")

Estatísticas descritivas

	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

3 Regressão

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

reg1 <- lm(y ~ x)  # regressao pela funcao lm nativa do R 
summary(reg1)  # ver o sumario de resultados do objeto reg1: saída da regressao

Call:
lm(formula = y ~ x)

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

Para obter uma saída mais organizada, pode-se utilizar o pacote stargazer, que fornece várias opções para a saída em formato científico.

# criterios de informacao
reg1$AIC <- AIC(reg1)
reg1$BIC <- BIC(reg1)
suppressMessages(library(stargazer))
stargazer(reg1, title = "Título: Resultado da Regressão OLS ou Mínimos Quadrados Ordinários", 
    align = TRUE, type = "text", style = "all", keep.stat = c("AIC", "BIC", 
        "rsq", "adj.rsq", "n"))

Título: Resultado da Regressão OLS ou Mínimos Quadrados Ordinários
===============================================
                        Dependent variable:    
                    ---------------------------
                                 y             
-----------------------------------------------
x                            0.509***          
                              (0.036)          
                            t = 14.243         
                            p = 0.00000        
Constant                     24.455***         
                              (6.414)          
                             t = 3.813         
                             p = 0.006         
-----------------------------------------------
Observations                    10             
R2                             0.962           
Adjusted R2                    0.957           
Akaike Inf. Crit.             69.562           
Bayesian Inf. Crit.           70.470           
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

3.2 Por matrizes

# 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 <- cbind(1, x)
print(X)

          x
 [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 se utiliza o mesmo codigo:
y1 <- cbind(y)
print(y1)

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

# 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

            x
    10   1700
x 1700 322000

# Para obter a inversa (X'X)-1 se deve primeiro ativar o pacote
# library(MASS), e usar ginv() para inverter a matriz determinante de (X'X)
# para ver se é inversível (precisa ser diferente de zero)
det(X_X)

[1] 330000

library(MASS)
invX_X <- (ginv(X_X))
invX_X  # (X'X)-1

             [,1]          [,2]
[1,]  0.975757576 -5.151515e-03
[2,] -0.005151515  3.030303e-05

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

       y
    1110
x 205500

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

              y
[1,] 24.4545455
[2,]  0.5090909

summary(reg1)  # comparar beta com coefficients do summary

Call:
lm(formula = y ~ x)

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

# 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

         y
y 337.2727

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

                y
 [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

# 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

[1] 42.15909

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

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

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

[1] 6.41381730 0.03574281

# estatistica t dos parametros
beta

              y
[1,] 24.4545455
[2,]  0.5090909

tbeta <- beta/sebeta
tbeta  # t dos parametros

             y
[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

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

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

     [,1]
[1,] 8890

sqres <- ulu
sqres  # soma dos quadrados dos residuos

         y
y 337.2727

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

         y
y 8552.727

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

          y
y 0.9620616

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

          y
y 0.9573193

F <- (sqreg/(k - 1))/(sqres/(n - k))
F  # coeficiente F de significação

         y
y 202.8679

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

             y
y 5.752746e-07

Referências

COLONESCU, Constantin. Principles of Econometrics with R. 2016. https://bookdown.org/ccolonescu/RPoE4/ .
FIGUEIREDO, A.M.R. Econometria. Campo Grande: UFMS, 2018. 225p. (Notas de aula do professor) (pp.1-52)
GUJARATI, Damodar N.; PORTER, Dawn C. Econometria básica. 5.ed. Porto Alegre: AMGH/Bookman/McGraw-Hill do Brasil, 2011.
ROMERO, Luis Quintana; GONZÁLEZ, Miguel Ángel Mendoza (coords.). Econometría aplicada utilizando R. Mexico, Universidad Nacional Autónoma de México/Facultad de Estudios Superiores Acatlán. Primera edición, marzo 2016. Disponible en: http://saree.com.mx/econometriaR/
WOOLDRIDGE, J.M. Introdução à econometria. São Paulo: CENGAGE, 4.ed. 2011. (cap. 1,2,3,4,7)