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).

Classes Sociais, Renda e Clusters

Breve introdução

Se tem um negócio mal explicado é a questão da divisão em classes sociais no Brasil. A impressão que tenho é que cada um usa a metodologia que mais lhe interessa para explicar aquilo que quer.

O IBGE

O IBGE tem sua metodologia que divide as classes baseados no número de salários mínimos de um domicílio. Ok, bacana, se a renda é maior do que ‘X’ a família é classe C, se é maior do que ‘Y’ classe B e por aí vai. Até que é legal, mas utilizar a renda de uma família e não a per capita nos deixa com o problema Homer Simpson:

Problema Homer

ABEP

Esse aqui é bacana!!! Sabemos que a pesquisa do IBGE tem vários problemas, entre eles a imprecisão das suas medidas, então alguém teve a ideia de contar os bens de uma pessoa para definir sua classe social, a ideia até que é interessante, mas ao tentar resolver o problema criamos tantos outros que vai ser difícil descrever em um post só…

A tabela é mais ou menos assim:

ABEP

Terceirizando o problema

É óbvio que não temos uma solução para isso, então vamos tentar encarar esse problema de outra forma.

Imagine um alienígena que veio para nosso país, começou a observar a população brasileira e tomar notas para depois responder a seguinte pergunta: ‘Como subdividir a população brasileira em classes sociais’?

Em ‘machine learning’ de modo bem grosseiro podemos dividir os algoritmos em dois:

Supervisionados e não supervisionados.

No caso da tarefa que passamos para nosso colega de outro mundo, estamos falando de um problema não supervisionado, pois ninguém definiu para o alien as classes e como subdividi-las, nem o que é ou não mais importante nesse momento, simplesmente demos para o coitado o fardo de realizar tal classificação.

O modelo

Em poucas palavras, vamos tentar separar a população brasileira em grupos baseados simplesmente na distância entre essas pessoas. Então, imagine uma linha reta:

X——————————————–X

Agora, imagine que temos algumas pessoas distribuídas ao longo da mesma:

X—OOO——OOOO——OOO——X

Só de olhar para essa linha nós institivamente conseguimos separar em 3 grupos de pessoas. Agora imagine que ao invés de uma linha medida em metros ou quilômetros, estamos falando da Renda das famílias. Para famílias que ganham X, Y ou Z, a separação em grupos seria praticamente a mesma.

A diferença com o nosso modelo é que ao invés de calcularmos a distância em um espaço unidimensional (uma linha) ou bidimensional, calculamos a distância em um espaço n-dimensional.

Os dados

Para simular o problema do alienígena, vamos utilizar a PNAD de 2013, porém o triste é que com ela não é possível reproduzir o método da ABEP, já que não temos o número de vários itens, só uma variável binária dizendo se a família o possui ou não, então vamos jogar a PNAD praticamente inteira no modelo e criar nossa própria classificação híbrida, parte ABEP e parte IBGE.

Ao invés de contar pontos de item por item, criamos uma variável chamada BENS que soma 18 questões diferentes da PNAD para dar o número de bens das famílias, como TV, DVD, Fogão, Celular, 3/4G e etc. (A ABEP atribui pesos para cada item, nós não o faremos por uma questão de simplificação, mas um refinamento deste modelo poderá ter pesos diferentes.)

Além desta, utilizamos mais 8 variáveis da PNAD domicílios, e 03 variáveis construídas da PNAD pessoa que são ‘MAIS_INSTR’ que é igual ao número de anos de estudo da pessoa mais instruída do domicílio, número de filhos e número de empregados domésticos.

##  [1] "Total_de_moradores" "Comodos"            "Banheiros"         
##  [4] "n_TV_telafina"      "n_TV_Tubo"          "Automovel"         
##  [7] "Renda"              "N_FILHO"            "N_DOMESTICO"       
## [10] "INSTR"              "BENS"

Resultados – Iteração #1

Os resultados ficaram assim (PART é a participação da CLASSE no total):

## Source: local data frame [6 x 3]
## 
##   CLASSE PART      RENDA
## 1      A  0.2 84225.6966
## 2     B1  0.9 23496.5027
## 3     B2  3.0 12404.7848
## 4     C1  8.5  6536.8275
## 5     C2 21.2  3394.0439
## 6    D-E 66.1   899.4745

plot of chunk unnamed-chunk-4

Ficou bem diferente do que estamos acostumados a ver, acho que nosso amigo alienígena errou feio!

