Dados coletados - Preço por barril do petróleo bruto Brent (FOB) (EIA366_PBRENT366)

Frequência: Diária de 04/01/1986 até 18/02/2025 Fonte: Energy Information Administration (EIA) Unidade: US$ Comentário: Preço por barril do petróleo bruto tipo Brent. Produzido no Mar do Norte (Europa), Brent é uma classe de petróleo bruto que serve como benchmark para o preço internacional de diferentes tipos de petróleo. Neste caso, é valorado no chamado preço FOB (free on board), que não inclui despesa de frete e seguro no preço. Mais informações: https://www.eia.gov/dnav/pet/TblDefs/pet_pri_spt_tbldef2.asp Atualizado em: 27/02/2025

datatable(preco_bar_petroleo, extensions = "Buttons", 
          options = list(dom = "Bftrp", 
                         buttons = c("excel"), pageLength = 10))

Estatísticas descritivas

Podemos perceber que no 3º quartil temos 75% dos preços abaixo de 76,86 USD e na mediana temos o valor central sugere que metado dos preços estão abaixo de 49,19 USD

summary(preco_bar_petroleo$Preco_barril)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.10   20.68   49.19   53.46   76.86  143.95
ggplot(preco_bar_petroleo, aes(y = Preco_barril))+
  geom_boxplot(fill = "lightblue", outlier.colour = "red")+
  labs(title = "Boxplot do Preço do Barril de Petróleo", y = "Preço (USD")+
  theme_minimal()

Criando um gráfico de Densidade para observar a distribuição dos dados ajudando a ver os quartis

ggplot(preco_bar_petroleo, aes(x = Preco_barril)) +
  geom_density(fill = "lightblue", alpha = 0.5) +
  geom_vline(xintercept = quantile(preco_bar_petroleo$Preco_barril, 0.75), 
             color = "orange", linetype = "dashed") +
  labs(title = "Densidade do Preço do Barril de Petróleo", 
       x = "Preço (USD)", 
       y = "Densidade") +
  theme_minimal() +
  geom_text(aes(x = quantile(preco_bar_petroleo$Preco_barril, 0.75), 
                y = 0, label = "3º Quartil"), vjust = 1, color = "orange")

Podemos perceber que os dados estão a grande maior parte no 3º quartil, 75% dos preços estão abaixo de 76,86 USD

Criando a variável log(preco)

preco_bar_petroleo$log_preco_petro <- log(preco_bar_petroleo$Preco_barril)

datatable(preco_bar_petroleo, extensions = "Buttons", 
          options = list(dom = "Bftrp", 
                         buttons = c("excel"), pageLength = 10))

Criando e a primeira diferença do preço do Petróleo

preco_bar_petroleo$dlog_preco_petro <- c(NA, diff(preco_bar_petroleo$log_preco_petro))

datatable(preco_bar_petroleo, extensions = "Buttons", 
          options = list(dom = "Bftrp", 
                         buttons = c("excel"), pageLength = 10))

Criando um gráfico Original Série Temporal do Preço de Petróleo

Criando um gráfico da primeira diferença do logarítmo do Preço de Petróleo

Função de autocorrelação - Abaixo podemos explicar que:

ACF do log do preço de petróleo ao longo de diferentes lags indica uma correlação em diferentes períodos do tempo, podemos perceber que todas as linhas ultrapassam o zero nos diferentes lags e que com base dos preços do passado tem influencia positiva nos valores atuais.Por exemplo, se o preço do petróleo estava alto em um período, as chances de que continue alto é alta e em períodos subsequentes ou vice-versa, indica que a série não é estacionária, sugere que a série pode ter tendência ou sazonalidade Em resumo, indica que a série temporal tem componentes de tendência ou ciclos que se repetem, isso abre uma necessidade de tratamento da série para garantir a estacionariedade antes de aplicar modelos de previsão.

