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:
Muito bom… quem me dera ter tido um professor que fizesse isso. Esta é de longe a melhor forma de aprender o que está por trás de cada estimação.
e os alunos? e que acharam? foram bem?
Acho que gostaram (preciso perguntar!). No geral, foram bem sim.. O tempo foi muito curto, então eles não podiam enrolar. Projetei um timer na lousa e dei play.
Acho que agora eles entenderam que o Stata não é uma caixa preta que joga resultados.