#######################
# AULA AMOSTRAS GRANDES #
#######################


# Para deletar tudo que está aberto
rm(list=ls())

# Definir diretório

getwd()

setwd("C:/Users/Wander Plassa/Desktop/Pos_Graduacao/Doutorado/R/Aula/Multi regressão")

# chamar o comando "foreign" para obter bases de dados de outros programas, por exemplo Stata

library(foreign)

# chamar a base em formato dta pelo comando read.dta ou TXT (formato raw)


#######################
# Introdução ##########
#######################

# Considerando os pressuposto de Gauss-Markov.

# Até agora lidavamos com amostras finitas ou pequenas amostras.
# No capítulo 4 foi adicionado uma sexta propriedade para se fazer inferência
# A normalidade: Afirma que o termo de erro u é normalmente distribuído e 
# independente das variáveis explicativas.

# Se o erro não é normalmente distribuído em amostras pequenas a distribuição de uma estatística t 
# Não é exatamente t, e uma estatística F não tem uma distribuição exata  F 
# para qualquer tamanho de amostra.

# Portanto é importante conhecer as propriedades assintóticas
# Ou propriedades de grandes amostras de estimadores e estatísticas de teste.


# Um achado importante é que mesmo sem a  Hipótese de normalidade
# (Pressuposto 6), estatísticas t e F têm aproximadamente
# distribuições t e F, pelo menos em  amostra grandes.


#######################
# Consistência #
#######################

# Nem sempre não viés é alcançado, contudo consistência é um requisito 
# mínimo que um estimado deve ter.

# A consistência implica que o plim â >>> a, isto é,
# O valor estimado tende em probabilidade para o valor verdadeiro
# conforme o n aumenta.

# Isso significa que podemos fazer nosso estimador
# Arbitrariamente perto de "a" se pudermos aumentar nossa amostra

# Naturalmente, no mundo real temos amostras fixas,
# Consistência envolve em: o que aconteceria com nosso estimador
# se a amostra fosse maior.

#######################
# Distribuição Normal #
#######################

# Além de consistência, é necessário que as distribuições amostrais,
# sejam normalmente distribuídas para fazer inferência.

# Uma variável aleatória normalmente distribuída é distribuída simetricamente
# em torno de sua média, e mais de 95% da área sob a distribuição está
# dentro de dois desvios=padrão.


# Analisando o exemplo 4.6 do Livro texto,
# Verificamos em primeiro lugar a normalidade da variável "prates"
# que é a porcentagem de participação em planos de pensão do tipo 401(k) , que se 
# Configura na variável dependente


Plansd = read.dta("401Kold.dta")   

# Como não observamos a variável u, via de regra focamos na normalidade
# da variável dependente

# Realizando um histograma dessa variável ou frequência de distribuição 
# Observamos que ela Fortemente concentrada na direita, em vez de estar normalmente
# distribuída

hist(Plansd$prate)


# Cerca de 40% dos valores estão concentrados no valor 100

# A normalidade é importante para fazermos inferência baseadas
# Nas estatísticas de t e F. Sem essa premissa (6 pressuposto), 
# A análise de inferência fica comprometida.


# Mesmo que o a var. dependente não tenha uma distribuição normal, podemos usar o Teorema Central do Limite (TCL)
# para concluir que os estimadores de MQO satisfazem ou não as hipóteses Normalidade, 
# o que significa que eles são aproximadamente normalmente distribuídos em amostra grandes o suficiente.


#### Teste t ####

# Lembramos que o t calculado é obtido a partir da seguinte formula:

# t = (â - a)/EP (â) ~ assintoticamente para t (n-k-1)

# A partir de uma amostra grande o suficiente de indivíduos podemos usar a apromixação 
# do teorema central do limite Usando a distribuição normal padrão. Esse termo "grande suficiente"
# não é bem definido. Além disso, a qualidade da aproximação não depende
# Apenas do n mas também dos graus de liberdade Na tabela do teste de t, com 
# graus de liberdade acima de 120 já podemos usar a normal padrão.

