Pacote: universo tidyverts

Eu nunca tinha ouvido falar desse universo de pacotes que o tidyverts se propõe trazer para Análise de Séries Temporais, mas provavelmente isso seja por eu estar imersa na área da saúde ou também por quase ninguém do meu círculo não conhecer ou escrever sobre esse universo de pacotes.

E você, conhecia esse universo?

O nome do pacote é bem sugestivo, mas não vou deixar de dizer que ele foi desenhado para compilar lindamente com o universo tidyverse.

Maravilhoso, né?

E como aprender nunca é demais e já que fiquei sabendo desse universo, vamos hablar um pouco sobre ele.

Para compreender as funções desse pacote, primeiro precisamos entender um pouco mais sobre Análise de Séries Temporais.

O conceito de Séries Temporais é chatinho, mas com os nossos colegas Morretin e Toloi (2011) fica um pouquinho mais difícil (hehehe). Eles dizem que os modelos utilizados para descrever séries temporais são processos estocásticos, isto é, processos controlados por leis probabilísticas, ou seja, uma série temporal é um processo estocástico Z = Z(t) : t ∈ T, onde T é um conjunto de índices que representa o tempo (cada dia, cada mês, cada ano).

Então, é possível que dois pesquisadores obtenham diferentes modelos para uma mesma série?

SIM! Os modelos dependem de vários fatores, incluindo estacionariedade, tendência e sazonalidade e isso resulta em modelos diferentes para a mesma série.

Mas o que é tendência e sazonalidade?

É, realmente não tem como falar de série temporal sem falar desses dois queridinhos. Então vamos lá...

A tendência e a sazonalidade são componentes não observáveis e que fazem parte da decomposição de uma série temporal e, outra componente importante e que não mencionei ali encima, é a componente aleatória de média zero e variância constante. Mas, comumente falando, o movimento sistemático de aumento ou decréscimo, linear ou não linear, de uma série temporal ao longo do tempo é chamado de tendência. A tendência representa o efeito de longo prazo ao redor da média, ou seja, o movimento dominante em uma série temporal e, em relação a sazonalidade, entende-se por sazonalidade, por sua vez, flutuações periódicas recorrentes que podem aparecer em períodos específicos do ano, quando as observações são intra anuais, ou seja, registradas mensalmente, trimestralmente ou semanalmente, por exemplo.

E o que isso tem a ver com a decomposição de uma série temporal?

TUDO! Uma série temporal é descrita conforme a equação abaixo, em que o Z de t seja as observações de uma série temporal no tempo t igual a um, dois, três até N e que essa é composta pela Tendência, Sazonalidade e a como ruído branco.

Pensa assim...

Se hoje (t) está calor a tendência (T) é que amanhã (t+1) também esteja fazendo calor, talvez um pouco menos ou um pouco mais e, essa dúvida, pode ser também devido a estação do ano (S) em que nós estamos. É esperado que nessa época seja um período mais quente, com chuvas isoladas e, dependendo do estado do Brasil, o mesmo período é de chuvas fortes e temperaturas baixas (S). E veja que os dados de séries temporais geralmente não são independentes, especialmente se os intervalos da amostra são curtos como coleta de temperatura diária. As observações próximas costumam ser mais parecidas que as mais distantes, mas pode haver algum erro (a).

Conseguiu sentir com esse exemplo o que é essa decomposição de uma série temporal por tendência e sazonalidade?

Já a estacionariedade faz-se respeito a uma série temporal que se desenvolve no tempo aleatoriamente ao redor de uma média constante, refletindo alguma forma de equilíbrio estável.

É importante ressaltar que para cada composição da série temporal, há um teste para validar a sua existência ou não, visto que é de interesse do Estatístico ou Analista modelá-los.

Mas aí vem uma dúvida:

E se descobrimos que uma série temporal é não estacionária, há uma forma de podermos torná-la estacionária?

Tem sim! E para você ir compreendendo essa dinâmica, abre o RStudio aí. Não tem instalado na sua máquina? Vem aqui que tenho uma forma de resolver o seu BO.

Abriu? Letisbora.

Vamos utilizar o pacote tsibbledata que faz parte do universo tidyverts e tem o intuito de trazer bancos de dados com características de diversos tipos de séries temporais. Vou escolher o banco de dados global_economy que contém indicadores macroeconômicos anuais para 263 países para dar uma graça na nossa brincadeira.

Os dados são provenientes do Banco Mundial e este é composto por variáveis que demostram o PIB, Crescimento percentual anual do PIB, IPC, Importações, Exportações e População, cuja disponibilidade varia entre os países.

