Exemplo Séries Temporais: consumo Morettin e Toloi: hierárquico - formato ts

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

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Exemplo Séries Temporais: consumo Morettin e Toloi: hierárquico - formatos ts. Campo Grande-MS, Brasil: RStudio/Rpubs, 2023. Disponível em https://rpubs.com/amrofi/morettin_toloi_hierarquico.

1 Introdução

Os dados vem do livro de Morettin e Toloi, Análise de Séries Temporais (https://www.ime.usp.br/~pam/ST.html), consumo no varejo de São Paulo, mensais de Jan/1984 a Out/1996, em https://www.ime.usp.br/~pam/CONSUMO.XLS.

Inicialmente, trabalharemos com o formato ts, conforme o livro de Hyndman (2018), segunda edição. O artigo segue a rotina de Hyndman (2016).

2 Dados em formato ts

# vou puxar os libraries que tenho costume de usar: destaque para o readxl e
# tseries
library(readxl)  # puxar Excel
library(fpp2)  # livro Hyndman formato ts - já carrega ggplot2 e forecast
# codigo para puxar do Excel dados <- read_excel('CONSUMO morettin R.xlsx',
# sheet = 'dados')
dados <- structure(list(obs = structure(c(441763200, 444441600, 446947200, 449625600,
    452217600, 454896000, 457488000, 460166400, 462844800, 465436800, 468115200,
    470707200, 473385600, 476064000, 478483200, 481161600, 483753600, 486432000,
    489024000, 491702400, 494380800, 496972800, 499651200, 502243200, 504921600,
    507600000, 510019200, 512697600, 515289600, 517968000, 520560000, 523238400,
    525916800, 528508800, 531187200, 533779200, 536457600, 539136000, 541555200,
    544233600, 546825600, 549504000, 552096000, 554774400, 557452800, 560044800,
    562723200, 565315200, 567993600, 570672000, 573177600, 575856000, 578448000,
    581126400, 583718400, 586396800, 589075200, 591667200, 594345600, 596937600,
    599616000, 602294400, 604713600, 607392000, 609984000, 612662400, 615254400,
    617932800, 620611200, 623203200, 625881600, 628473600, 631152000, 633830400,
    636249600, 638928000, 641520000, 644198400, 646790400, 649468800, 652147200,
    654739200, 657417600, 660009600, 662688000, 665366400, 667785600, 670464000,
    673056000, 675734400, 678326400, 681004800, 683683200, 686275200, 688953600,
    691545600, 694224000, 696902400, 699408000, 702086400, 704678400, 707356800,
    709948800, 712627200, 715305600, 717897600, 720576000, 723168000, 725846400,
    728524800, 730944000, 733622400, 736214400, 738892800, 741484800, 744163200,
    746841600, 749433600, 752112000, 754704000, 757382400, 760060800, 762480000,
    765158400, 767750400, 770428800, 773020800, 775699200, 778377600, 780969600,
    783648000, 786240000, 788918400, 791596800, 794016000, 796694400, 799286400,
    801964800, 804556800, 807235200, 809913600, 812505600, 815184000, 817776000,
    820454400, 823132800, 825638400, 828316800, 830908800, 833587200, 836179200,
    838857600, 841536000, 844128000), class = c("POSIXct", "POSIXt"), tzone = "UTC"),
    t = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
        21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
        40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58,
        59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
        78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96,
        97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112,
        113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127,
        128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
        143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154), consumo = c(114.13,
        110.79, 116.46, 111.57, 120.66, 121.15, 121.27, 127.02, 129.04, 133.3, 130.6,
        179.39, 120.64, 114.05, 130.6, 118.26, 145.54, 135.13, 153.35, 159.95, 150.01,
        164.93, 170.37, 220.96, 134.26, 133.11, 147.84, 164.46, 181.86, 170.44, 186.64,
        174.21, 181.62, 194.16, 181.9, 232.01, 140.16, 130.78, 119.04, 120.73, 129.81,
        111.04, 122.75, 133.95, 125.41, 132.05, 129.54, 176.37, 110.09, 113.25, 124.03,
        110.63, 116.72, 124.63, 124.38, 130.27, 119.87, 115.75, 122.44, 162.43, 105.89,
        115.59, 147, 131.7, 131.32, 136.66, 126.43, 134.88, 128.26, 125.32, 124.61,
        166.11, 116.25, 96.93, 89.27, 101.87, 125.57, 113.31, 109.39, 127.33, 120.56,
        117.73, 113.81, 147.25, 100.15, 95.11, 112.26, 109.39, 114.2, 113.8, 126.47,
        128.36, 115.71, 116.09, 99.53, 127.27, 87.08, 85.67, 82.02, 98.2, 96.44,
        90.23, 97.15, 95.08, 94, 93, 96.09, 129.21, 75.39, 77.7, 97.34, 84.97, 87.55,
        86.64, 90.52, 95.4, 95.2, 95.8, 101.23, 128.49, 85.63, 82.77, 96.55, 81.33,
        96.91, 83.76, 90.19, 114.84, 108.4, 106.05, 109.71, 143.86, 99.12, 99.28,
        114.75, 106.13, 110.02, 108.07, 112.52, 113.87, 107.84, 112.12, 112.03, 139.37,
        92.24, 93.56, 107.37, 102.89, 114.78, 102.88, 118.41, 119.23, 117.36, 122.06)),
    row.names = c(NA, -154L), class = c("tbl_df", "tbl", "data.frame"))