# Se o tamanho da amostra não é muito grande então a distribuição t pode ser uma má aproximação
# Para a distribuição das estatísticas t quando u não é normalmente distribuído.

# Na presença de heteroscedasticidade o teorema central do limite não é válido,
# Não imporando o tamanho da amostra. Mas essa é uma discussão para o capítulo 8.

## ERRO PADRÃO ASSINTÓTICO ##

# Pode esperar que os erros padrão diminuam a uma taxa  que é o inverso 
# da raiz quadrada do tamanho da amostra.

regplan = lm( Plansd$prate ~ Plansd$totpart + Plansd$totelg + Plansd$totemp, data=Plansd)

summary(regplan)

co = coef(summary(regplan))

Plansdmenor =  Plansd[sample(nrow(Plansd), 750), ]

regplanmenor = lm(Plansdmenor$prate ~ Plansdmenor$totpart + Plansdmenor$totelg + Plansdmenor$totemp, data=Plansdmenor)

summary(regplanmenor)

co_m = coef(summary(regplanmenor))

co
co_m

# Dividindo o erro padrão de um estimador qualquer da reg com n maior
# e n menor

3.946683/6.738165

# Obtemos 0.58. O livro texto mostra que esse calculo é aproximadamente
# igual a dividir n menor por n maior e tirar a raiz.

sqrt(750/1534)

# O erro padrão do estimador "totpart", ao usar a amostra maior deve ser cerca
# de 70% do erro padrão obtido ao usar a amostra menor.


#### Teste LM ####

# A forma da estatística de LM que derivamos aqui se baseia nas premissas de Gauss-Markov,
# Os mesmos pressupostos que justificam a estatística F em grandes amostras. Não precisamos
# do pressuposto da normalidade.

# Gostaríamos de testar se q variáveis têm parâmetros 
# populacionais iguais a zero: a hipótese nula é

# a1 = 0 e a2 = 0 ....

# A hipótese alternativa é que pelo menos um desses parâmetros
# é diferente de zero

# Como faz o teste de F

# Como exemplo pegamos a base de crimes

Crime = read.dta("CRIME1OLD.dta") 

# Iremos testar se duas variáveis (portanto q=2) são irrelevantes no modelo:
# "avgsen" e "tottime".

# PASSO A PASSO

# Primeiro rodamos a regressão sem essas duas variáveis (regressão restrita) e coletamos o resíduo

regcrime = lm( Crime$narr86 ~ Crime$pcnv + Crime$ptime86 + Crime$qemp86, data=Crime)

residcrime = residuals(regcrime)

# A intuição é se essas variáveis tem coeficientes populacionais iguais a zero,
# O resíduo da regressão restrita não deve ser correlacionada com as variáveis
# que permaneceram no modelo.

# Uma etapa seguinte é rodar uma regressão dos resíduos da regressão restrita
# contra todas as variáveis do modelo (isto é, com as variáveis excluídas também)

regrescrime = lm(residcrime ~ Crime$pcnv + Crime$ptime86 + Crime$qemp86
                 + Crime$avgsen + Crime$tottime, data=Crime)

summary(regrescrime)


# O R² do resíduo contra as variáveis indenpendentes (todas elas) foi de 0.0015.
# Ao multiplicar por n (tamanho da amostra) para obter a estatística de LM resulta em 

LMsta = 0.0015*2725

qchisq(.870458, df=2)

1 - pchisq(LMsta,2) 


# Que nos dá um valor de 4.08. Olhando os valores críticos da distribuição 
# Qui-quadrada com 2 graus de liberdade (número de restrições sendo testadas)
# e com 10% de significância é 4.61. Portanto como calculado < tabelado,
# Não conseguimos rejeitar a H0 de que os parâmetros propulacionais retirados
# são iguais a zero a 10%.

# A diferença do teste de F que não utilizamos os graus de liberdade 
# da regressão não restrita.

# Se fizermos o teste de F

# Que lembrando:

