BETS package

O pacote BETS (sigla para Brazilian Economic Time Series) para o R fornece, de maneira descomplicada, as mais relevantes séries temporais econômicas do Brasil e diversas ferramentas para analisá-las. O BETS preenche uma lacuna no processo de obtenção de dados no Brasil, na medida em que unifica os pontos de acesso às inúmeras séries e oferece uma interface bastante simples, flexível e robusta.

As séries presentes na base de dados do pacote são produzidas por três importantes e respeitados centros: o Banco central do Brasil (BACEN), o Instituto Brasileiro de Geografia e Estatística (IBGE) e o Instituto Brasileiro de Economia da Fundação Getúlio Vargas (FGV/IBRE e FGV dados). Em sua concepção original, o BETS apenas deveria reunir em um só lugar o maior número possível de séries destes centros, dada a dificuldade com que o pesquisador poderia se defrontar para obtê-las. Este objetivo foi cumprido e hoje o pacote disponibiliza mais de 18 mil séries econômicas brasileiras.

Por conta do vasto tamanho das bases de dados, o BETS foi expandido para prover mecanismos que auxiliam o analista na pesquisa, extração e exportação das séries. A partir daí, a formulação do pacote se transformou e o BETS caminhou para se tornar um ambiente integrado de análise e aprendizado. Foram incorporadas diversas funcionalidades ausentes do universo R, com a intenção de facilitar ainda mais a modelagem e a interpretação das séries fornecidas. Ademais, algumas funções já incluem a opção de gerar saídas inteiramente didáticas, incentivando e auxiliando o usuário a ampliar seus conhecimentos

Neste documento, discutiremos algumas formas básicas de uso do pacote, mostrando algumas de suas principais funções.

Instalando o BETS

versão CRAN

install.packages("BETS")
library(BETS)

versão Github

install.packages("devtools")
require(devtools)
install_github("pedrocostaferreira/BETS",force = TRUE)
require(BETS)

BETS.search

Devido ao tamanho considerável da base de dados, foi necessário criar um modo de pesquisar por séries a partir de seus metadados, isto é, uma ferramenta de busca que utilizasse uma ou mais informações das séries como palavras-chave. A função BETS.search realiza as pesquisas por cada campo da tabela de metadados.

O protótipo da BETS.search tem a forma:

BETS.search(description, src, periodicity, unit, code, view = TRUE, lang = "en")

Onde os argumentos recebem, respectivamente:

  • description: um character. String de busca com os termos que devem ou não estar presentes na descrição da série desejada.
  • src: um . A fonte dos dados.
  • periodicity um character. A frequência na qual a série é observada.
  • unit: um character. A unidade na qual os dados foram medidos.
  • code: um integer. O código único da série na base do BETS.
  • view: um boolean. Por padrão, TRUE. Se FALSE, os resultados serão mostrados direto no console do R.
  • lang: um character. Idioma da pesquisa. Por padrão, 'en', para inglês. Também é possivel fazer a pesquisa em português, bastando mudar o valor para 'pt'.

Para refinar as buscas, há regras de sintaxe para o parâmetro description:

  1. Para procurar palavras alternativas, separe-as por espaços em branco. Exemplo: description = 'núcleo ipca' significa que a descrição da série deve conter 'ipca' e 'núcleo'.

  2. Para procurar expressões inteiras, basta cercá-las com ' '. Exemplo: description = 'índice de 'núcleo ipca'' significa que deve conter na descrição da série 'núcleo ipca' e 'índice'.

  3. Para excluir palavras da busca, insira um \(\sim\) antes de cada um delas. Exemplo: description = 'ipca ~ núcleo' significa que a descrição da série deve conter 'ipca' e não pode conter 'núcleo'.

  4. Para excluir todas as expressões da busca, de forma semelhante ao item anteiror, basta cercá-los com ' ' e inserir um \(\sim\) antes de cada uma delas. Exemplo: description = '~ índice 'núcleo ipca'' significa que a descrição da série deve conter 'índice' e não pode conter 'núcleo ipca'.

  5. É possível pesquisar ou negar várias palavras ou expressões, desde que sejam respeitadas as regras precedentes.

  6. O espaço em branco após o sinal de negação (\(\sim\)) não é necessário. Mas os espaços em branco depois de expressões ou palavras são necessários.

Alguns exemplos de uso podem ser visto abaixo:

require(BETS)
BETS.search(description = "sales")
BETS.search(description = "sales ~ retail") # quero excluir retail (varejo)
BETS.search(description = "'sales volume index' ~ vehicles")
BETS.search(description = "distrito federal", periodicity = 'A', src = 'IBGE')
library(BETS)
BETS.search(description = "gdp accumulated", unit = "US", view = F)
## 
## BETS-package: Found 2 out of 16854 time series.
##   code                                            description
## 1 4192 GDP accumulated in the last 12 months - in US$ million
## 2 4386           GDP accumulated in the year - in US$ million
##            unit periodicity      start last_value    source
## 1 US$ (million)           M 31/01/1990   may/2017 BCB-Depec
## 2 US$ (million)           M 31/01/1990   may/2017 BCB-Depec
results = BETS.search(description = "consumption ~ 'seasonally adjusted' ~ private", view = F)
## 
## BETS-package: Found 56 out of 16854 time series.
head(results)
##   code                                           description
## 1 1393          Petroleum derivatives consumption - Gasoline
## 2 1394               Petroleum derivatives consumption - GLP
## 3 1395          Petroleum derivatives consumption - Fuel oil
## 4 1396        Petroleum derivatives consumption - Diesel oil
## 5 1397 Petroleum derivatives consumption - Other derivatives
## 6 1398 Petroleum derivatives consumption - Total derivatives
##                     unit periodicity      start last_value source
## 1 Barrels/day (thousand)           M 31/01/1979   apr/2017    ANP
## 2 Barrels/day (thousand)           M 31/01/1979   apr/2017    ANP
## 3 Barrels/day (thousand)           M 31/01/1979   apr/2017    ANP
## 4 Barrels/day (thousand)           M 31/01/1979   apr/2017    ANP
## 5 Barrels/day (thousand)           M 31/01/1979   apr/2017    ANP
## 6 Barrels/day (thousand)           M 31/01/1979   apr/2017    ANP

Para mais informações sobre a BETS.search, incluindo os valores válidos em cada campo, consulte o manual de referência, digitando ?BETS.search no console do R.

BETS.get

A BETS.get funciona unicamente através do código de referência da série, obtido com as consultas feita com a . Sua assinatura é:

BETS.get(code, data.frame = FALSE)

O parâmetro code é, obviamente, obrigatório. O argumento opcional data.frame representa o tipo do objeto que será retornado. Por padrão, seu valor é FALSE, indicando que o objeto devolvido pela função será um ts (time series). Caso data.frame = TRUE, a série será armazenada em um objeto do tipo data.frame.

Vamos extrair duas das séries pesquisadas anteriormente.

# Obter a serie do PIB acumulado em 12 meses, em dolares
library(BETS)
gdp_accum = BETS.get(4192)
window(gdp_accum, start = c(2014,1))
##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug
## 2014 2472429 2481077 2481608 2479220 2480118 2475031 2471294 2465255
## 2015 2403557 2348867 2300552 2245795 2189218 2140418 2084812 2028876
## 2016 1791221 1790785 1787186 1785780 1784608 1786830 1785172 1788233
## 2017 1816633 1834842 1857456 1875310 1898474 1917529 1939361 1962439
##          Sep     Oct     Nov     Dec
## 2014 2465302 2462106 2457281 2454846
## 2015 1969996 1910297 1853223 1796168
## 2016 1790145 1789188 1792866 1797234
## 2017 1983900 2007129 2030697 2054244
# Obter a serie do PIB do Distrito Federal, a precos de mercado
gdp_df = BETS.get(23992, data.frame = T)
head(gdp_df)
##         date       value
## 1 2002-01-01 53902199799
## 2 2003-01-01 58456124319
## 3 2004-01-01 67076505202
## 4 2005-01-01 75732681210
## 5 2006-01-01 84661405538
## 6 2007-01-01 93404000766

BETS.save

Para conferir versatilidade às formas de armazenamento das séries do , há a possibilidade de criar arquivos com as séries em formatos proprietários, isto é, formatos que pertencem a softwares pagos.

A BETS.save extrai a série temporal da base de dados do pacote na forma de um data.frame e cria um arquivo no formato especificado. No arquivo, há uma tabela onde a primeira coluna conterá as datas e a segunda, os dados.

A função possui três variações:

BETS.save.sas(code, data = NULL, file.name = "series")
BETS.save.spss(code, data = NULL, file.name = "series")
BETS.save.stata(code, data = NULL, file.name = "series")

Novamente, o parâmetro code recebe o código da série. O usuário pode fornecer sua própria série através do argumento data, que pode ser um data.frame ou um ts. Não é necessário acrescentar a extensão ao nome do arquivo no parãmetro file.name.

Alguns exemplos típicos de uso seriam:

# Salvar uma série qualquer no formato do SPSS
myseries <- BETS.get(code = 2078)
BETS.save.spss(data = myseries, file.name = "series_spss")

BETS.chart

A BETS.chart foi inicialmente projetada para ser uma função privada, auxiliar da BETS.dashboard. No entanto, pensamos ser de grande valia para o usuário dispor de um meio de obter os gráficos dos dashboards separadamente, de modo a poder incorporá-los em seus trabalhos.