Por outro lado, uma informação bem interessante que é possível derivar deste gráfico é que o brasil é um país extremamente desigual e a renda média da classe ‘A1’ acabou ficando em torno de $84k por mês!

We are the 99% – Iteração #2

A desigualdade no Brasil é tão grande, mas tão grande que enviesou completamente o modelo, porém até o documento da ABEP alerta sobre essa desigualdade indicando que para casos onde a renda é superior a R$30.000 o modelo de classes não ajuda muito.

Para resolver isso, vamos adotar o parâmetro ‘occupy Wall Street’ no nosso modelo e simplesmente ignorar o 1% mais rico da população, afinal, eles já venceram a corrida dos ratos.

Essas famílias são automaticamente classe ‘A’, da mesma forma que o pessoal de renda zero é classe ‘E’, porém como a PNAD contém várias omissões e erros, vamos ignorar o pessoal de renda Zero também.

## Source: local data frame [6 x 3]
## 
##   CLASSE  PART     RENDA
## 1      A  4.14 12657.533
## 2     B1  5.09  6973.986
## 3     B2  8.94  4720.220
## 4     C1 12.61  3274.278
## 5     C2 19.16  2245.416
## 6    D-E 50.06  1060.438

plot of chunk unnamed-chunk-6

Melhorou bastante! Agora está até lembrando as classificações mais tradicionais, apesar que ainda nosso modelo coloca mais pessoas nas classes D e E do que o que estamos acostumados, isso pode ser um problema da amostra ou político, afinal é bem conveniente classificar a maioria da população como vivendo no ‘meio do bolo’ para dar uma perpétua sensação de mobilidade social, porém os dados pintam um quadro diferente.

Ficando mais técnico – Controle da variância

Sabemos que as variáveis utilizadas no modelo podem ter impactos diferentes na classe social, por isso na hora de montar o modelo de clusterização (k-means) decidimos deixar as variáveis no estado mais puro o possível (direto da base), porém algumas variáveis podem roubar o show em um modelo de clusterização.

Vamos pegar um exemplo de livro texto.

plot of chunk unnamed-chunk-7

Lindo, certo? Em um espaço bidimensional, a relação entre as variáveis X e Y praticamente implora por uma clusterização, note que o modelo até comete alguns erros classificando os ‘O’s e os ‘+’s mas ele acabou fazendo um ótimo trabalho na clusterização no final das contas.

Agora vamos olhar outro exemplo

plot of chunk unnamed-chunk-8

Ao olhar para este segundo exemplo conseguimos até separar em dois clusters, mas essa separação visual é feita praticamente só com a informação do eixo ‘Y’, isso faz com que Y enviese o resultado por uma condição matemática e não real.

Apesar de controlarmos por várias variáveis no nosso modelo, isso pode acontecer, por exemplo com o número de empregados domésticos que vai de 0 a 2 no máximo, versus a renda que vai de 0 a mais de 150.000.

Pré-processamento

Para tentar tratar esse problema vamos utilizar técnicas básicas de pré-processamento centralizando e colocando as variáveis na mesma escala, dessa forma todas as variáveis do nosso modelo ficarão com média ‘0’ e variância ‘1’.

##                   VAR MEDIA VARIANCIA
## 1  Total_de_moradores     0         1
## 2             Comodos     0         1
## 3           Banheiros     0         1
## 4       n_TV_telafina     0         1
## 5           n_TV_Tubo     0         1
## 6           Automovel     0         1
## 7               Renda     0         1
## 8             N_FILHO     0         1
## 9         N_DOMESTICO     0         1
## 10              INSTR     0         1
## 11               BENS     0         1

Variáveis ajustadas – Iteração #3

## Source: local data frame [6 x 3]
## 
##   CLASSE  PART    RENDA
## 1      A  6.71 7656.665
## 2     B1 10.29 3998.186
## 3     B2 18.88 3000.812
## 4     C1 12.63 1976.285
## 5     C2 19.08 1849.001
## 6    D-E 32.40 1483.811

plot of chunk unnamed-chunk-11

Esse novo modelo ficou com uma distribuição melhor ainda que o anterior, mas o gráfico de Renda x Classe acabou ficando com vários outliers, indicando que agora com as variáveis ajustadas a renda não tem uma importância tão grande quanto tinha antes. É bem possível que a influência das outras variáveis na classe social seja maior. Ótimo, era mais ou menos essa a intenção quando misturamos as metodologias do IBGE e ABEP.