# H0 é que um conjunto de variáveis não tem nenhum efeito sobre a var. dependente, 
# Por exemplo H0: a1 =0, a2 = 0.

# com HA:  H0 é falsa

# A fórmula do teste de F é dado por:
# F = [(SQRm - SQR)/q]/[SQR/(n-k-1)]


# Para isso a H0: "avgsen"  = 0 ou "tottime" = 0

regcrma = lm( residcrime ~ Crime$pcnv + Crime$ptime86 + Crime$qemp86 + Crime$avgsen + Crime$tottime, data=Crime)

summary(regcrma)

residcrmaior <-(regcrma$residuals)

residcrmaior_sqr<-residcrmaior*residcrmaior

SQRmaior<-sum(residcrmaior_sqr)

summary(SQRmaior)



regcrme = lm( Crime$narr86 ~ Crime$pcnv + Crime$ptime86 + Crime$qemp86, data=Crime)

summary(regcrme)

residcrmenor <-(regcrme$residuals)

residcrmenor_sqr<-residcrmenor*residcrmenor

SQRmenor<-sum(residcrmenor_sqr)

summary(SQRmenor)

((SQRmenor - SQRmaior)/2) /( SQRmaior/(2725 - 6 - 1))

# F calculado = 2.03, F tabelado com gl>120 e q=2 = 2.30
# Portanto F cal< F tab. Mesma direção do resultado LM
# de que não conseguimos rejeitar H0 a 10%

#######################
##### EXEMPLOS ########
#######################

# 1) Exemplo salários

# Rodamos uma regressão multipla

Wage = read.dta("WAGE1OLD.dta")   

reg1 = lm( Wage$wage ~ Wage$educ + Wage$exper + Wage$tenure, data=Wage)

summary(reg1)

# Salvamos os resíduos e apresentamos em forma de Histograma.

resid1 = residuals(reg1)

hist(resid1, breaks = 100)


hist(resid1, prob=TRUE , breaks = 50)
curve(dnorm(x, mean=mean(resid1), sd=sd(resid1)), col="darkblue", lwd=2, add=TRUE, yaxt="n")


# Agora rodamos uma regressão com o salário em formato de log

reg2 = lm( log(Wage$wage) ~ Wage$educ + Wage$exper + Wage$tenure, data=Wage)

summary(reg2)

resid2 = residuals(reg2)

# O histograma do resíduo agora fica da seguinte forma

hist(resid2, breaks = 100)

hist(resid2, prob=TRUE , breaks = 50)
curve(dnorm(x, mean=mean(resid2), sd=sd(resid2)), col="darkblue", lwd=2, add=TRUE, yaxt="n")


# Aparentemente os resíduos da regressão considerando o log do salário estão mais normalmente,
# distribuídos quando comparado aos resíduos da regressão do salário em nível.
# Os resíduos da regressão do salário em nível estão mais concentrados à esquerda.

#######################
# TRANSFORMAÇÕES #
#######################

# Transformações de variáveis ao quadrado.

# Por exemplo, temos novamente a equação de salários,
# em que queremos ver o efeito de experiência e experiência ao quadrado.

regwage = lm( Wage$wage ~  Wage$exper + Wage$expersq, data=Wage)

summary(regwage)

0.2981001 / ( 2*(0.0061299))

# Como interpretar? Um ano adicional de experiência equivale a 30 centavos
# no salário hora. Mas, conforme os anos passam esse aumento é decrescente
# Se no primeiro ano é 30, no segundo reduz 0,6 (x) 2 centavos. E assim por diante.

regwage1 = lm( Wage$wage ~  Wage$exper + Wage$expersq + Wage$educ, data=Wage)

summary(regwage1)

0.268287 / (2*(0.004612))

# Adicionando educação, evitamos um possível viés. Verificamos que tanto o coeficiente 
# de experiência como experiência ao quadrado se reduzem. Antes com 24 anos o retorno
# da experiência se tornava negativo sobre o salário hora, agora


# Fim