############################################################################################
######################### Monitoria 04 de Econometria I ####################################
############################# Data: 10/10/2018 #############################################
#################### Curso: Matemática Aplicada a Negócios #################################
########################## Monitora: Laura Ogando ##########################################
######################## e-mail: lauraogando@usp.br ########################################
#################### Professor: Luiz Guilherme Scorzafave ##################################
############################################################################################

# Nesse script daremos continuidade a leitura da PNAD contínua e trabalheremos um pouco com regres-
# são múltipla.
# Para limpar todo o seu ambiente utiliza-se o comando abaixo:
rm(list = ls())

# Precisaremos instalar alguns pacotes para ler os dados
install.packages('SAScii')
install.packages('readr')

# Carregando os pacotes a serem utilizados
library('SAScii')
library('readr')

# Definindo o diretorio de trabalho onde se encontram os microdrados da PNAD e o dicionarios em SAS
# Você deve colocar o diretório que está definido no seu computador entre aspas. Não se esqueça de 
# inverter as barras.
setwd("C:/Users/Laura/Desktop/Mestrado/2. PAE - Econometria I/Monitorias de R/Monitoria 4")

# Os arquivos necessários já se encontram no e-disciplinas, mas caso um dia você precise baixar outros dados
# da PNAD continua, você pode utilizar os links a seguir.
# Para baixar os dados https://downloads.ibge.gov.br/downloads_estatisticas.htm?caminho=Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_continua/Anual/Microdados/Dados
# Dicionários https://downloads.ibge.gov.br/downloads_estatisticas.htm?caminho=Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_continua/Anual/Microdados/Documentacao

# Nesse script vamos usar somente os dados da primeira visita de 2016

# A função parse.SAScii toma como argumento o arquivo .txt que contem o dicionários para ler os dados em SAS
# Geralmente, o IBGE disponibiliza dicionários neste formato
# Lendo o dicionário:
dic16.1 <- parse.SAScii('Input_PNADC_1_visita_2016_20180816.txt')

# Os comandos abaixo irão montar as coordenadas de cada variavel dentro do arquivo de microdados  
# Basicamente, o parâmetro width, que já vem no dicionários, indica o tamanho da variável dentro da base de dados.
# Sendo assim, basta definir o ponto onde ela começa na base de dados e onde ela termina.
dic16.1$start <- cumsum(dic16.1$width)-dic16.1$width+1
dic16.1$end <- cumsum(dic16.1$width)

# Uma vez indicadas as coordenadas, basta ler os dados da PNAD através da função read_fwf, que 
# toma como argumentos o .txt com os dados e a posição das variáveis, definidas no dataframe do dicionário
PNADc_16.1 <- read_fwf('PNADC_2016_visita1_20180816.txt', 
                        fwf_positions(dic16.1$start, 
                                      dic16.1$end, 
                                      dic16.1$varname))

# Agora, vamos selecionar somente algumas variáveis de interesse.
# Por meio do arquivo em excel você pode escolher algumas variáveis que sejam interessantes para a sua pesquisa.
variaveis <-          c('ANO', 'TRIMESTRE', 'UF', 'CAPITAL', 
                        'RM_RIDE', 'UPA', 'ESTRATO', 'V1008', 
                        'V1014', 'V1022', 'V1023', 'V1030', 
                        'V1031', 'V1032', 'POSEST', 'V2001', 
                        'V2003', 'V2005', 'V2007', 'V2008', 
                        'V20081', 'V20082', 'V2009', 'V2010', 
                        'V3001', 'V3002', 'V3009A', 'V3014',
                        'V4012', 'V4019', 'V4029', 'V4048',
                        'VD2002', 'VD2003', 'VD2004', 'VD3001', 
                        'VD3002', 'VD4001', 'VD4002', 'VD4003', 
                        'VD4004A', 'VD4005', 'VD4007', 'VD4008', 
                        'VD4009', 'VD4010', 'VD4011', 'VD4012', 
                        'VD4013', 'VD4014', 'VD4015', 'VD4016', 
                        'VD4017', 'VD4018', 'VD4019', 'VD4020', 
                        'VD4031', 'VD4035', 'VD4036', 'VD4037',
                        'VD5005')

# Agora vamos manter somente essas variáveis que selecionamos na lista acima
PNADc_16 <- PNADc_16.1[variaveis]
  
# Para ficar mais organizado, podemos apagar as bases auxiliares e salvar a base da PNAD com o ano,
# Isso pode ajudar se você estiver trabalhando com PNADs contínuas de diferentes peíodos.
rm(dic16.1, PNADc_16.1, variaveis)

# Gerando o código identificador de pessoa
PNADc_16$PESSOA <- paste(PNADc_16$UPA, PNADc_16$V1008, PNADc_16$V1014, PNADc_16$V2003, sep="")

# Agora vamos supor que queremos retirar a parte categórica da variável VD3002, anos de estudo.
# Sendo assim, vamos retirar os valores de 15 anos os mais de estudo.
# Primeiro, vamos transformar essa variável em uma variável numérica
PNADc_16$VD3002 <- as.numeric(PNADc_16$VD3002)
summary(PNADc_16$VD3002)

# Agora, vamos retirar essa última categoria:
PNADc_16$VD3002[PNADc_16$VD3002 == 15] <- NA

# Note que o máximo deixa de ser 15, passando a ser 14
# Além disso, estamos atribuindo a essas variáveis que tem valor 15, o valor NA. Podemos notar que o
# número de NA's sobe:
summary(PNADc_16$VD3002)

# Vamos começar a analisar regressões múltiplas.
# O comando no R é o mesmo, apenas temos que incluir mais variáveis no nosso modelo
# Vamos começar criando um subset dos indivíduos da base que estaria na População Economicamente
# Ativa (PEA). Seriam pessoas com mais de 14 anos e menos de 65 anos. Além disso, iremos selecionar
# pessoas com rendimento do trabalho principal maior que zero.
# Vamos transformar as variáveis em numéricas
PNADc_16$V2009 <- as.numeric(PNADc_16$V2009)
PNADc_16$VD4017 <- as.numeric(PNADc_16$VD4017)

PEA <- subset(PNADc_16, subset = V2009>=14 & V2009<=65 & VD4017>0)

# Vamos começar a olhar pra equação minceriana. Portanto, vamos criar a variável ln do salário.
PEA$lnsal <- log(PEA$VD4017, exp(1))

# Criando uma proxy da experiência por meio da idade do indivíduo
PEA$exper <- PEA$V2009 - 14
summary(PEA$exper)
summary(PEA$lnsal)

# Fazendo um modelo de regressão linear simples só com anos de estudo
reg1 <- lm(lnsal ~ VD3002, PEA)
summary(reg1)

# Fazendo um modelo de regressão linear múltipla com anos de estudo e experiência
reg2 <-lm(lnsal ~ VD3002+exper, PEA)
summary(reg2)

# Criando uma variável de experiência ao quadrado
PEA$exper_2 <- PEA$exper^2
# Rodando a regressão com experiência ao quadrado
reg3 <- lm(VD4017 ~ VD3002+exper+exper_2, PEA)
summary(reg3)

# Podemos juntar todos os nossos resultados em uma única tabela:
library(stargazer)
stargazer(reg1, reg2, reg3, type = 'text', digits=2)