Já dá um CTRL+C CTRL_V do código abaixo para instalar todos os pacotes que vamos utilizar. Não vai esquecer de carregá-los.

install.packages("fable")
install.packages("fabletools")
install.packages("tsibble")
install.packages("tsibbledata")
install.packages("feasts")
install.packages("dplyr")

Este banco de dados foi coletado entre os anos de 1960 a 2017, totalizando 15150 observações em que para cada ano temos 260, 262 ou 263 medidas de acordo com o país. Como este post é voltado a fins didáticos, vamos filtrar a base de dados para o nosso Brasilsão.

dados_exemplo <- tsibbledata::global_economy %>% 
  dplyr::filter(Country == "Brazil") 

Também vamos utilizar a variável Growth que é o crescimento percentual anual do PIB do país. Quando você abrir o banco de dados e inspecionar a variável, você verá que no ano de 1960 o Brasil não possui observação, então vamos filtrar a base para retirar esse NA. Poderíamos tratar esse NA? Poderíamos e existem vários métodos para tal, mas esse não é o nosso objetivo aqui.

dados_exemplo <- tsibbledata::global_economy %>% 
  dplyr::filter(Country == "Brazil") %>% 
  dplyr::filter(!(is.na(Growth)))	

Com a base pronta...

Preciso contar um segredinho pra você:

Se você pegar qualquer banco de dados que seja uma série temporal devido as características dos dados, jogar no RStudio e não dizer que ele é uma série temporal, vai dar erro quando você utilizar qualquer função voltada a séries temporais.

Ficou chocade? Ninguém te falou isso?

Pois é, isso acontece (e te juro que nas melhores famílias também) e uma forma (antiga) de fazer essa declaração é utilizando a função ts do pacote stats. Vou te explicar como é a partir do exemplo abaixo.

stats::ts(1:10, frequency = 4, start = c(1959, 3))

O resultado do código acima é esse:

     Qtr1 Qtr2 Qtr3 Qtr4
1959              1    2
1960    3    4    5    6
1961    7    8    9   10

O parâmetro frequency dita o número de observações por unidade de tempo, ou seja, no ano de 1959 vamos ter 4 observações em que a partir do ano 1959 essa unidade de tempo vai iniciar na terceira casa (parâmetro start). Agora vai brincando aí para entender essa dinâmica alterando os parâmetros.

E você deve estar se perguntando indignade aí...

E que diferença faz Fê?

Essa diferença aqui:

freq_ts <- ts(1:10, frequency = 4, start = c(1959, 2))
class(freq_ts)
[1] "ts"

A classe do objeto passa a ser ts e isso faz TODA diferença. Mas, como eu disse, essa é a forma antiga de se declarar uma ts e é claro que há uma forma mais atual para ditar quando um dado é uma série temporal a partir de um dos pacotes do universo tidyverts.

A função as_tsibble é do pacote tsibble e faz brilhantemente o papel da função ts de uma forma até mais intuitiva. Veja seu preâmbulo abaixo.

tsibble::as_tsibble(
  x,
  key = NULL,
  index,
  regular = TRUE,
  validate = TRUE,
  .drop = TRUE,
  ...
)

Vamos de exemplo desse queridinho e observe que já estamos utilizando um pacote do universo tidyverse para a confecção do banco de dados.

tbl1 <- tibble::tibble(
  date = as.Date("2017-01-01") + 0:9,
  value = rnorm(10))

O resultado dele é:

> tbl1
# A tibble: 10 × 2
   date          value
   <date>        <dbl> 
   1 2017-01-01 -0.777 
   2 2017-01-02  0.766  
   3 2017-01-03  1.75   
   4 2017-01-04 -0.0665 
   5 2017-01-05 -0.313 
   6 2017-01-06 -0.397 
   7 2017-01-07 -0.235 
   8 2017-01-08 -0.0879 
   9 2017-01-09 -0.00182
  10 2017-01-10 -0.784 

Observou que ele é um tibble? Agora solicita a classe desse objeto tbl1 e guarda aí na sua mente. Abaixo vamos utilizar a transformação do pacote tsibble.

tbl1 <- tsibble::as_tsibble(tbl1, index = date)

A classe de tbl1 agora é:

> class(tbl1)
[1] "tbl_ts"     "tbl_df"     "tbl"        "data.frame"

Aí foi um detalhizinho que é importante e não podia passar despercebido por esse post, mas agora...

Bora brincar com essa série temporal com visualizações, testes e insigths?

A primeira missão que nós temos é a visualização da série e obter insigths de tendência, sazonalidade e estacionariedade. Lá na minha graduação nós utilizávamos o pacote stats, mas aqui vamos utilizar o pacote fabletools que fornece ferramentas para construção de pacotes de modelagem, com foco na previsão de séries temporais.

