Irei apresentar neste site os processos realizados para modelagem da volatilidade do café arábica, assim como os scripts utilizados.

library(tseries)  # Para o teste PP
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(readxl) # Para carregar o arquivo em xlsx 
library(fracdiff) # Para estimar o valor de d
library(forecast)
library(tseries)
library(randtests) # Para realizar o teste de Cox e Stuart
## 
## Anexando pacote: 'randtests'
## O seguinte objeto é mascarado por 'package:tseries':
## 
##     runs.test
library(descomponer) # Função periodograma    
library(TSA) # Função perodogram
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Anexando pacote: 'TSA'
## Os seguintes objetos são mascarados por 'package:stats':
## 
##     acf, arima
## O seguinte objeto é mascarado por 'package:utils':
## 
##     tar
library(boot) #estimativa bootstrap
library(lmtest)
## Carregando pacotes exigidos: zoo
## 
## Anexando pacote: 'zoo'
## Os seguintes objetos são mascarados por 'package:base':
## 
##     as.Date, as.Date.numeric
library(fGarch) #Modelo Garch
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")
coffe <- read_excel("coffe.xlsx")

dados <- coffe
head(dados); tail(dados)
## # A tibble: 6 × 1
##   Reais
##   <dbl>
## 1  141.
## 2  136.
## 3  151.
## 4  146.
## 5  144.
## 6  133.
## # A tibble: 6 × 1
##   Reais
##   <dbl>
## 1  878.
## 2  868.
## 3  803.
## 4  801.
## 5  846.
## 6  925.

Inicialmente plotamos a série de preços

y<-dados$Reais 
rm(y)

y = ts(dados$Reais, start = c(1995,1), end = c(2023,12), freq = 12)
plot(y, xlab= c("Tempo"),ylab = c("Preço (R$)"))

A série log de retornos foi plotada logo em sequência

y1 = ts(diff(log(y)), start = c(1995,1), end = c(2023,12), freq = 12) # r = ln[y_t/y_(t-1)]
plot(y1, xlab= c("Tempo"),ylab = c("Série de log retornos"))

Verificando ACF e PACF

n = length(y) #tamanho da amostra
par(mfrow = c(2,2))
acf(y, n)
acf(y1, n)
pacf(y, n)
pacf(y1, n)

O teste de Cox-Stuart verifica se a série possui alguma tendencia significativa.

cox.stuart.test(y1)
## 
##  Cox Stuart test
## 
## data:  y1
## statistic = 82, n = 174, p-value = 0.4952
## alternative hypothesis: non randomness

O Periodograma verifica se a série possui algum periodo e frequência significativos

pegram<-periodogram(y1)

Encontrando a frequência dominante

freq_dominant <- pegram$freq[which.max(pegram$spec)]
period <- 1 / freq_dominant
cat("Frequência dominante:", freq_dominant, "\n")
## Frequência dominante: 0.08333333
cat("Período dominante:", period, "\n")
## Período dominante: 12

Utilizarei o método de reamostragem bootstrap para verificar o intervalo de significância da frequência.

# Função de estatística de interesse
stat_function <- function(data, indices) {
  sample_data <- data[indices]
  pgram <- spectrum(sample_data, plot = FALSE)
  max(pgram$spec)
}

# Aplicando o bootstrap
boot_result <- boot(data = y1, statistic = stat_function, R = 1000)

# Calculando o intervalo de confiança
boot_ci <- boot.ci(boot_result, type = "perc")

# Visualizando os resultados
print(boot_ci)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_result, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 0.0211,  0.0515 )  
## Calculations and Intervals on Original Scale

Ajustando componentes seno e cosseno para modelar a periodicidade encontrada.

(wjj = 2*pi*freq_dominant) 
## [1] 0.5235988
t = seq(1,n)
sen = sin(wjj*t)
cos = cos(wjj*t)
Xs = cbind(cos)

Ajustando um modelo ARIMA (1,0,0) com uma componente cosseno.

arima = Arima(y1, order = c(1,0,0), xreg = Xs, include.mean = F); summary(arima); coeftest(arima)
## Series: y1 
## Regression with ARIMA(1,0,0) errors 
## 
## Coefficients:
##          ar1     cos
##       0.2472  0.0254
## s.e.  0.0519  0.0067
## 
## sigma^2 = 0.00503:  log likelihood = 428.05
## AIC=-850.09   AICc=-850.02   BIC=-838.53
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE     MAPE      MASE
## Training set 0.003963162 0.07071869 0.05442427 -182.6881 436.0142 0.7122079
##                    ACF1
## Training set -0.0010983
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 0.2472137  0.0519374  4.7598 1.937e-06 ***
## cos 0.0254105  0.0067294  3.7761 0.0001593 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
y_hat = arima$fitted
plot(y1, col = "black", type = "l")
lines(y_hat, col = "red", add = T)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" não é um parâmetro
## gráfico

