Introdução

Este material é dedicado a apresentação do modelo SARIMA. Para tal, fez-se uso da série temporal (ST) de vendas de passagens aéreas, mais conhecida como Air Passengers. Trata-se de uma série temporal mensal que registra o total de passageiros internacionais (em milhares) da linha aérea (Pan Am) no período de janeiro de 1949 a dezembro 1960, nos EUA.

A ST de vendas de passagens aéreas é um exemplo clássico de representação da modelagem de Box & Jenkins e a estrutura que “melhor” representa essa ST é um modelo SARIMA(0,1,1)(0,1,1). Dado a fama obtida por esse modelo, é equivalente dizer que uma ST segue um SARIMA(0,1,1)(0,1,1) ou um modelo Airline.

A abordagem de Box and Jenkins (1970) permite que valores futuros de uma série sejam previstos tomando por base apenas seus valores presentes e passados. Tais modelos são chamados de modelos autorregressivos integrados de médias móveis ou, simplesmente, ARIMA. Ao considerar relações sazonais, o modelo é nomeado SARIMA. Um modelo SARIMA\((p,d,q)(P,D,Q)_s\) é representado da seguinte forma:

\[\begin{equation} \phi(L)\Phi(L)\Delta^d\Delta^Dy_t = \theta(L)\Theta(L)\varepsilon_t \end{equation}\]

onde:
\(p\) é a ordem do polinômio autoregressivo não sazonal \(\phi(L)\);
\(P\) é a ordem do polinômio autoregressivo sazonal \(\Phi(L)\);
\(q\) é a ordem do polinômio de médias móveis não sazonal \(\theta(L)\);
\(Q\) é a ordem do polinômio de médias móveis sazonal \(\Theta(L)\);
\(d\) é a ordem de diferença não sazonal;
\(D\) é a ordem de diferença sazonal;
\(\phi(L) = (1 - \phi_1L - \phi_2L^2 -... - \phi_pL^p)\);
\(\Phi(L) = (1 - \Phi_1L^{s} - \Phi_2L^{2s} -... - \Phi_pL^{Ps})\);
\(\theta(L) = (1 - \theta_1L - \theta_2L^2 - ... - \theta_pL^q)\);
\(\Theta(L) = (1 - \Theta_1L^s - \Theta_2L^{2s} - ... - \Theta_pL^{Qs})\);
\(\Delta = 1 - L\);\ \(L\) é o operador de defasagem tal que \(L^ny_t = y_{t-n}\).

Ao longo dessa aula discutiremos as características dessa ST e os passos para modelá-la utilizando o software R, discutindo quais são os possíveis pacotes disponibilizados pelo programa. Conforme observaremos, essa é uma série temporal não estacionária nas partes sazonal e não sazonal e na variância. Aprenderemos a identificar essas características e qual é a maneira adequada de corrigí-las para fazermos uso da metodologia proposta por Box & Jenkins (Figura 1).

Figura 1: Ciclo Iterativo de Box & Jenkins

Ao ler esse material, pretende-se que o leitor esteja apto a modelar uma ST “não complexa”, seguindo a proposta de Box & Jenkins, utilizando o software R. Para atingir tal objetivo, nossa jornada pelo mundo dos modelos SARIMA terã as seguintes fases:

Preliminares: vamos definir o nosso diretório de trabalho e comentar sobre os pacotes que precisamos instalar para estimarmos e analisarmos de maneira correta o modelo SARIMA;

Análise Exploratória da ST de vendas de passagens aéreas: vamos explorar a ST de vendas de passagens aéreas observando sua tendência, variância e padrão sazonal;

Identificação: vamos aprofundar o nosso conhecimento sobre a ST que estamos trabalhando e discutir quais são os procedimentos que devemos adotar para modelá-la baseando-se no ciclo iterativo proposto por Box & Jenkins, iremos modelar a nossa ST e prevê-la 12 passos à frente

salvando a nossa previsão: mostraremos como exportar as previsões para arquivos .xlsx e .csv.

Preliminares

Instalando o pacote BETS

O pacote BETS (sigla para Brazilian Economic Time Series) para o R (R Core Team, 2012) 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 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). Em sua concepção original, o BETS pretendia 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 39 mil séries econômicas brasileiras.

Por conta do vasto tamanho da base 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 tornar-se 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.

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

## Loading required package: BETS

Demais pacotes de Séries Temporais

Além do BETS instalaremos os seguintes pacotes:

install.packages(“urca”): Unit root and cointegration tests for time series data (Pfaff, Zivot, and Stigler (2016));
install.packages(“TSA”): Time Series Analysis (Chan and Ripley (2012));
install.packages(“forecast”): Forecasting Functions for Time Series and Linear Models (Hyndman, Razbash, and Schmidt (2012));
install.packages(“lmtest”): Testing Linear Regression Models (Zeileis and Hothorn (2002));
install.packages(“normtest”): Tests for Normality (gavrilov2014package);
install.packages(“FinTS”): Companion to Tsay (2005) Analysis of Financial Time Series (Graves (2014));
install.packages(“xlsx”): Read, write, format Excel 2007 and Excel 97/2000/XP/2003 files (Dragulescu (2014)).