dados_exemplo %>% 
  fabletools::autoplot(Growth) +
  ggplot2::labs(title = "Crescimento do PIB - Brasil",
   y = "%", x = "Anos") +
  ggplot2::scale_x_continuous(breaks = c(1960, 1970, 1980,
                                         1990, 2000, 2010,2017), 
                              limits = c(1960, 2017))

E o produto do código acima é esse lindo gráfico com uma linha indo pra baixo e pra cima ao longo dos anos.

O que você vê nessa série temporal? Será que tem tendência? Ou alguma sazonalidade? Ou uma história por trás desses picos?

A história dos altos e baixos durante todos esses anos a gente tem certeza que tem, mas tendência, sazonalidade e estacionariedade vamos ter que testar.

Para testar a tendência vamos utilizar o teste de sequências de Wald-Wolfowitz (H0: Não há tendência na série), para sazonalidade o teste não paramétrico de Kruskal-Wallis (H0: Não há sazonalidade na série) e teste de Dickey-Fuller (H0: A série não é estacionária) Aumentado para testar a estacionariedade.

Os testes e o respectivo resultado para a série temporal se encontram abaixo.

• Estacionariedade

tseries::adf.test(dados_exemplo$Growth)
	Augmented Dickey-Fuller Test

data:  dados_exemplo$Growth
Dickey-Fuller = -3.135, Lag order = 3, p-value = 0.1167
alternative hypothesis: stationary

• Sazonalidade

stats::kruskal.test(dados_exemplo$Growth~dados_exemplo$Year)
	Kruskal-Wallis rank sum test

data:  dados_exemplo$Growth by dados_exemplo$Year
Kruskal-Wallis chi-squared = 56, df = 56, p-value = 0.4749

• Tendência

randtests::runs.test(dados_exemplo$Growth)
	Runs Test

data:  dados_exemplo$Growth
statistic = -2.967, runs = 18, n1 = 28, n2 = 28, n = 56, p-value = 0.003008
alternative hypothesis: nonrandomness

Logo, de acordo com os testes acima e assumindo o nível de 5% significância, a série temporal não é estacionária, tem tendência e não possui sazonalidade.

Em um nível informal, a estacionariedade fraca pode ser testada pelo correlograma de uma série temporal, que é um gráfico de autocorrelação em várias defasagens. Para a série temporal estacionária, o correlograma enfraquece rapidamente, enquanto para a série temporal não estacionária, ele enfraquece gradualmente. Para uma série puramente aleatória, as autocorrelações em todas as defasagens, 1 e superiores, são zero.

Então bora ser um pouquinho informal?

Essa outra forma citada acima e que podemos ter insigths, é através da estrutura da Função de Autocorrelação (ACF) da série. Nesse caso nós vamos utilizar o pacote feasts que fornece uma coleção de ferramentas para a análise de dados de séries temporais. O nome do pacote é um acrônimo que compreende seus principais recursos: Extração de recursos e estatísticas para séries temporais.

Mais um CTRL+C CTRL_V aí e letsbora...

dados_exemplo %>% 
 feasts::ACF(Growth, lag_max = 10)

Como resultado vamos ter o lag de 10 anos atrás.

# A tsibble: 10 x 3 [1Y]# Key:       Country [1]
   Country      lag     acf
   <fct>   <cf_lag>   <dbl> 
 1 Brazil        1Y  0.518 
 2 Brazil        2Y  0.368 
 3 Brazil        3Y  0.226 
 4 Brazil        4Y  0.163 
 5 Brazil        5Y  0.0656 
 6 Brazil        6Y  0.138 
 7 Brazil        7Y  0.168 
 8 Brazil        8Y  0.0600 
 9 Brazil        9Y  0.0504
 10 Brazil       10Y -0.0451

Veja que do 1Y ao 4Y temos correlações a princípio baixas e vão decaindo ao longo do tempo até chegar o momento do 5Y e ser extremamente baixo (0.0656). Um ponto a ser observado é que no 10Y a ACF é negativa (Dica: Faz para 20 anos e observe o comportamento). E lembre-se: a correlação se situa entre -1 a +1, como qualquer correlação coeficiente faz.

Vamos fazer um correlograma?

Vamos lá...

dados_exemplo %>% 
  feasts::ACF(Growth, lag_max = 20) %>% 
  feasts::autoplot()

O produto do código acima é esse lindo gráfico ACF.

Esse gráfico demostra através de seu caimento (um tanto brusco) entre as lags que essa série não é estacionária, tem tendência e, aparentemente, não possui sazonalidade.