O protótipo da BETS.chart é o que se segue:

BETS.chart(ts, file = NULL, open = TRUE, params = NULL)

O parâmetro ts pode receber uma dentre as várias opções pré-definidas de gráficos ou uma série do usuário. Há também a opção de salvar a saída no working directory, definindo o nome do arquivo file. Caso o arquivo deva ser aberto após a criação, open deve ser mantido como TRUE. O parâmetro params é reservado para gráficos de séries do usuário. É uma lista que pode conter o campo codace, que recebe um booleano e indica se devem ser desenhadas áreas sombreadas representando recessões datadas pelo CODACE (FGV/IBRE), e o campo start, que especifica qual deve ser a data de início da série. Uma vez que se tratam de gráficos de conjuntura, a data de fim não pode ser alterada e é sempre o último dado disponível.

Vejamos dois exemplos de uso da BETS.chart.

BETS.chart(ts = 'iie_br', file = "iie_br", open = TRUE)
library(BETS)
?BETS.chart
# misery_index = 13522 plus 24369
BETS.chart(ts = "misery_index", file = "misery_index.png", open = TRUE)
infla <- window(BETS.get(13522), start= c(2000,1))
BETS.chart(ts = infla)

Para uma lista completa dos gráficos diponíveis, consulte o manual de referência da BETS.chart.

BETS.dashboard

O BETS possui uma poderosa ferramenta de análise de conjuntura, os dashboards. Para criar um dashboard, chamamos a BETS.dashboard, cuja assinatura é:

BETS.dashboard(type = "business_cycle", saveas = NA)

Assim, obtemos um arquivo .pdf. No exemplo, o usuário escolhe salvar o arquivo com o nome de survey.pdf. Os gráficos que são exibidos também são de escolha do usuário, selecionados através do parâmetro charts (por default, seu valor é 'all'). O manual de referência possui uma lista completa dos gráficos disponíveis.

BETS.dashboard(type = "business_cycle", saveas = "survey.pdf")
Terceira página do Business Cycle Dashboard



Modelando a Produção de Bens Intermediários, IBGE

Preliminares

O primeiro passo é encontrar a série PBI na base de dados do BETS. Isso pode ser feito com a função BETS.search. O comando e sua saída são mostrados abaixo.

# Busca em português pela série de produção de bens intermediários
results = BETS.search(description = "'bens intermediarios'", lang = "pt", view = F)
## 
## BETS-package: Found 5 out of 21953 time series.
results
##    code
## 1  1334
## 2 11068
## 3 21864
## 4 25302
## 5 25328
##                                                                       description
## 1 Indicadores da produção (1991=100) - Por categoria de uso - Bens intermediários
## 2 Indicadores da produção (2002=100) - Por categoria de uso - Bens intermediários
## 3                        Indicadores da produção (2012=100) - Bens intermediários
## 4                                        Importações - Bens intermediários (CGCE)
## 5                                   Importações (kg) - Bens intermediários (CGCE)
##     unit periodicity      start     source
## 1 Índice           M 31/01/1975       IBGE
## 2 Índice           M 31/01/1991       IBGE
## 3 Índice           M 01/01/2002       IBGE
## 4    US$           M 01/01/1997 MDIC/Secex
## 5     kg           M 01/01/1997 MDIC/Secex

Agora, carregamos a série através da função BETS.get e guardamos alguns valores para, posteriormente, comparar com as previsões do modelo que será estabelecido. Também criaremos um gráfico, pois ele ajuda a formular hipóteses sobre o comportamento do processo estocástico subjacente.

# Obtenção da série de código 21864 (Produção de Bens Intermediários, IBGE)
library(BETS)
data <- BETS.get(21864)

# Guardar últimos valores para comparar com as previsões
data_test <- window(data, start = c(2015,11), frequency = 12)
data_training <- window(data, start = c(2002,1), end = c(2015,10), frequency = 12)
# Gráfico da série
plot(data_training, main = "", col = "royalblue", ylab = "PBI (Número Índice)")
abline(v = seq(2002,2016,1), col = "gray60", lty = 3)
Gráfico da série de produção de bens intermediários no Brasil.

Gráfico da série de produção de bens intermediários no Brasil.

Quatro características ficam evidentes. Primeiramente, trata-se de uma série homocedástica e sazonal na frequência mensal. Este último fato é corroborado pelo gráfico mensal da série, que mostra o nível de produção por mês (a média é a linha tracejada).

# Gráfico mensal da série
monthplot(data_training, labels = month.abb, lty.base = 2, col = "red", 
                          ylab = "PBI (Número Índice)", xlab = "Month")
Gráfico mensal da série em estudo.

Gráfico mensal da série em estudo.

