Já pensou em fazer uma Regressão na mão?

Já pensou em gerar todos os principais resultados de uma regressão sem usar o comando regress? Apliquei este desafio na última semana a meus pupilos na monitoria de Econometria I (graduação FEA-USP). Monitor malvado, hein? Nem tanto.. Já vínhamos trabalhando bastante nas aulas e os códigos já não eram novidade.

Não queria chamar a atividade de “prova” pois os alunos já tem muitas. Também não queria chamar de “atividade avaliativa” pois soaria pedagógico demais. Ficou como Econometrics Challenge — Mission 1: Regression by hand. Eles teriam 69 minutos para gerar cada saída do comando regress (beta chapéu, erro-padrão, t, p-valor, ANOVA, R², F…), tudo usando comandos de matriz e somatório.

Aqui estão o enunciado e a base de dados. Se você nunca fez, sinta-se desafiado!

Reproduzo abaixo uma possível solução (aqui o do-file e o log-file), usando um modelo mais geral. O código foi construído em Stata, mas imagino que a lógica é muito semelhante no R e Matlab.

Tentei fazer o código de modo que apenas os comandos iniciais (passos 1 e 2) são específicos do modelo usado, os passos seguintes são comuns a toda regressão.

use C:FEAadritEconometricsChallengeMission1.dta

**************************************************************
*** PASSO 1: Declarando o número de variáveis ****************
**************************************************************

matrix n_variaveis=(6)

**************************************************************
*** PASSO 2: Definindo Y, X e Ybarra *************************
**************************************************************

* Definindo o vetor Y
mkmat log_salario, matrix(y)

* Definindo a matriz X
gen um = 1
mkmat anosdeestudo exper idade sexo branco casado um, matrix(x)

* Gerando Ybarra (será útil para calcular R2 logo mais)
egen Ybarra= mean(log_salario)

**************************************************************
*** PASSO 3: Encontrando o número de observações *************
**************************************************************

egen n_obs= total(um)
mkmat n_obs, matrix(n_obs)

**************************************************************
*** PASSO 4: Calculando os betas chapéus *********************
**************************************************************

* Definindo a matriz X'X
matrix xlinhax=x'*x

* Definindo a inversa da matriz X'X
matrix inversaxlinhax=syminv(xlinhax)

* Definindo a matriz X'Y
matrix xlinhay=x'*y

* Definindo o vetor bchapeu
matrix bchapeu=inversaxlinhax*xlinhay

* Listando bchapeu
mat list bchapeu

**************************************************************
*** PASSO 5: Calculando os erros-padrão **********************
**************************************************************

* Definindo o vetor de resíduos
matrix res=y-x*bchapeu