PARF do log(preco petroleo no lag 1,2,3 mostra um valor significativo acima do interlvalo de confiança indicando que esses lags demonstram contribuição para a modelagem do preço do barril de petróleo

ACF da 1ª diferença do Log(preço do petróleo) Podemos ver valores significativos positivos nos primeiros lags, Lag 1 e Lag 2, sugerindo existe uma relação com as variações anteriores inferindo nas variações recentes nos preços do barril de petróleo, o que pode indicar um comportamento persistente nas mudanças de preço.

PACF da 1ª differença do log(preço do petróleo) lags decai rapidamente para zero,indicando que as relações autorregressivas são rapidamente inadequadas após um certo período. APACF mostra significância até o lag 2 e se torna insignificante a partir do lag 3, isso pode indicar que um modelo AR(2), Esse comportamento pode indicar tendências nos preços que precisam ser modelados possivelmente utilizando um modelo ARIMA que inclua termos autorregressivos até a ordem determinada pelas lags significativas.

Obs: Existem fatores externos que podem influenciar nos preços do petroleo, demanda mundial, fotores políticos, pandemias, etc. Que podem impactar diretamente nos preços. Então além das análises Estatisticas, também temos que considerar esses fatores.

par(mfrow =c(2,2)) # arganiza os gráficos lado a lado
acf(preco_bar_petroleo$log_preco_petro, main = "ACF do log(preço petroleo)")
pacf(preco_bar_petroleo$log_preco_petro, main = "PARF do log(preco petroleo)")

# Agora fazendo a 1ª diferença do Log do preço do barril de petróleo
acf(diff(preco_bar_petroleo$log_preco_petro), main = "ACF da 1ª diferença do Log(preço do petróleo)")
pacf(diff(preco_bar_petroleo$log_preco_petro), main = "PACF da 1ª differença do log(preço do petróleo)")

par(mfrow =c(1,1))

Estimação modelo ARMA (1,1) - Ajustando dois modelos ARMA, um objeto para método da soma dos quadrados condicionais e outro para o método que maximiza a função de verossimilhança

soma_quad_cond <- arima(diff(preco_bar_petroleo$log_preco_petro), order = c(1,0,1), method = "CSS")
print(soma_quad_cond)
## 
## Call:
## arima(x = diff(preco_bar_petroleo$log_preco_petro), order = c(1, 0, 1), method = "CSS")
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.4090  -0.5379      1e-04
## s.e.  0.0349   0.0317      2e-04
## 
## sigma^2 estimated as 0.0007048:  part log likelihood = 25110.62
Max_verossim <- arima(diff(preco_bar_petroleo$log_preco_petro), order = c(1,0,1), method = "ML")
print(Max_verossim)
## 
## Call:
## arima(x = diff(preco_bar_petroleo$log_preco_petro), order = c(1, 0, 1), method = "ML")
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.4110  -0.5397      1e-04
## s.e.  0.0349   0.0317      2e-04
## 
## sigma^2 estimated as 0.0007047:  log likelihood = 25111.1,  aic = -50214.21

Explicando: O metodo da soma dos quadrados condicional “CSS” ele estima os parametros minimizando a soma dos quadrados condicionais, mas pode não ser tão precisa quanto a outros métodos. Por outro lado o método do modelo de Máxima Verossimilhança busca parametros para maximizar a função de Verossimilhança, com isso acabam sendo mais robusto e precisos em comparação com o “CSS”

Vamos fazer uma comparação no gráfico para melhor visualização

css_coef <- coef(soma_quad_cond)
ml_coef <- coef(Max_verossim)
erro_css <- sqrt(diag(soma_quad_cond$var.coef)) # erro padrão extraídos soma dos quadrados condicionais
erro_ml <- sqrt(diag(Max_verossim$var.coef))

result_compara <- data.frame(
  Coeficiente = rep(names(css_coef), 2),
  Valor = c(css_coef, ml_coef),
  Erro_padrao = c(erro_css, erro_ml),
  Metodo = rep(c("CSS", "ML"), each = length(css_coef))
)


ggplot(result_compara, aes(x = Coeficiente, y = Valor, fill = Metodo))+
  geom_bar(stat = "Identity", position = "dodge", width = 0.7)+
  geom_errorbar(aes(ymin = Valor - Erro_padrao, ymax = Valor + Erro_padrao),
                width = 0.2, position = position_dodge(0.7))+
  geom_text(aes(label = round(Valor, 5)), 
            position = position_dodge(0.7), 
            vjust = -0.9,  # Ajusta a posição do texto acima da barra
            color = "black") +
    labs(title = "Comparação de Coeficientes dos Modelos ARMA",
       y = "VAlor do Coeficiente",
       x = "Coeficientes")+
  theme_minimal()

Tabela das comparações soma dos quadrados condicional “CSS” e método do modelo de Máxima Verossimilhança “ML”

datatable(result_compara, extensions = "Buttons", 
          options = list("Bftrp",pageLength = 10))

Mostrando as informações das estimações ARMA método da soma dos quadrados condicionais e método maximiza a função de verossimilhança

info_soma_quad_cond <- tidy(soma_quad_cond) %>%
  mutate(model = "soma_quad_cond")


info_Max_verossim <- tidy(Max_verossim) %>%
  mutate(model = "Max_verossim")


resultados <- bind_rows(info_soma_quad_cond, info_Max_verossim)%>%
  rename(Coeficiente = estimate, 
         Erro_padrão = std.error,
         Modelo = model)

Tabela dos dados soma dos quadrados condicionais e método maximiza a função de verossimilhança

datatable(resultados, extensions = "Buttons",
          options = list(dom = "Bftrp",
                         buttons =c("excel"), pageLength = 10))

Fazendo o Diagnóstico dos resíduos - Verificando o tamanho da amostra dos resíduos

Obs: Caso as amostras tenham mais de 5000 registros ou menor que 3, o teste de Shapiro-Wilk. não é recomendado para amostras grandes, pois tende a detectar pequenas desvios da normalidade.

resid_diag <- soma_quad_cond$residuals

#length(resid_diag)

print(length(resid_diag))
## [1] 11363
print("Não podemos usar o teste de Shapiro-Wilk, amostras muito grande")
## [1] "Não podemos usar o teste de Shapiro-Wilk, amostras muito grande"

Fazendo os testes de normalidade

Anderson-Darling Normality Test - O p-value muito baixo em 2.2e-16, o que indica que há evidências forte para rejeitar H0. Os resíduos seguem uma distribuição normal, logo os residuos sugere que não são normalmente distribuídos

print(ad.test(resid_diag))
## 
##  Anderson-Darling normality test
## 
## data:  resid_diag
## A = 169.32, p-value < 2.2e-16

Jarque-Bera Test - Verifica a normalidade com base na assimetria e na curtose Podemos perceber que o valor p-value ficou em 2.2e-16 confirmando o do Anderson-Darling, logo, os residuos sugere que não são normalmente distribuídos.Abaixo plotamos um gráfico para melhor visualização e confirmação

print(jarque.bera.test(resid_diag))
## 
##  Jarque Bera Test
## 
## data:  resid_diag
## X-squared = 1029117, df = 2, p-value < 2.2e-16
library(moments)
assimetria <- skewness(resid_diag)
curtose <- kurtosis(resid_diag)

result_assim_curt <- data.frame(
  Medida = c("Assimetria", "Curtose"),
  Valor = c(assimetria, curtose)
)


ggplot(result_assim_curt, aes(x = Medida, y = Valor, fill = Medida))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = round(Valor, 2)), vjust = -0.5)+
  labs(title = "Assimetria e Curtose dos Resíduos",
       y = "valor",
       x = "")+
  theme_minimal()

