varejoMS_arima

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

---
title: "Séries Temporais com R: Análise ARIMA do Consumo do Varejo em MS"
author: "Adriano Marcos Rodrigues Figueiredo, *e-mail: adriano.figueiredo@ufms.br*"
linkcolor: blue
abstract: 
 This is an exercise for class use. We analyse data on Retail consumption, from January/2000 to nowadays. 
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:
 html_document:
 code_download: true
 theme: default
 number_sections: true
 toc: yes
 toc_float: yes
 df_print: paged
 fig_caption: true
 pdf_document:
 toc: yes
---
```{r knitr_init, echo=FALSE, cache=FALSE}
library(knitr)
library(rmarkdown)
library(rmdformats)
## Global options
options(max.print="100")
opts_chunk$set(echo=TRUE,
	 cache=TRUE,
 prompt=FALSE,
 tidy=TRUE,
 comment=NA,
 message=FALSE,
 warning=FALSE)
opts_knit$set(width=100)
```
Licença {-#Licença}
===================
This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit <http://creativecommons.org/licenses/by-sa/4.0/> or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
![License: CC BY-SA 4.0](https://mirrors.creativecommons.org/presskit/buttons/88x31/png/by-sa.png){ width=25% }
Citação {-#Citação}
===================================
Sugestão de citação:
FIGUEIREDO, Adriano Marcos Rodrigues. Séries Temporais com R: Análise do Consumo do Varejo em MS. Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em <http://rpubs.com/amrofi/arima_varejoms>. 
Introdução
===================
Neste arquivo utilizo a série do Índice de volume de vendas no varejo Total de Mato Grosso do Sul, série mensal a partir de jan/2000 até jul/2019 obtida com o pacote `BETS` e importada do Banco Central do Brasil. Portanto, são 235 observações mensais.
Conforme Hyndman e Athanasopoulos (2018), se combinarmos a diferenciação com autoregressão e um modelo de média móvel, obteremos um modelo ARIMA não sazonal. ARIMA é um acrônimo para Média Móvel Integrada AutoRegressiva (neste contexto, “integração” é o inverso da diferenciação). O modelo completo pode ser escrito como
\[
\begin{equation}
 y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p}
 + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t}, \tag{1.1}
\end{equation}
\]
 
onde $y'_{t}$ é a série diferenciada (pode ter sido diferenciada mais de uma vez). Os "preditores" no lado direito incluem valores defasados de $y'_{t}$ assim como erros defasados. Chamamos isso de modelo ARIMA(p,d,q), em que p é a ordem da parte autoregressiva (AR(p)); d é a ordem de diferenciação envolvida (I(d)); q é a ordem da parte da média móvel (MA(q)).
Dados
===================
Farei de duas formas para o leitor. Uma carrega direto do site do Banco Central do Brasil com o pacote `BETS` (FERREIRA, SPERANZA e COSTA, 2018) e a outra eu gerei a estrutura idêntica pela função `dput()` para os leitores que não conseguirem por qualquer motivo o acesso ao site do Banco Central (as vezes vejo isso ocorrer dependendo dos bloqueios da sua rede de internet). A forma pelo dput assume o nome varejoms2 enquanto a extraída pelo BETS tem nome varejoms. Esclareço ao leitor que após baixar a série pelo BETS, fiz o dput e a partir de então, desabilitei o bloco (`Chunk`) que acessa o BETS apenas para agilizar os cálculos.
```
library(BETS)
# Pegando as séries a partir do site do Banco Central do Brasil
# Índice de volume de vendas no varejo Total de Mato Grosso do Sul
# mensal a partir de jan/2000 até jul/2019 
# 235 observações mensais
varejoms <- BETSget(1479) 
print(varejoms)
class(varejoms)
dput(varejoms) # opção para ter os dados como na structure abaixo
```
```{r warning=FALSE, comment=FALSE}
suppressMessages(library(readxl));suppressMessages(library(foreign))
suppressMessages(library(dynlm));suppressMessages(library(car))
suppressMessages(library(lmtest));suppressMessages(library(sandwich))
suppressMessages(library(fpp2));suppressMessages(library(tseries))
suppressMessages(library(zoo));suppressMessages(library(forecast))
suppressMessages(library(ggplot2))
# Pegando as séries a partir do site do Banco Central do Brasil
# Índice de volume de vendas no varejo Total de Mato Grosso do Sul
# mensal a partir de jan/2000 até jul/2019 
# Loading the dataset
# library(readxl)
# dados <- read_excel("dados.xlsx",sheet = "dados")
# attach(dados)
library(BETS)
# Pegando as séries a partir do site do Banco Central do Brasil
# Índice de volume de vendas no varejo Total de Mato Grosso do Sul
# mensal a partir de jan/2000 até jul/2019 (em 12.09.2019)
#varejoms <- BETSget(1479) 
varejoms<- structure(c(35.2, 35.6, 39.2, 40.5, 41.6, 40.4, 40.8, 38.7, 37.3, 
37.6, 35.6, 47.4, 34.2, 32.2, 38, 37.5, 38.8, 35, 38.4, 40.2, 
38.2, 39.3, 36, 46.3, 36.4, 34, 39, 37.8, 38.9, 35.3, 37.2, 38.1, 
35.6, 38.3, 35.6, 45.7, 32.2, 31.6, 35.2, 36.8, 37.5, 34.8, 38.4, 
38.1, 37, 39, 37.3, 49, 35.8, 35.3, 40.2, 41.3, 43.9, 41.6, 45.9, 
42.3, 42.2, 44, 41.3, 56.9, 38.5, 38.5, 45.3, 43.6, 46.2, 44.3, 
47.5, 46.5, 46.4, 46, 44.1, 61, 42, 40.2, 44.6, 44.5, 47.8, 45.3, 
46.5, 48.5, 47.7, 50.2, 49.3, 64.5, 47, 46.8, 51, 50.5, 55, 51.3, 
52.8, 55.3, 54.8, 55.6, 55.3, 72.2, 54.5, 52.1, 56.2, 57.2, 60.8, 
56.1, 61.8, 61.6, 59.8, 63.3, 57.7, 77.4, 61.4, 51.9, 57.3, 57.9, 
61.9, 57.3, 61.1, 61.1, 60.6, 65.6, 63.5, 83.1, 64.1, 60.2, 67.8, 
67.1, 72.8, 68.5, 71.1, 69.3, 69.9, 71.1, 67.9, 92.7, 67.5, 64.8, 
69.1, 69.4, 79.6, 70.2, 73.8, 72.5, 71.3, 75.6, 74.7, 100.8, 
79.5, 75.7, 82.4, 78, 84.8, 83.2, 84.8, 88.5, 86.3, 91.7, 92.8, 
111.4, 92.8, 83.7, 92.5, 88.3, 93.9, 88.8, 96, 95.9, 93.2, 98.3, 
100.5, 128.8, 97.2, 90.2, 94.3, 94.4, 101.1, 92, 96.4, 98.2, 
97.6, 105.8, 103.1, 129.6, 99.6, 87.8, 97, 94.8, 98.6, 93.4, 
98.4, 96.4, 92.5, 100.6, 97.2, 124.5, 91.5, 85.1, 91.6, 88.5, 
92.2, 87.4, 90.5, 88.1, 85.2, 89.4, 93.4, 116.9, 90.8, 84, 89.7, 
86.3, 90, 87.3, 90.8, 93.5, 93.7, 91.4, 93.5, 114.1, 87.8, 81.1, 
94.5, 83.2, 89.9, 88.8, 89.3, 93.7, 93.5, 96.3, 101.3, 118.3, 
93.8, 85.2, 90, 86.6, 90, 85.2, 90.9), .Tsp = c(2000, 2019.5, 
12), class = "ts")
```
A rotina de dados obtidos pelo BETS já retorna a série em formato `ts`, ou seja, série temporal. Farei então a criação de uma série em diferenças para observar o comportamento da série em nível e em diferenças. 
Inicialmente olharei as estatísticas descritivas da série. Em seguida farei um plot básico da série e o plot pelo pacote `dygraphs`, útil para ver os pontos de picos e momentos específicos.
```{r}
dvarejo<-diff(varejoms)
# estatisticas basicas
summary(varejoms)
# Min. 1st Qu. Median Mean 3rd Qu. Max. 
# 31.60 44.05 64.50 67.28 90.00 129.60 
# plot basico
# lembrar que em class(), ele já indicou que era ts = serie temporal
plot(varejoms) 
# pelo pacote dygraph dá mais opções
library(dygraphs)
dygraph(varejoms, main = "Índice de volume de vendas no varejo total de Mato Grosso do Sul <br> (Mensal) (2011=100) BCB 1479") %>%
 dyAxis("x", drawGrid = TRUE) %>%
 dyEvent("2005-1-01", "2005", labelLoc = "bottom") %>%
 dyEvent("2015-1-01", "2015", labelLoc = "bottom") %>%
 dyEvent("2018-1-01", "2018", labelLoc = "bottom") %>%
 dyEvent("2019-1-01", "2019", labelLoc = "bottom") %>%
 dyOptions(drawPoints = TRUE, pointSize = 2)
```
É possivel visualizar nos plots acima: sazonalidade (por exemplo, picos em dezembro de cada ano); a tendência aparentemente crescente até 2014 e decresce com a "crise" brasileira; e uma aparente não-estacionariedade (média e variância mudam no tempo). Mais a frente, aplicarei o teste de raiz unitária na série para avaliar a estacionariedade de modo mais explícito.
Análise da série
======================
Uma ressalva deve ser feita, que no presente exercício, não farei a divisão entre amostra teste e amostra
treino, de modo que usarei a série toda para os ajustes. O leitor deve em geral fazer estas divisões para certificar de que o modelo é um bom preditor. 
Geralmente, não é possível dizer olhando o gráfico de tempo, quais valores de lags p e q serão apropriados para os dados. Às vezes, é possível usar o gráfico ACF e o gráfico PACF para determinar valores apropriados para p e q. A primeira análise indica olhar os gráficos (correlogramas) da Função de Autocorrelação (FAC, em português, ou ACF em inglês) e Autocorrelação parcial (FACp, em português, ou PACF em inglês). 
um gráfico FAC mostra as autocorrelações que medem a relação entre 
$y_t$ e $y_{t-k}$ para diferentes valores de k. Agora, se $y_t$ e $y_{t-1}$ estão correlacionados, então $y_{t}$ e $y_{t-2}$ também deve ser correlacionado. Mas isso pode ser porque ambos estão correlacionados com $y_{t-1}$. Para evitar essas "interferências" é que se faz a FACp. As autocorrelações parciais na FACp medem a relação entre $y_t$ e $y_{t-k}$ depois de remover os efeitos dos lags 1, 2, 3, ... , (k-1). Portanto, a primeira autocorrelação parcial é idêntica à primeira autocorrelação, porque não há nada entre elas para remover. Cada autocorrelação parcial pode ser estimada como o último coeficiente em um modelo autoregressivo. 
Isto pode ser obtido com as funções `acf` e `pacf`, ou usando as funções `ggAcf` e `ggPacf` do pacote `forecast`. Faremos o mesmo procedimento para a série em nível e para a série em primeira diferença.
## Função de Autocorrelação (FAC) e Autocorrelação parcial (FACp) com defasagem 36
### Série em nível
Usarei a rotina do Hyndman e Athanasopoulos (2018). 
```{r}
varejo<-varejoms
varejo %>% ggtsdisplay(main="")
```
Neste caso, o `ggtsdisplay` do pacote `forecast` já retorna os gráficos da série, e respectivas ACF e PACF.
### Série em primeira diferença
Como a série apresenta variações sazonais importantes, assim como tendência importante, claramente não-estacionárias, vou olhar também em primeira diferença.
```{r}
varejo %>% diff() %>% ggtsdisplay(main="Série varejo de MS em primeira diferença")
```
### Primeira diferença sazonal e ACF e PACF
Farei agora a diferença sazonal.
```{r}
varejo %>% diff(lag=12) %>% ggtsdisplay(main="Série varejo de MS em primeira diferença sazonal") 
```
### Ajuste sazonal e ACF e PACF
Aplicaremos o ajuste sazonal tipo STL (Seasonal Decomposition of Time Series by Loess) aos dados.
```{r}
library(fpp2)
varejo %>% stl(s.window='periodic') %>% seasadj() -> varejoadj
autoplot(varejoadj)
```
E agora os plots de ACF e PACF na série ajustada sazonalmente.
```{r}
varejoadj %>% diff() %>% ggtsdisplay(main="Série varejo de MS em primeira diferença e ajuste sazonal")
```
Parece indicar para algum AR de ordem 4 na PACF, alguma sazonalidade marcante na ACF (indicando um termo SMA), e indicacao de um ARIMA (4,1,1)(0,1,1)[12].
Uma estimação inicial desse modelo daria algo como 
```{r}
(fit <- Arima(varejoadj, order=c(4,1,1),seasonal = c(0, 1, 1)))
```
É possível observar os AICc = 1046.03, suspeitar da não significância de ar3 e olhar os resíduos. A função `autoplot` retorna o círculo unitário e as raízes unitárias desta estimação, enquanto a `checkresiduals` fornece os correlogramas residuais.
O círculo unitário indica estabilidade das raízes unitárias do modelo, enquanto o Ljung-Box retorna a não-rejeição de H0 de que os parâmetros de autocorrelação residual são nulos (retorna um grande valor p, sugerindo também que os resíduos são ruído branco). O gráfico ACF dos resíduos do modelo ARIMA (4,1,1)(0,1,1)[12] mostra que todas as autocorrelações estão dentro dos limites à exceção de um lag 23, indicando que os resíduos estão se comportando como ruído branco.
```{r}
autoplot(fit)
checkresiduals(fit)
library(FitAR) # pratico para plotar graficos de Q
#LBQPlot(res, lag.max = 30)
### res é a série que temos interesse de realizar os testes
### lag.max é o número de lags que se deseja imprimir o gráfico
LBQPlot(residuals(fit),36)
# checar normalidade dos resíduos
library(tseries)
jarque.bera.test(residuals(fit))
```
ARIMA 
======================
Uma estimação automatizada pela função `auto.arima` {forecast} padrão indica um modelo diferente.
```{r}
fit2 <- auto.arima(varejo)
summary(fit2)
#p-values dos coeficientes
pvalues<-(1-pnorm(abs(fit2$coef)/sqrt(diag(fit2$var.coef))))*2
pvalues
library(FitAR) # pratico para plotar graficos de Q
LBQPlot(residuals(fit2),36)
# checar normalidade dos resíduos
library(tseries)
jarque.bera.test(residuals(fit2))
```
Neste caso, o modelo teve a indicação de ARIMA(0,1,1)(0,1,1)[12], diferente do anterior (ARIMA(4,1,1)(0,1,1)[12] - definido pela intuição dos correlogramas), como reflexo do default da função admitir lags máximos e procedimentos que restringem a otimização da função verossimilhança. O AICc foi de 1046.34, pouco acima do anterior, indicando ajuste inferior que o primeiro.
Modelo ARIMA automático mais geral
================================
A recomendação do Hyndman e Athanasopoulos (2018), nestas situações, é fazer um modelo `auto.arima` mais geral, que investigue mais possibilidades de modelos,fazendo os argumentos da `auto.arima` da forma `stepwise=FALSE` e `approximation=FALSE`. Se estiver usando a série com ajuste sazonal (dessazonalizada), será interessante usar tambem o argumento `seasonal=FALSE`.
```{r}
require(fpp2)
fit3 <- auto.arima(varejoms,stepwise=FALSE, approximation=FALSE)
summary(fit3)
#p-values dos coeficientes
pvalues<-(1-pnorm(abs(fit3$coef)/sqrt(diag(fit3$var.coef))))*2
pvalues
``` 
Neste cenário, a rotina `auto.arima` encontra o ARIMA(4,1,0)(0,1,1)[12], com AICc de 1044.14, portanto, melhor (menor) que as duas primeiras estimações. como intuido anteriormente, ele detecta o AR(4), a Integração (I(1)), a integração sazonal e o SMA(1). O termo ar3 não foi significativo, mas os demais foram.
### Checando resíduos
Os modelos podem ser considerados white noise, mas não normais.
```{r}
checkresiduals(fit3)
require(FitAR) # pratico para plotar graficos de Q
LBQPlot(residuals(fit3),36)
# checar normalidade dos resíduos
require(tseries)
jarque.bera.test(residuals(fit3))
```
Pode-se fazer o gráfico dos forecasts deste modelo:
```{r}
autoplot(forecast(fit3,h=24),
 title ="Forecasts de volume de vendas no varejo total de Mato Grosso do Sul",
 xlab = "Ano", ylab = "Índice (2011=100)"
 )
```
### Teste de raiz unitária
O teste de raiz unitária deveria ser feito antes de rodar o ARIMA, a fim de investigar e/ou confirmar a não-estacionariedade da série, que foi intuida a partir da ACF e PACF. A rotina `auto.arima` já faz esse teste em seu processo. Existem diferentes testes de raiz unitária possíveis (ADF, PP, KPSS) assim como diferentes testes para a sazonalidade.
#### Teste ADF
O teste de Dickey-Fuller aumentado (ADF) é nossa primeira escolha. O padrão é que H0: série não estacionária ou série tem raiz unitária.
```{r}
ADF.test<-adf.test(varejoms)
ADF.test
```
#### Teste ADF invertendo a hipótese nula
Neste caso, o padrão é que H0: a série é estacionária ou a série não tem raiz unitária. Neste caso, a hipótese alternativa é que a série é "explosiva".
```{r}
ADF.test12<-adf.test(varejoms,k=12,alternative = c("explosive"))
ADF.test12
```
ARIMA - Codigo do Hyndman e Athanasopoulos, secao 8.1
======================
A função `ndiffs` usa um teste de raiz unitária para determinar o número de diferenças necessárias para que séries temporais xsejam feitas estacionárias. 
Se `test="kpss"`, o teste KPSS é usado com a hipótese nula que x possui uma raiz estacionária contra uma alternativa de raiz unitária. Em seguida, o teste retorna o menor número de diferenças necessárias para passar no teste ao nível alpha. 
Se `test="adf"`, o teste ADF é utilizado e se `test="pp"`, o teste de Phillips-Perron é usado. Em ambos os casos, a
hipótese nula é que x tem uma raiz unitária contra uma alternativa de série estacionária. Em seguida, o teste retorna o menor número de diferenças necessárias para falhar o teste no nível alpha.
A função `nsdiffs` usa testes de raiz unitária sazonal para determinar o número de diferenças sazonais necessárias para que as séries temporais x sejam feitas estacionárias (possivelmente com algum diferencial de atraso também). Se `test="ch"`, o teste de Canova-Hansen (1995) é usado (com hipótese nula de sazonalidade determinística) e se `test="ocsb"` o teste de Osborn-Chui-Smith-Birchenhall (1988) é usado (com hipótese nula de existir uma raiz unitária sazonal).
```{r}
x<-varejoms
ns <- nsdiffs(x) #numero de diferencas sazonais
ns
if(ns > 0) {
 xstar <- diff(x,lag=frequency(x),differences=ns)
} else {
 xstar <- x
} 
```
Portanto, temos uma diferença sazonal necessária (ns=1) e a série na diferença sazonal está armazenada em xstar.
```{r}
nd <- ndiffs(xstar) # numero de diferencas
nd
if(nd > 0) {
 xstar <- diff(xstar,differences=nd)
}
# xstar será a série devidamente em diferenças sazonais e diferenças normais
autoplot(xstar)
auto.arima(xstar) # ver que o resultado é igual, mas agora nao tem diferencas
modelo<-auto.arima(xstar,stepwise=FALSE, approximation=FALSE)
summary(modelo)
autoplot(modelo) 
```
Comentário: uma vez que a série utilizada em `auto.arima` é a resultante xstar, o mmodelo não identifica a necessidade de diferenciação nem diferenciação sazonal. Com esse procedimento, ele encontra o ARIMA(4,0,0)(0,01)[12] para xstar que já teve as diferenças aplicadas sobre varejoms. Da mesma forma que em fit e fit3, o modelo com argumentos `stepwise=FALSE, approximation=FALSE` retorna os melhores ajustes.
Referências {-#Referências}
========================
FERREIRA, Pedro Costa; SPERANZA, Talitha; COSTA, Jonatha (2018). BETS: Brazilian Economic Time Series. R package version 0.4.9. Disponível em: <https://CRAN.R-project.org/package=BETS>. 
HYNDMAN, Rob. (2018). fpp2: Data for "Forecasting: Principles and Practice" (2nd Edition). R package version 2.3. Disponível em: <https://CRAN.R-project.org/package=fpp2>. 
HYNDMAN, R.J., & ATHANASOPOULOS, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. Disponível em: <https://otexts.com/fpp2/>. Accessed on 12 Set 2019.