Um terceiro aspecto marcante da série é a quebra estrutural em novembro de 2008, quando ocorreu a crise financeira internacional e a confiança dos investidores despencou. A quebra impactou diretamente na quarta característica importante da série: a tendência. Incialmente, a tendência era claramente crescente, mas não explosiva. A partir de novembro de 2008, porém, parece que o nível da série se manteve constante ou até mesmo descresceu. Em um primeiro momento, a quebra estrutural será desconsiderada na estimação dos modelos, mas logo o benefício de levá-la em conta explicitamente ficará claro.

A seguir, criaremos um modelo para a série escolhida de acordo com os passos definidos anteriormente.

Identificação

Testes para Estacionariedade

Esta subseção trata de um passo crucial na abordagem de Box & Jenkins: a determinação da existência e da quantidade total de raízes unitárias no polinômio autorregressivo não-sazonal e sazonal do modelo. De posse desses resultados, obtemos uma série estacionária através da diferenciação da série original. Assim, poderemos identificar a ordem dos parâmetros através da FAC e FACP, pois isso deve feito através de séries estacionárias de segunda ordem.

A função BETS.ur_test executa o teste Augmented Dickey Fuller (ADF). Ela foi construída em cima da função ur.df do pacote urca, que é instalado juntamente com o BETS. A vantagem da BETS.ur_test é a saída, desenhada para que o usuário visualize rapidamente o resultado do teste e tenha todas as informações de que realmente necessita. Trata-se de um objeto com dois campos: uma tabela mostrando as estatísticas de teste, os valores críticos e se a hipótese nula [presença de raíz unitária - ST não estacionária] é rejeitada ou não, e um vetor contendo os resíduos da equação do teste. Esta equação é mostrada abaixo.

\[\begin{equation} \Delta y_t = \phi + \tau_{1} t + \tau_{2} y_{t-1} + \delta_1 \Delta y_{t-1} + \cdots +\delta_{p-1} \Delta y_{t-p+1} + \varepsilon_t \end{equation}\]

As estatísticas de teste da tabela do objeto de saída se referem aos coeficientes \(\phi\) (média ou drift), \(\tau_{1}\) (tendência determinística) e \(\tau_{2}\) (raiz unitária). A inclusão da média e da tendência determinística é opcional. Para controlar os parâmetros do teste, a BETS.ur_test aceita os mesmos parâmetros da ur.df, além do nível de significância desejado.

df = BETS.ur_test(y = data_training, type = "none", lags = 11, 
                                  selectlags = "BIC", level = "1pct")

# Exibir resultado dos testes
df$results
##      statistic crit.val rej.H0
## tau1  0.951481    -2.58     no

Portanto, para a série em nível, observa-se que não se pode rejeitar a hipotése nula de existência de uma raiz unitária ao nível de confiança de 99%, pois a estatística de teste é maior do que o valor crítico. Agora, iremos aplicar a função diff à série repetidas vezes e verificar se a série diferenciada possui uma raiz unitária.

ns_roots = 0
d_ts = diff(data_training,lag = 1,differences = 1) #ST com 1 diferença na parte não sazonal

# Loop de testes de Dickey-Fuller. 
# A execução é interrompida quando não for possível rejeitar a hipótese nula
while(df$results[1,"statistic"]> df$results[1,"crit.val"]){
    ns_roots = ns_roots + 1
    d_ts = diff(d_ts)
    df = BETS.ur_test(y = d_ts, type = "none", lags = 11, 
                                selectlags = "BIC", level = "1pct")
 }

ns_roots
## [1] 1

Logo, para a série em primeira diferença, rejeita-se a hipótese nula de que há raiz unitária a 5% de significância. A FAC dos resíduos da equação do teste evidencia que ele foi bem especificado, pois a autocorrelação não é significativa até a décima primeira defasagem.

# Fazer FAC dos resíduos, com intervalo de confiança de 99%
BETS.corrgram(df$residuals,ci=0.99,style="normal",lag.max = 11)

Um outro pacote bastante útil que é instalado com o BETS é o forecast. Usaremos a função nsdiffs deste pacote para realizar o teste de Osborn-Chui-Smith-Birchenhall e identificar raízes unitárias na frequência sazonal (em nosso caso, mensal).

library(forecast)
# Testes OCSB para raízes unitárias na frequencia sazonal
d_ts <- diff(data_training,lag = 1,differences = 1)
nsdiffs(d_ts, test = "ocsb")
## [1] 0

Infelizmente, a nsdiffs não fornece nenhuma outra informação sobre o resultado do teste além do número de diferenças sazonais que devem ser tiradas para eliminar as raizes unitárias. Para o caso da série em análise, o programa indica que não há raiz unitária mensal, não sendo necessária, portanto, diferenças nesta frequência.

Funções de Autocorrelação