Kolmogorov-Smirnov Test - Compara a distribuição empírica dos resíduos com a distribuição normal teórica, usando a média e o desvio padrão dos resíduos, apesar do valor D ficar com maior diferença entre os dois testes anteriores podendo não ser significativo , pelo p-value muito baixo de 2.2e-16, podemos concluir que existe diferença significativa devemos rejeitar H0, que os resíduos seguem uma distribuição normal

print(ks.test(resid_diag, "pnorm", mean(resid_diag), sd(resid_diag)))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resid_diag
## D = 0.075123, p-value < 2.2e-16
## alternative hypothesis: two-sided
print("Conclusão: Todos os três testes de Normalidade dos Resíduos indicam que os resíduos não seguem uma Distribuição Normal")
## [1] "Conclusão: Todos os três testes de Normalidade dos Resíduos indicam que os resíduos não seguem uma Distribuição Normal"

Para encontrar lambda ideal, usamos os dados em nível (preço)

result_boxcox <- boxcox(lm(preco_bar_petroleo$Preco_barril ~ 1, data = preco_bar_petroleo), lambda = seq(-2, 2, by = 0.1))

Extrair o lambda ótimo

otimo_lambda <- result_boxcox$x[which.max(result_boxcox$y)]
print(paste("Lambda ótimo:", otimo_lambda))
## [1] "Lambda ótimo: 0.181818181818182"