Análise de resíduo

rm = residuals(arima)

hist(rm, 
     main = "Histograma dos Resíduos do Modelo ARIMA", # Título do histograma
     xlab = "Valores dos Resíduos",                    # Legenda do eixo X
     ylab = "Frequência",                              # Legenda do eixo Y
     col = "lightblue")

# Iniciando a Análise Residual 

shapiro.test(rm) # Não rejeitar H0 indica normalidade
## 
##  Shapiro-Wilk normality test
## 
## data:  rm
## W = 0.98722, p-value = 0.003688
Box.test(rm, 15, type = c("Box-Pierce")) # Não Rejeitar H0 indica ausência de correlação
## 
##  Box-Pierce test
## 
## data:  rm
## X-squared = 19.427, df = 15, p-value = 0.195
Box.test(rm^2, 15, type = c("Box-Pierce")) # Não Rejeitar H0 indica variância Constante
## 
##  Box-Pierce test
## 
## data:  rm^2
## X-squared = 36.455, df = 15, p-value = 0.001519
par(mfrow = c(2,1))
acf(rm, n); pacf(rm, n)

acf(rm^2, n); pacf(rm^2, n)

qqnorm(rm); qqline(rm, col = "red"); hist(rm, prob = F)

Para contornar a autocorrelação ainda presente no quadrado dos resíduos, ajustamos um modelo GARCH(1,1) juntamente com uma distribuição Skew_t para capturar a leptocurtica presente nos dados.

arch = garchFit(~ garch(1,1), include.mean = T, data = rm, trace = F, cond.dist = c("sstd"))
summary(arch)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = rm, cond.dist = c("sstd"), 
##     include.mean = T, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x0000020a8a22e9b0>
##  [data = rm]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##         mu       omega      alpha1       beta1        skew       shape  
## 4.4408e-03  7.4475e-05  4.6890e-02  9.3856e-01  1.1754e+00  1.0000e+01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     4.441e-03   3.467e-03    1.281  0.20024    
## omega  7.447e-05   6.583e-05    1.131  0.25789    
## alpha1 4.689e-02   2.099e-02    2.234  0.02549 *  
## beta1  9.386e-01   2.521e-02   37.228  < 2e-16 ***
## skew   1.175e+00   9.297e-02   12.644  < 2e-16 ***
## shape  1.000e+01   3.863e+00    2.589  0.00963 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  443.3588    normalized:  1.274019 
## 
## Description:
##  Sat Oct 26 10:43:43 2024 by user: Lucas 
## 
## 
## Standardised Residuals Tests:
##                                  Statistic      p-Value
##  Jarque-Bera Test   R    Chi^2  19.3181702 0.0000638429
##  Shapiro-Wilk Test  R    W       0.9879044 0.0053916406
##  Ljung-Box Test     R    Q(10)  10.8914341 0.3660386716
##  Ljung-Box Test     R    Q(15)  17.1744736 0.3085418293
##  Ljung-Box Test     R    Q(20)  23.1125016 0.2833042954
##  Ljung-Box Test     R^2  Q(10)  11.9281697 0.2898926130
##  Ljung-Box Test     R^2  Q(15)  16.8549609 0.3276085662
##  Ljung-Box Test     R^2  Q(20)  22.2484952 0.3271802939
##  LM Arch Test       R    TR^2   15.4299134 0.2187632180
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -2.513556 -2.447139 -2.514137 -2.487114
#plot(arch) # 10, 11, 1337.903.873 26.594.653

#cond.dist = c("norm", "snorm", "ged", "sged",
#"std", "sstd", "snig", "QMLE"),

Ao fim, plotamos a previsão de preços para os próximos mêses

prev = predict(arch, h = 3, plot = T, trace = F)

prev
##    meanForecast  meanError standardDeviation lowerInterval upperInterval
## 1   0.004440783 0.06676355        0.06676355    -0.1177099     0.1472168
## 2   0.004440783 0.06683545        0.06683545    -0.1178415     0.1473706
## 3   0.004440783 0.06690623        0.06690623    -0.1179710     0.1475220
## 4   0.004440783 0.06697591        0.06697591    -0.1180985     0.1476710
## 5   0.004440783 0.06704451        0.06704451    -0.1182240     0.1478177
## 6   0.004440783 0.06711203        0.06711203    -0.1183475     0.1479621
## 7   0.004440783 0.06717851        0.06717851    -0.1184691     0.1481043
## 8   0.004440783 0.06724395        0.06724395    -0.1185889     0.1482442
## 9   0.004440783 0.06730838        0.06730838    -0.1187067     0.1483820
## 10  0.004440783 0.06737182        0.06737182    -0.1188228     0.1485176