# atribuir objeto ts
dados.ts <- ts(dados$consumo, start = c(1984, 1), frequency = 12)
# plot inicial
plot(dados.ts, main = "Consumo do varejo de Sao Paulo, Morettin e Toloi (2006)",
    xlab = "Ano", ylab = "Indice")

# plot no ggplot
library(ggplot2)
datas <- seq(as.Date(paste(c(start(dados.ts), 1), collapse = "/")), by = "month",
    length.out = length(dados$consumo))
dados.df <- data.frame(date = datas, value = dados$consumo)
ggplot(dados.df) + aes(x = date, y = value) + geom_line(colour = "#112446") + labs(x = "Mês/Ano",
    y = "Índice (base 100)", title = "Série de consumo:", subtitle = "varejo de São Paulo",
    caption = "Fonte: Morettin e Toloi (2006). Disponível em: <https://www.ime.usp.br/~pam/CONSUMO.XLS>") +
    ggthemes::theme_economist_white()

2.0.1 Estatísticas descritivas

Olhando as estatísticas descritivas, o pesquisador percebe a amplitude de variação, e o gráfico do dygraphs permite “caminhar” com o mouse sobre a série. Perceberás que a sazonalidade está bem marcada com picos em dezembro de cada ano.

# estatisticas basicas
summary(dados.ts)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  75.39  102.12  116.36  120.94  130.52  232.01 

# Min. 1st Qu.  Median Mean 3rd Qu.  Max.  75.39 102.12 116.36 120.94 130.52
# 232.01 pelo pacote dygraph dá mais opções
library(dygraphs)
dygraph(dados.ts)

É possível visualizar nos plots acima: sazonalidade (por exemplo, picos em dezembro de cada ano); a tendência aparentemente crescente e depois decresce com a “crise” brasileira; e uma aparente não-estacionariedade (média e variância mudam no tempo). Em outra postagem aplicarei o teste de raiz unitária na série para avaliar a estacionariedade de modo mais explícito.

3 Hierarquico - thief

O pacote thief é um acrônimo para Temporal HIErarchical Forecasting. A ideia é pegar uma série temporal sazonal e calcular todas as agregações temporais possíveis que resultem em um número inteiro de observações por ano. Por exemplo, uma série temporal trimestral é agregada a bianual e anual; enquanto uma série temporal mensal é agregada a bimestral, trimestral, quadrimestral, semestral e anual. Cada uma das séries temporais resultantes é prevista e, em seguida, as previsões são reconciliadas usando o algoritmo de reconciliação hierárquica descrito por Hyndman (2016). Isso tende a fornecer melhores previsões, mesmo que nenhuma informação nova tenha sido adicionada, especialmente para séries temporais com longos períodos sazonais. Também permite que diferentes tipos de previsões para diferentes horizontes de previsão sejam combinados de maneira consistente.

A função principal é thief que pode ser usada como qualquer uma das funções de previsão do pacote forecast.

# Hierarquico
library(thief)
library(ggplot2)
fc <- thief(dados.ts)
autoplot(fc)

summary(fc)

