###################
# AULA INFERENCIA #
###################


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

###################
# BASE WOOLDRIDGE #
###################

# 1) Definir diretório

getwd()

setwd("C:/Users/Wander Plassa/Desktop/Pos_Graduacao/Doutorado/R/Aula/Stata")

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

library(foreign)


# 3) Obter bases do Livro texto (Wooldridge)

# http://www.cengage.com/cgi-wadsworth/course_products_wp.pl?fid=M20b&product_isbn_issn=9780324581621&discipline_number=413


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

CEOSAL = read.dta("CEOSAL2.dta")

CEOSALTXT = read.table("CEOSAL2.raw") 

## caso de erro de versão mude a versão do arquivo stata para uma mais antiga:
## use o comando "saveold nomedonovobanco" no stata
## Agora é só abrir o banco de dados modificado da seguinte forma:

CEOSAL <- read.dta("CEOSAL2OLD.dta")

## Caso queira salvar o banco no formato do R use o seguinte comando:

save(CEOSAL,file="CEOSAL2OLD.Rda")

## Para carregar esse banco de dados no novo formato em uma nova ocasião basta usar:

load("CEOSAL2OLD.Rda")

# 4) Realizando regressões, em que o salário é a variável dependente.

# A primeira é uma regressão simples, com apenas uma variável independente (vendas).

reg1 = lm(log(CEOSAL$salary) ~ log(CEOSAL$sales), data=CEOSAL)

summary(reg1)


## A segunda é a regressão adicionando a variável "mktval" 
## Que é o valor de mercado da firma e "profmarg" que é o lucro como porcentagem das vendas

reg2 = lm(log(CEOSAL$salary) ~ log(CEOSAL$sales) + log(CEOSAL$mktval) + CEOSAL$profmarg, data=CEOSAL)

summary(reg2)

## A terceira é a regressão completa com "ceoten" que é anos como CEO na firma
## e "comten" é o total de anos com a firma

reg3 = lm(log(CEOSAL$salary) ~ log(CEOSAL$sales) + log(CEOSAL$mktval) + CEOSAL$profmarg +
CEOSAL$ceoten + CEOSAL$comten, data=CEOSAL)

summary(reg3)

######################
##### TESTE T ########
######################

# Teste para UM [UNICO] parâmetro populacional

# 1) Calculo teste de T

# No resultado da regressão o teste de t já aparece para cada parâmetro,
# contudo caso queria calcular manualmente basta fazer a seguinte conta

#  t = â - a / sd(â)

# onde "â" é o parametro estimado, "a" o valor a ser testado e sd(â) é
# o desvio padrão estimado do parâmetro estimado. Por exemplo, para Log(sales)
# se o meu objetivo é testar a hipótese nula de que H0: a = 0, podemos fazer
# a seguinte conta:

(0.187787 - 0)/0.040003

# Valores de t calculado suficientemente longe de zero resultará na rejeição de H0.
# Isto é, rejeitará a hipótese que a variável independente "log(sales)" não tem efeito sobre a
# variável depedente em questão. Mas quão distante de zero ele deve estar?
) Nivel de significância:
# 2

# Para tal, precisamos definir níveis de significância,
# Isto é, a probabilidade de reijeitar a H0 quando ela é verdadeira.
# Portanto, quanto menor o nível de significância mais robustos
# serão resultados obtidos, uma vez que a chance de falarmos
# que a covariada tem efeito sobre a variável dependente
# quando na verdade ela não tem, se torna pequena.

# Valores comuns do nível de significância: 1%, 5% e 10%.
# Com valor t calculado em mãos e o nível de significância,
# podemos encontrar o t tabelado (t com distribuição n-k-1).
# Onde n é o tamanho da amostra e k o número de parâmetros na regressão.
# E -1 refere-se ao intercepto

# No nosso caso, a 5%, considerando uma comparação bilateral (usado quando
# o sinal de "a" não é bem determinado pela teoria - mais prudente), teremos
# um t tabelado com 177 - 5 - 1 (s graus de liberdade) = 1,96.
# Isto é, t calculado  acima de 1,96, nesse exemplo, nos indica que,
# com um nível de significância a 5%, as covariadas afetam o log(salário)

# Por exemplo, o t calculado da covariada "profmarg" é de -1.050. Em módulo,
# ela não ultrapassa o valor de 1,96. Portanto não podemos rejeitar a 
# hipótese nula a 5%.

# 3) Testando outras Hipóteses sobre "a"

# Podemos testar outros valores específicos para "a" além de zero.
# Por exemplo, fazer uma H0 : "a" = 1. Para isso, basta substituir
# 1 no lugar de "a" e proceder da mesma forma, escolhendo os 
# níves de significância.

# (â - 1)/sd(â)

# Por exemplo, queremos testar se a variável "ceoten" é estatisticamente
# diferente de 1 a 5%.

(-0.009238 - 1)/0.003337

# O valor em módulo é 302.4387. Bem acima de 1.96 (valor de t tabelado a 5%).
# Podemos rejeitar a H0 que "comten" é igual a 1.


# 4) Computando o valor de "p" ou p-valor para o teste de t

# outro valor que o programa R já nos apresenta no sumário da regressão
# é o p-value. Ele contorna o problema de arbitrariedade
# que fazemos ao escolher o valor de signficância.

