Machine Learning e a Copa do Mundo

Não adianta esbravejar pelo fato de que economistas, matemáticos e estatísticos gastam horas desenvolvendo modelos matemáticos sofisticados que apontam, basicamente, para o mesmo resultado quando se trata de estimar os resultados de um evento esportivo como a Copa do Mundo.

Adivinhar o campeão é apenas a cereja do bolo neste processo. Muito mais divertido é poder dar um palpite elegante sobre quem será o vencedor num confronto do tipo Bélgica e Inglaterra, principalmente para alguém que manja tanto de futebol quanto uma criança de 7 anos de idade – que é o caso dos autores desse estudo/blog; ou mesmo tentar prever quais seleções avançarão para cada etapa da competição, marcar alguns pontos a mais no bolão da firma e, com sorte, ganhar alguns trocados. Dito isso, aqui vai um alerta de spoiler: no exercício que se segue, Brasil e Alemanha têm chances muito próximas de ganhar a Copa do Mundo da Rússia, com ligeira vantagem para a seleção canarinho.

Mas antes de apresentarmos o modelo e seus resultados, o que explica essa concentração de apostas do “mercado” em torno de poucos países, como Brasil, Alemanha e Espanha?

O desempenho corrente e a experiência passada importam consideravelmente na avaliação do destino de cada seleção participante de uma Copa do Mundo. Dos 20 torneios mundiais disputados desde 1930, o Brasil ergueu cinco taças, a Itália e a Alemanha ganharam quatro vezes, a Argentina e o Uruguai saíram vitoriosos duas vezes cada, ao passo que Inglaterra, França e Espanha ganharam uma edição. Portanto, é impossível negar que a escolha de qualquer um desses times como campeão é uma aposta segura (com o cuidado de se observar o ranking atual e, obviamente, a presença deles no evento). A taxa de presença desses times nas Copas do Mundo ao longo da história e suas campanhas passadas reforçam esse tipo de aposta: o Brasil é o único time que participou de todas as edições, ao passo que a Alemanha não esteve presente em apenas dois torneios e a Itália deixou de ir ao maior evento esportivo do mundo em três ocasiões (já contabilizando 2018).

O número de vezes em que estes times disputaram as finais também é algo digno de nota: oito no caso da Alemanha, sete do Brasil e seis da Itália. É ainda mais impressionante o fato de que, até hoje, apenas oito times conseguiram chegar até a etapa final da competição, tendo como base as últimas 13 edições: além dos três destacados acima, foram os casos de Argentina, Espanha, França, Inglaterra e Holanda.

As semifinais são eventos um pouco mais imprevisíveis do que as finais. Vinte e quatro times chegaram até esta etapa ao longo da história do futebol. No entanto, uma vez mais, apenas cinco deles (Alemanha, Argentina, Brasil, França e Itália) responderam por mais de 50% das vagas nesta fase da competição. Vale destacar que, ao que parece, em todo torneio há pelo menos um azarão nas semifinais. Em 2002 foram dois países: Turquia e Coréia do Sul, apesar deste último ter sido uma das seleções da casa (o que confirma a relevância do fator-casa em modelos da Copa); em 2006, muitos apontam que foi o caso de Portugal; e em 2010 essa alcunha de azarão foi dada ao Uruguai.

Como destacado acima, ser seleção do país-sede certamente é uma vantagem. Os times da casa ergueram a taça em pouco menos de 1/3 das Copas organizadas até agora, chegaram às semifinais em 50% das vezes e estiveram presentes nas oitavas de final em todos os eventos desde 1986, exceto na África do Sul.

Neste momento o leitor já deve ter se convencido de que, embora diversos fatores (muitas vezes imprevisíveis) influenciem os resultados jogo-a-jogo, nenhum modelo estatístico bem-especificado poderia sugerir uma semifinal sem a presença de duas ou três das seleções destacadas acima. No mesmo sentido, vale destacar também que nem abastecido com todos os dados históricos do mundo um modelo poderia transformar um time pouco competitivo em um campeão. Exemplo disso é que algumas seleções reconhecidas por sua tradição no futebol ficaram de fora da Copa do Mundo de 2018, caso da Itália e da Holanda.

Mas é possível avaliar de forma objetiva o desempenho passado e corrente de uma seleção de modo a termos um bom palpite para os jogos da Copa da Rússia?

Os sites de aposta (como o Bet365 e Betfair) podem ser um bom ponto de partida como uma medida subjetiva de força de cada seleção. Também é possível basear-se no ranking oficial da FIFA, apesar de a instituição ter declarado recentemente que a fórmula de cálculo do ranking sofrerá alteração. Nesta simulação nos baseamos no rating ELO, justamente por ser uma base padronizada e bem organizada de todas as partidas oficiais (desde 1872!) envolvendo times nacionais. A partir destes dados temos algo em torno de 45 mil jogos para alimentar nossos modelos de previsão.

Aqui vale uma pequena pausa para explicar o que é o ELO rating e como este indicador é calculado. O ELO foi desenvolvido por um físico Húngaro-Americano chamado Arpad Elo (1903-1992), com o objetivo de medir e ranquear a habilidade de jogadores de xadrez. Os ratings ELO têm sido utilizados em diversos esportes, como o tênis, mas a sua aplicação no mundo do futebol ganhou relevância nos últimos anos por se tratar de uma métrica de habilidade (ou força) dos times com qualidade superior ao ranking da FIFA.

O método ELO de ranqueamento leva em conta não só o número de vitórias, derrotas e empates de cada time, mas também as condições sob as quais o evento ocorreu. Como resultado, derrotar um time forte como o Brasil ou a Alemanha aumentará o ELO do time vencedor numa proporção muito maior do que uma vitória sobre Andorra ou Zâmbia. Além disso, uma vitória fora de casa gera um incremento maior sobre o ELO do que uma vitória dentro do próprio país; ganhar uma qualificatória é melhor que ganhar um amistoso; ou ainda, uma goleada de 7 a 0 gera mais pontos que uma vitória de 2 a 1. A partir do ELO também é possível comparar a força de cada seleção ao longo do tempo. Por exemplo, o time alemão que venceu a Copa do Brasil de 2014 era mais forte (no início do torneio) do que as seleções germânicas que venceram a competição em 1954, 1974 e 1990.