As conclusões anteriores são corroboradas pela função de autocorrelação da série em primeira diferença. Ela mostra que autocorrelações estatisticamente significativas, isto é, fora do intervalo de confiança, não são persistentes para defasagens múltiplas de 12 e no entorno destas, indicado a ausência da raiz unitária sazonal.

As funções do BETS que utilizamos para desenhar correlogramas é a BETS.corrgram. Diferentemente de sua principal alternativa, a ACF do pacote forecast, a BETS.corrgram retorna uma gráfico atraente e oferece a opção de calcular os intervalos de confiança de acordo com a fórmula proposta por Bartlett. Sua maior vantagem, contudo, não pôde ser exibida aqui, pois depende de recursos em flash. Caso o parâmetro style seja definido como 'plotly', o gráfico torna-se interativo e mostra todos os valores de interesse (autocorrelações, defasagens e intervalos de confiança) com a passagem do mouse, além de oferecer opções de zoom, pan e para salvar o gráfico no formato png.

# Correlograma de diff(data)
BETS.corrgram(d_ts, lag.max = 48, mode = "bartlett", style="plotly", knit = T)

Função de Autocorrelação de \(\nabla Z_t\)

O correlograma acima ainda não é suficiente para determinamos um modelo para a série. Faremos, então, o gráfico da função de autocorrelação parcial (FACP) de \(\nabla Z_t\). A BETS.corrgram também pode ser utilizada para este fim.

# Função de autocorrelação parcial de diff(data)
BETS.corrgram(d_ts, lag.max = 36, type = "partial", style="plotly", knit = T)

Função de Autocorrelação Parcial de \(\nabla Z_t\)

A FAC e a FACP podem ter sido geradas por um processo SARIMA(2,1,2) (1,0,2). Minha conjectura:

1. FAC (parte não sazonal): corte brusco a partir do lag 2 (inclusive), indica uma parte MA(2);
2. FAC (parte sazonal): corte brusco no lag 24 (inclusive), indica SMA(2);
3. FACP (parte não sazonal): corte brusco a partir do lag 2 (inclusive), indica uma parte AR(2);
4. FACP (parte sazonal): corte brusco no lag 12 (inclusive), indica SAR(1);

Seguindo a lógica retratada acima, estimaremos, primeiramente, um SARIMA(2,1,2) (1,0,2)[12]. Obviamente que iremos aumentar e diminuir o número de parâmetros para ver se esse é realmente o modelo que melhor descreve a série temporal de bens intermediários.

Estimação

Para estimar os coeficientes do modelo SARIMA(2,1,2)(1,0,2)[12], será aplicada a função Arima do pacote forecast. Os testes t serão feitos através da função BETS.t_test do BETS, que recebe um objeto do tipo arima ou Arima, o número de variáveis exógenas do modelo e o nível de significância desejado, devolvendo um data.frame contendo as informações do teste e do modelo (coeficientes estimados, erros padrão, estatísticas de teste, valores críticos e resultados dos testes).

# Estimacao dos parâmetros do modelo
library(forecast)
model1 = Arima(data_training, order = c(2,1,2), seasonal = c(1,0,2), lambda = 0) # lambda = 0 -> transf log

# Teste t com os coeficientes estimados
# Nível de significância de 5%
BETS.t_test(model1, alpha = 0.05)
##          Coeffs  Std.Errors          t Crit.Values Rej.H0
## ar1  -1.1540895 0.005584597 206.655819    1.975092   TRUE
## ar2  -0.9995708 0.001153483 866.567245    1.975092   TRUE
## ma1   1.1403843 0.024074586  47.368801    1.975092   TRUE
## ma2   0.9992786 0.030151110  33.142348    1.975092   TRUE
## sar1  0.9974117 0.002720581 366.617113    1.975092   TRUE
## sma1 -0.9338877 0.096264508   9.701268    1.975092   TRUE
## sma2  0.1069412 0.091467703   1.169169    1.975092  FALSE

Duas considerações importante:

1. rejeitamos H0 do SMA2
2. os nossos coeficientes estão muito próximos de 1, precisamos reduzir o modelo.
# Estimacao dos parâmetros do novo modelo
library(forecast)
model2 = Arima(data_training, order = c(2,1,2), seasonal = c(1,0,1), lambda = 0) # lambda = 0 -> transf log

# Teste t com os coeficientes estimados
# Nível de significância de 5%
BETS.t_test(model2, alpha = 0.05)
##          Coeffs  Std.Errors         t Crit.Values Rej.H0
## ar1  -1.1531750 0.005909417 195.14194    1.974996   TRUE
## ar2  -0.9994390 0.001446706 690.83764    1.974996   TRUE
## ma1   1.1389815 0.023878833  47.69838    1.974996   TRUE
## ma2   0.9995717 0.028426457  35.16343    1.974996   TRUE
## sar1  0.9987779 0.002087312 478.49962    1.974996   TRUE
## sma1 -0.8804195 0.095863133   9.18413    1.974996   TRUE