Agora lascou...

Nem tudo está lascado, tá? Existe uma coisinha comum que costumamos fazer ou aplicar na série quando a série "não passa" em algum dos testes. Essa coisinha se chama tomar diferenças da série, mas é importante ressaltar que há outras formas de se tratar esse problema e você deve dar um Google aí e dedicar um tempo de estudo para essas outras opções.

E agora você deve estar se perguntando:

O que é isso? Diferenciação? Como assim Fê?

Eu te entendo. Então vamos lá...

Se uma ou mais séries temporais são não estacionárias, temos de transformar as séries temporais não estacionárias para torná-las estacionárias. Uma propriedade importante dos processos estacionários é que os processos obtidos através das combinações lineares dos processos estacionários são também estacionários, ou seja, se Z de t é estacionário então o processo W de t = Z de t - Z de t-1 é também estacionário. Quando aplicamos o método da diferenciação, nós costumamos dizer que a série é estacionária em diferenças.

A forma que fazemos essa diferença no R está abaixo.

dados_exemploDIFF <- base::diff(dados_exemplo$Growth)

E quando você rodar...

Quantas observações você tem após a diferenciação?

Neste caso nós passamos a ter 56 observações da variável Growth. E aqui eu ficava com muitas dúvidas e pode ser que você também e, por isso, vou explicar de outra forma essa diferenciação.

Pensa assim...

Fomos ao supermercado comprar azeite durante o período de um ano. Se o preço de 1 mês atrás estiver a R$ 20,00 e variou para R$ 30,00, o valor da diferenciação será de R$ 10,00 e se o valor de hoje estiver a R$ 40,00 então a diferença será de -R$ 10,00.

Conseguiu pegar a ideia?

Eu espero que sim, mas se não conseguiu, escreve aí nos comentários a sua dúvida que eu ou comunidade te dá um help.

Outro ponto importante é falar sobre o impacto dessa diferenciação. Em geral, se uma série temporal (não estacionária) precisa ser diferenciada d vezes para tornar-se estacionária, denominamos essa série temporal integrada de ordem d. Se uma série é estacionária, ela será integrada de ordem 0. Veja que os termos “série temporal estacionária” e “série temporal integrada de ordem zero” significam a mesma coisa. E a cada diferenciação nós teremos uma série integrada de ordem 1, 2, 3 e quantas forem necessárias, mas normalmente a primeira diferenciação já resolve o nosso problema. Essa integração é denominada por "I(número de diferenciações)" e você provavelmente já viu muitos ARIMA por aí, né? O I vem dessa integração.

Mas Fê, e se a série não tiver tendência e sazonalidade?

Daí nós vamos falar do famigerado Ruído Branco. O Ruído Branco surge quando os dados não estão correlacionados ao longo do tempo, em que estes possuem média zero e variância constante, ou seja, independência e, consequentemente, sem padrão.

E você faz o quê com isso, né?

Faz nada. Ou até faz, mas vai ser um caminho dolorido e longo. Os modelos mais conhecidos como AR, MA e ARIMA não performariam muito bem nesse caso e você teria que dar uma estudada a mais para modelar uma série temporal aleatória.

Voltando ao exemplo do R...

Após toda diferenciação você deve avaliar através dos testes a estacionariedade, tendência e sazonalidade novamente. Nesse caso vou deixar como exercício fazer esses testes após a aplicação do método de diferenciação e, se você vir a fazer, comenta aí os seus resultados.

E agora chegou aquele momento...

Faço o quê agora?

Bem, agora podemos pensar que a qualidade de um experimento ou análise pode ser avaliada pela magnitude do erro experimental e pelo atendimento das pressuposições do modelo estatístico que está sendo utilizado. Na Análise de Séries Temporais, os erros devem ser independentes e essa independência dos erros existe quando os erros controlados de um fator de variação não interferem com os erros controlados de outros fatores.

O que isso significa?

Significa que a probabilidade de erro de uma observação ter certo valor, não depende dos valores dos erros de outras observações. E é importante apontar que essa abordagem é clássica e, se você tiver interesse, a abordagem bayesiana em séries temporais é útil e trata as preposições de uma outra forma.

E, na minha opinião, esse post já está ENORME e iniciar pressupostos, qualidade do modelo estimado e previsão utilizando séries temporais não seria o ideal.

Vamos combinar (jurado de dedinho) de você voltar aqui semana que vem para ler o outro post sobre o universo tidyverts?

Espero que você tenha gostado e que passe a utilizar e estudar mais sobre universo de séries temporais. Até mais!

Fernanda Kelly R. Silva | Estatística