Forecast method: THieF-ETS

Model Information:
NULL

Error measures:
                     ME     RMSE      MAE        MPE     MAPE      MASE
Training set -0.5387926 11.85398 8.645059 -0.1400359 7.037349 0.5160809
                  ACF1
Training set 0.4713074

Forecasts:
           Jan       Feb       Mar       Apr       May       Jun       Jul
1996                                                                      
1997 100.83958  98.39163 109.67194 106.59647 115.93014 109.41115 115.86202
1998 100.89544  98.46595 109.78544 106.67949 116.05506 109.48834 115.95766
           Aug       Sep       Oct       Nov       Dec
1996                               118.62168 155.55950
1997 121.89536 117.11705 119.26108 118.69737 155.83767
1998 122.00948 117.20760 119.34892                    

A função thief calcula os agregados temporais, depois calcula todas as previsões (usando ets por padrão) e, finalmente, reconcilia as previsões usando a abordagem descrita em Hyndman (2016). As previsões reconciliadas correspondentes à série original são retornadas.

Existem opções para calcular previsões usando outros métodos, como auto.arima, o método Theta (thetaf) ou uma função de previsão definida pelo usuário.

Existem também funções para fazer cada parte separadamente: tsaggregates calculará todas as séries temporais agregadas, reconcilethief reconciliará uma lista de séries temporais (ou uma lista de objetos de previsão) de diferentes níveis de agregação temporal.

# Construct all temporal aggregates
totalagg <- tsaggregates(dados.ts)
autoplot(totalagg, main = "Varejo Morettin - tsaggregates")

4 Cálculo dos Forecasts

Faremos com ets.

# Compute base forecasts
base <- list()
for (i in seq_along(totalagg)) base[[i]] <- forecast(ets(totalagg[[i]]), h = 2 *
    frequency(totalagg[[i]]), level = 80)

# Reconcile forecasts
reconciled <- reconcilethief(base)

# Plot original and reconciled forecasts
par(mfrow = c(2, 3), mar = c(3, 3, 1, 0))
for (i in 6:1) {
    plot(reconciled[[i]], main = names(totalagg)[i], xlim = c(1984.1, 1998.1), flwd = 1)
    lines(base[[i]]$mean, col = "red")
}

A linha vermelha representa as previsões de base originais e a linha azul as previsões reconciliadas. Observe como as previsões básicas são diferentes nos níveis de agregação. No nível anual, não há tendência capturada na previsão devido ao ajuste limitado da amostra, e há uma tendência fraca capturada nas previsões semanais devido aos dados ruidosos. Mas a tendência é bastante forte nos níveis intermediários de agregação temporal. Da mesma forma, a sazonalidade é capturada com relativa precisão nos níveis trimestral, mensal e quinzenal, mas não no nível semanal. Uma vez conciliados, a tendência nos níveis anual e semanal é mais forte, e a sazonalidade semanal parece mais razoável.

Faremos com auto.arima.

# Compute base forecasts
baseaa <- list()
for (i in seq_along(totalagg)) baseaa[[i]] <- forecast(auto.arima(totalagg[[i]],
    stepwise = F, approximation = F), h = 2 * frequency(totalagg[[i]]), level = 80)

# Reconcile forecasts
reconciledaa <- reconcilethief(baseaa)

# Plot original and reconciled forecasts
par(mfrow = c(2, 3), mar = c(3, 3, 1, 0))
for (i in 6:1) {
    plot(reconciledaa[[i]], main = names(totalagg)[i], xlim = c(1984.1, 1998.1),
        flwd = 1)
    lines(baseaa[[i]]$mean, col = "red")
}

5 Referências

HYNDMAN, Rob J. fpp2: Data for “Forecasting: Principles and Practice” (2nd Edition). CRAN, R package version 2.3. 2018. Disponível em: https://CRAN.R-project.org/package=fpp2.

HYNDMAN, R.J. The thief package for R: Temporal HIErarchical Forecasting. Hyndsight blog. 22 August 2016. Disponível em: https://robjhyndman.com/hyndsight/thief/. Acesso em: 20.04.2023.

MORETTIN, Pedro A.; TOLOI, Clélia M.C. Análise de Séries Temporais. São paulo: Editora Edgard Blücher, 2004. https://www.ime.usp.br/~pam/st/