rm(list=ls()) # Remoção de arquivos no espaço de trabalho.
library(tidyverse) # Data science.
library(forecast) # Previsão de séries temporais no R.
library(sidrar) # Raspagem de dados do SIDRA-IBGE.
library(zoo) # Tratamento de séries temporais no R.
library(quantmod) # Dados financeiros em R.LDE-Séries Temporais
Análise de Séries Temporais com ARIMA no R
Primeiro, vamos acionar os pacotes que serão utilizados neste estudo. Vejam a utilidade de cada pacote pesquisando na internet. Para abrir chunks: ctrl+Alt+I.
Séries temporais são dados organizados de forma sequencial no tempo. São dados muito comuns em Economia.
Vamos pegar, como exemplo, dados sobre a Pesquisa Industrial Mensal - Produção Física (PIM-PF) do IBGE, Tabela 8888, no Sidra.
pim <- get_sidra(api="/t/8888/n1/all/v/11601/p/all/c544/129314/d/v11601%201")Agora vamos limpar nossa série. Você pode usar uma LLM para explicar os detalhes do código.
pim_pf <- pim %>% # Vamos trabalhar com pipes, para canalizar o que queremos do df salvo.
select(Mes = `Mês (Código)`, # Usamos essa função para selecionar as colunas de nosso interesse.
Valor = Valor) %>%
mutate(Mes = as.numeric(Mes), # Estamos transformando as colunas selecionadas.
Valor = as.numeric(Valor)) %>% # Queremos que elas represente números.
arrange(Mes) # Essa função é usada para ordenar os dados.A série começa em janeiro de 2002. Vamos montar agora o objeto série temporal (ts) propriamente dito.
pim_ts <- ts(pim_pf$Valor,
start = c(2002, 01), # Queremos que a série começe nessa data.
frequency = 12) # Frequência, no caso, 12 meses em um ano.
pim_ts Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2002 0.0 2.2 -0.3 0.2 0.0 0.3 0.4 -0.9 0.4 0.4 -0.1 -0.3
2003 0.4 -0.4 -0.2 -0.2 -1.5 0.6 -0.9 1.4 2.9 0.5 1.0 -1.5
2004 2.0 1.5 0.6 0.8 0.9 0.4 0.7 1.1 0.4 -0.6 0.0 0.6
2005 0.5 -1.4 1.0 0.0 3.5 -2.1 -0.6 0.1 -0.7 -0.5 0.9 1.7
2006 0.4 0.8 -0.6 0.5 0.7 -1.2 0.8 -0.1 -1.2 0.9 1.3 0.6
2007 -0.2 1.3 1.3 1.0 0.2 0.6 -0.3 1.4 -1.1 2.8 -0.9 0.3
2008 1.7 -0.3 0.2 0.1 0.8 1.3 0.5 -1.1 1.1 -2.7 -6.1 -11.1
2009 2.1 2.3 0.5 1.3 1.2 1.4 1.3 1.6 1.4 2.2 1.3 0.4
2010 1.5 -0.2 2.3 0.4 -0.3 -0.4 -0.6 -0.2 0.2 -0.3 0.7 0.5
2011 0.3 1.2 0.3 -2.6 3.1 -1.8 -0.2 -0.8 -1.6 -1.1 0.5 2.6
2012 -4.2 0.8 -0.4 -0.1 0.5 0.2 1.4 2.2 -2.0 1.5 -1.3 0.6
2013 1.0 -2.9 2.0 1.6 -0.2 0.4 -0.6 0.3 0.5 -1.1 0.5 -3.5
2014 1.8 0.4 -0.4 -0.8 -1.0 -2.8 2.5 0.4 0.2 0.1 -0.9 -2.8
2015 0.5 -0.9 -0.4 -1.9 -1.1 -0.8 -1.8 0.4 -2.3 -0.1 -1.7 -1.5
2016 0.9 -1.1 0.6 0.0 0.4 -1.4 0.7 -1.9 1.1 -1.4 1.0 1.7
2017 0.7 1.2 -1.6 0.0 0.4 1.0 0.2 0.1 0.3 1.0 0.4 2.7
2018 -2.3 -0.6 1.4 0.4 -10.9 12.4 -0.2 -0.9 -2.9 1.4 -0.1 -0.2
2019 -0.5 0.3 -0.6 0.5 0.0 -0.9 -0.5 1.0 0.0 2.3 -2.4 -1.4
2020 1.4 0.6 -7.7 -19.8 7.8 10.0 9.2 2.9 2.7 1.5 0.4 -0.5
2021 0.7 -1.5 -2.1 -1.8 1.1 -0.4 -1.3 -0.1 -0.3 -0.5 0.0 1.5
2022 -0.1 0.9 -0.8 0.6 -0.5 -0.1 1.2 -1.6 -1.5 1.5 0.0 0.0
2023 0.0 -0.2 0.8 -0.4 0.6 -0.4 -0.5 0.4 -0.4 0.4 1.1 1.5
2024 -1.2 0.2 0.5 -0.4 -1.2 4.4 -1.4 -0.3 1.0 0.0 -0.7 -0.4
2025 0.2 0.0 1.8 -0.7 -0.5 0.1 -0.1 0.7 -0.4 Gráfico da série com o ggplot2
Para desenhar o gráfico da série temporal usando o pacote ggplot precisamos da função autoplot. Mas para que isso funcione, precisamos do pacote forecast, desenvolvido para se trabalhar séries temporais.
autoplot(pim_ts, colour = "red", size = 1.2) + # Essa função no ggplot permite fazer gráficos de objetos como séries temporais.
ggtitle("Produção Física - 2002-2025") +
ylab("Variação Percentual") +
xlab("Mês/Ano") +
theme_minimal()
Analisando a série
Vamos começar decompondo a série em seus elementos, principais, a saber, tendência (trend), componente sazonal (seazonal) e ruído (random). Normalmente, a série temporal é a soma desses três componentes.
decomposicao <- decompose(pim_ts)
plot(decomposicao)
Vamos usar o pacote monthplot para visualizar a sazonalidade da série de maneira mais nítida. Esse pacote calcula a médias de cada mês.
monthplot(pim_ts, col = "darkgreen", lwd = 2,
main = "Sazonalidade da PIM-PF", ylab = "Variação percentual")
Pelo gráfico, notamos que a menor média é em abril de cada ano e a maior em julho.
Previsão (modelo ARIMA)
ARIMA é a sigla para:
AutoRegressive Integrated Moving Average, ou seja, modelo autorregresivo, integrado e com média móvel.
É um modelo usado para fazer previsões sobre valores futuros de uma série. O modelo é composto de três fatores:
O fator autorregressivo, AR(p).
O fator de integração, I(d).
O fator da média móvel, MA(q).
Os parâmetros p, d & q representam as dimensões de cada parte do modelo ARIMA. Vamos procurar entender o que são essas partes.
Componente autorregressivo, AR(p)
A autorregressão quer dizer que procuramos entender como uma série se comporta a partir de seus valores passados. Em símbolos,
\(y_t = f(y_{t-1})\)
Se a série é explicada, com duas defasagens, ou seja, \(p = 2\), o componente autorregressivo do modelo pode ser escrito da seguinte maneira:
\(y_t = \phi_1y_{t-1} + \phi_2y_{t-2} + \epsilon_t\)
Onde, c é um parâmetro constante que dá conta dos fatores que afetam a série e não tem nada a ver com seus componentes autorregressivos. Seria o intercepto de uma equação linear.
Os parâmetros \(\phi\) revelam como os componentes defasados da série a afetam. Por exemplo, digamos que temos os seguintes valores para cada item da série:
\(\phi_1 = 0,7\)
\(\phi_2 = 0,2\)
\(y_{t-1} = 4\)
\(y_{t-2} = 6\)
\(\epsilon_t = -0.17\)
Então,
\(y_t = 3,83\)
Que tal nós simularmos no R um série com padrão autorregressivo de \(p = 1\)? Na simulação vamos usar a função arima.sim do pacote stats, do R-base.
# Simulação de um processo AR(1)
set.seed(123) # Valores aleatórios, mas fixos.
phi <- 0.8 # Esamos determinando o coeficiente autorregressivo.
n <- 100 # Tamanho da série
y <- arima.sim(model = list(ar = phi), n = n)
# Plotar a série simulada
plot(y, type = "l", col = "blue", lwd = 2,
main = paste("Série AR(1) Simulada (φ =", phi, ")"),
xlab = "Tempo", ylab = expression(y[t]))
grid()
Esclarecendo alguns pontos do código acima. Usamos a função arima.sim para simular um modelo autorregressivo, com lag = 1. Como argumento da função arima.sim a função model pede que um conjunto de argumentos (uma lista, na verdade, daí o uso da função list). E que argumentos são esses em nosso exemplo? O termo de defasagem e o número de observações. A função cria, então o conjunto de dados, gerados aleatoriamente, com um erro normalmente distribuído, ou seja,
\(\epsilon_t \sim \mathcal{N}(0,1)\)
Ou seja, o erro segue um padrão de ruído branco.
A rigor, portanto, cada ponto da linha azul foi gerado levando em conta a seguinte equação:
\(y_t = 0,8y_{t-1} + \epsilon_{t}\)
É importante esclarecer que nós não geramos um modelo Arima de séries temporais. Geramos, artificialmente, um série temporal: valores ao longo do tempo. Esses valores são unidos pela linha azul, o que dá origem ao gráfico acima.
Quando queremos estimar os parâmetros \(\phi\) do modelo AR que representa a série temporal, os pacotes utilizam o método da máxima verossimilhança: o software encontra os parâmetros que maximizam a verossimilhança, ou seja, que tornam os dados observados mais prováveis.
Vamos visualizar as consequências de escolhas diferentes para \(\phi\). Para isso, utilizamos o seguinte código:
set.seed(123)
n <- 100
# Três séries com diferentes níveis de autorregressão
y1 <- arima.sim(model = list(ar = 0.2), n = n)
y2 <- arima.sim(model = list(ar = 0.8), n = n)
y3 <- arima.sim(model = list(ar = -0.8), n = n)
#png("graficoAR.png", width = 1000, height = 800) # Para exportar o gráfico em PNG.
par(mfrow = c(3,1))
plot(y1, type = "l", col = "darkgreen", lwd = 2,
main = "AR(1) com φ = 0.2 (fraca dependência)",
xlab = "Tempo", ylab = expression(y[t]))
grid()
plot(y2, type = "l", col = "blue", lwd = 2,
main = "AR(1) com φ = 0.8 (forte dependência positiva)",
xlab = "Tempo", ylab = expression(y[t]))
grid()
plot(y3, type = "l", col = "red", lwd = 2,
main = "AR(1) com φ = -0.8 (forte dependência negativa)",
xlab = "Tempo", ylab = expression(y[t]))
grid()
Analisando esses resultados:
- \(\phi = 0,2\) : a influência de valores passados é pequena. A série é quase um ruído branco.
- \(\phi = 0,8\) : a série fica menos errática, se move com mais lentidão. Há memória persistente na sua composição.
- \(\phi = -0,8\) : valores positivos e negativos se alterna.
A intuição por trás do modelo AR: As séries econômicas geralmente são do tipo AR, ou seja, os resultados ou números de hoje são reflexo dos de ontem. É como se em Economia, a Lei de Brandford é bem presente.
Componente de média móvel, MA(q)
O termo média móvel tem dois significados distintos em estátística das séries temporais ou econometria.
O primeiro significado é o mais comum: média móvel é um processo de suavização de uma série temporal. Por exemplo, durante a pandemia do Covid-19, era sempre informada a média móvel de casos registrados. Média móvel aqui é a média de \(k\) observações da série. Então, cada \(k\) observações eram substituídas por por sua média, recriando uma série suavizada e com menos valores. O outro sentido é o de uma estimativa dos erros usada em modelos do tipo ARIMA.
Vamos ilustrar a média móvel como a suavização de uma série usando dados sobre a Covid-19 no Brasil. No site https://covid.saude.gov.br/, clique em Painel Interativo. Em seguida, clique em Baixar Dados de Casos novos por semana. Infelizmente o arquivo vem em formato proprietário. Eu abri o arquivo, exclui a coluna de datas, nomeie a coluna de valores como “casos” e o salvei com o nome “covid.csv”.
Vamos carregar esses dados:
covid <-read.csv("covid.csv", header = T)Vamos acrescentar uma coluna de datas semanais ao nosso data frame.
covid$semana <- seq(
as.Date("2025-01-01"),
by = "week",
length.out = nrow(covid)
)Vamos acrescentar a média móvel de 3 semanas ao nosso data frame.
covid_mm <- covid %>%
mutate(media_movel = rollmean(casos, k = 3, fill = NA, align = "right"))Explicando o código acima: Usamos um dos mais antigos e importantes pacotes do R que trabalha com séries temporais, o pacote zoo. A função é a rollmean que produz médias móveis. Escolhemos fazer médias de 3 semanas. Quando não houver valor, a informação fica com NA. E de fato, como precisamos de 3 observações para se calcular a média móvel, as duas primeiras ficam com NA. Calculamos uma média e passamos para a seguinte, avançando uma semana (daí o termo móvel). O alinhamento à direita significa que estamos representando os dados da última semana sendo considerada.
Vamos fazer o gráfico da série original com sua média móvel.
ggplot(covid_mm, aes(x = semana)) +
geom_line(aes(y = casos), color = "darkgreen", linewidth = 1.2) +
geom_line(aes(y = media_movel), color = "blue", linewidth = 1.2) +
labs(
title = "Série com Média Móvel",
x = "Data",
y = "Valor"
) +
theme_minimal()Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_line()`).
Notem o aviso de que dois valores foram desconsiderados.
O outro tipo de erro é o estocástico e será tratado no contexto de um modelo de séries temporais ARMA e ARIMA.
Modelos ARMA, prelúdio para os modelos ARIMA
Modelos ARMA são usados para representar séries temporais consideradas estacionárias (mais sobre isso depois). São dois os componentes de um modelo ARMA:
Um componente autorregressivo (AR) — que, como vimos, usa os valores passados da própria série.
Um componente de média móvel (MA) — que usa os erros (ou choques) passados do modelo.
A estrutura do modelo é a seguinte:
\(y_t = \mu + \phi_1y_{t-1} + ... + \phi_py_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + ... + \theta_q\epsilon_{t-q}\)
Onde:
\(y_t\) é o valor da série no momento t.
\(\mu\) é o valor médio da série.
\(\phi_1\) é o parâmetro estimado relacionado à observação da série no tempo t-1.
\(\epsilon_t\) o erro da série no momento t.
\(\theta_1\) o parâmetro associado ao erro defasado de um período.
O componente MA(1)
O fator de média móvel se refere à importância de estimativas dos erros para o valor de uma série. O valor atual é resultado de uma combinação do erro atual e do erro passado:
\(y_t = \mu + \epsilon_t + \theta_1\epsilon_{t-1}\)
Portanto, o termo média móvel quer dizer aqui média de erros e não média de observações, como no caso anterior.
No caso de um ARMA(1,1) estamos considerando que a série pode ser representada pelo seguinte modelo:
\(y_t = \mu + \phi_1y_{t-1} + \epsilon_t + \theta_1\epsilon_{t-1}\)
Intuitivamente, estamos afirmando que um bom modelo a representar uma série temporal é o leva em conta o passado recente da série e um componente de correção do erros de previsão cometidos no passado também recente.
Se \(\theta_1 > 0\) erro de ontem “reforça” o erro de hoje. Se \(\theta < 0\), os erros de ontem amenizam a ocorrência do erro de hoje. É como aprender com seus próprios erros. Vamos estudar esse assunto com ajuda de simulações feitas no R, mais uma vez.
Intuição por trás de um modelo com coeficiente de MA negativo (\(\theta<0\)): aqui está a crença de que os agentes econômicos aprendem com seus erros. Ou seja, suas expectativas estão sendo adaptadas conforme erros de previsão aconteçam. É um caso de expectativas adaptativas.
Simulando modelos de média móvel no R.
set.seed(123)
y1 <- arima.sim(model = list(ma = 0.8), n = 100)
y2 <- arima.sim(model = list(ma = -0.8), n = 100)
#png("graficoMA.png", width = 1000, height = 800) # Para exportar o gráfico em PNG.
par(mfrow = c(2,1))
plot(y1, type = "l", col = "blue", main = "MA(1) com θ = 0.8")
plot(y2, type = "l", col = "red", main = "MA(1) com θ = -0.8")
Os dois gráficos mostram o seguinte:
\(\theta = 0,8\) : o erro passado reforça o erro atual.
\(\theta = -0,8\) : o erro passado corrige o erro atual.
Para estimar os parâmetros \(\theta\) com base em um conjunto de dados, usamos o método da máxima verossimilhança mais uma vez.
Condições para aplicação do modelo ARMA:
A série precisa ser estacionária (mais sobre isso depois).
O erro é do tipo erro branco.
Procedimentos estatísticos são usados para definir as ordens dos componentes AR(p) e MA(q).
Simulando um modelo ARMA
Exemplo no R. Acompanhe o seguinte código, usando, mais uma vez, a função arima.sim. Agora temos dois argumentos, dando mais sentido à lista no argumento da função.
set.seed(123)
# parâmetros
phi <- 0.6 # AR(1)
theta <- 0.8 # MA(1)
n <- 100 # tamanho da série
# simulação
y <- arima.sim(model = list(ar = phi, ma = theta), n = n)
# gráfico
plot(y, type = "l", col = "purple", lwd = 2,
main = paste("Série ARMA(1,1): φ =", phi, ", θ =", theta),
xlab = "Tempo", ylab = expression(y[t]))
grid()
Vamos comparar os três modelos que estudamos até aqui: AR, MA e ARMA.
set.seed(123)
# AR(1)
y_ar <- arima.sim(model = list(ar = 0.6), n = 100)
# MA(1)
y_ma <- arima.sim(model = list(ma = 0.8), n = 100)
# ARMA(1,1)
y_arma <- arima.sim(model = list(ar = 0.6, ma = 0.8), n = 100)
#png("graficoARMA.png", width = 1000, height = 800) # Para exportar o gráfico em PNG.
par(mfrow = c(3,1))
plot(y_ar, type = "l", col = "blue", main = "AR(1) φ=0.6")
plot(y_ma, type = "l", col = "red", main = "MA(1) θ=0.8")
plot(y_arma, type = "l", col = "purple", main = "ARMA(1,1) φ=0.6 θ=0.8")
Resultados: no primeiro gráfico vê-se uma forte persistência. No segundo gráfico note a elevada flutuação de valores. E no terceiro gráfico temos a combinação dos dois valores.
Estimando um modelo ARMA(1,1)
Digamos que tenhamos escolhidos os parâmetros \(p = 1, q = 1\) para nosso modelo ARMA, ou seja, acreditamos (por razões estatísticas) que nossa série é razoavelmente estacionária e que informações recentes sobre os valores da variável e sobre erros (choques) a explicam bem.
O que significa, então, estimar um modelo ARMA?
O modelo a ser estimado é o seguinte:
\(y_t = \mu + \phi_1y_{t-1} + \epsilon_t + \theta_1\epsilon_{t-1}\)
Precisamos estimar os parâmetros: \(\mu\), \(\phi_1\) e \(\theta_1\).
Uma dificuldade prática com essa estimação é que, ao contrário do que acontece com \(y_t\), não conhecemos previamente os valores do erro. Então não podemos simplesmente aplicar um método de regressão linear padrão.
Então a ideia é escolher os parâmetros que tornam o modelo o mais próximo possível do conjunto de dados originais. Ou seja, parâmetros que têm a maior chance de representar de forma satisfatória os dados disponíveis. Esse é sentido de um modelo de máxima verossimilhança.
Trecho técnico: o componente de erro
No estimação do modelo pelo método da máxima verossimilhança (ML), supomos que os erros são normalmente distribuídos.
\(\epsilon_t \sim \mathcal{N}(0,1)\)
Portanto, a função densidade de probabilidade (FDP) dos erros é uma típica função normal. Uma FDP estima a probabilidade de ocorrência de eventos para o caso de variáveis contínuas.
\(f(\epsilon_t) = \frac{1}{\sqrt{2\pi\sigma^2}} exp^{(-\frac{\epsilon^2_i}{2\sigma^2})}\)
A fórmula acima expressa a probabilidade de ocorrência de um dado erro (\(\epsilon_t\)). Levando em conta todos os erros previstos, a probabilidade total será dada pelo produto dessas probabilidades individuais, considerando que os erros são independentemente distribuídos.
\(L(μ,ϕ,θ,σ^2)= \prod_{t=1}^{T}(2\pi\sigma^2)^{-1/2}exp(-\frac{\epsilon^2_t}{2\sigma^2})\)
Cada \(\epsilon_t\) é resultado da diferença entre o valor verdadeiro da variável e o valor estimado pelo modelo (por isso, os parâmetros estão implícitos na função \(L(.)\) acima.
Maximizar a verossimilhança consiste em determinar os valores dos parâmetros (\(\mu\), \(\phi\), \(\theta\) e \(\sigma^2\)) que maximizam a função \(L(.)\).
Intuitivamente, estamos dizendo que a função \(L(.)\) expressa a probabilidade de ocorrência dos valores observados, dado os conjunto de parâmetros, ou seja:
\(L(\mu, \phi, \theta, \sigma^2) = P(y_1, ...y_T| \mu, \phi, \theta, \sigma^2)\)
Se o modelo está bem ajustados, os \(\epsilon_t\) serão pequenos. (Observe que na fórmula do \(L(.)\) o erro é um expoente com sinal negativo. Portanto, minimizar o erro é o mesmo que maximizar a função de máxima verossimilhança.
Abrindo \(\epsilon_t\)
\(\epsilon_t = \hat{y}_t - y_t\)
onde:
\(y_t = \mu + \phi_1y_{t-1} + \epsilon_t + \theta_1\epsilon_{t-1}\)
e
\(\hat{y}_t = \mu + \phi_1y_{t-1} + \theta_1\epsilon_{t-1}\)
Método interativo
A escolha dos parâmetros do modelo, ou seja, aqueles que maximizam a função \(L(.)\) e minimizam os \(\epsilon_t\) se dá por interação.
Tecnicamente falando, o algoritmo que, por exemplo, o R usará para definir os parâmetros do modelo é um método numérico interativo de otimização*
O processo incia com o algoritmo lançando um conjunto de valores inciais para os parâmeros, digamos:
\(\mu^{(0)}\), \(\phi^{(0)}\), \(\theta^{(0)}\), \(\sigma^{2(0)}\)
Os resíduos dessa primeira rodada são calculados:
\(\epsilon^{(0)}_t = y_t - \mu^{(0)} - \phi^{(0)}y_{t-1} - \theta^{(0)}\epsilon^{(0)}_{t-1}\)
A verossimilhança é calculada com esses resíduos (na verdade, por facilidade, o log da verossimilhança. Os parâmetros são levemente ajustados e a verossimilhança é calculada outra vez. E assim por diante. Até o ponto em que um novo valor de parâmetros não produza resultados significativos.
Em cada interação é calculado o gradiente da função no ponto. Esse gradiente aponta na direção para que os novos parâmetros sejam escolhidos de forma a aumentar o valor de \(L(.)\). Isso é feito até que não haja novos aumentos significativos, ou seja, se chegou à cominação ótima de parâmetros que otimiza \(L(.)\) e o gradiente, nesse ponto máximo, é aproximadamente igual a zero. Intuitivamente, é como se o algorítimo usasse o gradiente com bússola, ela sinaliza para qual direção a função deve crescer mais.
Componente de integração (d) ou diferenciação.
O que são séries não estacionárias? Tecnicamente, dizemos que uma série é estacionária quando atende aos seguintes pré-requisito:
A média da série é constante.
A variância da série é contante.
A autocovariância da série depende apenas do intervalo de tempo considerado e não de uma “memória” entre valores passados e presente.
Trecho técnico: autocovariância
A covariância entre duas variáveis é uma medida de quão relacionadas estão. Por exemplo, se uma variável \(X\) cresce no momento em que outra variável, \(Y\), podemos dizer que essas duas séries apresentam covariância positiva. A fórmula da covariância é a seguinte:
\(Cov(X, Y) = E[(X - E[Y])(Y - E[Y])]\)
Onde E(.) significa esperança matemática ou uma média de valores. Se os valores de X são maiores que sua média no momento em que os valores de Y são, também, maiores que sua média; os valores de X são menores que sua média exatamente quando os valores de Y são, também, menores que sua média, o resultado da multiplicação será positivo. Isso evidencia que os dados variam na mesma direção.
No caso da autocovariância, comparamos não os dados de duas séries, mas dados de uma mesma série em momentos diferentes no tempo. Considere uma série \(y_t\). Vamos então apresentar a fórmula da covariância para um intervalo (lag) de tempo igual a \(k\):
\(\gamma_k = Cov(y_t, y_{t-k}) = E[(y_t - \mu)(y_{t-k} - \mu)]\)
Em que \(\mu\) representa o valor médio da série.
Para entender melhor, digamos que a série em questão é o PIB per capita (PIBpc), de 2001 até 2024. Se \(k=3\), então estamos comparando os seguintes valores: PIBpc de 2024 com o de 2021; PIBpc de 2023 com o de 2020 e assim por diante. Ou seja, criamos duas séries. uma de valores “atuais” e outra com valores defasado. Mas a média dessas “duas” séries é uma só, como dissemos, \(\mu\).
Então, para que uma série seja considerada estacionária, é suposto que a autocovariância depende mais das defasagens, por exemplo o valor, o valor do PIBpc de 2024 está relacionado (varia no mesmo sentido) com o PIBpc de 2021, assim como o de 2023 está relacionado (mesma covariância) com o de 2020 e assim por diante. Ou seja, a relação está determinada pelo tamanho da defasagem e a autocovariância é constante.
O que seria um indicativo de não estacionariedade? Se a autocovariância não se mantêm constante (a cada \(k\) períodos de defasagem). Ou seja, a autocovariância se desvincula do lag e passa a ser determinada pelo tempo, t, mesmo.
Vamos usar uma analogia com temperaturas para entender melhor esse ponto. Considere que temos uma série de temperaturas diárias e queremos calcular a autocovariância para um lag de 90 dias ou três meses.
Se entre os dias 21 e 22 de setembro de 2020 a temperatura subiu 3 graus e o mesmo aconteceu entre os dias 21 e 22 de junho dos mesmo ano, podemos considerar que essa série é estacionária. Agora se a temperatura aumentou 1 grau por dia em um mês de verão e caiu 1 grau por dia em um mês de inverno (o lag de três meses os separando), dizemos que a série é não estacionária.
White noise & Random walk
Para visualizar a questão, vamos comparar o gráfico de uma série estacionária (ruído branco) com uma série não estacionária (passeio aleatório).
Ruído branco:
set.seed(123)
x <- rnorm(100) # 100 valores aleatórios N(0,1)
# Visualização
plot(x, type="l", main="Série Estacionária: Ruído Branco", col="blue", lwd = 2)
Passeio aletatório:
# Série não estacionária: passeio aleatório
set.seed(123)
y <- cumsum(rnorm(100)) # soma cumulativa
# Visualização
plot(y, type="l", main="Série Não Estacionária: Passeio Aleatório", col="red", lwd = 2)
Questão: para você, qual gráfico representa melhor uma série com preços de uma ação?
Diferenciação
Considere a seguinte série de PIBpc a preços do ano anterior:
### Pegar a série.
pibpc <- get_sidra(api = "/t/6784/n1/all/v/9809/p/all")All others arguments are desconsidered when 'api' is informed### Limpar a série.
pib_pc <- pibpc %>% # Vamos trabalhar com pipes, para canalizar o que queremos do df salvo.
select(Ano = `Ano (Código)`, # Usamos essa função para selecionar as colunas de nosso interesse.
Valor = Valor) %>%
mutate(Ano = as.numeric(Ano), # Estamos transformando as colunas selecionadas.
Valor = as.numeric(Valor)) %>% # Queremos que elas represente números.
arrange(Ano) # Essa função é usada para ordenar os dados.Criando uma série temporal.
pib_pc_ts <- ts(pib_pc$Valor,
start = c(1996), # Queremos que a série comece nessa data.
frequency = 1) # Frequência, no caso, 12 meses em um ano.Gráfico:
plot(pib_pc_ts,
main = "PIB per capita Brasil (1996-2022)",
xlab = "Ano",
ylab = "PIB per capita",
col = "blue",
lwd = 2)
grid(col = "darkgreen", lty = "dotted") # grade leve
Claramente, a série não é estacionária, apresenta uma forte tendência de crescimento.
Vamos diferenciar a série para tentar torná-la mais estacionária. Em seguida, vamos ao seu gráfico.
dpib_pc_ts <- diff(pib_pc_ts)plot(dpib_pc_ts,
main = "PIB per capita Brasil - série diferenciada",
xlab = "Ano",
ylab = "PIB per capita",
col = "blue",
lwd = 2)
grid(col = "darkgreen", lty = "dotted") # grade leve
Vamos diferenciar duas vezes a série.
dpib_pc_ts2 <- diff(pib_pc_ts, differences = 2)plot(dpib_pc_ts2,
main = "PIB per capita Brasil - série diferenciada",
xlab = "Ano",
ylab = "PIB per capita",
col = "blue",
lwd = 2)
grid(col = "darkgreen", lty = "dotted") # grade leve
A série acima está mais estacionária.
Série de preços de ações
Baixar dados da ação PETR4 (Petrobras - preferencial nominativa é a ação preferencial da Petrobras, com prioridade em dividendos, mas sem direito a voto. O pacote utilizado será o quantmod.
getSymbols("PETR4.SA", src = "yahoo", from = "2020-01-01", to = Sys.Date())[1] "PETR4.SA"class(PETR4.SA)[1] "xts" "zoo"No código acima, a fonte dos dados é o site “Yahoo” e Sys.Date() representa a data de hoje. No resultado, xts significa eXtensible Time Series, ou seja, um objeto de tipo série temporal que é extensível. É uma estrutura de dados especial usada no R para lidar com informações que variam ao longo do tempo — como preços diários de ações, taxas de câmbio, PIB mensal, temperaturas, etc.
Vamos transforma o objeto de tipo xts em data frame.
petra_df <- data.frame(data = index(PETR4.SA), fechamento = as.numeric(Cl(PETR4.SA)))
class(petra_df)[1] "data.frame"No código acima, Cl(.) é uma função do pacote quantmod para pegar o preço de fechamento das ações. Vamos ao gráfico dessa ação.
ggplot(petra_df, aes(x = data, y = fechamento)) +
geom_line(color = "blue", size = 1.1) +
labs(title = "Preço de Fechamento - PETR4 (Petrobras)",
x = "Data",
y = "Preço (R$)") +
theme_minimal()Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Vamos calcular uma média móvel de 15 dias para inserir no gráfico. Vamos usar a função SMA do pacote TTR. Esse pacote já vem junto com o quantmod e serverllll .
mm_15 <- SMA(Cl(PETR4.SA), n = 15)Agora vamos criar um data frame com média móvel de 15 dias.
petra_dfm <- data.frame(
data = index(PETR4.SA),
fechamento = as.numeric(Cl(PETR4.SA)),
mm15 = as.numeric(mm_15)
)O gráfico mais completo fica, então na seguinte forma.
ggplot(petra_dfm, aes(x = data)) +
geom_line(aes(y = fechamento), color = "blue", size = 1.2) +
geom_line(aes(y = mm15), color = "red", size = 1) +
labs(title = "PETR4 - Preço de Fechamento e Média Móvel de 15 dias",
x = "Data",
y = "Preço (R$)") +
theme_minimal()Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).
Diferenciando a série
Vamos diferenciar a série original com os preços de fechamento, segundo a fórmula abaixo:
\(variacaopct= \frac{P_t - P_{t-1}}{P_{t-1}}\)
Para isso, criamos um data frame só com os preços de fechamento. e, em seguida, usamos a função Delt() do pacote quantmod.
fechamento <- Cl(PETR4.SA)Em seguida aplicamos ao fechamento a função Delt(), do pacote quantmod.
variacao <- Delt(fechamento) * 100Agora estruturamos um data frame com o preço de fechamento e sua variação percentual.
petra_dfpct <- data.frame(
data = index(PETR4.SA),
fechamento = as.numeric(fechamento),
variacao_pct = as.numeric(variacao)
)E, finalmente, o seu gráfico.
ggplot(petra_dfpct, aes(x = data, y = variacao_pct)) +
geom_line(color = "darkgreen", size = 1) +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
labs(title = "Variação Percentual Diária - PETR4",
x = "Data",
y = "Variação (%)") +
theme_minimal()Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).
Compare os dois gráficos: os dados originais que esboçam um Random walk e os dados da série diferenciada, mais estacionária.