| data | total |
|---|---|
| 2003-09-01 | 98.472 |
| 2002-07-01 | 96.199 |
| 2016-06-01 | 139.567 |
| 2023-08-01 | 153.478 |
| 2003-06-01 | 76.772 |
| 2007-02-01 | 118.806 |
| 1994-10-01 | 102.604 |
| 2005-07-01 | 110.037 |
| 1994-02-01 | 73.365 |
| 2020-09-01 | 161.071 |
| 2017-05-01 | 163.284 |
| 1991-03-01 | 50.278 |
| 2017-03-01 | 158.055 |
| 2001-08-01 | 107.488 |
| 2016-11-01 | 148.675 |
| 1994-12-01 | 115.876 |
| 1999-02-01 | 44.579 |
| 2010-06-01 | 192.457 |
| 1998-06-01 | 99.928 |
| 2009-05-01 | 195.720 |
EST562 - Análise de Séries Temporais
Especialização em Estatística Computacional Aplicada - Prof. Thiago Rezende
Introdução
O relatório apresentado a seguir foi produzido para fins acadêmicos. Trata-se do trabalho prático final da disciplina Análise de Séries Temporais. O objetivo é aplicar os conceitos vistos em sala de aula. Toda a análise foi realizada usando a linguagem R (v. 4.3.1)1.
Dados
Para o trabalho, foi escolhida a série temporal mensal de vendas de veículos pelas concessionárias - automóveis disponibilizada pelo portal de dados abertos do Banco Central do Brasil e cuja a fonte é a Federação Nacional da Distribuição de Veículos Automotores. Atualmente, o período da série é de 01/01/1990 a 01/02/2024.
Conceito: unidades comercializadas de automóveis, geradas pelos emplacamentos. Todo emplacamento de veículo origina emissão de uma nota fiscal de venda. Os emplacamentos englobam as vendas diretas e vendas a varejo. O que diferencia uma venda direta de uma venda a varejo é a identificação do emissor dessa Nota Fiscal. Se o CNPJ do emissor for de um fabricante de veículos, esta será uma venda direta. Já as demais, serão consideradas vendas a varejo2.
Sobre a série temporal, é interessante observar as mudanças abruptas de nível em dezembro de 1997 e em fevereiro de 2015. Essa constatação influenciou o ajuste de um dos modelos a seguir. Entre 2020 e 2021 também houve uma nova mudança, provavelmente devido à pandemia de COVID-19. No entanto, nesse caso, a mudança de nível foi breve.
| Mínimo | 1º Quartil | Mediana | Média | 3º Quartil | Máximo | |
|---|---|---|---|---|---|---|
| total | 23.705 | 96.042 | 127.773 | 137.072 | 176.562 | 325.722 |
Os dados têm um coeficente de variação de 0,422. A assimetria é de 0,431 enquanto que a curtose é de 2,66, indicando que a distribuição a) tende aos maiores valores e b) é platicúrtica (< 3), o que significa uma menor propenção a outliers.
O teste de Shapiro-Wilk evidencia ao nível de 95% de confiança que os dados não seguem uma distribuição Normal. Dada essa característica, foi aplicado uma transformação de Box-Cox (\(\lambda\): 0,697). A transformação usando raíz quadrada não resolveu o problema inteiramente mas foi útil para estabilizar os dados.
Finalmente, o teste de Ljung-Box evidencia com 95% de confiança que há correlação nos resíduos até, no mínimo, o lag 120.
Shapiro-Wilk normality test
data: dados$total
W = 0.97826, p-value = 7.973e-06
Shapiro-Wilk normality test
data: dados$total_sqrt
W = 0.99258, p-value = 0.03971
Box-Ljung test
data: dados_treino_ts
X-squared = 9433.2, df = 120, p-value < 2.2e-16Modelos
Regressão linear
Modelo 1
Tempo (\(t\)) e \(t^2\), em uma tentativa de captura de uma tendência quadrática.
Call:
lm(formula = total ~ tempo + tempo_quad, data = dados_treino_lm_01)
Residuals:
Min 1Q Median 3Q Max
-204.253 -40.348 0.106 38.708 153.852
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.958e+02 7.671e+00 25.52 <2e-16 ***
tempo 1.533e+00 8.923e-02 17.18 <2e-16 ***
tempo_quad -2.648e-03 2.177e-04 -12.17 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 50.63 on 393 degrees of freedom
Multiple R-squared: 0.6108, Adjusted R-squared: 0.6088
F-statistic: 308.3 on 2 and 393 DF, p-value: < 2.2e-16Análise dos resíduos
A análise dos resíduos deixa claro que há um componente de sazonalidade que não está sendo capturado pelo modelo.
Teste de Durbin-Watson lag Autocorrelation D-W Statistic p-value
1 0.8045471 0.3896325 0
Alternative hypothesis: rho != 0Modelo 2
Adicionando o componente de sazonalidade.
Call:
lm(formula = total ~ tempo + tempo_quad + sazonalidade, data = dados_treino_lm_02)
Residuals:
Min 1Q Median 3Q Max
-192.038 -39.015 2.097 35.717 138.125
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.775e+02 1.108e+01 16.016 < 2e-16 ***
tempo 1.530e+00 8.708e-02 17.569 < 2e-16 ***
tempo_quad -2.648e-03 2.124e-04 -12.467 < 2e-16 ***
sazonalidade2 -9.674e+00 1.216e+01 -0.795 0.426866
sazonalidade3 1.738e+01 1.216e+01 1.429 0.153816
sazonalidade4 7.249e+00 1.216e+01 0.596 0.551532
sazonalidade5 1.583e+01 1.216e+01 1.301 0.193888
sazonalidade6 1.621e+01 1.216e+01 1.333 0.183367
sazonalidade7 2.682e+01 1.216e+01 2.205 0.028025 *
sazonalidade8 3.490e+01 1.216e+01 2.870 0.004337 **
sazonalidade9 2.395e+01 1.216e+01 1.969 0.049700 *
sazonalidade10 2.870e+01 1.216e+01 2.360 0.018791 *
sazonalidade11 2.555e+01 1.216e+01 2.101 0.036318 *
sazonalidade12 4.039e+01 1.216e+01 3.321 0.000984 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 49.4 on 382 degrees of freedom
Multiple R-squared: 0.6398, Adjusted R-squared: 0.6275
F-statistic: 52.18 on 13 and 382 DF, p-value: < 2.2e-16Análise dos resíduos
A incorporação de um componente de sazonalidade não foi o suficiente para a melhoria da qualidade do modelo. Os resíduos permanecem exibindo um padrão sazonal e o teste de Durbin-Watson indica que há correlação de ordem 1 nos resíduos.
Teste de Durbin-Watson lag Autocorrelation D-W Statistic p-value
1 0.8459556 0.3054439 0
Alternative hypothesis: rho != 0Modelo 3
Inclusão de duas intervenções do tipo passo (contínua) e de um termo autoregressivo (resíduos defasados).
Call:
lm(formula = total ~ tempo + tempo_quad + sazonalidade + passo_t1 +
passo_t2 + residuo_def, data = dados_treino_lm_03)
Residuals:
Min 1Q Median 3Q Max
-161.812 -12.315 1.063 13.641 86.122
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.727e+02 6.140e+00 28.129 < 2e-16 ***
tempo 1.763e+00 9.324e-02 18.903 < 2e-16 ***
tempo_quad -2.791e-03 2.097e-04 -13.313 < 2e-16 ***
sazonalidade2 -9.779e+00 6.296e+00 -1.553 0.121219
sazonalidade3 1.794e+01 6.299e+00 2.848 0.004639 **
sazonalidade4 7.631e+00 6.299e+00 1.212 0.226412
sazonalidade5 1.604e+01 6.298e+00 2.546 0.011294 *
sazonalidade6 1.624e+01 6.298e+00 2.579 0.010294 *
sazonalidade7 2.668e+01 6.299e+00 4.236 2.87e-05 ***
sazonalidade8 3.458e+01 6.299e+00 5.490 7.36e-08 ***
sazonalidade9 2.345e+01 6.300e+00 3.723 0.000227 ***
sazonalidade10 2.803e+01 6.300e+00 4.449 1.13e-05 ***
sazonalidade11 2.471e+01 6.301e+00 3.921 0.000105 ***
sazonalidade12 3.937e+01 6.303e+00 6.247 1.12e-09 ***
passo_t1 -3.578e+01 8.270e+00 -4.326 1.94e-05 ***
passo_t2 -2.778e+01 7.564e+00 -3.673 0.000275 ***
residuo_def 7.171e-01 3.653e-02 19.630 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25.58 on 379 degrees of freedom
Multiple R-squared: 0.9042, Adjusted R-squared: 0.9002
F-statistic: 223.6 on 16 and 379 DF, p-value: < 2.2e-16Análise dos resíduos
No modelo 3, a inclusão das intervenções e dos resíduos defasados eliminou por completo a correlação de ordem 1 nos resíduos.
Teste de Durbin-Watson lag Autocorrelation D-W Statistic p-value
1 -0.08026602 2.148552 0.204
Alternative hypothesis: rho != 0Comparação dos modelos de regressão linear
# Comparison of Model Performance Indices
Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | R2 | R2 (adj.) | RMSE | Sigma
--------------------------------------------------------------------------------------------------------------------
modelo_reglinear_01 | lm | 4237.0 (<.001) | 4237.1 (<.001) | 4252.9 (<.001) | 0.611 | 0.609 | 50.436 | 50.628
modelo_reglinear_02 | lm | 4228.3 (<.001) | 4229.6 (<.001) | 4288.1 (<.001) | 0.640 | 0.627 | 48.521 | 49.402
modelo_reglinear_03 | lm | 3709.8 (>.999) | 3711.6 (>.999) | 3781.5 (>.999) | 0.904 | 0.900 | 25.020 | 25.575Analisando a comparação acima, confirma-se que o modelo 3 é o melhor dentre os modelos ajustados anteriormente, uma vez que ele consegue explicar 90% da variação nos dados e possui também o menor AIC (medida de comparação de modelos - quanto menor melhor).
SARIMA
Como ficou explícito durante a análise dos resíduos dos modelos anteriores a presença de sazonalidade na série temporal, o modelo SARIMA foi o segundo escolhido para a análise dos dados. Mas antes da aplicação do modelo, é necessário primeiro a realização de um teste de raíz unitária para a verificação de estacionaridade (ou falta de) da série.
Phillips-Perron Unit Root Test
data: dados_treino_ts
Dickey-Fuller = -4.7628, Truncation lag parameter = 5, p-value = 0.01
Augmented Dickey-Fuller Test
data: dados_treino_ts
Dickey-Fuller = -2.126, Lag order = 7, p-value = 0.524
alternative hypothesis: stationaryCuriosamente, os testes mostram cenários opostos. O teste de Phillips-Perron diz que a série é estacionária. Já o teste de Dickey-Fuller mostra que a série não possui raiz unitária, o que indica uma série não estacionária e, portanto, a necessidade da aplicação de uma diferença. Analisando os gráficos da série temporal durante a fase de exploração dos dados, é nítido que a série não mantém seu nível durante o período de análise, o que configura uma série não estacionária.
Identificação
Para o processo de identificação do modelo, foi utilizado o método de Box-Jenkins. Utilizar o método mencionado significa fazer uma análise das autocorrelações (FAC) e das autocorrelações parciais (FACP) para a identificação dos parâmetros do modelo.
Mesmo aplicando uma diferenciação nos dados, a identificação de um padrão na FAC e na FACP se tornou uma tarefa complexa. Uma saída para o problema foi o uso de um método empírico de definição dos parâmetros. A função forecast::auto.arima(y, seasonal = T) sugeriu o modelo ARIMA(1,1,2)(1,0,0)[12]. Como é evidente a não estacionaridade da série, o parâmetro de diferenciação sazonal \(D\) foi alterado, de 0 para 1.
Estimação
Verificação
$fit
Call:
arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))
Coefficients:
ar1 ma1 ma2 sar1
-0.5743 0.2608 -0.5165 -0.3979
s.e. 0.1243 0.1199 0.0527 0.0474
sigma^2 estimated as 867.7: log likelihood = -1840.51, aic = 3691.02
$degrees_of_freedom
[1] 379
$ttable
Estimate SE t.value p.value
ar1 -0.5743 0.1243 -4.6206 0.0000
ma1 0.2608 0.1199 2.1761 0.0302
ma2 -0.5165 0.0527 -9.8087 0.0000
sar1 -0.3979 0.0474 -8.3853 0.0000
$ICs
AIC AICc BIC
9.637133 9.637410 9.688674 
Shapiro-Wilk normality test
data: modelo_sarima$fit$residuals
W = 0.95954, p-value = 5.602e-09Uma análise dos resultados confirma que ainda há correlação e não Normalidade dos resíduos. Em simulações realizadas, uma diferenciação \(D = 2\) não resolveu o problema. Essa quebra dos pressupostos impactará significamente na confiabilidade das predições.
Alisamento exponencial
Antes do ajuste do terceiro e último modelo, é interessante apresentar a decomposição da série temporal em estudo.
Após a análise do gráfico acima, o modelo de Holt-Winters foi ajustado considerando um nível, uma tendência e uma sazonalidade aditiva.
Sumário ajuste Holt-Winters Length Class Mode
fitted 1536 mts numeric
x 396 ts numeric
alpha 1 -none- numeric
beta 1 -none- numeric
gamma 1 -none- numeric
coefficients 14 -none- numeric
seasonal 1 -none- character
SSE 1 -none- numeric
call 3 -none- call Análise dos resíduos
Comparação
A métrica escolhida para ser um critério de qualidade entre os modelos foi o Erro Quadrático Médio (EQM).
Regressão linear: 626,014;
SARIMA: 839,232;
Alisamento exponencial: 763,407.
O menor EQM mostra que a regressão linear possui o menor erro de predição.
Predição
Para finalizar, os modelos ajustados foram utilizados para a realização de \(h = 14\) predições, que foram comparadas aos dados reais observados.
O gráfico com as predições deixa claro que, ao mesmo tempo que a regressão linear possui as melhores métricas de ajuste, o alisamento exponencial gerou as melhores predições. O pior modelo nesse cenário é o SARIMA, mas com as importantes ressalvas que o desempenho ruim é devido ao fato da não Normalidade dos dados (astsa::sarima() não oferece a opção de ajuste pelo método de mínimos quadrados - cuja suposição de Normalidade é opcional) e da dificuldade em encontrar os melhores parâmetros através do método de Box-Jenkins. Sugestões de melhorias nos comentários são sempre bem-vindas.
Notas de rodapé
R Core Team (2023). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/↩︎
Banco Central do Brasil. Portal de Dados Abertos. Vendas de veículos pelas concessionárias - Automóveis. Disponível em: https://dadosabertos.bcb.gov.br/dataset/7384-vendas-de-veiculos-pelas-concessionarias---automoveis. Acesso em: 26/04/2023.↩︎