Category

Dicas de R

Modelos VAR no R

By | Dicas de R

Quem trabalha com modelagem e previsão provavelmente já teve que construir modelos baseados em vetores autorregressivos. Nessa Dicas de R - disponível toda quarta-feira aqui no blog da AM -vamos mostrar como é possível implementar esse tipo de modelo no R e utilizá-lo para fins de previsão. Essa aula, inclusive, faz parte do nosso novíssimo Curso de Previsão Macroeconométrica usando o R.

Em análises de séries temporais, é bastante comum considerar simultaneamente duas ou mais séries. Por exemplo, o nível de ociosidade da economia certamente tem influência sobre a taxa de crescimento dos preços, assim como um aumento da taxa de juros tem relação com a taxa de desemprego. Ademais, pode ser necessário em algumas ocasiões avaliar o quanto um choque em X_{t} afeta Y_{t}. Essas questões, bastante pertinentes no dia a dia, são tratadas dentro do que chamamos de análise multivariada de séries temporais. São objetivos básicos da análise multivariada de séries temporais:

  • Estudar as relações dinâmicas entre séries diversas;
  • Melhorar as previsões sobre uma variável específica.

Seja z_{t} = (z_{1t}, z_{2t}, ..., z_{kt})^{'} um vetor de dimensão k que contém séries temporais observadas em um período de tempo comum. Por exemplo, z_{1t} é o PIB trimestral brasileiro e z_{2t} é a taxa de desemprego também trimestral. Acaso, estudemos z_{1t} e z_{2t} conjuntamente poderemos verificar a dependência contenporânea e passada que existe entre essas duas variáveis. Poderemos verificar o quanto um choque no PIB afeta a taxa de desemprego e quanto tempo isso tende a durar. Analogamente, podemos estar interessados na relação entre a taxa de desemprego e a inflação.


desemprego = Quandl('BCB/24369', start_date='2012-03-01',
end_date='2018-12-01', order='asc', type='ts')
inflacao = Quandl('BCB/13522', start_date='2013-03-01',
end_date='2019-12-01', order='asc', type='ts')
data = data.frame(desemprego, inflacao)

A correlação entre as variáveis é -0,81, bastante elevada, mas isso, claro, não é o bastante para identificarmos a natureza da relação entre elas. Algumas questões devem ser levadas em consideração. Em primeiro lugar, para uma análise econométrica, precisamos definir se as séries são estacionárias. Essa primeira (grande) questão definirá o tipo de análise que faremos entre elas. Vamos supor que elas sejam estacionárias. Podemos, nesse caso, estimar um \textbf{Vetor Autoregressivo} de primeira ordem, de modo que:

(1)   \begin{align*} Unemployment_{t} = \delta_{1} + \gamma_{11}Unemployment_{t-1} + \gamma_{12}Inflation_{t-1} + \varepsilon_{1t} \\ Inflation_{t} = \delta_{2} + \gamma_{21}Inflation_{t-1} + \gamma_{22}Unemployment_{t-1} + \varepsilon_{2t} \end{align*}

O VAR(1) vai descrever a evolução dinâmica da interação entre a inflação e desemprego. Uma vez estimado esse modelo, poderemos nos perguntar se existe causalidade nessa relação, isto é, se inflação de fato ajuda a prever o desemprego, se o contrário ocorre ou se há uma simultaneidade. Nesse último caso, dizemos que existe uma causalidade bidirecional.

Ademais, esse tipo de análise também irá nos permitir estimar funções de impulso-resposta, onde analisamos a resposta a impulsos em uma das variáveis. Por exemplo, um choque no desemprego tem qual efeito sobre a inflação? E o caso contrário? Uma análise como essa pode ser muito interessante para avaliar a relação que existe entre duas ou mais variáveis, não é mesmo?

Um dos pressupostos, entretanto, do modelo acima é que as variáveis sejam estacionárias, no sentido que discutimos anteriormente. Isso significa que se essa hipótese for violada, teremos que partir para outro tipo de análise. No nosso Curso de Previsão Macroeconométrica usando o R, discutiremos a existência de cointegração entre duas ou mais variáveis, isto é, para o caso em que lidamos com séries não estacionárias, pode ser o caso de nos perguntarmos se existe uma relação entre de longo-prazo entre elas. Para o exemplo acima, se as séries não forem estacionárias, como de fato aparentam não ser, as mesmas podem estar relacionadas ao longo do tempo. Se for esse o caso, dizemos que as séries possuem uma tendência comum, o que nos permite avaliar a relação entre elas ao longo do tempo.

Para além disso, vamos agora ilustrar a aplicação de um modelo VAR a dados reais brasileiros. Estamos, por suposto, interessados em construir um modelo VAR que seja utilizado para fins de previsão da inflação mensal medida pelo IPCA. Para isso, vamos começar, como de praxe, carregando alguns pacotes.


### Pacotes
library(ggplot2)
library(forecast)
library(mFilter)
library(BETS)
library(vars)
library(scales)
library(gridExtra)
library(stargazer)

Os dados que precisaremos são coletados com uso do pacote \texttt{BETS} abaixo. Pegamos as séries de inflação, taxa Selic, taxa de câmbio e a produção industrial. Ademais, já tratamos as séries, preparando-as para o modelo VAR.


### Coletar os dados
inflacao <- BETSget(433)
selic <- BETSget(4189)
cambio <- BETSget(3697)
industria <- BETSget(21940)
# Transformações
dselic <- diff(selic)
dcambio <- diff(cambio)
# Criar e projetar hiato t+1
hp <- hpfilter(industria, freq=14400, type='lambda')
hiato <- ts(hp$cycle, start=start(industria), freq=12)
hiatof <- forecast(auto.arima(hiato, max.p=4, max.q=4,
seasonal = F), h=1, level=40)$mean
hiato <- ts(c(hiato, hiatof), start=start(industria), freq=12)
# Criar dummies sazonais
dummies <- window(ts(seasonaldummy(inflacao),
start=start(inflacao), freq=12), start=c(2007,01))
# Juntar séries
data <- window(ts.intersect(inflacao, dselic, dcambio, hiato),
start=c(2007,01))
colnames(data) <- c('inflacao', 'dselic', 'dcambio', 'hiato')

Um gráfico com as séries que utilizaremos é colocado abaixo.

Coletados e tratados os dados, o modelo VAR é então estimado com o código abaixo.


lag <- VARselect(data, lag.max=12, type='trend', season = 12)
var <- VAR(data, min(lag$selection), type='both', exogen = dummies)

Com o modelo em mãos, podemos gerar uma previsão como abaixo.


# Previsão
h <- 7
fvar <- stats::predict(var, n.ahead=h, ci=.4, dumvar=head(dummies,h))
fvar$fcst$inflacao

O output da previsão pode ser visto abaixo.

fcst lower upper CI
1 0.38 0.25 0.51 0.13
2 0.45 0.30 0.59 0.15
3 0.41 0.26 0.56 0.15
4 0.40 0.25 0.55 0.15
5 0.36 0.20 0.51 0.15
6 0.30 0.14 0.45 0.15
7 0.24 0.09 0.39 0.15

Esse tipo de modelagem é bastante utilizado entre profissionais que precisam de previsões acuradas no curto prazo para diversos tipos de variáveis. Ela é vista em detalhes nos nossos Cursos avançados de Macro Aplicada e de Análise de Séries Temporais.
____________________

equatiomatic: transformando modelos em equações LaTeX

By | Dicas de R

Quem trabalha com modelagem costuma ter que escrever equações de forma bastante rotineira. Para exercícios simples, é bastante tranquilo escrever uma ou outra equação em \LaTeX. O problema é quando você tem muitas equações no mesmo documento. Nessa Dicas de R - disponível toda quarta-feira aqui no blog da AM -vamos divulgar um novo pacote, o equatiomatic, que trata justamente desse problema.

library(equatiomatic)
library(palmerpenguins)
library(ggplot2)
library(latex2exp)

Usei a própria documentação do pacote para exemplificar o seu uso. Primeiro, rodamos um lm qualquer como o abaixo.


m <- lm(bill_length_mm ~ bill_depth_mm + flipper_length_mm, penguins)

Agora, basta usar a função extract_eq para que tenhamos acesso à equação.


extract_eq(m)

(1)   \begin{equation*} \operatorname{body\_mass\_g} = \alpha + \beta_{1}(\operatorname{bill\_length\_mm}) + \epsilon \end{equation*}

A equação extract_eq contém, inclusive, alguns argumentos que permitem a customização da equação a ser exibida. Para além disso, outra coisa legal do pacote é poder plotar gráficos com equações, como abaixo.


# Fit an lm model
m <- lm(body_mass_g ~ bill_length_mm, penguins)
# extract equation with `ital_vars = TRUE` to avoid the use of `\operatorname`
m_eq <- extract_eq(m, use_coef = TRUE, ital_vars = TRUE)
# swap escaped underscores for dashes
prep_eq <- gsub("\\\\_", "-", m_eq)
# swap display-style 
prep_eq <- paste("$", as.character(prep_eq), "$", sep = "")
# Plot
ggplot(penguins, aes(x = bill_length_mm, y = body_mass_g)) +
geom_point() +
geom_smooth(method = "lm") +
labs(title = "Relation between bill length and body mass",
subtitle = TeX(prep_eq))

____________

Um pdf e um script com todo o código desse exercício está disponível para os membros do Clube AM.

 

Análise de microdados da PNAD Covid com o R

By | Dicas de R

Com o R, é possível acessar diversas bases de dados e baixar o que precisa diretamente para o RStudio. Um exemplo disso são os microdados da PNAD Covid. Nessa Dicas de R - disponível toda quarta-feira aqui no blog da AM -vamos mostrar como pegar esses dados diretamente do ftp do IBGEComo de praxe, o código começa carregando alguns pacotes que utilizaremos.

library(tidyverse)

A seguir, nós podemos fazer o download dos dados diretamente do FTP do IBGE. Escolhemos os últimos microdados disponíveis da PNAD Covid.


## Baixar os dados via ftp
url = "ftp://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_PNAD_COVID19/Microdados/Dados/PNAD_COVID_112020.zip"
download.file(url, destfile='PNAD_COVID_112020.zip', mode='wb')
unzip('PNAD_COVID_112020.zip')
pnad_covid <- read_csv("PNAD_COVID_112020.csv",
col_types = cols(.default = "d"))

No código acima, é preciso observar que estamos baixando os dados para a pasta previamente setada no RStudio. Baixado os dados e importado para o RStudio, podemos agora baixar também o dicionário dos dados.


## Baixar o dicionário dos dados
url_2 = "ftp://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_PNAD_COVID19/Microdados/Documentacao/Dicionario_PNAD_COVID_112020.xls"
download.file(url_2, destfile='Dicionario_PNAD_COVID_112020.xls', mode='wb')

De posse desses dados, podemos agora explorá-los. Vamos supor que estamos interessados em saber quem recebeu o auxílio emergencial do governo por tipo de trabalho que a pessoa ocupa. Esse foi o questionamento do meu aluno Rodrigo Ashikawa, da nossa Formação Data Science para Economistas, que acabou me dando a ideia desse Dicas de R de hoje.

Com o código a seguir, nós criamos uma nova variável com o tipo de trabalho que a pessoa ocupa, de acordo com o dicionário dos microdados da pesquisa.


pnad_covid <- mutate(pnad_covid,
tipo_emprego = case_when(
C007 == 1 ~ "Trabalhador Doméstico",
C007 == 2 ~ "Militar",
C007 == 3 ~ "Policial ou Bombeiro",
C007 == 4 ~ "Setor Privado",
C007 == 5 ~ "Setor Público",
C007 == 6 ~ "Empregador",
C007 == 7 ~ "Conta Própria",
C007 == 8 ~ "Trabalhador Familiar",
C007 == 9 ~ "Fora do mercado",
is.na(C007) ~ 'Não informado'))

Uma vez feito isso, nós podemos agora selecionar as variáveis que queremos dos dados e também criar uma nova variável dummy que diz se a pessoa recebeu ou não o auxílio emergencial. O código a seguir implementa.


auxilio_emprego <- select(pnad_covid, D0051, tipo_emprego) %>%
mutate(auxilio_emergencial = ifelse(D0051 == 1, "Sim", "Não"))

Por fim, nós podemos agrupar os nossos dados, de modo a saber quantas pessoas da nossa base receberam o auxílio emergencial de acordo com o trabalho que ocupam.


auxilio_emprego <- auxilio_emprego %>%
group_by(tipo_emprego) %>%
summarise(
count_total = n(),
recebeu = sum(D0051 == 1),
n_recebeu = sum(D0051 == 2)) %>%
mutate("Auxilio Emergencial (%)" = (recebeu / count_total)*100)

O gráfico a seguir ilustra.

O gráfico resume a pesquisa que fizemos nos microdados da PNAD Covid. 58% dos que se identificaram como "Conta Própria" receberam o auxílio emergencial, enquanto 21,4% que se identificaram como "Militar" tiveram acesso ao benefício. Dos que colocaram "Não informado", o maior contingente, 52% tiveram acesso ao auxílio.

Um pdf e um script com todo o código desse exercício está disponível para os membros do Clube AM.

Caso você tenha interesse em microdados, conheça nosso Curso de Microeconometria usando o R.

____________

Baixando dados do Banco Mundial com o R

By | Dicas de R

Com o R, é possível acessar diversas bases de dados e baixar o que precisa diretamente para o RStudio. Um exemplo disso é a base de dados do Banco Mundial. Nessa Dica de R - sim, volto a publicar toda quarta-feira uma dica de R aqui no Blog - vamos mostrar como pegar os dados sobre poupança e taxa de juros de diversos países com o pacote WDI. Como de praxe, o código começa carregando alguns pacotes que utilizaremos.


library(WDI)
library(ggplot2)
library(ggrepel)
library(png)
library(grid)

A seguir, podemos pegar os dados que precisamos.


interest = WDI(country='all',
indicator = 'FR.INR.RINR',
start=2019, end=2019)

saving = WDI(country = 'all',
indicator = 'NY.GNS.ICTR.ZS',
start=2019, end=2019)

data = cbind(interest, saving)
data = data[complete.cases(data),]
data$iso2c = data$iso2c = data$country = data$year = data$year = NULL
colnames(data) = c('interest', 'country', 'saving')
data = subset(data, interest>0 & saving>0)

Um gráfico com os dados é posto abaixo.

 

______________

Para acessar os códigos completos desse exercício, é preciso fazer parte do Clube AM.

Google Trends no R: o pacote gtrendsR

By | Dicas de R

Com o avanço da pandemia do coronavírus, muitas consultorias e departamentos de research têm avançado na busca de dados de "alta frequência" para quantificar os seus efeitos sobre a economia. Os dados do google como o Google Trends e de geolocalização têm sido cada vez mais utilizados de forma a quantificar os efeitos da peste sobre o nível de atividade.

Como já tratei várias vezes nesse espaço, os dados do google podem ser inclusive utilizados para forecasting de variáveis econômicas. Um exemplo dessa abordagem é visto na edição 68 do Clube do Código, que busca replicar o paper The predictive power of google search in forecasting US unemployment, publicado no International Journal of Forecasting, para o Brasil.

Nesse paper e no exercícioutilizamos a pesquisa pela palavra "emprego" como uma das variáveis que explicariam o avanço da taxa de desemprego ao longo do tempo.

Na situação atual, contudo, talvez seja interessante pesquisar por outros termos, como, por exemplo, "seguro desemprego". Podemos para isso utilizar o pacote gtrendsR para fazer a pesquisa e os pacotes tidyverse para tratar e visualizar os dados.

Uma dica aqui é que a versão disponível no CRAN não rodou para mim. Tive que instalar a versão disponível no github. Para isso, você pode rodar a linha de comando abaixo.


if (!require("devtools")) install.packages("devtools")
devtools::install_github("PMassicotte/gtrendsR")

Uma vez instalado o pacote, podemos pegar tanto as buscas por "emprego" quanto "seguro desemprego", como no código abaixo.


data_gtrends = gtrends(keyword = c("seguro desemprego", 'emprego'),
geo = "BR", time='all', onlyInterest=TRUE)

De posse dos dados, nós selecionamos e mensalizamos as buscas por "seguro desemprego".


seguro_desemprego = data_gtrends$interest_over_time %>%
filter(keyword == 'seguro desemprego') %>%
mutate(mes = floor_date(date, "month")) %>%
group_by(mes) %>%
summarize(interesse = mean(hits)) %>%
mutate(date = as.Date(mes)) %>%
select(date, interesse)

Por fim, podemos gerar um gráfico com o ggplot2 como abaixo.

Como esperado, há um forte aumento em abril nas pesquisas por "seguro desemprego".

__________________

(*) Aprenda R em nosso Curso de Introdução ao R para Análise de Dados.

Cadastre-se em nossa Lista VIP para receber conteúdos exclusivos!

Fazer Inscrição
{"cart_token":"","hash":"","cart_data":""}