Após a instalação, é preciso usar a função library() para carregar os pacotes, mas faremos isso ao longo do texto para que fique claro para o leitor em quais pontos estamos usando os pacotes.

Análise Exploratória da ST

Por ser uma ST conhecida, o R já a disponibiliza na sua base de dados, tornando-se muito fácil a sua leitura. Basta executar os seguintes comandos:

ts.plot(AirPassengers)

Vamos vizualiza-la utilizando o dygraph?!

install.packages("dygraphs")
library(dygraphs)
dygraph(AirPassengers)

Observando o gráfico da ST de vendas de passagens aéreas podemos perceber que há uma tendência crescente do número de passageiros. As oscilações de picos e vales podem ser relacionadas às estações do ano, nas quais, mais especificamente, temos períodos de férias, feriados, etc. Essas oscilações, como observadas, acontecem anualmente, o que nos faz acreditar que há presença de sazonalidade. Do começo do ano a outubro percebemos um comportamento crescente, seguido de um comportamento decrescente da série. Isso se repete todo ano.

Nesse sentido, apenas observando o gráfico podemos “levantar” as seguintes hipóteses sobre essa ST:

Tendência: parece haver aumento do número de passageiros transportados ao longo dos anos pela Pam Am. Coerente com a teoria econômica pois espera-se que ao longo do tempo a empresa cresça e, consequentemente aumente as vendas de passagens aéreas.

Variância: observa-se que, além do aumento do número de passagens vendidas, a distância entre os meses com maiores e menores vendas também está aumentando, indicando aumento da variância. Fato este também coerente com a teoria econômica pois ao aumentar o volume de vendas, espera-se maiores oscilações em relação ao valor médio.

Sazonalidade: verifica-se um comportamento sazonal das vendas de passagens aéreas. Isto é, nos meses de março (feriado de Páscoa) e julho (Dia da Independência e férias escolares) há um aumento nas vendas quando comparado com os meses anteriores. Além disso, parece que a sazonalidade é crescente ao longo do tempo.

Observe que a análise gráfica nos permitiu conhecer a nossa ST e é uma fase muito importante para esse tipo de modelagem. Obviamente, como bons econometristas que somos, iremos testar estatisticamente todas os pontos levantados anteriormente. Antes disso, vamos tentar entender um pouco mais o comportamento sazonal da nossa ST.

Uma análise um pouco mais profunda da sazonalidade

Os gráficos monthplot e ggfortify ajudam a detectar visualmente a presença de sazonalidade na ST. Como se pode verificar, esta ST apresenta a média e a variância não constantes, indícios de não estacionariedade na parte sazonal da ST.

Observando o gráfico, podemos ver que o número médio de passageiros (traços horizontais) aumenta nos meses de férias (indício de sazonalidade). Analisando os traços verticais, verifica-se um aumento contínuo na venda de passagens aéreas ano a ano, indício de não estacionariedade na parte sazonal da Série Temporal.

monthplot(AirPassengers)

library(fpp2)
ggseasonplot(AirPassengers, polar = T)

Decomposição da ST

Por fim, seguindo a decomposição clássica, decompõe-se a ST utilizando filtros de médias móveis em três componentes principais:

tendência + ciclo
sazonalidade
resíduo (componente irregular, inovação)

plot(decompose(AirPassengers))

Conforme observado no gráfico da ST, verifica-se que o número mínimo de passagens vendidas foi de 104 (Nov-1949) e o máximo de 622 (Jul-1960). Ao observar o componente de tendência, observa-se que a ST é fortemente afetada por esse componente (em torno de 85%). Com relação ao componente sazonal, verifica-se que o mesmo também está presente nessa ST e gira em torno de 10%. Sobrando uma pequena parte de componente irregular, o qual é levemente “contaminado” pela parte sazonal, mostrando que método de decomposição utilizado não é muito eficiente.

Essa análise é interessante pois mostra que, basicamente, precisamos modelar as componentes de tendência e sazonalidade (em torno de 95% da ST), componentes “bem” modelados pelos modelos SARIMA(p,d,q)(P,D,Q). Este fato mostra porque essa ST é tão utilizada para exemplificar o uso dessa metodologia.

Prosseguindo, nosso próximo passo será testar estatisticamente as percepções levantadas anteriomente. Isto é, a ST de vendas de passagens aéreas é realmente não estacionária na parte não sazonal? A ST de vendas de passagens aéreas é realmente não estacionária na parte sazonal? Como faremos para “corrigir” esses “problemas”?

Conhecendo a ST antes de iniciar a modelagem BJ

Para responder aos questionamentos feitos na seção anterior, iremos abordar os dois tópicos:

Testar a estacionariedade da parte não sazonal
Testar a estacionariedade da parte sazonal

Testando a estacionariedade da parte não sazonal

Há, basicamente, quatro maneiras de observar se a ST em estudo é ou não estacionária:

Análise gráfica;
Comparar a média e a variância em diferentes períodos de tempo da ST;
Observar a FAC (Função de Autocorrelação);
Testes de Raiz Unitária.

Já vimos que a análise gráfica nos mostrou indícios de não estacionariedade. Fica claro também que, se “fatiássemos” a ST e calculássemos as médias de cada ano, observaríamos uma tendência de alta nas médias, indicando não estacionariedade das mesmas.

Outra maneira de ver a não estacionariedade da ST é visualizando o gráfico da FAC. A figura a seguir mostra que as autocorrelações plotadas pela FAC não decrescem exponencialmente ou de forma senoidal conforme descrito pela teoria de Box & Jenkins}. Esse é mais um indicativo de que a ST é não estacionária

BETS.corrgram(AirPassengers,lag.max = 36)

Nesse momento, o leitor atento pode estar se fazendo a seguinte pergunta: Para que tantas maneiras de se observar a estacionáriedade, se ao observar o gráfico da ST já está claro que a mesma é não estacionária? A resposta a esse questionamento é que nenhuma das maneiras, vistas até o momento, de se verificar se a ST é ou não estacionária nos dá uma resposta “clara” (com significância estatística) se a ST é ou não estacionária. Mais ainda, tais métodos, não nos dizem quantas diferenciações precisaremos fazer na ST em estudo para torná-la estacionária e qual é o tipo de não estacionariedade (determinística ou estocástica). Para obter essas respostas precisamos testar a estacionariedade através dos testes de Raiz Unitária.

Os testes de Raíz Unitária (RU) foram uma grande revolução na Econometria na década de 1980. Existe uma grande quantidade de testes e, basicamente, todos têm a mesma ideia, isto é, a hipótese nula é que a série temporal possui uma raiz unitária (a ST é não-estacionária) e a hipótese alternativa a que a série é estacionária, com exceção do teste KPSS que tem as hipóteses alternadas. Abaixo podemos ver alguns exemplos de testes de Raiz Unitária:

Augmented Dickey Fuller (ADF)
Phillips-Perron (PP)
Kwiatkowski-Phillips-Schmidt-Shin (KPSS)
Dickey Fuller GLS (DF-GLS)
Elliott–Rothenberg-Stock point optimal (ERS)

Apesar de haver uma grande quantidade de testes, nesse livro abordaremos apenas o teste de Dickey Fuller Aumentado (ADF), que tem as seguintes hipóteses:

H0: a ST possui uma RU. A série é não estacionária.
Ha: a ST não possui RU. A série é estacionária

A regra de rejeição da hipótese nula funciona da seguinte forma: se o valor observado para a estatística de teste for inferior ao valor crítico, rejeitamos a hipótese nula e, portanto, concluímos que a ST é estacionária de acordo com o nível de confiança adotado previamente. Caso contrário, a ST será não estacionária.

Como existem várias especificações consistentes com a não-estacionariedade, irão existir várias formas de testá-la. Na prática, a questão importante é escolher a forma do teste de RU adequada para a ST em questão. O teste ADF apresenta as seguintes formas:

Raiz unitária + constante + tendência determinística (R: type = "trend")
Raiz unitária + constante (R: type = "drift")
Raiz unitária (R: type = "none")

Para executar o teste no R, usaremos a função ur.df() do pacote urca. Os principais argumentos dessa função são:

ur.df(y, type = c("none", "drift", "trend"), lags = 1,
      selectlags = c("Fixed", "AIC", "BIC"))

    y: ST em que será testada a raiz unitária;
    type: tipo da especificação do teste que o usuário deseja realizar;
    lags: número de defasagens a serem usadas para captar o comportamento da ST e, consequentemente,   corrigir o problema da autocorrelação residual;
    selectlags: a função pode definir automaticamente, baseada em um critério de informação, o número de lags a serem inclusos dado um valor máximo no argumento lags.

Antes de iniciar o teste é importante observar que o número de lags que serão incluídos na equação do teste ADF será definido com base na análise dos resíduos da regressão e não somente nos critérios de informação.

Dando início aos testes, vamos testar a estacionariedade da ST considerando a equação sem tendência e com constante. É importante você saber que testamos a equação com tendência determinística antes, porém o parâmetro dessa variável não foi significativo. Nessa fase o parâmetro mais difícil e importante de definir é o lag, isto é, você precisa encontrar um número de lags que corrija a autocorrelação dos resíduos e seja parcimonioso com relação ao número de parâmetros da equação do modelo.

Estipulamos, inicialmente, o lag máximo como 24 e o critério de informação a minimizar sendo o AIC. A seguir observamos o resultado do teste ADF e a FAC do resíduos, a qual mostra que não há presença de autocorrelação.

require(urca)
adf.drift <- ur.df(y = AirPassengers, type = c("drift"),lags = 24, selectlags = "AIC")
summary(adf.drift)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.763  -5.824  -1.537   5.107  31.972 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.91071    3.45201   2.002 0.048143 *  
## z.lag.1       0.02250    0.01211   1.858 0.066240 .  
## z.diff.lag1  -0.52049    0.10323  -5.042 2.20e-06 ***
## z.diff.lag2  -0.30032    0.10929  -2.748 0.007177 ** 
## z.diff.lag3  -0.28030    0.11026  -2.542 0.012634 *  
## z.diff.lag4  -0.44588    0.11167  -3.993 0.000129 ***
## z.diff.lag5  -0.12369    0.11733  -1.054 0.294456    
## z.diff.lag6  -0.15506    0.11581  -1.339 0.183764    
## z.diff.lag7  -0.24401    0.11651  -2.094 0.038888 *  
## z.diff.lag8  -0.28984    0.11884  -2.439 0.016591 *  
## z.diff.lag9   0.09797    0.12135   0.807 0.421488    
## z.diff.lag10 -0.34144    0.12258  -2.785 0.006453 ** 
## z.diff.lag11 -0.28896    0.11426  -2.529 0.013090 *  
## z.diff.lag12  0.52135    0.11359   4.590 1.36e-05 ***
## z.diff.lag13  0.26911    0.12106   2.223 0.028590 *  
## z.diff.lag14 -0.23993    0.12398  -1.935 0.055944 .  
## z.diff.lag15 -0.04974    0.12420  -0.400 0.689730    
## z.diff.lag16 -0.04429    0.12377  -0.358 0.721283    
## z.diff.lag17 -0.19230    0.12548  -1.533 0.128715    
## z.diff.lag18 -0.41198    0.12676  -3.250 0.001596 ** 
## z.diff.lag19 -0.10697    0.12564  -0.851 0.396697    
## z.diff.lag20 -0.30006    0.12190  -2.462 0.015639 *  
## z.diff.lag21 -0.46157    0.12496  -3.694 0.000369 ***
## z.diff.lag22 -0.23550    0.12723  -1.851 0.067281 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.55 on 95 degrees of freedom
## Multiple R-squared:  0.9326, Adjusted R-squared:  0.9163 
## F-statistic: 57.18 on 23 and 95 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: 1.8582 7.9144 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

require(TSA)
BETS.corrgram(adf.drift@res,lag.max = 36)

Ao analisar a estatística de teste (tau2 = 1,8582) notamos que seu valor é superior ao valor crítico associado ao nível de confiança de 95% (-2,88). Dessa forma, conclui-se que a ST não é estacionária (não rejeição da hipótese nula).

O leitor pode visualizar mais informações sobre o teste de RU, como a equação ajustada por exemplo, usando a função summary(adf.drift).

Ao concluir que a ST tem raiz unitária, precisamos descobrir o número de diferenciações necessárias para torná-la estacionária. É importante observar que esse é apenas um exercício para que o leitor observe o comportamento da ST e da FAC antes e após a diferenciação, pois, como veremos nas próximas seções, faremos as “correções” de não estacionariedade da ST na própria função que estimará o modelo SARIMA.

Dado que a nossa ST é não estacionária, vamos tentar torná-la estacionária fazendo uma diferenciação e vamos observar o gráfico e a FAC novamente.

ts.plot(diff(AirPassengers, lag = 1, differences = 1))
require(BETS)
BETS.corrgram(diff(AirPassengers, lag = 1, differences = 1),lag.max = 36)

Observe que ao aplicar a diferenciação, a ST aparenta estar estacionária na média, mas a variância é crescente ao longo do tempo. Como sabemos, um dos pressupostos da teoria Box & Jenkins é que a ST seja também estacionária na variância. Para tal, iremos passar o log na ST em questão.

ts.plot(diff(log(AirPassengers),lag = 1,differences = 1))
require(BETS)
BETS.corrgram(diff(log(AirPassengers), lag = 1, differences = 1),lag.max=48)

Note agora que temos uma série temporal estacionária tanto na média quanto na variância. Ao analisarmos a FAC, a pergunta que fica é: essa FAC é adequada para idenficarmos a estrutura do nosso modelo SARIMA?

Avaliando a estacionariedade da parte sazonal

Com relação a pergunta feita na seção anterior, o leitor atento já observou que nos lags sazonais a função de autocorrelação também apresenta um decrescimento lento, indicando que a ST é não estacionária na parte sazonal*.

(*) Observe que na FAC `tradicional` gerada pelo R os lags sazonais são 1(=12), 2(=24), 3(=36), etc
Existem testes estatísticos para avaliar a presença de não estacionariedade sazonal, um dos mais conhecidos é o teste de HEGY.

Para corrigir esse problema precisamos diferenciar a parte sazonal, para isso diferenciaremos a ST já diferenciada na parte não sazonal. Tal procedimento é feito mudando o parâmetro lag da função diff() para 12, conforme pode ser observado abaixo:

require(BETS)
BETS.corrgram(diff(diff(log(AirPassengers), lag = 1, differences = 1),lag = 12, differences = 1), lag.max = 48)

Observe que agora a FAC apresenta cortes bruscos nos lags 1 e 12 e não apresenta mais decrescimento lento nem na parte sazonal nem na não sazonal.

Vamos refazer o teste de RU para confirmar a estacionariedade da ST após aplicar as transformações anteriores. O valor da estatística de teste (tau2 = -4,0398) é inferior ao valor crítico (-2,88) ao nível de significância de 95%. Assim, podemos concluir que a série é estacionária.

adf.drift2 <- ur.df(y = diff(diff(log(AirPassengers), lag = 1), lag = 12), 
                    type = "drift", lags = 24, selectlags = "AIC")
require(BETS)
BETS.corrgram(adf.drift2@res, lag.max = 36)

Ao analisar a FAC para os resíduos do teste ADF, o leitor pode notar que alguns lags aparecem significativos, porém não são relevantes (apresentam correlação muito baixa). Dessa forma, confirmamos a validade do teste e podemos começar a nossa modelagem.

Modelando a ST

A metodologia Box & Jenkins para séries temporais estacionárias e construção dos modelos ARIMA segue um ciclo iterativo composto por cinco partes:

Especificação: a classe geral de estruturas SARIMA(p,d,q)(P,D,Q) é analisada. Identificação: com base na FAC e FACP amostrais e outros critérios são definidos os valores de p,q e P,Q.

Estimação: os parâmetros do modelo identificado são estimados e testados estatisticamente sobre sua significância.

Diagnóstico: faz-se uma análise dos resíduos (devem ser ruído branco) e testes de verificação (Ljung-Box) para ver se o modelo sugerido é adequado. Em seguida, verifica-s os modelos que apresentam menores valores para os critérios AIC e BIC. Caso haja problemas no diagnóstico, volta-se à identificação.

Modelo definitivo: para previsão ou controle. Verificar quais modelos têm as melhores medidas RMSE (Root Mean Square Error) e MAPE (Mean Absolute Percentage Error) (este não vale para dados próximos de zero, sendo preferível a utilização de outro método para a análise dos erros), por exemplo.

Um processo ARIMA(p,d,q) é um ARMA diferenciado d vezes até estar estacionário. Os modelos SARIMA são usados para séries temporais que apresentam comportamento periódico em s espaços de tempo, isto é, quando ocorrem desempenhos semelhantes a cada intervalo de tempo. Este é o caso da série a ser trabalhada neste capítulo.

Identificação

Como sabemos, o primeiro passo para identificar o nosso modelo SARIMA é observando a FAC e a FACP. Como os modelos propostos por Box & jenkins são da década de 1970, o esforço computacional para estimar o modelo era muito grande, portanto essa fase era fundamental para se ter um modelo adequado à ST em análise. Atualmente, graças aos avanços computacionais, observar a FAC e a FACP é útil, principalmente, para se ter uma ideia inicial do modelo a ser testado, pois, como veremos mais adiante, o ideal é escolher um modelo que minimize os critérios de informação.

Assim, vamos observar a FAC e a FACP (função de autocorrelação parcial) da ST de vendas de passangens aéreas diferenciada na parte sazonal e não sazonal e com transformação logarítmica.

BETS.corrgram( diff(diff(log(AirPassengers), lag = 1, differences = 1),lag = 12, differences = 1), lag.max = 48)


BETS.corrgram(diff(diff(log(AirPassengers), lag = 1, differences = 1),lag = 12, differences = 1), lag.max = 48,type= "partial")

Observando os gráficos com um pouco de boa vontade podemos pensar nos seguintes modelos:

- SARIMA(1,1,1)(1,1,1) - corte brusco na FAC e na FACP nas partes sazonais e não sazonais;
- SARIMA(0,1,1)(0,1,1) - corte brusco na FAC e decrescimento das partes sazonais e não sazonais.

Uma vez identificado os possíveis modelos, passa-se para o próximo passo: a estimação.

Estimação

Para estimar o modelo, deve-se testar as possibilidades dos SARIMAs que idealizamos a partir da visualização da FAC e da FACP. Para tanto, utilizaremos a função Arima() do pacote forecast. Com relação ao método de estimação dos parâmetros neste trabalho, usaremos o default do R, que utiliza o método de Máxima Verossimilhança, denotado como ML (Maximum Likelihood).

Dessa forma, o primeiro modelo a ser estimado será uma SARIMA(1,1,1)(1,1,1)\(_{12}\). Observe que na função Arima() a variável de entrada é a ST original, mas ajustar o argumento lambda em zero permite que seja feita a tranformação logarítmica na série. Também não é necessário diferenciar a ST antecipadamente pois o própria função faz isso.

library("forecast")
fit.air <- Arima(AirPassengers, order = c(1,1,1), seasonal = c(1,1,1),
                  method = "ML", lambda = 0)
fit.air

## Series: AirPassengers 
## ARIMA(1,1,1)(1,1,1)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1      ma1     sar1    sma1
##       0.1668  -0.5616  -0.0994  -0.497
## s.e.  0.2458   0.2114   0.1540   0.136
## 
## sigma^2 estimated as 0.00138:  log likelihood=245.16
## AIC=-480.31   AICc=-479.83   BIC=-465.93

Para verificar, de forma rápida, se os parâmetros do modelo são significativos, desenvolvemos uma função no pacote BETS chamada t.test(), o código da função está disponibilizado a seguir. Consideramos nessa função o nível de confiança de 95%.

BETS.t_test(fit.air)

##           Coeffs Std.Errors         t Crit.Values Rej.H0
## ar1   0.16679124  0.2457980 0.6785705    1.977304  FALSE
## ma1  -0.56163441  0.2114211 2.6564723    1.977304   TRUE
## sar1 -0.09938487  0.1539918 0.6453907    1.977304  FALSE
## sma1 -0.49700743  0.1360485 3.6531644    1.977304   TRUE

Como podemos observar, os parâmetros da parte AR não sazonal e sazonal são não significativos, logo, tais parâmetros não devem permanecer no modelo. Então, estes foram retirados e o modelo foi reestimado.

fit.air <- Arima(AirPassengers, order = c(0,1,1), seasonal = c(0,1,1),
                 method = "ML", lambda=0)
fit.air

## Series: AirPassengers 
## ARIMA(0,1,1)(0,1,1)[12] 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 estimated as 0.001371:  log likelihood=244.7
## AIC=-483.4   AICc=-483.21   BIC=-474.77

BETS.t_test(fit.air)

##          Coeffs Std.Errors        t Crit.Values Rej.H0
## ma1  -0.4018268 0.08964405 4.482470    1.977054   TRUE
## sma1 -0.5569466 0.07309948 7.619023    1.977054   TRUE

Conforme pode ser observado, temos um modelo SARIMA(0,1,1)(0,1,1)\(_{12}\) onde todos os parâmetros são significativos e que minimiza os critérios de informação (BIC,AIC e AICc).

Diagnóstico

Após definir a “melhor” estrutura e estimar os parâmetros do modelo, outra etapa fundamental é a fase de diagnóstico do modelo. Nesta fase as seguintes características dos resíduos precisam ser analisadas e confirmadas:

Ausência de autocorrelação linear;
Ausência de heterocedasticidade condicional;
Normalidade.

Para uma visão geral dos resíduos, utiliza-se a função tsdiag(). Esta disponibiliza a distribuição dos resíduos padronizados, a função de autocorrelação dos resíduos e os p-valores da estatística do teste Ljung-Box. Conforme podemos observar no primeiro gráfico a seguir, os dados aparentam estar distribuídos simetricamente em torno da média zero, indicação de distribuição normal. Observe também que não temos nenhuma informação discrepante (muito fora do intervalo [-3,3]). A única exceção é o resíduo de janeiro de 1954, neste caso, poderíamos testar se a venda de passagens aéreas nesse mês é um outlier ou não.

O segundo gráfico disponibilizado pela função tsdiag() é a FAC dos resíduos. Este gráfico é extremamente útil para observar se há a presença de autocorrelação linear nos resíduos. Conforme verificamos, não há nenhum lag significante, logo, toda a parte linear da ST de vendas de passagens aéreas foi capturada pelo modelo SARIMA(0,1,1)(0,1,1)\(_{12}\).

O terceiro gráfico mostra o p-valor da estatística Ljung-Box para diferentes defasagens após a defasagem 14. De acordo com o gráfico, verificamos que não rejeitamos a hipótese nula da não existência de dependência serial para todas as defasagens. Tal resultado está em consonância com a análise feita anteriormente, isto é, não há dependência linear nos resíduos. Contudo, este gráfico não é confiável uma vez que os graus de liberdade usados para calcular os p-valores são baseados nos lags e não em (lags - (p+q)). Isto é, o processo usado para calcular os p-valores não leva em conta o fato de os resíduos terem sido gerados a partir de um modelo ajustado. Portanto, precisamos tomar cuidado ao observar esse gráfico.

diag <- tsdiag(fit.air, gof.lag = 20)

Bem, conforme observamos a função tsdiag() já nos deu bastante informação sobre os nossos resíduos, entretanto, vamos realizar alguns testes estatísticos específicos para cada uma das características em decorrência do problema da estatística Ljung-Box e da necessidade de testarmos a normalidade e a homocedasticidade dos resíduos.

Primeiramente, vamos testar a autocorrelação linear dos resíduos através do teste de Ljung Box. Como sabemos o teste de Ljung Box nos dá a presença ou não de autocorrelação serial dos resíduos para o “L” primeiros lags. Outro teste de autocorrelação residual muito conhecido é o teste de Durbin & Watson.

Box.test(x = fit.air$residuals, lag = 24,type = "Ljung-Box", fitdf = 2)

Conforme podemos observar, o resultado do teste de Ljung Box mostra que a 95% de confiança não rejeitamos a hipótese nula de não existência de autocorrelação serial até o lag 24. É importante observar o argumento fitdf, neste caso igual a 2 (p+q), pois o teste é feito nos resíduos de um modelo SARIMA com dois parâmetros.

Confirmada a ausência de autocorrelação linear nos resíduos, vamos testar a estacionariedade da variância. Para tal, faremos o teste Multiplicador de Lagrange para heterocedasticidade condicional autorregressiva (ARCH LM) disponível no pacote FinTS.

require(FinTS)
ArchTest(fit.air$residuals,lags = 12)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  fit.air$residuals
## Chi-squared = 14.859, df = 12, p-value = 0.2493

Conforme mostrado pelo teste, a hipótese nula é que não há presença de efeito ARCH. Dessa forma, dado o valor do p-valor, não rejeitamos a hipótese nula a 95% de confiança, logo, a variância é estacionária.

Por fim, precisamos testar a normalidade do nosso resíduo. Para tal, faremos o teste de Jarque Bera baseando-se no pacote normtest.

require(normtest)
jb.norm.test(fit.air$residuals, nrepl=2000)

## 
##  Jarque-Bera test for normality
## 
## data:  fit.air$residuals
## JB = 5.2265, p-value = 0.049

Os resultados mostram que a 95% de confiança não rejeitamos a hipótese nula de normalidade. Feito o diagnóstico dos resíduos, o nosso próximo passo será fazer as previsões.

Uma maneira mais direta de fazer o diagnóstico dos resíduos é utilizando a função checkresiduals do pacote forecast

library(forecast)
checkresiduals(fit.air)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 26.446, df = 22, p-value = 0.233
## 
## Model df: 2.   Total lags used: 24

Previsão

Após fazermos o diagnóstico dos resíduos e concluirmos que estamos modelando toda a parte linear da ST de vendas de passagens aéreas, o nosso próximo passo é fazer a previsão. Afinal de contas, esse é nosso objetivo final. Nessa etapa, basicamente, queremos conhecer a nossa previsão, saber qual é o intervalo de confiança (neste caso, 95%) e analisar as métricas de desempenho do modelo.

Para a previsão utilizaremos o pacote forecast e a função com o mesmo nome. Observe que, ao usar esta função, precisamos definir os seguintes parâmetros: (a) object: output do modelo SARIMA estimado; (b) h: horizonte de previsão (quantos passos à frente queremos prever); e (c) level: é o nível de confiança que desejamos para o nosso intervalo de confiança.

require(forecast)
plot(forecast(object = fit.air, h=12, level = 0.95))

Observando o gráfico, parece que fizemos uma “boa” previsão. Porém, uma maneira mais adequada de certificar isso é analisando as métricas de previsão. Conforme podemos observar, as métricas a seguir confirmam a análise gráfica. Analisando o MAPE, por exemplo, que é uma medida percentual do módulo dos erros e que não é contaminada pela escala da ST, observamos que o erro de previsão está apenas em 2,62%, o que é muito bom!*

(*) Dois pontos que gostaria de destacar aqui: primeiro, que uma análise da previsão fora da amostra seria importante para corroborar a performance do nosso modelo. Segundo, essa ideia de bom ou ruim é muito relativa, isto é, é sempre bom termos um modelo benchmark para compararmos nossas previsões.}.

accuracy(fit.air)

##                      ME     RMSE      MAE          MPE     MAPE      MASE
## Training set 0.05140376 10.15504 7.357555 -0.004079321 2.623637 0.2297061
##                     ACF1
## Training set -0.03689736

Exportando as Previsões

Imagine o seguinte: você trabalha em uma empresa na área financeira e seu chefe lhe pede a previsão das vendas de um determinado produto para os próximos 12 meses. Ainda, imagine também que ninguém na sua empresa conheça o R (não é tão dificíl de imaginar isso, certo?!).

A solução para o primeiro problema você já tem e já aprendeu ao longo desse capítulo. A solução para o segundo problema pode ser treinar toda a equipe da área financeira para trabalhar com o R ou então extrair as previsões e os intervalos de confiança estimados para um programa mais conhecido como o Excel através de um arquivo .csv ou .xlsx.

Como veremos, essa tarefa é muito fácil de fazer no R e pode ser executada com apenas uma linha de comando.

Em formato .csv:

write.csv(data.frame(forecast(object = fit.air, h=12, level = 0.95)),"previsao.csv")

Em formato .xlsx:

#require(xlsx)
#write.xlsx(data.frame(forecast(object = fit.air, h=12, level = 0.95)),"previsao.xlsx")

Considerações finais

Neste capítulo aprendemos empiricamente como modelar uma série temporal mensal com base na metodologia proposta por Box & Jenkins utilizando o software R. Aprendemos como fazer uma análise exploratória de uma ST, quais são os possíveis “problemas” que ela pode ter para ser modelada utilizando o arcabouço proposto por Box & Jenkins e como “consertar” esses problemas, através, por exemplo, da diferenciação da ST.

Foram abordados também alguns “pacotes” úteis para esse tipo de modelagem, discutimos algumas funções e chamamos a atenção para algumas limitações das mesmas. Apesar de ter sido uma experiência interessante, sabemos que ainda ficaram faltando alguns pontos a serem abordados, como por exemplo, não tratamos da identificação e “correção” de possíveis outliers, não mostramos como “corrigir” a presença de heterocedasticidade condicional nos resíduos, quando ela existir, etc.