* Definindo o estimador de sigma quadrado
matrix sigma2chapeu=(res'*res)/(n_obs[1,1]-n_variaveis[1,1]-1)

* Listando sigma chapeu quadrado
matrix list sigma2chapeu

* Definindo a matriz de variância-covariância estimada
matrix var_cov=(inversaxlinhax)*sigma2chapeu

* Listando a matriz de variância-covariância estimada 
matrix list var_cov

* Definindo o vetor que contém as variâncias estimadas 
matrix var_betachapeu=vecdiag(var_cov)

* Listando o vetor que contém as variâncias estimadas 
mat list var_betachapeu

* Definindo o vetor que contém os erros-padrão

scalar s=n_variaveis[1,1]+1
local s=n_variaveis[1,1]+1
matrix ep_betachapeu = J(1,s,0)

forvalues i = 1/1 {
forvalues j = 1/`s' {
matrix ep_betachapeu[`i',`j']= sqrt(var_betachapeu[`i',`j'])
}
}

* Listando o vetor que contém os erros-padrão
mat list ep_betachapeu

**************************************************************
*** PASSO 6: Calculando as estatísticas t ********************
**************************************************************

* Listando bchapeu
mat list bchapeu

* Listando o vetor que contém os erros-padrão dos betas chapeus
mat list ep_betachapeu

* Veja, portanto, que bchapeu é vetor coluna e ep_betachapeu é vetor linha

* Definindo o vetor que contém as estatísticas t
matrix t_est = J(1,s,0)

forvalues i = 1/1 {
forvalues j = 1/`s' {
matrix t_est[`i',`j']= bchapeu[`j',`i'] / ep_betachapeu[`i',`j']
}
}

mat list t_est

**************************************************************
*** PASSO 7: Calculando os p-valores *************************
**************************************************************

* Definindo o vetor que contém as estatísticas t absolutas
matrix t_est_a = J(1,s,0)

forvalues i = 1/1 {
forvalues j = 1/`s' {
matrix t_est_a[`i',`j']= abs(t_est[`i',`j'])
}
}

mat list t_est_a

matrix p_value = J(1,s,0)

forvalues i = 1/1 {
forvalues j = 1/`s' {
matrix p_value[`i',`j']= 2*ttail(n_obs[1,1]-n_variaveis[1,1]-1, t_est_a[`i',`j'])
}
}

mat list p_value

**************************************************************
*** PASSO 8: Calculando a parte inferior do IC ***************
**************************************************************

matrix ic_inf = J(1,s,0)

forvalues i = 1/1 {
forvalues j = 1/`s' {
matrix ic_inf[`i',`j']= bchapeu[`j',`i'] - invttail(n_obs[1,1]-n_variaveis[1,1]-1,0.025)*ep_betachapeu[`i',`j']
}
}

mat list ic_inf

**************************************************************
*** PASSO 9: Calculando a parte superior do IC ***************
**************************************************************

matrix ic_sup = J(1,s,0)

forvalues i = 1/1 {
forvalues j = 1/`s' {
matrix ic_sup[`i',`j']= bchapeu[`j',`i'] + invttail(n_obs[1,1]-n_variaveis[1,1]-1,0.025)*ep_betachapeu[`i',`j']
}
}

mat list ic_sup

**************************************************************
*** PASSO 10: SQR ********************************************
**************************************************************

matrix y_chapeu = x*bchapeu
matrix list y_chapeu
matrix res = y-y_chapeu
matrix list res

** Sabemos que SQR=res'*res **

matrix SQR = res'*res
matrix list SQR

**************************************************************
*** PASSO 11: SQT ********************************************
**************************************************************

** Sabemos que SQT=y'y - n*(Ybarra)^2 **

matrix ylinhay = y'*y
matrix list ylinhay
mkmat Ybarra, matrix(Ybarraexp)
matrix Ybarra=( Ybarraexp[1,1])
matrix nYbarraQuad = n_obs[1,1]*Ybarra*Ybarra
matrix list nYbarraQuad
matrix SQT= ylinhay - nYbarraQuad
matrix list SQT

**************************************************************
*** PASSO 12: SQE ********************************************
**************************************************************

matrix SQE = SQT - SQR
matrix list SQE

**************************************************************
*** PASSO 13: Graus de liberdade *****************************
**************************************************************

** gl do modelo = k
matrix glmodelo = n_variaveis
matrix list glmodelo

** gl do residuo = n-k-1
matrix glresiduo = n_obs[1,1] - n_variaveis - 1
matrix list glresiduo

** gl total = n-1
matrix gltotal = n_obs[1,1] - 1
matrix list gltotal

**************************************************************
*** PASSO 14: MS *********************************************
**************************************************************

matrix MSmodelo = SQE[1,1]/glmodelo[1,1]
matrix list MSmodelo

matrix MSresiduo = SQR[1,1]/glresiduo[1,1]
matrix list MSresiduo

matrix MStotal = SQT[1,1]/gltotal[1,1]
matrix list MStotal

**************************************************************
*** PASSO 15: R2 *********************************************
**************************************************************

matrix R2 = SQE[1,1]/SQT[1,1]
matrix list R2

**************************************************************
*** PASSO 16: R2 ajustado ************************************
**************************************************************

matrix R2ajust = 1 - (1 - R2[1,1])*(n_obs[1,1] - 1)/(n_obs[1,1] - n_variaveis[1,1] - 1)
matrix list R2ajust

**************************************************************
*** PASSO 17: Root MSE ***************************************
**************************************************************

matrix RootMSE = sqrt(SQR[1,1]/(n_obs[1,1] - n_variaveis[1,1] - 1))
matrix list RootMSE

**************************************************************
*** PASSO 18: Teste F ****************************************
**************************************************************

* Estatística F *

matrix Fnumerador = (R2[1,1])/n_variaveis[1,1]
matrix list Fnumerador

matrix Fdenominador = (1 - R2[1,1])/(n_obs[1,1] - n_variaveis[1,1] - 1)
matrix list Fdenominador

matrix F= Fnumerador[1,1]/Fdenominador[1,1]
matrix list F

* P-valor *
display Ftail(n_variaveis[1,1],n_obs[1,1] - n_variaveis[1,1] - 1,F[1,1])

* Econometrics Challenge completed 🙂
**************************************************************

Divertido, né?

Agora podemos confirmar os resultados usando regress: