Lista_7

IMPORTAÇÃO DE DADOS

library(readxl)

nome_arquivo <- 'wage1.xls'
if(nome_arquivo %in% dir(getwd())){
wages_file <- read_excel(nome_arquivo)
print("ARQUIVO IMPORTADO")}else{ print("ARQUIVO NÃO DISPONÍVEL")}

## [1] "ARQUIVO IMPORTADO"

REGRESSÃO LINEAR SIMPLES

library(stargazer)

reg1 <- lm(formula = wage~educ,data = wages_file)
tabela_reg1 <- stargazer(reg1, type = 'html',style = 'all')

, , , , , , , , , , , , , , , , , ,


	Dependent variable:

	wage

educ	0.541^***
	(0.053)
	t = 10.167
	p = 0.000
Constant	-0.905
	(0.685)
	t = -1.321
	p = 0.188

Observations	526
R²	0.165
Adjusted R²	0.163
Residual Std. Error	3.378 (df = 524)
F Statistic	103.363^*** (df = 1; 524) (p = 0.000)

Note:	p<0.1; p<0.05; p<0.01

VETORES DA REGRESSÃO

coef1 <- coef(reg1)
resid1 <- resid(reg1)
vals1 <- predict(reg1)

print('Coeficientes:')
print(coef1)
print('-------------------------------------------------')
print('Resíduos (amostra de 30):')
print(resid1[1:30])
print('-------------------------------------------------')
print('Valores (amostra de 30):')
print(vals1[1:30])

## [1] "Coeficientes:"
## (Intercept)        educ 
##  -0.9048517   0.5413593 
## [1] "-------------------------------------------------"
## [1] "Resíduos (amostra de 30):"
##           1           2           3           4           5           6 
## -1.95010017 -2.35145943 -2.05010017  2.57397760 -0.29145943  0.99310354 
##           7           8           9          10          11          12 
##  2.41038502 -0.59145943 -1.99145943  9.88174428 -1.50689646  1.99718131 
##          13          14          15          16          17          18 
##  3.17854057 -0.09145943 16.60854057  9.57310354  1.90854057  4.49718131 
##          19          20          21          22          23          24 
## -1.99145943 -1.09145943  1.28854057  2.88854057 -1.42689646 -5.06145943 
##          25          26          27          28          29          30 
##  0.94989983  1.80310354  0.02310354  4.74310354  5.28446280 -0.17602240 
## [1] "-------------------------------------------------"
## [1] "Valores (amostra de 30):"
##        1        2        3        4        5        6        7        8 
## 5.050100 5.591459 5.050100 3.426022 5.591459 7.756896 8.839615 5.591459 
##        9       10       11       12       13       14       15       16 
## 5.591459 8.298256 7.756896 6.132819 5.591459 5.591459 5.591459 7.756896 
##       17       18       19       20       21       22       23       24 
## 5.591459 6.132819 5.591459 5.591459 5.591459 5.591459 7.756896 5.591459 
##       25       26       27       28       29       30 
## 5.050100 7.756896 7.756896 7.756896 7.215537 3.426022

COEFICIENTES PERCENTUAIS E ELASTICIDADE

wages_filt <- wages_file[which(wages_file$educ != 0),]

reg_pct <- lm(formula = log(wage) ~ educ,data = wages_file)
reg_el <- lm(formula = log(wage) ~ log(educ),data = wages_filt)

print('Coeficientes percentuais:')
print(coef(reg_pct))
print('-------------------------------------------------')
print('Coeficientes em elasticidade:')
print(coef(reg_el))

## [1] "Coeficientes percentuais:"
## (Intercept)        educ 
##  0.58377267  0.08274437 
## [1] "-------------------------------------------------"
## [1] "Coeficientes em elasticidade:"
## (Intercept)   log(educ) 
##  -0.4446768   0.8252071

COMPARAÇÃO ENTRE AS TRÊS REGRESSÕES

reg2 <- lm(formula = wage~ 0 + educ ,data = wages_file,)
reg3 <- lm(formula = wage~ 1 ,data = wages_file,)
tabela3 <- stargazer(reg1,reg2,reg3, type='html',style='all',
          column.labels = c("completo", "pela origem","constante") )

, , , , , , , , , , , , , , , , , , , ,


	Dependent variable:

	wage
	completo	pela origem	constante
	(1)	(2)	(3)

educ	0.541^***	0.473^***
	(0.053)	(0.011)
	t = 10.167	t = 41.247
	p = 0.000	p = 0.000
Constant	-0.905		5.896^***
	(0.685)		(0.161)
	t = -1.321		t = 36.616
	p = 0.188		p = 0.000

Observations	526	526	526
R²	0.165	0.764	0.000
Adjusted R²	0.163	0.764	0.000
Residual Std. Error	3.378 (df = 524)	3.381 (df = 525)	3.693 (df = 525)
F Statistic	103.363^*** (df = 1; 524) (p = 0.000)	1,701.328^*** (df = 1; 525) (p = 0.000)

Note:	p<0.1; p<0.05; p<0.01

par(bg = '#bad9ff')

plot(wages_file$educ,
     wages_file$wage,
     xlab = 'educ',
     ylab = 'wage',
     cex=0.5)

abline(reg1,col = 'blue',lwd=2,lty = 1)
abline(reg2,col = 'red',lwd=2,lty=2)
abline(reg3,col = 'dark green',lwd=2,lty=3)
legend("topleft",
       c("completa","pela origem","constante"),
       lwd=2,lty=1:3, 
       col = c('blue','red','dark green'))

AMOSTRAS ALEATÓRIAS

dens_medias <- function(n1,n2,n3){
medias1 <- numeric(1000)
medias2 <- numeric(1000)
medias3 <- numeric(1000)

for (i in 1:1000){
  medias1[i] <- mean(rnorm(n1,0,1))
  medias2[i] <- mean(rnorm(n2,0,1))
  medias3[i] <- mean(rnorm(n3,0,1))
}

par(new)

plot(density(medias3),col = 'red', main = paste('Densidade das médias de n amostras (média = 0, dp = 1)'))
lines(density(medias2), col = 'blue')
lines(density(medias1), col = 'dark green')
abline(v=0,lty=2)
legend('topleft', legend = c(n3,n2,n1),lwd=2, col = c('red','blue','dark green'))

}



dens_medias(10,50,100)

Verificamos portanto que quanto maior o tamanho da amostra na simulação, mais próximo da média populacional tenderá a observação

REGRESSÃO LINEAR MÚLTIPLA

library(readxl)
library(stargazer)

nome_arquivo <- 'wage1.xls'
if(nome_arquivo %in% dir(getwd())){
wages_file <- read_excel(nome_arquivo)
print("ARQUIVO IMPORTADO")}else{ print("ARQUIVO NÃO DISPONÍVEL")}

## [1] "ARQUIVO IMPORTADO"

reg <- lm(formula = wage ~ educ + exper,data = wages_file)
tabela_reg <- stargazer(reg, type = 'html',style = 'all')

, , , , , , , , , , , , , , , , , , , , , ,


	Dependent variable:

	wage

educ	0.644^***
	(0.054)
	t = 11.974
	p = 0.000
exper	0.070^***
	(0.011)
	t = 6.385
	p = 0.000
Constant	-3.391^***
	(0.767)
	t = -4.423
	p = 0.00002

Observations	526
R²	0.225
Adjusted R²	0.222
Residual Std. Error	3.257 (df = 523)
F Statistic	75.990^*** (df = 2; 523) (p = 0.000)

Note:	p<0.1; p<0.05; p<0.01

Observamos que o coeficiente para a variável educação e experiência são ambos positivos, sendo significativamente maior para educação. Isso sugere que um ano adicional de educação tem um efeito maior sobre o salário do que um ano adicional de experiência, em média.

Além disso, para ambos os coeficientes o valor zero (hipótese nula) está bastante distante, considerando os intervalos de confiança. Confirmamos isso verificando que o p-valor é baixo (p < 0.05), indicando que as correlações positivas para ambas as variáveis são significativas, podemos rejeitar a hipótese nula.

Nota-se também que o valor de R² é bastante baixo (mesmo o R² ajustado), indicando que essas duas variáveis independentes apenas não são suficientes para explicar satisfatoriamente toda a variação do salário, explicam apenas uma parte (22%).

SIMULAÇÃO MONTE CARLO -DISTRIBUIÇÃO NORMAL

dens_medias_normal <- function(n1,n2,n3,n4){
N_SORTEIOS <- 1000
medias1 <- numeric(N_SORTEIOS)
medias2 <- numeric(N_SORTEIOS)
medias3 <- numeric(N_SORTEIOS)
medias4 <- numeric(N_SORTEIOS)

for (i in 1:N_SORTEIOS){
  medias1[i] <- mean(rnorm(n1,0,1))
  medias2[i] <- mean(rnorm(n2,0,1))
  medias3[i] <- mean(rnorm(n3,0,1))
  medias4[i] <- mean(rnorm(n4,0,1))

}

par(new)

plot(density(medias4), xlim = c(-1,1),col = 'red', main = paste('Densidade das médias de n amostras - normal (média = 0, d.p. = 1)'))
lines(density(medias3), col = 'blue')
lines(density(medias2), col = 'dark green')
lines(density(medias1), col = 'orange')


abline(v=0,lty=2)
legend('topleft', legend = c(n4,n3,n2,n1),lwd=2, col = c('red','blue','dark green','orange'))

}



dens_medias_normal(10,100,1000,10000)

Verificamos portanto que a distribuição das médias das amostras aleatórias de distribuição normal, conforme era esperado pelo teorema do limite central, se aproxima do formato de uma distribuição normal, sendo que quanto maior o número de amostras de cada medida, menor a variância observada.

SIMULAÇÃO MONTE CARLO -DISTRIBUIÇÃO QUI-QUADRADO

dens_medias_chisq <- function(n1,n2,n3,n4){
N_SORTEIOS <- 1000
medias1 <- numeric(N_SORTEIOS)
medias2 <- numeric(N_SORTEIOS)
medias3 <- numeric(N_SORTEIOS)
medias4 <- numeric(N_SORTEIOS)

for (i in 1:N_SORTEIOS){
  medias1[i] <- mean(rchisq(n1,5,1))
  medias2[i] <- mean(rchisq(n2,5,1))
  medias3[i] <- mean(rchisq(n3,5,1))
  medias4[i] <- mean(rchisq(n4,5,1))

}

par(new)

plot(density(medias4), xlim = c(5,7),col = 'red', main = paste('Densidade das médias de n amostras - qui-quadrado (df = 5, ncp = 1)'))
lines(density(medias3), col = 'blue')
lines(density(medias2), col = 'dark green')
lines(density(medias1), col = 'orange')


abline(v=6,lty=2)
legend('topleft', legend = c(n4,n3,n2,n1),lwd=2, col = c('red','blue','dark green','orange'))

}



dens_medias_chisq(10,100,1000,10000)

Verificamos, como no exemplo anterior, que a distribuição das médias também se aproxima de uma distribuição normal, mesmo que a distribuição chi-quadrado não seja igual à normal, conforme prevê o teorema do limite central. Além disso, também como no exemplo anterior, quanto maior o número de amostras da observação, menor a variância.

MONTE CARLO - REGRESSÃO LINEAR

library(lambda.r)

N_SIMUL <- 30
N_AMOSTRA <- 1000

b0_estimado <- numeric(N_SIMUL)
b1_estimado <- numeric(N_SIMUL)

#Valores reais:
b0 <- 5
b1 <- 0.25
v <- 2

x <- rnorm(N_AMOSTRA,2,1)


for (i in 1:N_SIMUL){
  u <- rnorm(N_AMOSTRA,0,v)
  y <- b0 + b1*x + u

  estim <- coefficients(lm(y ~ x))
  b0_estimado[i] <- estim[1]
  b1_estimado[i] <- estim[2]
  }

par(bg = '#d4e5ff')
plot(density(b1_estimado), main = 'Distribuição dos b1 estimados', ylim = c(0,8))
abline(v = mean(b1_estimado),lty=2)
abline(v = b1,lwd = 3)
legend('topleft', 
       legend = c('B1 real','B1 estimado (média)'),
       lty = c(1,2))

GRÁFICOS DAS SIMULAÇÕES

par(bg='#f7d094')

plot(NULL,xlim = c(0,10),ylim = c(4,8), xlab = 'VAR IND',ylab = 'VAR DEP')

for (i in 1:N_SIMUL)
  {
  abline(b0_estimado[i],b1_estimado[i], col = 'orange',lwd=1)
  
}
abline(a = b0,b = b1,col='blue',lwd=3)

legend('topleft', 
       legend = c('Real','Estimado'),
       col = c('blue','orange'),
       lwd = 3)

VIOLAÇÃO DOS RESULTADOS - HIPÓTESE 4

b0_estimado <- numeric(N_SIMUL)
b1_estimado <- numeric(N_SIMUL)

for (i in 1:N_SIMUL){
  #Valor esperado de X não é 0 e depende de X
  u <- rnorm(N_AMOSTRA,x/4,v)
  y <- b0 + b1*x + u

  estim <- coefficients(lm(y ~ x))
  b0_estimado[i] <- estim[1]
  b1_estimado[i] <- estim[2]
  }

par(bg = '#d4e5ff')
plot(density(b1_estimado), 
     main = 'Distribuição dos b1 estimados', 
     ylim = c(0,8), 
     xlim = c(0,1))
abline(v = mean(b1_estimado),lty=2)
abline(v = b1,lwd = 3)
legend('topleft', 
       legend = c('B1 real','B1 estimado (média)'),
       lty = c(1,2))

par(bg='#f7d094')

plot(NULL,xlim = c(0,10),ylim = c(4,8), xlab = 'VAR IND',ylab = 'VAR DEP')

for (i in 1:N_SIMUL)
  {
  abline(b0_estimado[i],b1_estimado[i], col = 'orange',lwd=1)
  
}
abline(a = b0,b = b1,col='blue',lwd=3)

legend('topleft', 
       legend = c('Real','Estimado'),
       col = c('blue','orange'),
       lwd = 3)

VIOLAÇÃO DOS RESULTADOS - HIPÓTESE 5

b0_estimado <- numeric(N_SIMUL)
b1_estimado <- numeric(N_SIMUL)

for (i in 1:N_SIMUL){
  #Variância do erro não é fixa, depende de X (heteroscedasticidade)
  u_var <- ((x + 0.1) * 1.5)
  
  u <- rnorm(N_AMOSTRA,0,u_var)
  y <- b0 + b1*x + u

  estim <- coefficients(lm(y ~ x))
  b0_estimado[i] <- estim[1]
  b1_estimado[i] <- estim[2]
  }

par(bg = '#d4e5ff')
plot(density(b1_estimado), 
     main = 'Distribuição dos b1 estimados', 
     ylim = c(0,8),
     xlim = c(0,1))
abline(v = mean(b1_estimado),lty=2)
abline(v = b1,lwd = 3)
legend('topleft', 
       legend = c('B1 real','B1 estimado (média)'),
       lty = c(1,2))

par(bg='#f7d094')

plot(NULL,xlim = c(0,10),ylim = c(4,8), xlab = 'VAR IND',ylab = 'VAR DEP')

for (i in 1:N_SIMUL)
  {
  abline(b0_estimado[i],b1_estimado[i], col = 'orange',lwd=1)
  
}
abline(a = b0,b = b1,col='blue',lwd=3)

legend('topleft', 
       legend = c('Real','Estimado'),
       col = c('blue','orange'),
       lwd = 3)

Observamos que a violação de hipótese 4 tem efeito mais significativo do que a violação da hipótese 5.

INTERVALOS DE CONFIANÇA E TESTES T

#Using R for Introductory Econometrics - Florian Heiss
#1.9.3 Simulation of confidence intervals and t tests

N_SIMUL <- 50
N_AMOSTRA <- 100
mu <- 10
s <- 2

cilower <- numeric(N_SIMUL) ; ciupper <- numeric(N_SIMUL)
pvalue1 <- numeric(N_SIMUL) ; pvalue2 <- numeric(N_SIMUL)
reject1 <- logical(N_SIMUL) ; reject2 <- logical(N_SIMUL)
color1 <- character(N_SIMUL);color2 <- character(N_SIMUL)

for (i in 1:N_SIMUL){
  amostra <- rnorm(N_AMOSTRA,mu,s)
  
  teste1 <- t.test(amostra,mu = 10, alternative = 'two.sided')
  cilower[i] <- teste1$conf.int[1]
  ciupper[i] <- teste1$conf.int[2]
  reject1[i] <- teste1$p.value < 0.05
  
  teste2 <- t.test(amostra,mu = 9.5, alternative = 'two.sided')
  reject2[i] <- teste2$p.value < 0.05
 
  }

mean(reject1)
table(reject1)
mean(reject2)
table(reject2)


for (i in 1:N_SIMUL){
  color1[i] <- if (reject1[i]==TRUE){'red'}else{'blue'}
  color2[i] <- if (reject2[i]==TRUE){'red'}else{'blue'}
}

par(mfrow=c(1,2))


barplot(ciupper-cilower,
        offset = cilower,
        space = 1,
        col = color1,
        horiz=TRUE,
        xlim = c(9,11),
        main = 'Média = 10',
        )
abline(v=10, lty=2)

barplot(ciupper-cilower,
        offset = cilower,
        space = 1,
        col = color2,
        main = 'Média = 9.5',
        xlim = c(9,11),
        horiz=TRUE)
abline(v=9.5, , lty=2)

A variável aleatória possui média 10 e desvio-padrão 2. Quando testamos a hipótese de que a média é igual a 10, praticamente todos os intervalos de confiança simulados contêm a média 10 - logo, corretamente aceitamos a hipótese. No entanto, ao propor a hipótese de que a média é igual a 9.5 para os mesmos dados, menos de 95% dos intervalos não incluem 9.5, tal que podemos rejeitar esta hipótese.

Capítulo 2

Caio_Donalisio

15/10/2020