---
title: "Análise Box & Jenkins usando o BETS"
subtitle: "Minicurso do II Seminário Internacional de Estatística com R"
author: "Talitha Speranza, Pedro Ferreira e Jonatha Azevedo"
date: "28 de Maio de 2017"
output: html_document
---

```{r setup, echo = FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.width = 8, fig.height = 3.84, fig.align="center", fig.pos="H", prompt = TRUE)
```

```{r echo = FALSE}
library(BETS)
```

Vamos modelar a série de produção de bens intermediários (PBI) com o BETS usando a metodologia Box \& Jenkins. No final, seremos capazes de usar o modelo para fazer previsões.

# Metodologia Box \& Jenkins

O método de Box \& Jenkins permite que os valores futuros da série em estudo sejam previstos somente com base nos valores passados e presentes da mesma série, isto é, as previsões são feitas a partir de modelos univariados. 

Estes modelos são chamados SARIMA, uma sigla para o termo em inglês *Seasonal Auto-Regressive Integrated Moving Average*, e têm a forma:

\begin{equation}
\Phi_{P}(B)\phi_{p}(B)\nabla^{d}\nabla^{D}Z_t = \Theta_{Q}(B)\theta_{q}(B)a_t.
\end{equation}

onde

* $Z_t$ é a série em estudo 
* $a_t$ é um ruído branco 
* $\nabla^{d} = (1-B)^{d}$ é o operador de diferenças e $d$ o número de raizes unitárias 
* $\nabla^{D} = (1-B^{s})^{D}$ é o operador de diferenças na frequência sazonal $s$ e $D$ o número de raízes unitárias sazonais 
* $\phi_{p}(B)$ é o polinômio autorregressivo 
* $\Phi_{P}(B)$ é o polinômio autorregressivo sazonal 
* $\theta_{q}(B)$ é o polinômio de médias móveis 
* $\Theta_{Q}(B)$ é  o polinômio de médias móveis sazonal 

Em sua concepção original, que será adotada aqui, a metodologia de Box \& Jenkins se divide em três estágios iterativos:

(i) Identificação e seleção dos modelos: verificação da estacionariedade e da sazonalidade, com as devidas correções para os casos não estacionários e sazonais, e determinação das ordens dos polinômios descritos acima, fazendo uso das funções de autocorrelação (FAC) e autocorrelação parcial(FACP). 
(ii) Estimação dos parâmetros do modelo, usando máxima verossimilhança ou regressão dinâmica. 
(iii) Diagnóstico da conformidade do modelo, através de testes estatísticos.

Se o modelo não for aprovado na fase (iii), volta-se ao passo (i). Caso contrário, o modelo pode ser utilizado para fazer previsões. Na próxima seção, conforme o exemplo for evoluindo, cada um desses estágios será observado de perto e mais será dito sobre a metodologia.

# Modelagem

## Preliminares

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

```{r echo = c(1,2,4)}
# Busca em português pela série de produção de bens intermediários
results = BETSsearch(description = "'bens intermediarios'", lang = "pt", view = F)
results$last_value = NULL
results
```

Agora, carregamos a série através da função `BETSget` e guardamos alguns valores para, posteriormente, comparar com as previsões do modelo que será estabelecido. Também criaremos um gráfico (figura \ref{fig:chart}), pois ele ajuda a formular hipóteses sobre o comportamento do processo estocástico subjacente. 

```{r series}
# Obtenção da série de código 21864 (Produção de Bens Intermediários, IBGE)
data <- BETSget(21864)

# Guardar últimos valores para comparar com as previsões
data_test <- window(data, start = c(2015,11), end = c(2016,4), frequency = 12)
data <- window(data, start = c(2002,1), end = c(2015,10), frequency = 12)
```

```{r chart, fig.cap="Gráfico da série de produção de bens intermediários no Brasil.", echo = -1}
par(cex.axis = 0.75, cex.lab = 0.8, mar = c(5.1, 4.1, 0.5, 2.1))
# Gráfico da série
plot(data, main = "", col = "royalblue", ylab = "PBI (Número Índice)")
abline(v = seq(2002,2016,1), col = "gray60", lty = 3)
```

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 (figura \ref{fig:month}), que mostra o nível de produção por mês (a média é a linha tracejada). 

```{r month, fig.cap="Gráfico mensal da série em estudo.", echo = -1}
par(cex.axis = 0.75, cex.lab = 0.8, mar = c(5.1, 4.1, 0.5, 2.1))
# Gráfico mensal da série
monthplot(data, labels = month.abb, lty.base = 2, col = "red", 
                          ylab = "PBI (Número Índice)", xlab = "Month")
```

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.

## 1. Identificação

### 1.1. 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 `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 `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 é 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 `ur_test` aceita os mesmos parâmetros da `ur.df`, além do nível de significância desejado.

```{r adf, fig.cap= "Função de Autocorrelação dos Resíduos do Teste de Dickey-Fuller"}
df = ur_test(y = diff(data), type = "none", lags = 11, 
                                  selectlags = "BIC", level = "1pct")

# Exibir resultado dos testes
df$results
```

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 95%, 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.

```{r}
ns_roots = 0
d_ts = diff(data)

# 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
```

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.\footnote{A FAC foi feita com a função `BETS.corrgram`.}

```{r}
# Fazer FAC dos resíduos, com intervalo de confiança de 99%
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).

```{r}
library(forecast)

# Testes OCSB para raízes unitárias na frequencia sazonal
nsdiffs(data, test = "ocsb")
```

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.

### 1.2. 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 (figura \ref{fig:fac1}). 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 `corrgram`. Diferentemente de sua principal alternativa, a `Acf` do pacote `forecast`, a `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.

```{r fac1, fig.cap= "Função de Autocorrelação de $\\nabla Z_t$"}
# Correlograma de diff(data)
corrgram(diff(data), lag.max = 48, mode = "bartlett", style="plotly", knit = 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.

```{r fac2, fig.cap = "Função de Autocorrelação Parcial de $\\nabla Z_t$"}
# Função de autocorrelação parcial de diff(data)
corrgram(diff(data), lag.max = 36, type = "partial", style="plotly", knit = T)
```

A FAC da figura \ref{fig:fac1} e a FACP da figura \ref{fig:fac2} podem ter sido geradas por um processo `SARIMA(0,0,2) (1,0,0)`. Esta conjectura se baseia na observação de que as defasagens múltiplas de 12 parecem apresentar corte brusco na FACP a partir da segunda (isto é, a de número 24) e decaimento exponencial na FAC. Além disso, as duas primeiras defasagens da FAC parecem significativas, enquanto as demais, não. Há, ainda, alguma evidência de decaimento exponencial na FACP, exceto na frequência sazonal. Os dois últimos fatos indicam que o polinômio de médias móveis (não sazonal) pode ter ordem 2. Por estas razões, o primeiro modelo proposto para $Z_t$ será um `SARIMA(0,1,2)(1,0,0)[12]`. 

## 2. Estimação

Para estimar os coeficientes do modelo `SARIMA(0,1,2)(1,0,0)[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).

```{r}
# Estimacao dos parâmetros do modelo
model1 = Arima(data, order = c(0,1,2), seasonal = c(1,0,0))

# Teste t com os coeficientes estimados
# Nível de significância de 1%
t_test(model1, alpha = 0.01)
```

Concluímos pela coluna `Rej.H0` que os dois coeficientes do modelo, quando estimados por máxima verossimilhança, são estatisticamente significativos a 99% de confiança.

## 3. 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 (figura \ref{fig:stdr1}) será feito com o auxílio da função `BETS.std_resid`, que foi implementada especificamente para isso. 

```{r stdr1, echo = -2, fig.cap="Resíduos padronizados do primeiro modelo proposto", echo = -1}
par(cex.axis = 0.75, cex.lab = 0.8, mar = c(5.1, 4.1, 0.5, 2.1))
# Gráfico dos resíduos padronizados
resids = BETS.std_resid(model1, alpha = 0.01)

# Evidenciar outlier
points(2008 + 11/12, resids[84], col = "red")
```

Observamos que há um *outlier* proeminente e estatisticamente significativo em novembro de 2008. Este ponto corresponde à data da quebra estrutural que identificamos na figura \ref{fig:chart}. 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 `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.

```{r}
dummy = dummy(start = c(2002,1), end = c(2015,10), from = c(2008,9), to = c(2008,11))
dummy
```

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 (figura \ref{fig:stdr2}) também mostra que a inclusão de $D_t$ foi adequada, uma vez que não há mais evidência de quebra estrutural. 

```{r stdr2, echo = c(-1,-9), fig.cap="Resíduos padronizados do modelo proposto após a detecção de quebra estrutural"}
par(cex.axis = 0.75, cex.lab = 0.8, mar = c(5.1, 4.1, 0.5, 2.1))
# Estimacao dos parâmetros do modelo com a dummy 
model2 = Arima(data, order = c(0,1,2), seasonal = c(1,0,0), xreg = dummy)

# Teste t com os coeficientes estimados
# Nível de significância de 1%
t_test(model2, alpha = 0.01)

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

# Evidenciar novembro de 2008
points(2008 + 11/12, resids[84], col = "red")

# Mostrar BIC dos dois modelos estimados
model1$bic
model2$bic
```

Notamos, ainda, que o *Bayesian Information Criteria* (BIC) do modelo com a *dummy* é menor. Logo, também por este critério, o modelo com a *dummy* deve ser preferido ao anterior. 

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.

```{r lb, fig.cap = "Função de autocorrelação dos resíduos do modelo com a {\\it dummy}."}
# Teste de Ljung-Box nos resíduos do modelo com a dummy
boxt = Box.test(resid(model2), type = "Ljung-Box",lag = 11)
boxt

# Correlograma dos resíduos do modelo com a dummy
corrgram(resid(model2), lag.max = 20, mode = "bartlett", style = "normal")
```

O p-valor de `r  format(round(boxt$p.value,4),scientific = F)` 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 a figura \ref{fig:lb}. Concluímos, então, que o modelo foi bem especificado.

## 4. Previsões

O `BETS` fornece uma maneira conveniente para fazer previsões de modelos `SARIMA`. A função `predict` recebe os parâmetros da função `forecast` do pacote homônimo ou da `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 `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 `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.

```{r echo=-1, fig.cap="Gráfico das previsões do modelo SARIMA proposto."}
par(cex.axis = 0.75, cex.lab = 0.8, mar = c(5.1, 4.1, 0.5, 2.1))
new_dummy = dummy(start = start(data_test), end = end(data_test))

preds = predict(object = model2, xreg = new_dummy, 
                     actual = data_test, xlim = c(2012, 2016.2), ylab = "Milhões de Reais", 
                     style = "normal", legend.pos = "bottomleft")
```

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. 

```{r}
preds[['accuracy']]
```

o outro campo deste objeto, `'predictions'`, contém o objeto retornado pela `forecast` (ou pela `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.

# O uso da BETS.report para a modelagem SARIMA

A função `report` executa toda a modelagem Box \& Jenkins para qualquer conjunto de séries à escolha e gera relatórios com os resultados, como foi dito no início desta seção. Ela permite que as previsões feitas através do modelo sejam salvas em um arquivo de dados 

```{r eval = F}
parameters = list(
  lag.max = 48,
  n.ahead = 12 ) 

report(ts = 21864, parameters = parameters)
```

O resultado abre automaticamente, na forma de um arquivo `html`.


