Faz parte de todo grande projeto saber qual o melhor momento de lançar-se no mercado, publicar um lançamento ou entrar em competição. Para tal, é preciso saber se a economia da região em questão está aquecida de forma a aproveitar ao máximo a onda do crescimento local e conseguir se beneficiar e se preparar também para as sazonalidades adversas, próprias do ambiente.
Para cumprir esse objetivo, e como forma de praticar a teoria econométrica, farei a análise e modelagem do Índice de Volume de Vendas no Comércio Varejista, da Pesquisa Mensal de Comércio, PMC-IBGE, para o Brasil e Ceará, com a base do índice sendo 2014 (2014=100). Seguirei de perto a metodologia proposta por Silveira Bueno , em seu livro Econometria De Séries Temporais; o livro do Enders: Applied Econometric Time Series; e o livro de Hyndman & Athanasopoulos: Forecasting: Principles and Practice. Sobre os procedimentos usados, é importante ter em mente que a modelagem econométrica não é um algoritmo linear, como uma receita de bolo; está mais para um movimento helicoidal ascendente, onde você volta nas mesmas ferramentas em diferentes estágios de maturação do modelo. No geral, começamos a modelagem abertos para todas as possibilidades: processos autoregressivos, ou médias móveis, ou ambos; tendências; sazonalidades; mudança estrutural e instabilidade paramétrica. Enquanto avançamos nos testes de significância e diagnósticos de resíduos, cortamos e refinamos essas possibilidades até encontrar um modelo parcimonioso e que gera ruídos brancos. Acerca da rotina de modelagem e do perigo de cair no overfitting, recomendo fortemente o estudo da seção 8 do capítulo 2 do Enders.
Usarei os pacotes Tidyverse e, para séries temporais, Tidyverts que são a espinha dorsal da gramática de manipulação de dados no R, com vários outros pacotes agregáveis, vale a pena conhecê-los profundamente, juntamente com Stats e Forecast.
library(stats)
library(tseries)
library(forecast)
library(tidyverse)
library(plotly)
library(knitr)
library(kableExtra)
library(gridExtra)
library(ggpubr)
# tidyverts core
library(tsibble)
library(fable)
library(feasts)
# tidyverts improvements tools
library(tsibbletalk)
library(fable.prophet)
library(fabletools)Os dados são da Tabela 3416, precisam ser importados e confirmados como data de frequência mensal. Caso tenha alguma dúvida em alguma função usada, digite help(nome_da_função) no console R para conhecê-la melhor. Use também as referências no site do pacote.
# importando
pmc <- read_csv('pmc.csv')
# nova variável
pmc$dt <- seq(from = lubridate::ym('2000 Jan'),
to = lubridate::ym('2021 Aug'), # alternativa: length=
by = "month")
# forçar a nova variável a ser mensal
pmc <- pmc |>
dplyr::mutate(dt = tsibble::yearmonth(as.character(dt))) |>
as_tsibble(index = dt)
pmc$Data <- NULL # apagar a data antiga
pmc <- pmc |> # defini-la como dados tsibble
as_tsibble()
# confirma-se que a periodicidade é mensal, pela forçada q demos
pmc |> interval()<interval[1]>
[1] 1MVou criar uma função set_tab(dados_tab, cap_tab) que simplificará a tabulação dos nossos dados. Os argumentos da função são os dados a serem tabelados e o título da tabela :
Visão geral das variáveis do dado:
Rows: 260
Columns: 29
$ Brasil <dbl> 46.0, 46.5, 48.5, 48.5, 50.5, 48.1, 49.8, 49.7, …
$ Rondônia <dbl> 24.0, 25.3, 25.5, 26.1, 26.0, 25.2, 26.5, 26.4, …
$ Acre <dbl> 28.0, 27.8, 29.9, 32.3, 30.4, 33.5, 34.2, 32.2, …
$ Amazonas <dbl> 39.8, 42.4, 42.6, 41.6, 47.2, 46.2, 47.3, 47.9, …
$ Roraima <dbl> 32.5, 33.5, 38.3, 29.8, 33.9, 33.4, 33.2, 33.1, …
$ Pará <dbl> 39.1, 38.0, 38.9, 41.1, 43.3, 41.7, 43.8, 45.8, …
$ Amapá <dbl> 32.5, 32.2, 32.3, 34.8, 38.0, 38.9, 38.1, 40.0, …
$ Tocantins <dbl> 20.7, 19.2, 20.1, 20.9, 22.4, 18.6, 22.8, 22.6, …
$ Maranhão <dbl> 29.1, 28.0, 27.1, 26.4, 28.9, 28.5, 30.8, 29.4, …
$ Piauí <dbl> 47.6, 44.3, 36.7, 34.8, 39.0, 37.7, 39.0, 37.1, …
$ Ceará <dbl> 37.5, 35.3, 34.3, 35.5, 38.7, 37.0, 40.0, 40.9, …
$ `Rio Grande do Norte` <dbl> 38.7, 37.7, 38.4, 37.3, 39.1, 39.1, 38.9, 38.9, …
$ Paraíba <dbl> 34.9, 32.7, 30.7, 31.8, 35.0, 35.7, 35.6, 36.1, …
$ Pernambuco <dbl> 46.3, 45.0, 44.3, 44.2, 47.2, 47.2, 45.1, 47.3, …
$ Alagoas <dbl> 36.8, 35.7, 35.6, 36.0, 39.3, 39.8, 37.6, 38.4, …
$ Sergipe <dbl> 44.6, 41.6, 43.0, 44.5, 45.9, 44.9, 44.6, 44.7, …
$ Bahia <dbl> 44.8, 44.7, 45.4, 45.4, 46.6, 48.2, 46.9, 47.7, …
$ `Minas Gerais` <dbl> 45.2, 45.6, 46.9, 47.8, 48.5, 47.4, 48.8, 48.1, …
$ `Espírito Santo` <dbl> 45.2, 45.1, 45.6, 44.5, 45.7, 45.7, 47.1, 46.3, …
$ `Rio de Janeiro` <dbl> 50.0, 50.4, 50.7, 54.5, 56.9, 54.9, 58.4, 58.3, …
$ `São Paulo` <dbl> 45.4, 46.5, 49.1, 48.0, 50.2, 46.3, 48.9, 48.9, …
$ Paraná <dbl> 47.9, 48.9, 51.8, 52.9, 53.0, 51.0, 52.2, 51.8, …
$ `Santa Catarina` <dbl> 47.5, 47.8, 48.6, 48.5, 49.1, 47.3, 48.2, 47.2, …
$ `Rio Grande do Sul` <dbl> 57.6, 57.1, 62.1, 62.6, 64.9, 60.9, 62.1, 60.2, …
$ `Mato Grosso do Sul` <dbl> 35.2, 35.6, 39.2, 40.5, 41.6, 40.4, 40.8, 38.7, …
$ `Mato Grosso` <dbl> 43.9, 48.0, 50.5, 49.5, 53.2, 54.5, 53.3, 50.7, …
$ Goiás <dbl> 42.9, 42.0, 46.4, 47.0, 50.0, 48.1, 48.5, 47.4, …
$ `Distrito Federal` <dbl> 53.9, 54.9, 57.1, 57.9, 59.6, 58.7, 59.9, 59.8, …
$ dt <mth> 2000 jan, 2000 fev, 2000 mar, 2000 abr, 2000 mai…Selecionamos Brasil e Ceará:
Índice de Volume de Vendas no Comércio Varejista. Ceará e Brasil. Jan2000-Ago2021
Vemos que a série cearense parece ter uma maior volatilidade que a brasileira. Para enxergar o que acontece dentro de cada frequência (ano e mês), vamos fazer o plot agrupando essas frequências para o Ceará:
A primeira impressão das visualizações acima é que temos um período estacionário até o fim de 2004, quando começa uma tendência positiva até o fim de 2013, quando passa a cair levemente. Claramente, há sazonalidade dentro do ano em todos os períodos. Confirmaremos nossa intuição visual de estacionariedade (<2005), tendência (2005~2013), significância da queda (>2013) e sazonalidade com os testes de raíz unitária, FAC e FACP.
Somente para fins didáticos, vamos fazer uma decomposição clássica aditiva da série, repartindo série = tend + sazon + aleat, que nos ajuda a entender melhor a série. A tendência (\hat{T_t}) é feita com a média móvel simples de seis meses. A sazonalidade é construída subtraindo a série original da tendência (\ddot{y_t} = y_t - \hat{T_t} ), agrupando cada mês e subtraindo da sua média (\hat{S}_t = \ddot{y}_{mês} - \ddot{\overline{y}}_{mês}). O resíduo são as sobras desse processo (\hat{R}_t = y_t - \hat{T}_t - \hat{S}_t). Como vemos no gráfico abaixo, esse processo não torna o resíduo um Ruído Branco, além desse tipo de tendência não ser indicado em Passeios Aleatórios nem para previsões.
Como podemos ver pelo gráfico da tendência, ela pode ser repartida nas três fases comentadas. Não existe, para a série inteira, somente um componente gerador de tendência.
Os testes Kwiatkowski-Phillips-Schmidt-Shin (KPSS) e Phillips-Perron (PP) são usados para confirmar nossa intuição visual de tendência nos determinados períodos. Resumidamente, as raízes do polinômio real de operadores de defasagem devem estar dentro do ciclo unitário para que o processo AR(.) seja estacionário e para que o processo MA(.) seja invertível. Caso seja detectado raíz fora do ciclo, teremos que remover a tendência, diferenciando a série. Usaremos o teste PP para raíz unitária por ser robusto, se comparado a Dickey-Fuller (DF) na presença de tendência temporal, intercepto e correlação serial nos erros e o teste KPSS para estacionariedade pelo baixo poder do DF na presença de médias móveis perto do círculo unitário.
A hipótese nula do KPSS é do processo ser estacionário ao entorno de uma tendência determinística, simplificamos com H_0: y_t \tilde{} I(0) enquanto que a do PP é de raíz unitária. Então, para confirmarmos nossa intuição, o p-valor de um deve ser baixo e do outro deve ser alto.
Como vemos abaixo, os testes para o primeiro e segundo períodos foram conclusivos para o Ceará e o Brasil: até 2004 é estacionária; 2005~2012 é integrável. Temos um problema na série cearence 2013~2021: os p-valores estão sinalizando características diferentes, podemos culpar o tamanho dos dados e uma posível heterocedasticidade. Veremos qual processo consegue melhor produzir ruídos brancos e i.i.d. O mesmo não ocorre para o Brasil, que mostra estacionariedade para o último período.
# CE ~~~~~~~~~~~~~~~~~~~~~
ce <- ce_br |>
dplyr::filter(dt <= tsibble::yearmonth('2004 Dec')) |>
features(ce, c(unitroot_kpss, unitroot_pp))
# BR ~~~~~~~~~~~~~~~~~~~~~
br <- ce_br |>
dplyr::filter(dt <= tsibble::yearmonth('2004 Dec')) |>
features(br,c(unitroot_kpss, unitroot_pp))
dplyr::bind_rows(list(Ce = ce, Br = br), .id = 'Região') |>
set_tab(cap_tab = "Testes KPSS e PP entre 2000-2004.")| Região | kpss_stat | kpss_pvalue | pp_stat | pp_pvalue |
|---|---|---|---|---|
| Ce | 0.1497 | 0.1 | -5.3517 | 0.01 |
| Br | 0.1806 | 0.1 | -6.2344 | 0.01 |
# CE ~~~~~~~~~~~~~~~~~~~~~
ce <- ce_br |>
dplyr::filter(dt >= tsibble::yearmonth('2005 Jan') &
dt <= tsibble::yearmonth('2012 Dec')) |>
features(ce, c(unitroot_kpss, unitroot_pp))
# BR ~~~~~~~~~~~~~~~~~~~~~
br <- ce_br |>
dplyr::filter(dt >= tsibble::yearmonth('2005 Jan') &
dt <= tsibble::yearmonth('2012 Dec')) |>
features(br, c(unitroot_kpss, unitroot_pp))
dplyr::bind_rows(list(Ce = ce, Br = br), .id = 'Região') |>
set_tab(cap_tab = "Testes KPSS e PP entre 2005-2012.")| Região | kpss_stat | kpss_pvalue | pp_stat | pp_pvalue |
|---|---|---|---|---|
| Ce | 2.3817 | 0.01 | -1.7654 | 0.1 |
| Br | 2.3879 | 0.01 | -2.2605 | 0.1 |
# CE ~~~~~~~~~~~~~~~~~~~~~
ce <- ce_br |>
dplyr::filter(dt >= tsibble::yearmonth('2013 Jan')) |>
features(ce, c(unitroot_kpss, unitroot_pp))
# BR ~~~~~~~~~~~~~~~~~~~~~
br <- ce_br |>
dplyr::filter(dt >= tsibble::yearmonth('2013 Jan')) |>
features(br, c(unitroot_kpss, unitroot_pp))
dplyr::bind_rows(list(Ce = ce, Br = br), .id = 'Região') |>
set_tab(cap_tab = "Testes KPSS e PP entre 2013-2021.")| Região | kpss_stat | kpss_pvalue | pp_stat | pp_pvalue |
|---|---|---|---|---|
| Ce | 0.9405 | 0.01 | -5.5370 | 0.01 |
| Br | 0.1457 | 0.10 | -7.1568 | 0.01 |
Precisamos diferenciar a série, remover sua tendência. Surge então dois desafios: 1) identificar qual processo gera a tendência, se determinística ou estocática; e 2) se aplicamos diferenciação sazonal e em qual lag. Pela característica das tendências, suspeito que seja uma estocástica, um random walk, que trata-se aplicando diferenciação. Onde e quantas diferenciações aplicar, além de quais ordens AR() e MA() fazer, respondemos pelas autocorrelações (FAC & FACP) e testes de autocorrelações (Ljung-Box).
Relembrando que a FAC é usada para evidência de MA(q) e suas sazonalidades; a FACP é usada para identificação de AR(p) e suas sazonalidades. Existem dois tipos de sazonalidade: a aditiva e a multiplicativa. A primeira contém componentes AR(.) e MA(.) degenerados, o que causam significância estatística isoladas na FACP e FAC. Já a sazonalidade multiplicativa causa significância também nos períodos ao entorno do “coeficiente principal”, como mostro no exemplo abaixo. (Lembrando que L é um operador de lag, não um escalar, de forma que y_t·L^h = y_{t-h}):
AR(1)(2)_6 \ \ \therefore (1- \phi_1 L)(1- \phi_6 L^6 - \phi_{12} L^{12}) y_t = \epsilon_t \ \ \rightarrow y_t = \phi_{1}y_{t-1} + \phi_{6}y_{t-6} - \phi_{1}\phi_{6}y_{t-7} + \phi_{12}y_{t-12} - \phi_{1}\phi_{12}y_{t-13} + \epsilon_t
MA(2)(2)_{11} \ \ \therefore y_t = (1- \theta_{1}L^{} - \theta_{2}L^{2})(1- \theta_{11}L^{11} - \theta_{22}L^{22})\epsilon_t \ \ \rightarrow y_t = \epsilon_t - \theta_{1}\epsilon_{t-1} - \theta_{2}\epsilon_{t-2} - \theta_{11}\epsilon_{t-11} + \theta_{1}\theta_{11}\epsilon_{t-12} + \theta_{2}\theta_{11}\epsilon_{t-13} - \theta_{22}\epsilon_{t-22} + \theta_{1}\theta_{22}\epsilon_{t-23} + \theta_{2}\theta_{22}\epsilon_{t-24}
Nos exemplos acima, vemos que os modelos com sazonalidade multiplicativa são particularmente parcimoniosos por considerarem y_{t-7}, y_{t-13}, \epsilon_{t-12}, \epsilon_{t-13}, \epsilon_{t-23} e \epsilon_{t-24} sem a adição de novos coeficientes para isso, sem diminuir o grau de liberdade nas estatísticas dos coeficientes estimados para a série.
Vamos fazer a primeira diferença das séries e analisar seus comportamentos.
Vemos que a variância em ambas as diferenciações parecem aumentar com o tempo, se essa diferença for estatisticamente significante, teremos problemas por não considerar a heterocedasticidade.
Agora, faremos o plot da FAC e FACP para o Ceará e Brasil:
Perceba que, em ambas as séries, para a FAC, vemos picos nos múltiplos de 12 que demoram a cair e coeficientes significantes ao entorno desses “picos principais”, o que reforça o uso das sazonalidades multiplicativas e mais uma diferenciação no lag 12.
Os FACP’s mostram significância até a primeira repetição do ciclo de 12 meses. É bem possível que esse AR(p)(1)_{12} seja invertível e esteja causando esse MA(q)(Q)_{12} com altíssima persistência em sQ \ \forall \ s=12,24,36... como apontado na FAC. Vamos diferenciar sazonalmente novamente as séries, em 12: \Delta^{12} y = y_{t} - y_{t-12}.
Faremos os mesmos plots da 1 diferenciação, perceba que a heterocedasticidade permanece visível:
FAC e FACP para o Ceará e Brasil:
Pelas novas FAC & FACP cearenses, podemos suspeitar de um MA(q)(Q) com q=1,3 e Q=1_{t=12} e um AR(p)(P).
Pelas novas FAC & FACP brasileiras, podemos suspeitar de um MA(q)(Q) com q=1,3 e Q=1_{t=12} e um AR(p)(P) com p=1,2,4,5 e P=1_{t=12}.
Os decaimentos e “quase significância” de demais picos podem criar raízes polinomiais muito próximas do ciclo unitário nos modelos que estimarmos.
A seguir, faço um conjunto de modelos propostos para ambos os locais. Os últimos modelos, ce_auto e br_auto, são feitos automaticamente, com a interação combinatória de diferentes ordens, onde se escolhe o melhor modelo segundo o critério de BIC (Bayesian Information Criterion). Vamos fazer os modelos e escolher o mais parcimonioso, ou seja, com o menor BIC e diagnosticar os resíduos. Não esqueceremos de diferenciar em (y_t - y_{t-1}) - y_{t-12}.
# CE ~~~~~~~~~~~~~~~~~~~~~
ce_models <- ce_br |>
model(ce_1 = ARIMA(ce ~ 0 + pdq(3,1,3) + PDQ(1,1,0,
period = "1 year")),
ce_11= ARIMA(ce ~ 1 + pdq(3,1,3) + PDQ(1,1,0,
period = "1 year")),
ce_2 = ARIMA(ce ~ 0 + pdq(3,1,3) + PDQ(0,1,1,
period = "1 year")),
ce_3 = ARIMA(ce ~ 0 + pdq(3,1,3) + PDQ(1,1,1,
period = "1 year")),
ce_4 = ARIMA(ce ~ 0 + pdq(3,1,0) + PDQ(1,1,0,
period = "1 year")),
ce_auto = ARIMA(ce,
stepwise = F,
approximation = F,
order_constraint =p+q+P+Q<= 8 &
(constant + d + D <= 3)))
# BR ~~~~~~~~~~~~~~~~~~~~~
br_models <- ce_br |>
model(br_1 = ARIMA(br ~ 0 + pdq(1,1,1) + PDQ(1,1,1,
period = "1 year")),
br_11= ARIMA(br ~ 1 + pdq(1,1,1) + PDQ(1,1,1,
period = "1 year")),
br_2 = ARIMA(br ~ 0 + pdq(1,1,4) + PDQ(1,1,1,
period = "1 year")),
br_3 = ARIMA(br ~ 0 + pdq(4,1,1) + PDQ(1,1,1,
period = "1 year")),
br_4 = ARIMA(br ~ 0 + pdq(4,1,4) + PDQ(1,1,1,
period = "1 year")),
br_auto = ARIMA(br,
stepwise = F,
approximation = F,
order_constraint =p+q+P+Q<= 8 &
(constant + d + D <= 3)))Tabulando os modelos criados:
# A mable: 6 x 2
# Key: Modelos [6]
Modelos Ordens
<chr> <model>
1 ce_1 <ARIMA(3,1,3)(1,1,0)[12]>
2 ce_11 <ARIMA(3,1,3)(1,1,0)[12] w/ poly>
3 ce_2 <ARIMA(3,1,3)(0,1,1)[12]>
4 ce_3 <ARIMA(3,1,3)(1,1,1)[12]>
5 ce_4 <ARIMA(3,1,0)(1,1,0)[12]>
6 ce_auto <ARIMA(2,1,2)(1,1,1)[12]># A mable: 6 x 2
# Key: Modelos [6]
Modelos Ordens
<chr> <model>
1 br_1 <ARIMA(1,1,1)(1,1,1)[12]>
2 br_11 <ARIMA(1,1,1)(1,1,1)[12] w/ poly>
3 br_2 <ARIMA(1,1,4)(1,1,1)[12]>
4 br_3 <ARIMA(4,1,1)(1,1,1)[12]>
5 br_4 <ARIMA(4,1,4)(1,1,1)[12]>
6 br_auto <ARIMA(3,1,4)(0,1,1)[12]>Organizando pelos critérios de informação:
# A tibble: 6 × 6
.model sigma2 log_lik AIC AICc BIC
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ce_auto 8.99 -588. 1191. 1191. 1215.
2 ce_3 9.05 -588. 1194. 1195. 1225.
3 ce_1 9.17 -590. 1195. 1196. 1223.
4 ce_2 9.18 -590. 1196. 1196. 1223.
5 ce_11 9.21 -590. 1197. 1198. 1229.
6 ce_4 9.41 -594. 1198. 1199. 1216.# A tibble: 6 × 6
.model sigma2 log_lik AIC AICc BIC
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 br_auto 5.82 -536. 1091. 1092. 1122.
2 br_4 5.68 -535. 1092. 1093. 1130.
3 br_3 6.08 -542. 1099. 1100. 1127.
4 br_2 6.14 -543. 1101. 1102. 1129.
5 br_1 6.25 -546. 1103. 1103. 1120.
6 br_11 6.28 -546. 1105. 1105. 1125.Vamos ficar com os modelos com o menor BIC. Pegarei os dois melhores modelos de cada região, mas mostrarei apenas o diagnóstico dos resíduos do melhor modelo, no caso, os que foram gerados automaticamente.
# CE ~~~~~~~~~~~~~~~~~~~~~
ce_models <- ce_br |>
model(ce_1 = ARIMA(ce ~ 0 + pdq(3,1,3) + PDQ(1,1,1,
period = "1 year")),
ce_auto = ARIMA(ce ~ 0 + pdq(2,1,2) + PDQ(1,1,1,
period = "1 year")))
# BR ~~~~~~~~~~~~~~~~~~~~~
br_models <- ce_br |>
model(br_1 = ARIMA(br ~ 0 + pdq(4,1,4) + PDQ(1,1,1,
period = "1 year")),
br_auto = ARIMA(br ~ 0 + pdq(3,1,4) + PDQ(0,1,1,
period = "1 year")))Vamos visualizar os resíduos:
Série Real e Estimada para o Ceará
Série Real e Estimada para o Brasil
Agora procederemos com três testes para os resíduos: para constatar normalidade, usaremos o Teste de Jarque-Bera (H_0: E(\epsilon^s_t)^3 =0); para conferir se há autocorrelação: o Ljung-Box (H_0: E(\epsilon_t \epsilon_{t+h}) = 0); para a heterocedasticidade, que foi um problema que suspeitamos, o teste Arch-LM ($H_0:_{j} = 0 j em (1-_j L^j)var(_i)=u_i $).
Jarque Bera Test
data: resid$.resid
X-squared = 415.18, df = 2, p-value < 2.2e-16 lb_stat lb_pvalue
0.1051355 0.7457518 stat_arch_lm
0.2592619 
Os testes para o modelo cearense mostraram resíduos com assimetria, sem autocorrelação entre os erros e heterocedasticidade, com a variância autocorrelacionada. Vemos a necessidade de tratar a variância. O que faremos em trabalhos futuros. Como suspeitamos, as raízes polinomiais estão muito perto do limite do ciclo unitário, o que indica alta persistência dos choques e valores passados.
Jarque Bera Test
data: resid$.resid
X-squared = 535.68, df = 2, p-value < 2.2e-16 lb_stat lb_pvalue
0.01317495 0.90861778 stat_arch_lm
0.12099 
Para o modelo brasileiro, vemos também assimetria nos resíduos, baixíssimas chances de autocorrelação e também uma variância autocorrelacionada, apesar de um p-valor menor que o teste para o Ceará.