# Ele nos responde a seguinte pergunta: Dado o valor observado da estatística t, 
# qual o menor nível de significância para o qual a hipótese nula seria rejeitado?

# No nosso exemplo, para a variável "profmarg", para a hipótese nula (H0: a = 0) 
# ser rejeitada, deveríamos aceitar um nível de significância maior do que 29%.

# Isso significa que os valores menores de p são evidência contra H0; 
# p-valores grandes oferece pouca evidência contra H0.

# Ele é um teste bilateral que testa a hipótese do coeficiente ser diferente de 0.

# 5) Intervalo de Confiança

# fornecem uma gama de valores prováveis para o parâmetro populacional, 
# e não apenas uma estimativa pontual.

# O calculo do intervalo se dá pela seguinte conta:
# â +/- (t tabelado*sd(â))

# Onde â é o parametro estimado, t tabelado e o valor t (n-k-1) obtido na
# tabela de t para o nível de signifcância desejado e bilateral.
# Para a variável "log(sales)" do nosso exemplo, com 95% de confiança (ou 5% de sig), 
# o intervalo seria:


# Para 95% de confiança
0.187787 + (1.96*0.040003) 

0.187787 - (1.96*0.040003)

# Para 99% de confiança
0.187787 + (2.617*0.040003)

0.187787 - (2.617*0.040003)


# O intervalor seria ->  0.1085811 e 0.2669929 (zero está fora do intervalo)

# Há uma chance de 95% de que o parâmetro verdadeiro de "sales"
# esteja nesse intervalo (não há uma garantia).

# Para ver todos os intervalos, basta usar o seguinte comando:
confint(reg3, level = 0.95)
confint(reg3, level = 0.99)

##### TESTE F ########

# testar hipóteses que envolvem mais do que um parâmetro populacional

# Testaremos várias restrições

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

# HA é que H0 é falsa

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

# Onde SQR é a soma dos Quadrados residuais obtido na ANOVA
# SQRm é referente a soma dos quadrados residuais da regressão
# sem o conjunto de covariadas que queremos testar
# q é o número de variáveis retiradas do modelo

# Exemplo: suponha que queremos testar que "log(sales)" e "log(mktval)",
# tenham efeito conjunto sobre a nossa variável dependente "log(salary)".

# Para isso a H0: log(CEOSAL$sales) = 0, log(CEOSAL$mktval) = 0

# Controlando pelas outras covariadas

summary(reg3)
anova(reg3)
resid<-(reg3$residuals)
resid2<-resid*resid
SQRnamao<-sum(resid2)
summary(SQRnamao)

# SQR = 41.856

# A regressão menor seria:

regmenor3 = lm(log(CEOSAL$salary) ~  CEOSAL$profmarg +
            CEOSAL$ceoten + CEOSAL$comten, data=CEOSAL)
summary(regmenor3)
anova(regmenor3)

residmenor <-(regmenor3$residuals)

residmenor2<-residmenor*residmenor

SQRnamaomenor<-sum(residmenor2)

summary(SQRnamaomenor)

# SQRm = 63.427

# Portanto o F calculado seria:
# F calculado = [(63.427 - 41.856)/2] / [41.856/(177 - 5 - 1)]

((SQRnamaomenor - SQRnamao)/2) /( SQRnamao/(177 - 5 - 1))

# F calculado = 44.06347

# Ao analisar a tabela do teste de F, considerando os
# graus de liberdade, acima de 120, com 5% de significância
# e um q=2, obtermos o valor 3.0718.

# Como F calculado > F tabelado, podemos rejeitar a H0
# de que as variáveis retiradas do modelo não possuem
# efeito sobre a variáveil dependente


###################
# BASE PNAD2013 #
###################

# 1) chamar a base da PNAD em formato dta pelo comando read.dta

STATA <- read.dta("PNAD2013old.dta")

# separando as variáveis de interesse 

Sal = data.frame(STATA$V0302, STATA$V8005, STATA$V0404, STATA$V4011, STATA$V4803, STATA$V4718, STATA$V9891, STATA$V9891*STATA$V9891)

# Considerando apenas pessoas com renda positiva e menor de 100.000,00 

Salexist = subset(Sal,  (STATA.V4718>0) & (STATA.V4718<100000))

# Renomeando as variaveis

colnames(Salexist) =  c("Sexo", "Idade", "Cor", "EstCivil", "AnoEstud", "renda", "experiência", "experiência2")

logrenda = log(Salexist$renda)

# regressão: EQUAÇÃO DE SALÁRIOS MINCERIANA 


regsal = lm( logrenda ~  AnoEstud +
               experiência + experiência2 + Cor + EstCivil + Sexo  , data=Salexist)
summary(regsal)

# Todos os regressores foram estatisticamente significativos pelo teste de T e pelo
# p-valor.


# Pegando uma amostra aleatória de 500 pessoas 

salexistM = Salexist[sample(nrow(Salexist), 500), ]

# Rodando uma regressão nessa amostra menor

regsalm = lm( log(renda) ~  AnoEstud +
               experiência + experiência2 + Cor + EstCivil + Sexo  , data=salexistM)
summary(regsalm)

# Em uma amostra menor, no entanto, tanto experiência como experiência ao quadrado
# Não foram estatisticamente significativas a 10%, por exemplo.