Transformar a Série com Box-Cox

if (otimo_lambda == 0) {
  preco_bar_petroleo$Preco_bcox <- log(preco_bar_petroleo$Preco_barril)  # Se lambda for 0, aplicamos log
} else {
  preco_bar_petroleo$Preco_bcox <- (preco_bar_petroleo$Preco_barril^otimo_lambda - 1) / otimo_lambda
}

Criando a 1ª diferença da série transformada

preco_bar_petroleo$log_preco_bcox <- log(preco_bar_petroleo$Preco_bcox)
preco_bar_petroleo$dlog_preco_petro <- c(NA, diff(preco_bar_petroleo$log_preco_bcox))

Ajustar modelo ARMA(1,1) na série transformada

st_boxcox <- arima(preco_bar_petroleo$dlog_preco_petro, order = c(1,0,1), method = "ML")
print(st_boxcox)
## 
## Call:
## arima(x = preco_bar_petroleo$dlog_preco_petro, order = c(1, 0, 1), method = "ML")
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.4344  -0.5532      0e+00
## s.e.  0.0373   0.0341      1e-04
## 
## sigma^2 estimated as 0.00011:  log likelihood = 35665.93,  aic = -71323.85
print(ad.test(st_boxcox$residuals))
## 
##  Anderson-Darling normality test
## 
## data:  st_boxcox$residuals
## A = 226.07, p-value < 2.2e-16

Gráfico dos resíduos - Podemos ver resíduos não seguem uma Distribuição Normal

par(mfrow = c(1,2))
hist(resid_diag, breaks = 20, col = "blue", main = "Histograma dos Resíduos", xlab = "Resíduo")
qqnorm(resid_diag); qqline(resid_diag, col = "red")

par(mfrow = c(1,1))
acf(resid_diag, main = "ACF dos Resíduos")

pacf(resid_diag, main = "PACF dos Resíduos")

Teste ARCH-LM para heterocedasticidade condicional - verificando se os erros (ou “resíduos”) de um modelo de previsão variam muito ao longo do tempo.

print(ArchTest(resid_diag, lags = 12))
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  resid_diag
## Chi-squared = 1109.3, df = 12, p-value < 2.2e-16

Podemos perceber que os erros não são constantes e precisamos levar em conta para melhorar nossas análises e previsões. esses erros não são sempre iguais. Às vezes, eles são maiores e às vezes, menores.Isso é importante notar pois afeta nas previsões. Pelo resultado de Chi-squared = 1109.3 e um valor-p < 2.2e-16 é extremamente baixo, isso indica que rejeitamos a hipótese nula. Ou seja, há evidência de que a variância dos erros não é constante, ou seja os dados mudam ao longo do tempo. conclusão de que a variância dos resíduos é heterocedástica, ou seja, varia em vez de ser constante.