Análise dos outliers

Vamos olhar para a renda para tentar entender o comportamento desses outliers.

plot of chunk unnamed-chunk-12

Até que no final das contas eles estão bem-comportados.

Próximos passos

Até que para uma análise de fim de semana o modelo ficou bem interessante. Alguns caminhos que podem ser seguidos daqui para frente é tentar analisar com esse modelo a evolução do conceito de classe social ao longo dos anos e tentar derivar algum conceito de mobilidade social.

Um exemplo disso é que de acordo com a metodologia da ABEP, se uma família comprar um carro ela marca mais pontos e pode ‘subir de classe’, porém se todo mundo fizer isso as famílias continuam semelhantes entre elas, o que não justifica uma mudança de classe social.

Uma forma interessante de analisar a ‘Tragédia dos comuns’…

A PNAD, os DVDS e as Florestas (Aleatórias)

Motivação

Em jun/15, em uma conversa de ‘boteco’ alguém sugeriu que eu comprasse um DVD da Galinha pintadinha para minha filha recém-nascida e a primeira coisa que pensei era: – ‘Ótimo, aí eu penduro o DVD no móbile dela’, pois não tenho onde reproduzir esse troço…

Isso me lembrou a época de ouro das locadoras de vídeo, onde era praticamente um passeio de família ir a locadora pegar um ou dois filmes para assistir no fim de semana e gigantes como a Blockbuster quando chegaram ao Brasil mudaram a cara deste mercado para sempre, da locadora de porta de garagem para algo grandioso, cheio de opções e principalmente sem filas de espera…

Porém no final de 2013 foi oficialmente declarada a morte da Blockbuster link, uma morte mais do que esperada, mas ainda assim o fim de uma era…

…Com isso, vamos fazer uma viagem no tempo, pegar os dados da PNAD 2013 e tentar responder a seguinte pergunta:

“Com tantas maneiras de assistir filmes, e as ruínas das locadoras espalhadas pelo país, o que leva uma família a ter um aparelho de DVD (em 2013)?”

Vamos tentar montar um modelo para explicar isso…

A PNAD

Para quem gosta de trabalhar com dados, o IBGE é um prato cheio. A PNAD é uma pesquisa divertidíssima onde temos uma bela amostra de dados (longe de ser perfeita, mas boa para algumas análises) e perguntas o suficiente para a realização das mais variadas análises.

Estávamos passando pelos dicionários da pesquisa de 2013 recentemente e encontramos a seguinte questão:

  • ‘Tem aparelho de DVD’?

Os dados

Não vamos nos preocupar em explicar como extraímos os dados da PNAD, as instruções estão no site do IBGE neste link e como isso por si só já dá um belo trabalho, talvez um dia no futuro escreveremos um post sobre como fazer isso.

A última PNAD anual foi realizada em 2013, então é essa pesquisa que vamos utilizar. Atualmente existe a PNAD contínua, mas para o efeito que queremos capturar nesse exercício a de 2013 está mais do que ótima.

Factor engineering

Essa parte é relativamente trabalhosa, pois temos que:

  1. Transformar variáveis numéricas em fatores. Olhar para o dicionário da PNAD é um exercício de paciência. Por exemplo:

‘Tem televisão em cores?’

2 = SIM
4 = NÃO

‘Tem aparelho de DVD?’

1 = SIM
3 = NÃO

PORQUEEEEEEEEEEEE?????

É muito difícil normalizar isso? Que tal 1=sim e 0=não? Ou, para facilitar, 2 = sim, 4 = não….sempre!

  1. Criar algumas variáveis que vão ajudar na nossa modelagem, como por exemplo ‘Acesso_INTERNET’ que combina as questões de banda larga e discada em uma variável
  2. Vamos usar os dados do questionário dos domicílios, mas no dicionário das pessoas vamos criar variáveis como EH_FILHO que é igual ao número de filhos de uma família, EH_FILHO_6 que é igual a quantidade de filhos menores do que 6 anos e FAIXA_IDADE_FILHO que é a média de idade dos filhos de um domicílio

Os datasets

A primeira coisa a fazer é criar os datasets de treino e teste, só que temos um problema, considerando que a amostra da PNAD domicílios é de mais de 148.000 observações onde 116.000 são válidas e cada um tem seu peso, na teoria teríamos que repetir cada domicílio um número de vezes igual ao seu peso, ou seja, se o peso médio é 560, temos mais de 64.000.000 domicílios.