Nesse sentido, é importante que o leitor que estiver usando esse manual para construir o seu modelo SARIMA tenha ciência de suas limitações e busque, sempre que possível, aprofundar o seu conhecimento sobre o assunto.

Exercício Complementar

Faça a previsão da venda de passagens aéreas utilizando os modelos ARIMA, NAIVE e SNAIVE. Observe e comente os resultados.

Modelo Naive

str(AirPassengers)

##  Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...

training <- window(AirPassengers, end = c(1959,12))
test <- window(AirPassengers, start = c(1960,1))
require(fpp2)
fc <- naive(training, h = 10,level = .9)
checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 294.49, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

autoplot(fc) + autolayer(test, series = "conjunto de teste") + autolayer(fitted(fc), series = "valores ajustados")

## Warning: Removed 1 rows containing missing values (geom_path).

require(forecast)
accuracy(fc,test)

##                     ME      RMSE      MAE        MPE      MAPE     MASE
## Training set  2.236641  31.33213 24.08397  0.4168179  8.979488 0.790935
## Test set     84.200000 112.38149 87.00000 15.3758783 16.091991 2.857143
##                   ACF1 Theil's U
## Training set 0.2863378        NA
## Test set     0.7025242  2.074598

Modelo SNAIVE

str(AirPassengers)

##  Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...

training <- window(AirPassengers, end = c(1959,12))
test <- window(AirPassengers, start = c(1960,1))
require(fpp2)
fcs <- snaive(training, h = 10)
checkresiduals(fcs)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 228.34, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

autoplot(fcs) + autolayer(test, series = "conjunto de teste") + autolayer(fitted(fcs), series = "valores ajustados")

## Warning: Removed 12 rows containing missing values (geom_path).

require(forecast)
accuracy(fcs,test)

##                    ME     RMSE   MAE      MPE     MAPE     MASE       ACF1
## Training set 30.16667 34.54828 30.45 11.23808 11.37483 1.000000  0.7658250
## Test set     51.90000 54.16918 51.90 10.64209 10.64209 1.704433 -0.1050355
##              Theil's U
## Training set        NA
## Test set      1.058275

Modelo SARIMA

training <- window(AirPassengers, end = c(1959,12))
test <- window(AirPassengers, start = c(1960,1))
fit.Air <- Arima(training, order = c(0,1,1), seasonal = c(0,1,1),
                 method = "ML",lambda = 0)
require(forecast)
checkresiduals(fit.Air)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 21.012, df = 22, p-value = 0.52
## 
## Model df: 2.   Total lags used: 24

forc.Air <- forecast(object = fit.Air, h=12, level = 0.95)


autoplot(forc.Air) + autolayer(test, series = "conjunto de teste") + autolayer(fitted(forc.Air), series = "valores ajustados")

require(forecast)
accuracy(forc.Air,test)

##                       ME      RMSE       MAE         MPE     MAPE
## Training set   0.2583698  9.000709  6.753757  0.04572919 2.587482
## Test set     -12.1637992 18.594929 13.262663 -2.66648999 2.904855
##                   MASE        ACF1 Theil's U
## Training set 0.2217982  0.02120815        NA
## Test set     0.4355554 -0.30142251 0.4094696

Materiais de Apoio

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;
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;
Khan Academy: ideal para aprender conceitos básicos de matemática e estatística;
BETS (Brazilian Economic Time Series package) – este é o pacote R que estamos desenvolvendo, fiquem à vontade para instala-lo. Iremos utiliza-lo no nosso curso;

# Instalação do pacote pelo CRAN 
install.package("BETS")
                      
# Para obter sempre a versão em desenvolvimento e acompanhar o precessod de criação:
install.packages("devtools")
require(devtools)
install_github("pedrocostaferreira/BETS")
require(BETS)

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

Referências

Box, G. E. P., and G. M. Jenkins. 1970. “Time Series Analysis: Forecasting and Control.” San Francisco: Holden Day.

Chan, Kung-Sik, and Brian Ripley. 2012. TSA: Time Series Analysis. https://CRAN.R-project.org/package=TSA.

Dragulescu, Adrian A. 2014. Xlsx: Read, Write, Format Excel 2007 and Excel 97/2000/Xp/2003 Files. https://CRAN.R-project.org/package=xlsx.

Graves, Spencer. 2014. FinTS: Companion to Tsay (2005) Analysis of Financial Time Series. http://cran.r-project.org/package=FinTS.

Hyndman, Rob J, Slava Razbash, and Drew Schmidt. 2012. “Forecasting Functions for Time Series and Linear Models.” R Package Version (Http://Cran. R-Project. Org/Web/Packages/Forec Ast/).

Pfaff, Bernhard, Eric Zivot, and Matthieu Stigler. 2016. Urca: Unit Root and Cointegration Tests for Time Series Data. https://cran.r-project.org/web/packages/urca/index.html.

Zeileis, Achim, and Torsten Hothorn. 2002. “Diagnostic Checking in Regression Relationships.” R News 2 (3): 7–10.

Business forecasting: Previsão das vendas de passagens aéreas utilizando os modelos (S)ARIMA

Pedro Costa Ferreira