Olhando um pouco para o histórico do ELO, em nenhuma ocasião o time vencedor iniciou o torneio com rating abaixo de 2030 pontos, exceto na vitória surpresa do Uruguai sobre o Brasil na Copa de 1950. O time mais forte a ganhar uma Copa do Mundo foi o da Alemanha de 2014, com um ELO de partida de 2223 pontos (mais informações em https://www.eloratings.net/).

A Copa do Mundo de 2018, que possui um ELO médio de 1828 para todos os 32 times envolvidos, é o terceiro torneio mais fraco da história sob esta ótica, ficando à frente apenas dos mundiais de 1930 e 1938. No evento deste ano, apenas três seleções contam com rating superior a 2000, ao passo que em 2014 cinco times iniciaram o torneio com ELO acima desse patamar. Mas o outlier nesse caso foi a Copa passada. Desde 1994, apenas 3 seleções, em média, iniciaram a Copa com ratings superiores a 2000.

Com relação ao rating da seleção brasileira no evento deste ano, com 2142 pontos temos a segunda escalação mais forte a participar de uma Copa do Mundo, ficando atrás apenas daquela que entrou em campo em 2014.  No caso da Alemanha, a seleção de hoje é a terceira mais forte a participar de um mundial (ELO 2077), ficando atrás dos times de 1978 (ELO 2083) e 2014 (ELO 2109).

Metodologia

Como destacado acima, utilizamos resultados de cerca de 45 mil jogos extraídos do site ELO ratings. Neste conjunto de dados constam diversas informações referentes a cada partida, como data do jogo, nome e código dos países, número de gols marcados pelas seleções, local (país) do jogo e nome do campeonato, além, é claro, do rating ELO das respectivas seleções, calculado ao final de cada partida.

Para esta simulação decidimos não percorrer o perigoso caminho de prever o país que tem maior probabilidade de vencer a Copa do Mundo da Rússia. Ao invés disso, nos debruçamos em estimar probabilidades jogo-a-jogo, o que nos dá a composição dos grupos durante a primeira fase (a colocação de cada seleção ao final desta etapa) e os resultados das oitavas, quartas, semifinais e, por fim, as probabilidades para os dois times que disputam a final do campeonato. Nos próximos parágrafos esmiuçamos um pouco mais a metodologia. De antemão peço desculpas pelo trecho carregado de termos técnicos. Ao leitor que não interessar este tipo detalhamento, recomendamos pular direto para a seção de resultados.

Tivemos como objetivo principal gerar probabilidades de vitória, empate e derrota para cada jogo da Copa da Rússia. Para isso utilizamos a abordagem de regressão logística multinomial, pois desta forma pode-se ter mais de duas categorias em sua variável dependente.

Definimos, portanto, três classes para determinar as probabilidades de vitória do time A sobre o time B, do time B sobre o time A e de empate nessa partida, tendo como variáveis/features (i) o rating ELO corrente e suas defasagens, (ii) uma variável binária para descontar o “fator casa”, (iii) médias móveis de gols a favor e (iv) médias móveis de gols contra. Para agregar mais features às simulações, criamos transformações polinomiais de 2o e 3o grau dessas variáveis, totalizando cerca de 100 variáveis no modelo.

O próximo passo foi utilizar a técnica de validação cruzada k-fold para treinar e testar nosso modelo de Machine Learning com as variáveis mais importantes, selecionadas pelo método LASSO (Least Absolute Shrinkage and Selection Operator). Nesta etapa, o algoritmo divide os dados de entrada em subconjuntos de dados (chamados de folds); depois o modelo é treinado em todos os subconjuntos, exceto um (k-1), e avaliado neste subconjunto que não foi usado para o treinamento. Este processo é repetido k vezes, cada uma com um subconjunto diferente reservado para avaliação (e excluído do treinamento).

A figura abaixo mostra de forma mais clara como a metodologia de validação cruzada 4-fold foi aplicada ao nosso modelo, gerando subconjuntos de treinamento e avaliação. A primeira iteração usa os primeiros 25% dos dados para avaliação e os 75% restantes para treinamento. A iteração seguinte usa o segundo quartil (de 25% a 50%) para avaliação e os três quartis restantes para treinamento, e assim por diante. Na próxima seção apresentamos os resultados dessa simulação.

 Copa2018_01

Resultados

A tabela abaixo mostra as probabilidades de vitória, empate e derrota para cada partida durante a fase de grupos. A última coluna indica o time com maior probabilidade de sair vitorioso. Os dados estão ordenados por grupo.

Os resultados das simulações mostram que a Rússia deve se beneficiar do “fator-casa” no confronto contra a Arábia Saudita.  Além disso, o modelo sugere o empate como uma boa aposta nos confrontos entre Rússia e Egito, Peru e Dinamarca, Islândia e Croácia, dentre outros.

Com base nestas simulações, os primeiros e segundos colocados dos grupos seriam, respectivamente: Uruguai e Rússia (A), Espanha e Portugal (B), França e Peru (C), Argentina e Croácia (D), Brasil e Suíça (E), Alemanha e México (F), Inglaterra e Bélgica (G) e Colômbia e Polônia (H).

Copa2018_02

Partindo dos resultados da fase de grupos, os resultados simulados para as demais etapas são apresentados na figura abaixo. Como destacado no início deste artigo, Brasil e Alemanha devem disputar a partida final e as simulações apontam para ligeiro favoritismo do Brasil. Rumo ao Hexa!

Copa2018_03

Advertisements

Dinâmica da pirâmide etária entre 1990 e 2050 e a reforma da previdência

Recentemente a equipe econômica do presidente interino Michel Temer anunciou um conjunto de medidas para tentar conter a rápida deterioração das contas públicas brasileiras. A mais importante delas, que será apresentada a partir de uma Proposta  de Emenda Constitucional, tem como objetivo limitar o crescimento das despesas primárias à inflação do ano anterior. O governo também pretende que o BNDES antecipe a devolução de pelo menos R$ 100 bilhões em recursos repassados pelo Tesouro (de um total de mais de R$ 500 bilhões) e indicou, dentre outras coisas, que poderá extinguir o Fundo Soberano Nacional e utilizar os R$ 2 bilhões de seu patrimônio atual. Foi a raspa do tacho na tentativa de reduzir em alguma medida o enorme déficit primário herdado da administração Dilma, sem propor (por enquanto!) nova rodada de aumento de impostos, mudança nas regras de benefícios sociais e/ou venda de ativos da União.

Embora o pacote de medidas tenha agradado o mercado, o governo preferiu por ora não entrar no debate da questão previdenciária. Michel Temer apenas indicou que discutirá a reforma da previdência com centrais sindicais e classe política e que nenhuma medida econômica será tomada sem a “concordância” da sociedade. Fato é que nós brasileiros precisamos nos preocupar com essa questão, pois o Estado já gasta mais de 12% do PIB com previdência (incluindo BPC-LOAS, de uma arrecadação total de cerca de 25% do PIB) e, se nada for feito, a situação se tornará insustentável em pouquíssimo tempo, uma vez que o processo de envelhecimento da população brasileira será muito rápido. O gráfico abaixo mostra de forma mais clara a velocidade deste processo de envelhecimento da população entre 1990 e 2050, utilizando as projeções do U.S. Census Bureau.

brazil_pyramid

Pelo estreitamento da base da pirâmide é possível notar que a proporção de crianças e adolescentes com até 14 anos vem caindo rapidamente nos últimos anos, reflexo da menor taxa de fecundidade das mulheres. Segundo o IBGE, o número médio de filhos passou de 2,4 nos anos 2000 para 1,7 em 2015, uma queda de quase 30%. Com efeito, a expectativa é de um aumento expressivo da população de idosos nos próximos 30-40 anos, característica bastante distinta da observada no início da década de 90, quando predominavam crianças e jovens na pirâmide etária brasileira.

O próximo gráfico mostra, a partir da mesma base do U.S. Census Bureau, que a taxa de fecundidade brasileira também é a menor da América Latina. Para visualizar os valores, basta clicar no gráfico.Taxa de Fertilidade, 2016 <br> Elaboração: dadosdadosdados. Fonte: U.S. Census BureauSe as informações demográficas apresentadas acima ainda não forem suficientes para convencer o leitor de que o debate sobre a reforma da previdência deve ser tratado em regime de urgência, sugiro a leitura da última versão (de março de 2015) do relatório de Projeções Financeiras e Atuariais para o Regime Geral de Previdência Social, disponível neste link. Na página 28 encontram-se as projeções oficiais (ou seja, do próprio governo), que indicam que, nas condições atuais, a despesa do INSS e o déficit da Previdência como proporção do PIB seguirão em tendência de crescimento até 2060, assim como mostrado nos gráficos abaixo. Importante também destacar que essas projeções levavam em conta um PIB de -0,9% em 2015, de +1,30% em 2016 e de +1,9% em 2017. A atualização do modelo com o recuo de -3,8% observado em 2015 e com as projeções mais recentes do Boletim Focus (-3,8% em 2016 e +0,55% em 2017), além de inflação mais elevada e crescimento nominal da massa salarial bem mais fraco, influenciará negativamente a estimativa de arrecadação, aumentando a previsão para o déficit da Previdência.

Despesa (azul) e Déficit (laranja) da Previdência como % do PIB <br> Elaboração: dadosdadosdados. Fonte: Previdência SocialOBS: para elaborar os gráficos deste post utilizei diversas bibliotecas do software R, sendo que as principais foram a idbr (que faz conexão direta com o API do US Census Bureau para obter os dados demográficos), ggplot e plotly (para criar os gráficos) e a animation (responsável por criar a animação do gráfico em GIF).

Google Trends, concorrência e crise econômica

F1As pesquisas na web, que já fazem parte do cotidiano de diversas pessoas há muito tempo, estão se popularizando cada vez mais através do acesso à internet móvel por smartphones e tablets. Hoje em dia não é mais necessário memorizar e digitar diversos endereços para encontrar a resposta que procura. Basta digitar no Google uma palavra-chave, que o sistema de busca se encarrega de trazer os resultados mais relevantes (abaixo, é claro, dos anúncios pagos). Geralmente a resposta estará nos primeiros cinco ou dez links.

O divertido de tudo isso é que cada termo procurado é armazenado e ranqueado (a partir da quantidade de buscas) em um sistema chamado Google Trends (link). Ao final de cada ano são divulgadas listas das pesquisas mais populares, separando os termos inclusive por tema. No Brasil, por exemplo, as três palavras mais pesquisadas em 2015 foram “BBB15”, “Cristiano Araújo” e “Resultado do Enem”. Nada muito diferente de 2014, quando os três termos mais buscados pelos brasileiros foram “BBB14”, “Enem” e “Amor à Vida”.  Mas é só isso? A única vantagem do Trends é medir a popularidade de um artista, música ou novela? Não.

Para a felicidade dos Data Scientists, esta ferramenta vai além da simples ordenação dos termos mais pesquisados; ela disponibiliza desde 2004 a evolução semanal do número de buscas de cada palavra, indicando se determinado tema ganhou ou perdeu relevância ao longo do tempo; e dependendo do termo pesquisado, podemos extrair algumas conclusões bacanas sobre concorrência, crescimento econômico, confiança dos agentes e até mesmo identificar padrões sazonais.

O software R, a partir da biblioteca gtrendsR, nos permite baixar e trabalhar com os dados do Google Trends de uma maneira fácil e intuitiva. Basta entrar com os dados da sua conta Gmail e informar os termos da pesquisa, como destacado no exemplo abaixo para “bares” e “restaurantes” no Brasil.

library(gtrendsR)
user <- "usuario@gmail.com"
psw <- "senha"
gconnect(user, psw)
rm(psw)
exemplo <- gtrends(c("bares", "restaurantes"), geo = "BR", res="week")
plot(exemplo)

Os termos utilizados no código de exemplo não têm nada de emocionante, então economizarei espaço omitindo o gráfico. Mas e se fizermos uma pesquisa com alguns sites de e-commerce no Brasil, seria possível identificar os maiores players e o que tem acontecido em termos de concorrência neste segmento somente a partir das buscas no Google?

G1

Através do gráfico acima é possível perceber que a quantidade de buscas pelo site Americanas.com é mais do que o dobro das buscas pelo site do Submarino. Além disso, vemos que este último perdeu relevância ao longo do tempo – que pode ter sido “roubada” em parte pela loja virtual do Walmart. Também interessante notar como as buscas por estes sites de e-commerce cresce nos dias de black friday (última semana de novembro de cada ano) e como este evento importado dos EUA tem ganhado popularidade mais recentemente.

Olhando agora para a evolução das pesquisas por “pós-graduação” e “MBA” (gráfico abaixo), é possível perceber uma forte tendência de queda na procura destes cursos de especialização, o que pode estar em parte relacionado ao elevado comprometimento de renda das famílias e ao enfraquecimento da atividade econômica nos últimos anos.

G2

Por fim, algumas pesquisas do Google, além de terem boa aderência com a evolução da atividade econômica, sugerem um padrão sazonal interessante. É o caso do termo “emprego”. Pelo gráfico abaixo nota-se que a procura por “emprego” aumenta bastante no início de cada ano e depois cai gradualmente até atingir o mínimo anual no mês de dezembro. Os meses de fevereiro e março também mostram algumas quedas abruptas, provavelmente relacionadas ao carnaval. Também destaca-se a forte alta na procura por emprego neste início de 2016, coerente com a deterioração observada no mercado de trabalho ao longo de 2015.

G3

Trabalhando com séries temporais no R – Parte II (estatísticas fiscais)

No post anterior (link) demonstrei como importar e definir séries temporais no R quando o vetor de datas não está no formato adequado (no caso de estar formatado como texto, por exemplo). Esse é um problema bastante comum e precisa ser tratado independente do software que o analista está utilizando.

Também destaquei que os canais de acesso à informação (IBGE, Banco Central, Ipeadata, etc.) não seguem um padrão de armazenamento dos dados – geralmente disponibilizam as séries no formato que julgam ser mais conveniente. Mas em alguns casos a formatação dos dados não está relacionada à preferência de quem os divulga. Algumas bases têm de seguir certas normas e padrões internacionais. Esse é o caso das séries de finanças públicas do Banco Central do Brasil.

Devido ao grande número de indicadores, as estatísticas fiscais seguem um padrão de divulgação bastante diferente daquela estrutura que trabalhamos no nosso primeiro exemplo; os dados estão dispostos nas colunas e os rótulos nas linhas. Além disso, na maior parte dos casos não há sequer um vetor de datas e qualquer preocupação de atribuir códigos exclusivos para cada uma das séries, de modo a facilitar a vida do analista quando, em uma mesma base, existem séries com nomes (rótulos) idênticos.

Neste exemplo utilizaremos as séries históricas de Necessidade de Financiamento do Setor Público, disponível na seção de séries temporais do site do Bacen (link). Este conjunto de dados foi escolhido por dois motivos. Primeiro porque atende os três pré-requisitos de base pentelha: (i) dados transpostos, (ii) sem vetor de datas e (iii) ausência de códigos exclusivos em cada indicador. Em segundo lugar, porque este é o assunto do momento – o mundo todo está atento para a forte e rápida deterioração da situação fiscal brasileira, principalmente por conta da incapacidade do governo de cortar gastos, mas também pela queda abrupta da arrecadação federal (afetada sobremaneira pela desaceleração da atividade econômica).

Nosso objetivo é exatamente o mesmo do alcançado no post anterior: deixar a série de dados no formato mais adequado, para posteriormente podermos utilizar todo instrumental de data analysis do R. Também lançaremos mão de um recurso incrível do R: baixar os dados diretamente da fonte e extrair o “zip” para uma pasta temporária do computador. Isso será especialmente útil para aqueles que precisam atualizar mensalmente bases, apresentações gráficas ou modelos econométricos.

O código foi escrito em oito passos. Para não poluirmos os comandos com diversos comentários, descreverei abaixo o que está sendo feito em cada etapa.

  • 1o passo: cria pasta temporária; baixa o arquivo “zip” do site do Bacen e salva na pasta criada; extrai os arquivos para a pasta; cria base “original” com os dados da aba “$-12meses(corrente)”.
  • 2o passo: deleta as linhas e colunas da base “original” onde todas as células estão como “NA”; cria “base1” deixando apenas as linhas com as informações que interessam.
  • 3o passo: cria um vetor de datas baseado no número de colunas (exclusive a coluna A de nomes), definindo apenas a data de início da série; renomeia os rótulos das colunas da “base1” com o vetor de datas, e as linhas com o número de linhas.
  • 4o passo: cria um vetor simples de referência (REF) que “codifica” as linhas para facilitar a escolha dos indicadores; adiciona esse vetor na primeira coluna da “base1”.
  • 5o passo: cria a “base2” convertendo a “base1” para o formato long (mais adequado).
  • 6o passo: converte cada coluna da “base2” para o formato adequado. No vetor de nomes (coluna “Indicador”), remove espaçamento extra (“ Nominal” vira “Nominal”). O vetor de datas é convertido para “aaaa/mm/dd”, o padrão internacional. Já o vetor de Valor é convertido para o formato numérico.
  • 7o passo: cria a “base3”, a partir da “base2”, selecionando o indicador de interesse (REF: 21 = “Primário”), que representa o resultado primário do setor público consolidado no acumulado em 12 meses; divide a coluna de valor por 1.000 (para representar em R$ Bilhões) e inverte o sinal.
library(XLConnect)
library(lubridate)
library(ggplot2)
library(dplyr)
library(reshape2)
library(xlsx)
library(scales)

# 1o PASSO
temp <- tempfile()
download.file("http://www.bcb.gov.br/ftp/notaecon/nfspp.zip", temp)
unzip(temp)
original <- read.xlsx(temp, file = "Nfspp.xls", 
                       sheetName = "$-12meses(corrente)", as.data.frame = TRUE,  
                       encoding = "UTF-8")

unlink(temp)

# 2o PASSO
original <- original[rowSums(is.na(original)) != ncol(original),]
original <- original[,colSums(is.na(original)) != nrow(original)]
base1 <- original[rowSums(is.na(original)) == 0,]

# 3o PASSO
meses <- ncol(base1)-1
datas <- seq.Date(as.Date("2002-11-01"), by = "month", length.out = meses)
colnames(base1) <- c("Indicador", as.character.POSIXt(datas))
row.names(base1) <- 1:nrow(base1)

# 4o PASSO
REF <- row.names(base1)
base1 <- cbind(REF, base1)

# 5o PASSO
base2 <- melt(base1, id = c("REF", "Indicador"), variable.name = "Data", value.name = "Valor")

# 6o PASSO
base2$Indicador <- gsub("\\s+"," ", base2$Indicador) 
base2$Data <- as.POSIXct(parse_date_time(base2$Data, "%Y/%m/%d"))
base2$Valor <- as.numeric(base2$Valor)

# 7o PASSO
selecao <- "21"
base3 <- filter(base2, REF == selecao)
base3$Valor <- base3$Valor/1000*(-1)

8o passo: cria gráfico de linha a partir da “base3”; as opções colocam título, linha reta cruzando o eixo Y no ponto zero, modificam a escala do eixo X para meses e colocam títulos nos eixos.

Pelo gráfico abaixo é possível perceber a forte deterioração do resultado primário do setor público consolidado (que considera o governo federal, os estados e municípios, bem como as empresas estatais não financeiras, exceto a Petrobras e o grupo Eletrobras) desde o final do primeiro trimestre de 2014. É mais do que justificável, portanto, a desconfiança do mercado na capacidade do governo de produzir um superávit primário de R$ 8,747 bilhões (0,15% do PIB) em 2015, uma vez que o resultado acumulado até agosto está negativo em R$ 43,8 bilhões.

# 8o PASSO
graf1 <- ggplot(base3, aes(Data, Valor)) + 
  geom_line(colour = "red", size = .7) + 
  ggtitle("Resultado Primário do Setor Público \n Acum. 12 meses. Elaboração: dadosdadosdados") + 
  theme(plot.title = element_text(lineheight=.9, face="bold")) + 
  geom_hline(yintercept=0, col = "black", linetype = "dashed") +
  scale_x_datetime(breaks = date_breaks("8 months"), labels = date_format("%b/%y")) +
  theme(axis.text.x = element_text(angle = 90))
update_labels(graf1, list(y = "R$ Bilhões", x = ""))

plot of chunk unnamed-chunk-3

Por fim, selecionaremos mais de 1 indicador para colocar no gráfico, e também reduziremos a amostra de análise. Assim poderemos observar composição mais recente do resultado primário por esfera de governo (nível federal e regional).

O resultado final é o destacado no gráfico abaixo.

# Selecionando mais indicadores e reduzindo a amostra  
selecao <- c("22", "27", "30", "31")
base3 <- filter(base2, REF %in% selecao)
base3$Valor <- as.numeric(base3$Valor)
base3$Valor <- base3$Valor/1000*(-1)

amostra <- filter(base3, Data >= "2012-12-01" & Data <= "2015-08-01")
## Warning in data.matrix(data): NAs introduced by coercion
  graf2 <- ggplot(amostra, aes(x=Data, y=Valor, fill=Indicador)) + 
    geom_bar(stat="identity") + 
    ggtitle("Resultado Primário do Setor Público \n Acum. 12 meses. Elaboração: dadosdadosdados") + 
    theme(plot.title = element_text(lineheight=.9, face="bold")) + 
    geom_hline(yintercept=0, col = "black", linetype = "dashed") +
    scale_x_datetime(breaks = date_breaks("2 months"), labels = date_format("%b/%y")) +
    theme(axis.text.x = element_text(angle = 90), legend.position = "bottom") 
  update_labels(graf2, list(y = "R$ Bilhões", x = ""))
## Warning: Stacking not well defined when ymin != 0

plot of chunk unnamed-chunk-4

Trabalhando com séries temporais no R – Parte I

O principal objeto para armazenamento de dados no software R é o data.frame, que representa uma estrutura retangular onde cada coluna pode abrigar informações de diferentes tipos (números, textos, datas, etc.). Embora o R seja uma ferramenta poderosa para se trabalhar com qualquer tipo de dado, no seu estado natural (ou seja, sem a ajuda das bibliotecas criadas pela comunidade de usuários), manejar séries de tempo pode ser frustrante, sobretudo se os dados não seguirem um padrão fixo de espaçamento – como ocorre, por exemplo, com séries intraday do mercado financeiro.

O pacote básico do R inclui duas classes para operações com séries de tempo: a ts para dados univariados e a mts para dados multivariados, que são úteis para trabalhar com dados mensais e trimestrais (de espaçamento fixo), mas que contam com poucas funções (métodos) de manipulação e processamento dos dados, além de uma interface gráfica bem pobre.

A partir dessas limitações, diversos pacotes foram criados para tornar mais fácil a vida dos que querem brincar com séries temporais no R. Os mais conhecidos são o zoo e o xts, que basicamente armazenam os dados em dois componentes: (i) um vetor de datas e/ou tempo e (ii) um data.frame com as observações (dados). Os pacotes de séries de tempo servem para combinar estes dois elementos e possibilitar a aplicação de diversas funções.

Mas antes de começarmos a trabalhar com as séries de tempo, serão necessárias mais algumas bibliotecas (nesse caso, o XLConnect e o lubridate), visto que precisaremos importar os dados de interesse para o R e muito provavelmente padronizar o vetor de datas. Aliás, em um típico processo de análise, a qualquer momento as etapas podem incluir alguma transformação dos dados, imputação de valores inexistentes, adicionar/deletar variáveis e recomeçar o processo inteiro novamente (como mostrado na figura abaixo).

fig 1

Nesse momento você deve estar se perguntando: por que diabos preciso me preocupar em padronizar o vetor de datas? A resposta passa pelo fato de que, no Brasil, cada fonte divulga seus dados da forma como acha conveniente e, na maioria das vezes, este formato não é o mais apropriado. Um exemplo clássico é o Ipeadata (link) – repositório gratuito de dados macroeconômicos, regionais e sociais de diversas fontes domésticas e internacionais -, que divulga o vetor de datas como texto (COMO TEXTO!). Em séries mensais o formato utilizado pelo Ipea é o “YYYY.MM”. Ou seja, a referência para set/15, que deveria ser 01/09/2015, é armazenada na coluna A como “2015.09”. Aí dificulta a vida do cidadão! Mas não se preocupe… com uma linha de comando o R resolve isso.

Neste post mostrarei como importar uma base para o R, transformar dados em séries temporais e também os primeiros passos de como gerar gráficos mais elegantes e elaborados. Nas próximas publicações partiremos para análises de decomposição, construção/avaliação de modelos de previsão e projeções de valores futuros.

No nosso exemplo, utilizaremos os dados do IBC-Br (indicador proxy do PIB, criado pelo Banco Central do Brasil), baixados lá no site Ipeadata. Basta acessar o link e digitar no campo de pesquisa o termo “IBC-Br” (sem as aspas). A série que nos interessa é a “Índice de Atividade Econômica do Banco Central (IBC-Br) (2002=100)”. Agora é só exportar em formato “xls”. De preferência, nomeie a planilha como “IBCBR.xls”, pois é este o nome que utilizaremos no código.

A tabela abaixo resume as funções e pacotes que utilizaremos no código.

fig 2

library(zoo)
library(XLConnect)
library(lubridate)
library(ggplot2)

#Importando os dados e renomeando as colunas
setwd("C:\\Users\\Tho\\Google Drive\\Blog\\201509\\")
wb <- loadWorkbook("IBCBR.xls")
  Matriz <- readWorksheet(wb, sheet = "Séries", startRow = 0, startCol = 0)

#alterando o vetor de datas para o formato correto   
  names(Matriz) <- c("Data", "IBCBr")

  Matriz$Data <- parse_date_time(Matriz$Data, "%Y.%m")

#transformando dados em time-series     
  tsIBCBr <- zoo(Matriz$IBCBr, Matriz$Data)

O resultado final é o destacado no gráfico abaixo.

#gerando gráficos mais elaborados      
  lc <- qplot(Data, IBCBr, data = Matriz, geom = "line") + 
    geom_line(colour = "red", size = .7)
  lc + ggtitle("IBC-Br (2002=100, Banco Central do Brasil) \n Elaboração: dadosdadosdados") + 
    theme(plot.title = element_text(lineheight=.9, face="bold"))

plot of chunk unnamed-chunk-3

Dessazonalizando o PIB brasileiro pelo R – Parte II

No post anterior, abordamos rapidamente a importância de aplicar o ajustamento sazonal em séries temporais. Agora focaremos apenas nos passos para replicar a metodologia do IBGE para a dessazonalização da série trimestral do PIB brasileiro. Vamos lá?

  1. Baixando e extraindo o X-13-ARIMA-SEATS no site do U.S. Census Bureau
    • A versão do X-13 que precisamos é a “x13ashtmlall.zip” (no momento que escrevo a versão mais atual é “_V1.1_B19”), disponível neste link;
    • Depois de feito o download, extraia os arquivos em qualquer pasta do seu computador (OBS: a pasta será utilizada no código);
    • Para facilitar, manteremos no código do R o nome original da pasta do X-13: “x13ashtml”;
  2. Baixando a série histórica do PIB Total das Contas Nacionais Trimestrais, disponível neste link
    • No site do IBGE, precisamos das “Tabelas Completas”, com link disponível no canto esquerdo da página;
    • Aconselhamos salvar e extrair a planilha na mesma pasta do X-13, pois ela também servirá de output para os dados dessazonalizados pelo código R.
    • Os dados que nos interessam nesta planilha estão na aba “Série Encadeada”, coluna R (referente ao PIB Total, em número índice).
  3. Verificando os parâmetros de dessazonalização do IBGE, disponível neste link
    • O IBGE divulga os parâmetros de dessazonalização do PIB na “Publicação completa” das Contas Nacionais Trimestrais” (em formato PDF, também no canto esquerdo da página). Dentro do documento, procure pelo quadro de “modelos adotados no ajuste sazonal”, disponível no anexo de notas metodológicas.
    • Importante destacar que os parâmetros de ajustamento sazonal não são fixos, de modo que sempre que houver nova divulgação as Contas Nacionais Trimestrais, o IBGE estimará novamente o modelo ARIMA e gerará diferentes outliers (efeitos de intervenção).
      Pronto, já temos todos os elementos para executar a rotina de dessazonalização no R. Mas antes cabem mais algumas ressalvas sobre as linhas do código: (i) atenção para a AMOSTRA que se pretende dessazonalizar – no exemplo abaixo, “ ,2015.1” significa ajustar do início da amostra até o 1º trimestre de 2015 e (ii) atenção para as pastas onde foram salvos o programa X-13 e os dados do IBGE – no exemplo, os caminhos estão definidos como “C:\Users\Tho\Google Drive\Blog\201508\x13ashtml”

O código precisará de três bibliotecas: seasonal, XLConnect e dplyr. Se algum erro surgir na execução do código, tente reinstalar as bibliotecas utilizando o install.packages(“XLConnect”, “dplyr”, “seasonal”)

library(XLConnect)
library(dplyr)

# Especificar o caminho onde foi salvo o programa X-13
# Para rodar da primeira vez, necessita apontar o caminho do x13
Sys.setenv(X13_PATH = "C:\\Users\\Tho\\Google Drive\\Blog\\201508\\x13ashtml")
library(seasonal)

# Mudar a AMOSTRA de acordo com a base de dados
AMOSTRA <- ",2015.1"

#especificar o caminho onde foi salva a planilha "Tab_Compl_CNT_1T15.xls" com a série do PIB. Neste exemplo, utilizamos a mesma pasta do X-13
arq <- 'C:\\Users\\Tho\\Google Drive\\Blog\\201508\\x13ashtml\\Tab_Compl_CNT_1T15.xls'
Original = readWorksheetFromFile(arq, sheet = "Série Encadeada",
                                 startRow = 4, startCol = 18, endCol = 18, useCachedValues = TRUE)

PIBT <- ts(Original[,1], start = 1996, frequency = 4)

plot(PIBT)

plot of chunk unnamed-chunk-2

#abaixo são definidas a decomposição do modelo (mult para Multiplicativa ou add para Aditiva), o modelo ARIMA "(0 1 1)(0 1 1), por exemplo" e os efeitos de intervenção ("TC1996.3", "AO1998.1", ...). Se a decomposição for Aditiva (x11.mode = "add"), colocar transform.function = NULL. Se x11.mode = "mult", então transform.function = "log"

PIBTSA <- seas(x = PIBT,
               series.modelspan = AMOSTRA,
               arima.model = "(0 1 1)(0 1 1)",
               regression.variables = c("LS2008.4", "Easter[1]"),
               forecast.maxlead = 6,
               forecast.lognormal = "yes",
               regression.aictest = NULL,
               outlier = NULL,
               transform.function = "log",
               automdl = NULL,
               x11.mode = "mult")

O resultados finais estão destacados no gráfico e nas tabelas abaixo.

#testes, gráficos e resultados
plot(PIBTSA)

plot of chunk unnamed-chunk-3

summary(PIBTSA)
##
## Call:
## seas(x = PIBT, transform.function = "log", regression.aictest = NULL,
##     outlier = NULL, automdl = NULL, series.span = AMOSTRA, regression.variables = c("TC1996.3",
##         "AO1998.1", "LS2008.4", "TC2009.1", "Easter[15]"), arima.model = "(0 2 2)(0 1 1)",
##     x11.mode = "mult")
##
## Coefficients:
##                    Estimate Std. Error z value Pr(&amp;amp;amp;amp;gt;|z|)
## TC1996.3           0.036471   0.006370   5.726 1.03e-08 ***
## AO1998.1          -0.025803   0.005112  -5.047 4.49e-07 ***
## LS2008.4          -0.053941   0.008096  -6.663 2.69e-11 ***
## TC2009.1          -0.028206   0.006147  -4.589 4.46e-06 ***
## Easter[15]        -0.006106   0.001515  -4.029 5.60e-05 ***
## MA-Nonseasonal-01  0.635040   0.113684   5.586 2.32e-08 ***
## MA-Nonseasonal-02  0.237159   0.113764   2.085   0.0371 *
## MA-Seasonal-04     0.738180   0.092061   8.018 1.07e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## X11 adj.  ARIMA: (0 2 2)(0 1 1)  Obs.: 77  Transform: log
## AICc: 257.1, BIC: 274.6  QS (no seasonality in final):    0
## Box-Ljung (no autocorr.): 27.65   Shapiro (normality): 0.9656 *
final(PIBTSA)
##          Qtr1     Qtr2     Qtr3     Qtr4
## 1996 100.0209 100.3977 104.7880 103.3985
## 1997 103.8248 104.8589 106.6537 107.1998
## 1998 104.3627 106.8397 107.1789 105.7283
## 1999 105.5060 105.9662 106.7138 108.0361
## 2000 109.5622 110.6865 111.6985 112.9109
## 2001 113.2017 113.0359 112.1191 112.2626
## 2002 114.3034 115.0038 116.8060 118.2141
## 2003 116.6877 116.6858 117.6694 119.1089
## 2004 121.2424 123.6196 125.3305 126.5121
## 2005 126.5316 128.4487 128.1028 129.3265
## 2006 130.9385 132.2814 134.0533 135.5494
## 2007 138.0144 140.4844 141.9269 144.2904
## 2008 146.7550 149.0367 151.7199 145.5084
## 2009 142.3680 146.2933 149.7956 153.1179
## 2010 155.9768 158.3498 160.1333 162.0052
## 2011 163.3604 166.4313 165.5810 166.1209
## 2012 166.5143 167.2685 169.3465 169.9862
## 2013 171.2721 173.3229 173.4634 173.5229
## 2014 174.7563 172.3464 172.5598 173.0335
## 2015 172.7567
resid(PIBTSA)
##               Qtr1          Qtr2          Qtr3          Qtr4
## 1996 -7.631667e-04 -2.201111e-03  3.141728e-03  2.574970e-03
## 1997  3.174729e-03  6.212903e-03  9.204532e-03 -7.465839e-03
## 1998 -9.728653e-03 -6.811626e-03 -5.964617e-03 -1.765721e-02
## 1999  4.723447e-05  4.385868e-03 -3.390666e-03  1.649382e-02
## 2000  9.088939e-03  4.091440e-03 -5.358513e-03  1.031379e-02
## 2001 -5.678862e-03 -6.629682e-03 -1.881175e-02  4.415687e-03
## 2002  1.588737e-02 -2.421547e-03  9.973459e-03  4.824630e-03
## 2003 -2.221758e-02  1.505260e-03 -1.293940e-04  6.493953e-03
## 2004  1.705166e-02  9.154115e-03 -1.592008e-03 -1.634295e-03
## 2005 -3.946405e-03  4.692298e-03 -1.941391e-02  5.489389e-03
## 2006  1.143963e-02 -7.772136e-03  4.263996e-03  3.790999e-04
## 2007  1.489767e-02 -3.828167e-03 -2.912046e-03  7.137433e-03
## 2008  5.731716e-03 -7.618234e-03  6.515553e-03 -2.574858e-03
## 2009 -6.548068e-03 -1.033460e-02  6.591033e-03  6.130106e-03
## 2010  3.502380e-03 -5.098606e-03 -3.128757e-03 -3.658016e-03
## 2011 -4.866644e-03  3.833883e-03 -2.047175e-02 -5.243749e-03
## 2012 -5.751823e-03 -7.500150e-03  9.345213e-03 -8.699906e-03
## 2013  4.248774e-03  8.307355e-04 -6.875811e-03 -4.441129e-03
## 2014  4.045059e-03 -2.492166e-02  5.948884e-03  3.909982e-04
## 2015 -4.024602e-03
coef(PIBTSA)
##          TC1996.3          AO1998.1          LS2008.4          TC2009.1
##       0.036471448      -0.025802666      -0.053941263      -0.028206039
##        Easter[15] MA-Nonseasonal-01 MA-Nonseasonal-02    MA-Seasonal-04
##      -0.006105615       0.635040284       0.237159253       0.738180493
fivebestmdl(PIBTSA)
## Error in seas(x = PIBT, transform.function = "log", regression.aictest = NULL, : object 'AMOSTRA' not found
spc(PIBTSA)
## series{
##   title = "PIBT"
##   file = "C:\Users\Tho\AppData\Local\Temp\RtmpgZIDc2/x13/data.dta"
##   format = "datevalue"
##   period = 4
##   span = (,2015.1)
## }
##
## transform{
##   function = log
##   print = aictransform
## }
##
## regression{
##   variables = (TC1996.3 AO1998.1 LS2008.4 TC2009.1 Easter[15])
## }
##
## arima{
##   model = (0 2 2)(0 1 1)
## }
##
## x11{
##   mode = mult
##   save = (d10 d11 d12 d13 d16 e18)
## }
##
## estimate{
##   save = (model estimates lkstats residuals)
## }
##
## spectrum{
##   print = qs
## }

Construindo o output em Excel, sendo que a primeira linha define o vetor de datas

Data &amp;amp;amp;amp;lt;- seq(as.POSIXlt("1996-3-1"), as.POSIXlt("2015-3-1"), by = "quarter")
PIBTSA <- as.numeric(series(PIBTSA, "d11"))
PIBTSA <- as.data.frame(PIBTSA)
PIBT <- as.numeric(PIBT)
PIBT <- as.data.frame(PIBT)

Matriz <- cbind(Data, PIBT, PIBTSA)

# Tire o "#" de comentário na linha abaixo para gravar a série dessazonalizada na nova aba "Dessaz"
# writeWorksheetToFile("Tab_Compl_CNT_1T15.xls", Matriz, sheet = "Dessaz", startCol = 1, clearSheets = TRUE)

Dessazonalizando o PIB brasileiro pelo R – Parte I

A desciclopédia (desciclopedia.org), segundo própria definição do site, é uma enciclopédia cheia de desinformações e mentiras grotescas – uma espécie de paródia do portal colaborativo Wikipedia. Embora controverso por conter uma série de artigos politicamente incorretos, o site traz também algumas definições sarcásticas, mas que em alguns casos até se aproximam da realidade.
Quando pesquisamos pelo termo “Economista”, além de descobrirmos que Mick Jagger e Arnold Schwarzenegger são formados em economia (o que é verdade!), também nos deparamos com a seguinte definição sobre a rotina deste profissional: “Seu tempo de trabalho é dividido da seguinte forma: 50% do tempo eles passam fazendo modelos sobre o futuro, e a outra metade explicando porque suas previsões não aconteceram”.

Na minha visão, a realidade não está muito distante do definido acima. Dedicamos boa parte do nosso tempo tratando bases de dados e construindo modelos abstratos para ajudar a entender como funcionam as complexas relações do mundo real. Em outras palavras, economistas usam modelos para simplificar a realidade e permitir um melhor entendimento do mundo; e os erros derivam justamente do fato de os modelos não serem capazes de representar 100% da realidade. Tanto que se fosse possível prever com 100% de confiança o exato valor do dólar um ano à frente, eu já estaria milionário.

Neste sentido, ao tratar os dados e aprimorar seus modelos, o economista está sempre buscando minimizar o erro de suas previsões; e uma das técnicas de tratamento de séries temporais mais utilizadas pelos economistas é o ajustamento sazonal. Isso porque este método nos permite ter uma ideia mais precisa do comportamento tendencial da série, limpando as flutuações sazonais, que podem ser de ordem natural (clima), cultural (dia das mães), religiosa (Natal) e de diversos outros tipos. A sazonalidade quase sempre aparece nas séries temporais, sobretudo nas de frequência mensal e trimestral; e a ocorrência desses eventos pode levar a conclusões inadequadas a respeito da série em estudo – por isso deve ser tratada. Um exemplo clássico é o aumento da oferta de emprego no final do ano para atender a maior demanda do período de festas, e no momento seguinte o fechamento destes postos de trabalho (temporários), o que causa diminuição do pessoal ocupado. Nesse caso, para o economista interessa separar o que ocorre periodicamente do que de fato está acontecendo com os dados observados.

Já que é tão importante, é possível fazer o ajustamento sazonal pelo R? Sim!
Existem vários métodos de ajustamento sazonal, sendo o X-13-ARIMA-SEATS, do U.S. Census Bureau e Bank of Spain, o mais atual e também um dos mais utilizados. Não é nosso objetivo aqui explorar a diferença entre os diversos métodos de dessazonalização; para mais detalhes o leitor pode consultar nota técnica do IBRE/FGV neste link.

No nosso exemplo, que será explorado na 2ª parte desse artigo, utilizaremos as séries do PIB Total (índice encadeado) e os parâmetros de dessazonalização pré-definidos pelo IBGE. Portanto, para prepararmos nossa receita, precisaremos da panela (o software R) e três ingredientes:

  1. O programa do X-13-ARIMA-SEATS, que será “chamado” pelo R para fazer a dessazonalização;
  2. A série histórica do PIB Total em bases trimestrais;
  3. Os parâmetros de dessazonalização pré-definidos pelo IBGE;