O dataset do nosso modelo vai consistir em 100.000 observações selecionadas aleatoriamente do universo de 64.000.000

DICA DO DIA: Nunca tente criar um vetor de 64.000.000 de observações no R e selecionar 100.000 linhas aleatoriamente!!!

Depois de travar o PC mudamos de estratégia…

Como?

Vamos assumir um universo de 3 domicílios com pesos 1.000.000, 50.000.000 e 13.000.000, criamos um vetor de pesos acumulados, 1.000.000, 51.000.000 e 64.000.000, fazemos 100.000 sorteios aleatórios de 1 a 64.000.000 e simplesmente contamos o número de itens no vetor de pesos acumulados menores do que o sorteio mais 1, por exemplo:

Sorteio: 10; número < 0; posição = 1
Sorteio: 1.100.000; numero < 1; posição = 2
Sorteio: 64.000.000; numero < 2; posição = 3

Isso faz com que a probabilidade se selecionarmos um domicílio seja igual seu peso no total.

Treino e teste

Agora que temos as 100.000 observações aleatórias, vamos dividir em dois datasets, um para treinar o nosso modelo e outro para testar seus resultados, faremos isso com a razão de 70% treino e 30% teste

 

O modelo

Para o modelo temos um problema de classificação, ou seja, para um certo grupo de inputs queremos uma resposta 0 ou 1 para classificar um domicílio como ‘Tem DVD’ ou ‘Não tem DVD’, esses inputs são algumas variáveis socioeconômicas dos domicílios e dos seus habitantes.

library(randomForest)
modelFit <- randomForest(Tem_DVD ~.,
                         data=train,
                         do.trace=T, 
                         ntree=200,
                         importance=TRUE,
                         mtry=10)
plot(modelFit,main='randomForest error rate')

plot of chunk unnamed-chunk-4

O pessoal costuma ir para o Overkill rodando milhares de árvores, mas a verdade é que não precisamos de tantas para estabilizar o erro, além disso demora uma eternidade! É possível ver pelo gráfico acima que depois de 100 árvores o erro já está mais do que controlado e o ganho é marginal.

Mostrando o pau

Agora que matamos a cobra, vamos aplicar nosso modelo no dataset de testes e ver se o mesmo ficou aceitável ou se somente capturou o ruído do dataset de treino.

pred <- predict(modelFit,newdata=test)
CM <- confusionMatrix(test$Tem_DVD, pred)
CM
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     0     1
##          0  5573  2835
##          1  1578 20013
##                                           
##                Accuracy : 0.8529          
##                  95% CI : (0.8488, 0.8569)
##     No Information Rate : 0.7616          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6179          
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.7793          
##             Specificity : 0.8759          
##          Pos Pred Value : 0.6628          
##          Neg Pred Value : 0.9269          
##              Prevalence : 0.2384          
##          Detection Rate : 0.1858          
##    Detection Prevalence : 0.2803          
##       Balanced Accuracy : 0.8276          
##                                           
##        'Positive' Class : 0               
## 

TMI [Too much information] – Mais para frente faremos um post explicando detalhes do que é cada um desses resultados, por enquanto vamos focar em duas coisas:

Accuracy : 0.8529

O modelo tem uma precisão de 85%, ou seja, com base em variáveis socioeconômicas é possível prever com 85% de precisão se um domicílio tem ou não DVD. O intervalo de confiança é bem pequeno (95% CI : (0.8488, 0.8569)), o que é bom.

CM$table
##           Reference
## Prediction     0     1
##          0  5573  2835
##          1  1578 20013

Temos dois tipos de erro, quando a referência diz que um domicilio tem um DVD e o modelo diz que não e quando a referência = não e o modelo = sim. Como diria um colega, se o modelo diz algo e a realidade diz outro é a realidade que tem um problema…

Importância das variáveis

plot of chunk unnamed-chunk-7

Um negócio bacana das florestas é poder observar a importância das variáveis e para esse modelo Renda vence de longe como sendo a mais importante, depois a UF (afinal estamos no Brasil e Renda+UF explicam praticamente tudo!) e só depois a média de idade dos filhos.

Uma coisa interessante é que as variáveis de filhos que criamos para este modelo (EH_FILHO = número de filhos no domicílio e EH_FILHO_6 = Número de filhos com 6 anos ou menos) não fazem muita diferença no modelo…

Pelo visto ninguém vai comprar um aparelho de DVD só pela Galinha Pintadinha…

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;