Os parâmetros do modelo 2 são todos significantes, mas os parâmetros do modelo esão muito próximos de 1, podendo indicar não estacionariedade do modelo. Para não ficar rodando infinitos modelos, irei usar a função ‘auto.arima’ para me ajudar.
Observe o resíduo do ‘auto.arima’.

# Estimacao dos parâmetros do modelo
auto <- auto.arima(data_training,lambda = 0)
BETS.corrgram(resid(auto), lag.max = 20, mode = "bartlett", style = "normal")


Temos um parâmetro significante no lag 12. Precisamos aumentar a ordem do modelo para tentarmos corrigir esse problema. Irei inserir um lag sazonal na parte MA.

# Estimacao dos parâmetros do modelo

model3 = Arima(data_training, order = c(1,1,0), seasonal = c(1,0,1), lambda = 0) # lambda = 0 -> transf log
model3
## Series: data_training 
## ARIMA(1,1,0)(1,0,1)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1  sar1     sma1
##       -0.1379     1  -0.7906
## s.e.   0.0756     0   0.0010
## 
## sigma^2 estimated as 0.0006067:  log likelihood=365.95
## AIC=-723.9   AICc=-723.65   BIC=-711.48
# Teste t com os coeficientes estimados
# Nível de significância de 5%
BETS.t_test(model3, alpha = 0.05)
##          Coeffs   Std.Errors            t Crit.Values Rej.H0
## ar1  -0.1378790 7.561336e-02 1.823473e+00    1.974716  FALSE
## sar1  0.9999979 4.068826e-09 2.457706e+08    1.974716   TRUE
## sma1 -0.7905844 9.592666e-04 8.241550e+02    1.974716   TRUE
comp_CI <- data.frame(model1$aicc,model2$aicc,model3$aicc)
comp_CI
##   model1.aicc model2.aicc model3.aicc
## 1   -727.3294   -728.3703   -723.6543

observações:

a. olhando para os critérios de informação (AICc, neste caso) escolheríamos o modelo 2, contudo, como os parâmetros estão muito próximos de 1, por experiência, escolheria o modelo 3. Dessa forma, farei todo o diagnóstico e previsão a partir do modelo 3, mas no final iremos comparar as previsões;
b. o parâmetro AR(1) foi não significante, mas como temos alguns meses com outliers, irei solucionar esse problema e observarei se esse é realmente um parâmetro não significante.

Testes de Diagnóstico

O objetivo dos testes de diagnóstico é verificar se o modelo escolhido é adequado. Neste trabalho, duas conhecidas ferramentas serão empregadas: a análise dos resídos padronizados e o teste de Llung-Box.

O gráfico dos resíduos padronizados será feito com o auxílio da função BETS.std_resid, que foi implementada especificamente para isso.

# Gráfico dos resíduos padronizados
resids = BETS.std_resid(model3, alpha = 0.01)

# Evidenciar outlier
points(2008 + 11/12, resids[84], col = "red")
Resíduos padronizados do primeiro modelo proposto

Resíduos padronizados do primeiro modelo proposto

Observamos que há um outlier proeminente e estatisticamente significativo em novembro de 2008. Este ponto corresponde à data da quebra estrutural que identificamos no início do exercício. Portanto, foi proposto um segundo modelo, que inclui uma dummy definida como se segue:

\[\begin{equation} D_{t} = \begin{cases} & 0 \text{, } t < \text{setembro de 2008} \\ & 1 \text{, } \text{setembro de 2008} <= t <= \text{novembro de 2008} \\ & 0 \text{, } t > \text{novembro de 2008} \end{cases} \end{equation}\]

Esta dummy pode ser criada com a função BETS.dummy, como mostramos abaixo. Os parâmetros start e end indicam o início e o fim do período coberto pela dummy, que nada mais é que uma série temporal cujos valores podem ser apenas 0 ou 1. Os campos from e to indicam o intervalo em que a dummy deve assumir valor 1.

head(data_training)
##       Jan  Feb  Mar  Apr  May  Jun
## 2002 77.2 74.8 83.3 84.3 87.0 85.3
tail(data_training)
##        May   Jun   Jul   Aug   Sep   Oct
## 2015  96.3  97.1  99.0 102.2  96.0  98.9
dummy = BETS.dummy(start = c(2002,1), end = c(2015,10), from = c(2008,9), to = c(2008,11))
dummy
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2002   0   0   0   0   0   0   0   0   0   0   0   0
## 2003   0   0   0   0   0   0   0   0   0   0   0   0
## 2004   0   0   0   0   0   0   0   0   0   0   0   0
## 2005   0   0   0   0   0   0   0   0   0   0   0   0
## 2006   0   0   0   0   0   0   0   0   0   0   0   0
## 2007   0   0   0   0   0   0   0   0   0   0   0   0
## 2008   0   0   0   0   0   0   0   0   1   1   1   0
## 2009   0   0   0   0   0   0   0   0   0   0   0   0
## 2010   0   0   0   0   0   0   0   0   0   0   0   0
## 2011   0   0   0   0   0   0   0   0   0   0   0   0
## 2012   0   0   0   0   0   0   0   0   0   0   0   0
## 2013   0   0   0   0   0   0   0   0   0   0   0   0
## 2014   0   0   0   0   0   0   0   0   0   0   0   0
## 2015   0   0   0   0   0   0   0   0   0   0

Como podemos ver nos resultados da execução do trecho de código a seguir, a estimação deste modelo através de máxima verossimilhança resultou em coeficientes estatisticamente diferentes de 0 ao nível de significância de 5%, inclusive para a dummy. O gráfico dos resíduos padronizados dos valores ajustados pelo novo modelo também mostra que a inclusão de \(D_t\) foi adequada, uma vez que não há mais evidência de quebra estrutural.

# Estimacao dos parâmetros do modelo com a dummy 
model4 = Arima(data_training, order = c(1,1,0), seasonal = c(1,0,1),xreg = dummy, lambda = 0) # lambda = 0 -> transf log
model4
## Series: data_training 
## Regression with ARIMA(1,1,0)(1,0,1)[12] errors 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1    sar1     sma1   dummy
##       -0.2190  0.9989  -0.8985  0.0521
## s.e.   0.0762  0.0027   0.1217  0.0163
## 
## sigma^2 estimated as 0.0005747:  log likelihood=364.37
## AIC=-718.75   AICc=-718.37   BIC=-703.22
# Teste t com os coeficientes estimados
# Nível de significância de 5%
BETS.t_test(model4, alpha = 0.05)
##            Coeffs  Std.Errors          t Crit.Values Rej.H0
## ar1   -0.21895438 0.076192215   2.873711    1.974716   TRUE
## sar1   0.99892817 0.002672497 373.780875    1.974716   TRUE
## sma1  -0.89849281 0.121654525   7.385609    1.974716   TRUE
## dummy  0.05205229 0.016301160   3.193165    1.974716   TRUE
resids = BETS.std_resid(model4, alpha = 0.01)

# Evidenciar novembro de 2008
points(2008 + 11/12, resids[84], col = "red")
Resíduos padronizados do modelo proposto após a detecção de quebra estrutural

Resíduos padronizados do modelo proposto após a detecção de quebra estrutural 

comp_CI <- data.frame(model1$aicc,model2$aicc,model3$aicc,model4$aicc)
comp_CI
##   model1.aicc model2.aicc model3.aicc model4.aicc
## 1   -727.3294   -728.3703   -723.6543   -718.3691

Notamos, que pelo AICc escolheríamos o modelo 2, mas, por parcimônia, ficaremos com o modelo 4.

O teste de Ljung-Box para o modelo escolhido pode ser executado através da função Box.test do pacote stats. Para confirmar o resultado dos testes, fazemos os correlogramas dos resíduos e vemos se há algum padrão de autocorrelação.

# Teste de Ljung-Box nos resíduos do modelo com a dummy
boxt = Box.test(resid(model4), type = "Ljung-Box",lag = 11)
boxt
## 
##  Box-Ljung test
## 
## data:  resid(model4)
## X-squared = 30.765, df = 11, p-value = 0.0012
# Correlograma dos resíduos do modelo com a dummy
BETS.corrgram(resid(model4), lag.max = 20, mode = "bartlett", style = "normal")
Função de autocorrelação dos resíduos do modelo com a {\it dummy}.

Função de autocorrelação dos resíduos do modelo com a {}.

O p-valor de 0.0012 indica que há grande probabilidade de a hipótese nula de que não há autocorrelação nos resíduos não seja rejeitada. Parece ser o caso, como mostra acima. Concluímos, então, que o modelo foi bem especificado.

Previsões

O BETS fornece uma maneira conveniente para fazer previsões de modelos SARIMA. A função BETS.predict recebe os parâmetros da função forecast do pacote homônimo ou da BETS.grnn.test (a ser tratada adiante, no segundo estudo de caso) e devolve não apenas os objetos contendo as informações da previsão, mas também um gráfico da série com os valores preditos. Essa visualização é importante para que se tenha uma ideia mais completa da adequação do modelo. Opcionalmente, podem também ser mostrados os valores efetivos do período de previsão.

Chamaremos a BETS.predict para gerar as previsões do modelo proposto. Os parâmetros object (objeto do tipo arima ou Arima), h (horizonte de previsão) e xreg (a dummy para o período de previsão) são herdados da função forecast. Os demais são da própia BETS.predict, sendo todos parâmetros do gráfico, com exceção de actual, os valores efetivos da série no período de previsão.

new_dummy = BETS.dummy(start = start(data_test), end = end(data_test))

preds = BETS.predict(object = model4, xreg = new_dummy, 
                     actual = data_test, xlim = c(2012, 2016.2), ylab = "PBI (Número Índice)", 
                     style = "normal", legend.pos = "bottomleft")
Gráfico das previsões do modelo SARIMA proposto.

Gráfico das previsões do modelo SARIMA proposto.

As áreas em azul em torno da previsão são os intervalos de confiança de 85% (azul escuro) e 95% (azul claro). Parece que a aderência das previsões foi satisfatória. Para dar mais significado a esta afirmação, podemos verificar várias medidas de ajuste acessando o campo 'accuracy' do objeto retornado.

preds[['accuracy']]
##                ME     RMSE     MAE      MPE    MAPE      ACF1 Theil's U
## Test set -5.64955 5.924782 5.64955 -6.60687 6.60687 0.5034323  1.267422
preds2 = BETS.predict(object = model2, 
                     actual = data_test, xlim = c(2012, 2016.2), ylab = "PBI (Número Índice)", 
                     style = "normal", legend.pos = "bottomleft")

preds2[['accuracy']]
##                 ME     RMSE      MAE       MPE     MAPE      ACF1
## Test set -5.300856 5.616893 5.300856 -6.216425 6.216425 0.4403434
##          Theil's U
## Test set  1.206521

Observe agora o MAPE do modelo 2 que tem muito mais parâmetros. Veja que temos um erro de previsão muito similar mas utilizando um modelo muito mais parcimonioso.

preds2 = BETS.predict(object = model2, 
                     actual = data_test, xlim = c(2012, 2016.2), ylab = "PBI (Número Índice)", 
                     style = "normal", legend.pos = "bottomleft")

preds2[['accuracy']]
##                 ME     RMSE      MAE       MPE     MAPE      ACF1
## Test set -5.300856 5.616893 5.300856 -6.216425 6.216425 0.4403434
##          Theil's U
## Test set  1.206521

O outro campo deste objeto, 'predictions', contém o objeto retornado pela forecast (ou pela BETS.grnn.test, se for o caso). Na realidade, este campo ainda conta com um dado adicional: os erros de previsão, caso sejam fornecidos os valores efetivos da série no período de previsão.

Materiais de Apoio

  1. Introduction to R for Data Science: este é um curso introdutório e irá ajuda-lo a entender conceitos básicos de programação em R;

  2. Vídeos sobre o R (FGV/IBRE | NMEC): vídeos em português produzidos pelo nosso time da FGV que facilitarão o entendimento de conceitos básicos de programação em R;

  3. Khan Academy: ideal para aprender conceitos básicos de matemática e estatística;

  4. Forecasting using R: este é um excelente curso do Rob Hyndman que poderá lhe ajudar a compreender um pouco mais o mundo de séries temporais.

Prof. Dr. Pedro Costa Ferreira


Doutor em Engenharia Elétrica - (Decision Support Methods) e Mestre em Economia. Co-autor dos livros “Planejamento da Operação de Sistemas Hidrotérmicos no Brasil” e “Análise de Séries Temporais em R: curso introdutório”. É o primeiro e único pesquisador da América Latina a ser recomendado pela empresa RStudio Inc.


Atuou em projetos de Pesquisa e Desenvolvimento (P&D) no setor elétrico nas empresas Light S.A. (e.g. estudo de contingências judiciais), Cemig S.A, Duke Energy S.A, entre outras. Atuou como consultor em Big Data e Data Science nas empresas, Coca-Cola Brasil, Light SA, Duratex, ONS, entre outras. Ministrou cursos de estatística e séries temporais na PUC-Rio e IBMEC e em empresas como o Operador Nacional do Setor Elétrico (ONS), Petrobras e CPFL S.A.


Atualmente é professor de Econometria de Séries Temporais e Estatística, cientista chefe do Núcleo de Métodos Estatísticos e Computacionais (FGV|IBRE), coordenador do curso Big Data e Data Science (FGV|IDE) e sócio-diretor da empresa Model Thinking Br (MTBr). É também revisor de importantes journals, como Energy Policy e Journal of Applied Statistics. Principais estudos são em modelos Econométricos, Incerteza Econômica, Preços, R software e Business Analytics [e.g detecção de fraudes; HR analytics].


contatos:
email: pedro.guilherme@fgv.br
website: pedrocostaferreira.github.io
GitHub: github.com/pedrocostaferreira
Linkedin: linkedin.com/pedro-costa-ferreira