MODELO GARCH - CAFÉ ARÁBICA

 UNIVERSIDADE FEDERAL DA PARAÍBA

## UNIVERSIDADE FEDERAL DA PARAÍBA
## DEPARTAMENTO DE ECONOMIA
## CURSO DE GESTÃO DE RISCO
## PROFESSOR: SINÉZIO FERNANDES MAIA
## ALUNO: JOSUÉ DE MENESES LOPES

## Estimação de Modelos GARCH
## ============================================================================

## Limpeza do ambiente e configurações iniciais
remove(list=ls())
options(scipen=999)
options(max.print = 1000000)
par(mfrow=c(1,1))

## Carregamento das bibliotecas necessárias
library(tidyverse)
library(BatchGetSymbols)
library(rbcb)
library(readxl)
library(tseries)
library(PerformanceAnalytics)
library(deflateBR)
library(stats)
library(fBasics)
library(moments)
library(dynlm)
library(lmtest)
library(whitestrap)
library(FinTS)
library(fGarch)
library(rugarch)
library(urca)

Dados

## CAFÉ SPOT - MENSAL

cafe_mensal_spot <- read_excel("cafe_mensal_spot.xlsx");nrow(cafe_mensal_spot)
cafe_mensal_spot_ts <- ts(cafe_mensal_spot$fechamento_mensal_spot, frequency = 12, start = c(2014,01), end = c(2024,07))


## CAFÉ FUTURO - MENSAL

cafe_mensal_fut <- read_excel("cafe_mensal_fut.xlsx");nrow(cafe_mensal_fut)
cafe_mensal_fut_ts <- ts(cafe_mensal_fut$fechamento_mensal_fut, frequency = 12, start = c(2014,01), end = c(2024,07))


## ÍNDICE CAFÉ - MENSAL

indice_cafe_mensal <- get_series(29041, start_date = "2014-01-01", end_date = "2024-08-01");nrow(indice_cafe_mensal)
indice_cafe_mensal_ts <- ts(indice_cafe_mensal$"29041", frequency = 12, start = c(2014,01), end = c(2024,07))

Gráfico de Preços Spot e Futuro do Café Arábica

plot(cafe_mensal_spot_ts, col = "blue", lwd = 2, ylim = range(c(cafe_mensal_spot_ts, cafe_mensal_fut_ts)), xlab = "Ano", ylab = "Preço em Dólares", main = "Preço Spot e Futuro do Café Arábica - Nominal")
lines(cafe_mensal_fut_ts, col = "red", lwd = 2)
legend("topleft", legend = c("Preço Spot", "Preço Futuro"), col = c("blue", "red"), lwd = 2)
grid()

Gráfico de Preços Futuro e Índice de Commodities do Café Arábica

plot(cafe_mensal_fut_ts, col = "red", lwd = 2, ylim = range(c(cafe_mensal_fut_ts, indice_cafe_mensal_ts)), xlab = "Ano", ylab = "Preço em Dólares", main = "Preço Futuro do Café Arábica e Índice de Commodities - Nominal")
lines(indice_cafe_mensal_ts, col = "darkgreen", lwd = 2)
legend("topleft", legend = c("Preço Futuro", "Índice de Commodities"), col = c("red", "darkgreen"), lwd = 2)
grid()

Deflacionamento dos Preços

times <- seq(as.Date("2014/1/1"), by = "month", length.out = 127)
cafe_spot_R <- deflate(cafe_mensal_spot_ts, nominal_dates = times, real_date = "07/2024", index = "igpdi")
cafe_fut_R <- deflate(cafe_mensal_fut_ts, nominal_dates = times, real_date = "07/2024", index = "igpdi")
cafe_indice_R <- deflate(indice_cafe_mensal_ts, nominal_dates = times, real_date = "07/2024", index = "igpdi")

Gráfico de Preços Deflacionados

plot(cafe_spot_R, col = "blue", lwd = 2, ylim = range(c(cafe_spot_R, cafe_fut_R)), xlab = "Ano", ylab = "Preço em Dólares", main = "Preço Spot e Futuro do Café Arábica - Deflacionado")
lines(cafe_fut_R, col = "red", lwd = 2)
legend("topleft", legend = c("Preço Spot", "Preço Futuro"), col = c("blue", "red"), lwd = 2)
grid()

Análise dos Retornos do Preço Futuro - Nominal

par(mfrow=c(2,1))
chart.TimeSeries(cafe_mensal_fut_ts, main = "Preços Café Futuro - Nominal")
chart.TimeSeries(diff(log(cafe_mensal_fut_ts)), main = "Retorno Café Futuro - Nominal")

Análise dos Retornos do Preço Futuro - Deflacionado

par(mfrow=c(2,1))
chart.TimeSeries(cafe_fut_R, main = "Preços Café Futuro - Deflacionado")
chart.TimeSeries(diff(log(cafe_fut_R)), main = "Retorno Café Futuro - Deflacionado")

Análise dos Retornos do Preço Spot - Nominal

par(mfrow=c(2,1))
chart.TimeSeries(cafe_mensal_spot_ts, main = "Preços Café Spot - Nominal")
chart.TimeSeries(diff(log(cafe_mensal_spot_ts)), main = "Retorno Café Spot - Nominal")

Análise dos Retornos do Preço Spot - Deflacionado

par(mfrow=c(2,1))
chart.TimeSeries(cafe_spot_R, main = "Preços Café Spot - Deflacionado")
chart.TimeSeries(diff(log(cafe_spot_R)), main = "Retorno Café Spot - Deflacionado")

Análise dos Retornos do Índice de Commodities - Nominal

par(mfrow=c(2,1))
chart.TimeSeries(indice_cafe_mensal_ts, main = "Preços Índice de Commodities - Nominal")
chart.TimeSeries(diff(log(indice_cafe_mensal_ts)), main = "Retorno Índice de Commodities - Nominal")

Análise dos Retornos do Índice de Commodities - Deflacionado

par(mfrow=c(2,1))
chart.TimeSeries(cafe_indice_R, main = "Preços Índice de Commodities - Deflacionado")
chart.TimeSeries(diff(log(cafe_indice_R)), main = "Retorno Índice de Commodities - Deflacionado")

Média - Café Futuro

Media.cafe_fut_R=mean(diff(log(cafe_fut_R)));Media.cafe_fut_R

[1] -0.000599504

Desvio - Café Futuro

Desvio.cafe_fut_R=sd(diff(log(cafe_fut_R)));Desvio.cafe_fut_R

[1] 0.08733543

Histograma - Café Futuro

par(mfrow=c(1,1))
chart.Histogram(diff(log(cafe_fut_R)))

summary(diff(log(cafe_fut_R)))

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.2656255 -0.0586035 -0.0060793 -0.0005995  0.0479850  0.3430820 

Statcafe_fut_R<-basicStats(diff(log(cafe_fut_R)));Statcafe_fut_R

            diff.log.cafe_fut_R.
nobs                  126.000000
NAs                     0.000000
Minimum                -0.265625
Maximum                 0.343082
1. Quartile            -0.058603
3. Quartile             0.047985
Mean                   -0.000600
Median                 -0.006079
Sum                    -0.075538
SE Mean                 0.007780
LCL Mean               -0.015998
UCL Mean                0.014799
Variance                0.007627
Stdev                   0.087335
Skewness                0.323954
Kurtosis                1.362544

Skewness - Café Futuro

skewness(data.frame(diff(log(cafe_fut_R))))

diff.log.cafe_fut_R.. 
            0.3278489 

Curtose - Café Futuro

kurtosis(data.frame(diff(log(cafe_fut_R))))

diff.log.cafe_fut_R.. 
             4.432624 

Jarque Bera - Café Futuro

jarqueberaTest((diff(log(cafe_fut_R))))

Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 13.0323
  P VALUE:
    Asymptotic p Value: 0.001479 

ACF e PACF - Café Futuro

acf((diff(log(cafe_fut_R))))

pacf((diff(log(cafe_fut_R))))

Box-Pierce - Café Futuro

Box.test((diff(log(cafe_fut_R))), lag=1)

    Box-Pierce test

data:  (diff(log(cafe_fut_R)))
X-squared = 1.6917, df = 1, p-value = 0.1934

Box.test((diff(log(cafe_fut_R))), lag=6)

    Box-Pierce test

data:  (diff(log(cafe_fut_R)))
X-squared = 8.5738, df = 6, p-value = 0.199

Box.test((diff(log(cafe_fut_R))), lag=10)

    Box-Pierce test

data:  (diff(log(cafe_fut_R)))
X-squared = 14.41, df = 10, p-value = 0.1551

Box.test((diff(log(cafe_fut_R))), lag=24)

    Box-Pierce test

data:  (diff(log(cafe_fut_R)))
X-squared = 24.297, df = 24, p-value = 0.4447

CAFÉ SPOT

Média - Café Spot

Media.cafe_spot_R=mean(diff(log(cafe_spot_R)));Media.cafe_spot_R

[1] 0.00001325791

Desvio - Café Spot

Desvio.cafe_spot_R=sd(diff(log(cafe_spot_R)));Desvio.cafe_spot_R

[1] 0.06964253

Histograma - Café Spot

par(mfrow=c(1,1))
chart.Histogram(diff(log(cafe_spot_R)))

summary(diff(log(cafe_spot_R)))

       Min.     1st Qu.      Median        Mean     3rd Qu.        Max. 
-0.16597329 -0.04584737 -0.00552718  0.00001326  0.04353712  0.23316669 

Statcafe_spot_R<-basicStats(diff(log(cafe_spot_R)));Statcafe_spot_R

            diff.log.cafe_spot_R.
nobs                   126.000000
NAs                      0.000000
Minimum                 -0.165973
Maximum                  0.233167
1. Quartile             -0.045847
3. Quartile              0.043537
Mean                     0.000013
Median                  -0.005527
Sum                      0.001670
SE Mean                  0.006204
LCL Mean                -0.012266
UCL Mean                 0.012292
Variance                 0.004850
Stdev                    0.069643
Skewness                 0.344783
Kurtosis                 0.368316

Skewness - Café Spot

skewness(data.frame(diff(log(cafe_spot_R))))

diff.log.cafe_spot_R.. 
             0.3489284 

Curtose - Café Spot

kurtosis(data.frame(diff(log(cafe_spot_R))))

diff.log.cafe_spot_R.. 
              3.422424 

Jarque Bera - Café Spot

jarqueberaTest((diff(log(cafe_spot_R))))

Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 3.4936
  P VALUE:
    Asymptotic p Value: 0.1743 

ACF e PACF - Café Spot

acf((diff(log(cafe_spot_R))))

pacf((diff(log(cafe_spot_R))))

Box-Pierce - Café Spot

Box.test((diff(log(cafe_spot_R))), lag=1)

    Box-Pierce test

data:  (diff(log(cafe_spot_R)))
X-squared = 4.411, df = 1, p-value = 0.03571

Box.test((diff(log(cafe_spot_R))), lag=6)

    Box-Pierce test

data:  (diff(log(cafe_spot_R)))
X-squared = 7.5921, df = 6, p-value = 0.2695

Box.test((diff(log(cafe_spot_R))), lag=10)

    Box-Pierce test

data:  (diff(log(cafe_spot_R)))
X-squared = 11.803, df = 10, p-value = 0.2985

Box.test((diff(log(cafe_spot_R))), lag=24)

    Box-Pierce test

data:  (diff(log(cafe_spot_R)))
X-squared = 30.853, df = 24, p-value = 0.1581

CAFÉ ÍNDICE

Média - Café Índice

Media.cafe_indice_R=mean(diff(log(cafe_indice_R)));Media.cafe_indice_R

[1] -0.003658919

Desvio - Café Índice

Desvio.cafe_indice_R=sd(diff(log(cafe_indice_R)));Desvio.cafe_indice_R

[1] 0.02992322

Histograma - Café Índice

par(mfrow=c(1,1))
chart.Histogram(diff(log(cafe_indice_R)))

summary(diff(log(cafe_indice_R)))

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.132814 -0.020378 -0.001337 -0.003659  0.013966  0.061403 

Statcafe_indice_R<-basicStats(diff(log(cafe_indice_R)));Statcafe_indice_R

            diff.log.cafe_indice_R.
nobs                     126.000000
NAs                        0.000000
Minimum                   -0.132814
Maximum                    0.061403
1. Quartile               -0.020378
3. Quartile                0.013966
Mean                      -0.003659
Median                    -0.001337
Sum                       -0.461024
SE Mean                    0.002666
LCL Mean                  -0.008935
UCL Mean                   0.001617
Variance                   0.000895
Stdev                      0.029923
Skewness                  -0.825526
Kurtosis                   2.385806

Skewness - Café Índice

skewness(data.frame(diff(log(cafe_indice_R))))

diff.log.cafe_indice_R.. 
              -0.8354525 

Curtose - Café Índice

kurtosis(data.frame(diff(log(cafe_indice_R))))

diff.log.cafe_indice_R.. 
                5.472324 

Jarque Bera - Café Índice

jarqueberaTest((diff(log(cafe_indice_R))))

Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 46.7476
  P VALUE:
    Asymptotic p Value: 0.00000000007061 

ACF e PACF - Café Índice

acf((diff(log(cafe_indice_R))))

pacf((diff(log(cafe_indice_R))))

Box-Pierce - Café Índice

Box.test((diff(log(cafe_indice_R))), lag=1)

    Box-Pierce test

data:  (diff(log(cafe_indice_R)))
X-squared = 3.6329, df = 1, p-value = 0.05665

Box.test((diff(log(cafe_indice_R))), lag=6)

    Box-Pierce test

data:  (diff(log(cafe_indice_R)))
X-squared = 12.58, df = 6, p-value = 0.0502

Box.test((diff(log(cafe_indice_R))), lag=10)

    Box-Pierce test

data:  (diff(log(cafe_indice_R)))
X-squared = 16.856, df = 10, p-value = 0.07761

Box.test((diff(log(cafe_indice_R))), lag=24)

    Box-Pierce test

data:  (diff(log(cafe_indice_R)))
X-squared = 19.621, df = 24, p-value = 0.7181

Estimação dos Retornos por ARIMA

ARIMA para o Preço Futuro

ArimaFut111 <- arima(cafe_fut_R, order = c(1,1,1), method = c("CSS"))
ResiduosFUT <- residuals(ArimaFut111)
par(mfrow=c(1,1))
chart.TimeSeries(ResiduosFUT)

chart.ACF(ResiduosFUT)

tsdiag(ArimaFut111)

ARIMA para o Preço Spot

ArimaSpot111 <- arima(cafe_spot_R, order = c(1,1,0), method = c("CSS"))
ResiduosSPOT <- residuals(ArimaSpot111)
par(mfrow=c(1,1))
chart.TimeSeries(ResiduosSPOT)

chart.ACF(ResiduosSPOT)

tsdiag(ArimaSpot111)

ARIMA para o Índice de Commodities

ArimaIndice110 <- arima(cafe_indice_R, order = c(1,1,0), method = c("CSS"))
ResiduosINDICE <- residuals(ArimaIndice110)
par(mfrow=c(1,1))
chart.TimeSeries(ResiduosINDICE)

chart.ACF(ResiduosINDICE)

tsdiag(ArimaIndice110)

Dickey-Fuller

Dickey-Fuller - Café Futuro

ddickyFUT <- ur.df(ResiduosFUT, lags = 0, type="trend")
summary(ddickyFUT)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-54.729 -12.579  -1.589  12.291  65.931 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept) -8.56367    3.83833  -2.231              0.0275 *  
z.lag.1     -1.04595    0.09004 -11.617 <0.0000000000000002 ***
tt           0.10797    0.05231   2.064              0.0411 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.02 on 123 degrees of freedom
Multiple R-squared:  0.5232,    Adjusted R-squared:  0.5154 
F-statistic: 67.48 on 2 and 123 DF,  p-value: < 0.00000000000000022


Value of test-statistic is: -11.6168 44.9835 67.4752 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

Dickey-Fuller - Café Spot

ddickySPOT <- ur.df(ResiduosSPOT, lags = 0, type="trend")
summary(ddickySPOT)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.331  -8.896  -0.473  10.010  57.470 

Coefficients:
            Estimate Std. Error t value            Pr(>|t|)    
(Intercept) -3.67102    2.80916  -1.307               0.194    
z.lag.1     -1.13041    0.08939 -12.646 <0.0000000000000002 ***
tt           0.04843    0.03837   1.262               0.209    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.59 on 123 degrees of freedom
Multiple R-squared:  0.5653,    Adjusted R-squared:  0.5582 
F-statistic: 79.97 on 2 and 123 DF,  p-value: < 0.00000000000000022


Value of test-statistic is: -12.6465 53.311 79.9665 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

Dickey-Fuller - Café Índice

ddickyIND <- ur.df(ResiduosINDICE, lags = 0, type="trend")
summary(ddickyIND)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.512  -3.423   0.414   3.429  11.422 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept) -2.60658    0.99411  -2.622              0.00984 ** 
z.lag.1     -1.06519    0.09013 -11.818 < 0.0000000000000002 ***
tt           0.02818    0.01346   2.093              0.03842 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.406 on 123 degrees of freedom
Multiple R-squared:  0.5317,    Adjusted R-squared:  0.5241 
F-statistic: 69.84 on 2 and 123 DF,  p-value: < 0.00000000000000022


Value of test-statistic is: -11.8181 46.5577 69.835 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

Estimação do modelo MQO Dinâmico - CAFÉ FUTURO

dcafe_fut_R <- diff(log(cafe_fut_R))
cbind(cafe_fut_R, dcafe_fut_R)

Regressão Dinâmica usando a série diferenciada (MQO Dinâmico)

lag <- stats::lag
MQODin <- dynlm(dcafe_fut_R ~ lag(dcafe_fut_R, -1))
coef(MQODin)

         (Intercept) lag(dcafe_fut_R, -1) 
        -0.003412599         -0.115857422 

summary(MQODin)

Time series regression with "ts" data:
Start = 2014(3), End = 2024(7)

Call:
dynlm(formula = dcafe_fut_R ~ lag(dcafe_fut_R, -1))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.244870 -0.050422 -0.005169  0.053533  0.190029 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -0.003413   0.007310  -0.467    0.641
lag(dcafe_fut_R, -1) -0.115857   0.083702  -1.384    0.169

Residual standard error: 0.08173 on 123 degrees of freedom
Multiple R-squared:  0.01534,   Adjusted R-squared:  0.007332 
F-statistic: 1.916 on 1 and 123 DF,  p-value: 0.1688

Resíduos da regressão

residuosMQO_FUT <- residuals(MQODin)
residuosMQO_FUT
chart.TimeSeries(residuosMQO_FUT)

Histograma

chart.Histogram(residuosMQO_FUT)

## Gráfico de histograma com curva de distribuição normal sobreposta
par(mfrow = c(1, 1))
hist(residuosMQO_FUT, main = "", col = "cadetblue", prob = TRUE, 
     xlab = names(residuosMQO_FUT)[1], breaks = 30)
curve(expr = dnorm(x, mean = mean(residuosMQO_FUT), sd = sd(residuosMQO_FUT)), 
      col = "red", add = TRUE, lwd = 2)

Gráfico QQ-Norm (Normal Quantile-Quantile Plot) para verificar a normalidade dos resíduos

qqnorm(residuosMQO_FUT, col = "blue")
qqline(residuosMQO_FUT, col = "red")

Testes de Normalidade

Teste de Jarque-Bera

jarqueberaTest(residuosMQO_FUT)

Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 0.1218
  P VALUE:
    Asymptotic p Value: 0.9409 

Teste de Shapiro-Wilk

shapiro.test(residuosMQO_FUT)

    Shapiro-Wilk normality test

data:  residuosMQO_FUT
W = 0.99104, p-value = 0.6008

Testes de Heterocedasticidade

length(residuosMQO_FUT)

[1] 125

length(residuosMQO_FUT) * 0.15

[1] 18.75

Teste de Goldfeld-Quandt

gqtest(MQODin, fraction = 18.75, alternative = "greater")

    Goldfeld-Quandt test

data:  MQODin
GQ = 1.1831, df1 = 52, df2 = 51, p-value = 0.2746
alternative hypothesis: variance increases from segment 1 to 2

Teste de Breusch-Pagan

bptest(MQODin)

    studentized Breusch-Pagan test

data:  MQODin
BP = 0.68827, df = 1, p-value = 0.4068

Teste White

white_test(MQODin)

White's test results

Null hypothesis: Homoskedasticity of the residuals
Alternative hypothesis: Heteroskedasticity of the residuals
Test Statistic: 0.89
P-value: 0.639392

Teste ARCH

ArchTest(residuosMQO_FUT, lags = 4)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  residuosMQO_FUT
Chi-squared = 2.5216, df = 4, p-value = 0.6408

Testes de autocorrelação dos resíduos

Teste de Durbin-Watson

dwtest(MQODin)

    Durbin-Watson test

data:  MQODin
DW = 2.0632, p-value = 0.6409
alternative hypothesis: true autocorrelation is greater than 0

Teste de Breusch-Godfrey

bgtest(MQODin, order = 4)

    Breusch-Godfrey test for serial correlation of order up to 4

data:  MQODin
LM test = 7.1619, df = 4, p-value = 0.1276

Teste Box-Pierce

Box.test(residuosMQO_FUT, lag = 12, type = "Box-Pierce")

    Box-Pierce test

data:  residuosMQO_FUT
X-squared = 15.316, df = 12, p-value = 0.2246

Função de Autocorrelação dos resíduos com 36 defasagens

par(mfrow = c(1, 1))
acf(residuosMQO_FUT, lag = 36, col = "red", lwd = 2, 
    main = "ACF com 36 Defasagem", xlab = "Defasagem")

ESTIMAÇÃO GARCH / RETORNO E VARIÂNCIA

Estimação do modelo GARCH(1,1) com uma estrutura ARMA(1,1) na média

Garch11 <- garchFit(~ arma(1,1) + garch(1,1),
                    data = diff(log(cafe_fut_R)), 
                    include.mean = FALSE, 
                    trace = FALSE)


summary(Garch11)

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(1, 1) + garch(1, 1), data = diff(log(cafe_fut_R)), 
    include.mean = FALSE, trace = FALSE) 

Mean and Variance Equation:
 data ~ arma(1, 1) + garch(1, 1)
<environment: 0x0000017f73087488>
 [data = diff(log(cafe_fut_R))]

Conditional Distribution:
 norm 

Coefficient(s):
         ar1           ma1         omega        alpha1         beta1  
 0.113633405  -0.276333484   0.000038138   0.000000010   0.994522256  

Std. Errors:
 based on Hessian 

Error Analysis:
          Estimate  Std. Error  t value            Pr(>|t|)    
ar1     0.11363340  0.18431627    0.617               0.538    
ma1    -0.27633348  0.18724562   -1.476               0.140    
omega   0.00003814  0.00020673    0.184               0.854    
alpha1  0.00000001  0.05464259    0.000               1.000    
beta1   0.99452226  0.03343175   29.748 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
 138.8969    normalized:  1.102356 

Description:
 Tue Aug 20 13:41:54 2024 by user: josue 


Standardised Residuals Tests:
                                  Statistic   p-Value
 Jarque-Bera Test   R    Chi^2   0.01392616 0.9930611
 Shapiro-Wilk Test  R    W       0.99205198 0.6947955
 Ljung-Box Test     R    Q(10)  14.30175282 0.1596673
 Ljung-Box Test     R    Q(15)  19.80543025 0.1795244
 Ljung-Box Test     R    Q(20)  24.68802830 0.2136279
 Ljung-Box Test     R^2  Q(10)  10.21624187 0.4217310
 Ljung-Box Test     R^2  Q(15)  13.61402216 0.5549807
 Ljung-Box Test     R^2  Q(20)  17.35827965 0.6296015
 LM Arch Test       R    TR^2   10.44449670 0.5770260

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-2.125347 -2.012796 -2.128339 -2.079621 

round(Garch11@fit$matcoef, 6)

        Estimate  Std. Error   t value Pr(>|t|)
ar1     0.113633    0.184316  0.616513 0.537556
ma1    -0.276333    0.187246 -1.475781 0.140003
omega   0.000038    0.000207  0.184484 0.853634
alpha1  0.000000    0.054643  0.000000 1.000000
beta1   0.994522    0.033432 29.747840 0.000000

vconds1 = Garch11@h.t

Gráfico da série temporal da variância condicional GARCH(1,1)

plot.ts(vconds1, type = "l",
        main = "Variância Condicional GARCH(1,1)",
        xlab = "",
        ylab = "Variância Condicional")

plot(Garch11, which = 3)

plot.ts(vconds1)

Cálculo do desvio-padrão condicional

DP <- sqrt(vconds1)

Cálculo da volatilidade anualizada assumindo 252 dias úteis no ano

vol <- DP * sqrt(252)

plot.ts(vol)

Cálculo da média da volatilidade anualizada

mean(vol)

[1] 1.289633

Estimação do modelo MQO Dinâmico - CAFÉ SPOT

dcafe_spot_R <- diff(log(cafe_spot_R))
cbind(cafe_spot_R, dcafe_spot_R)

Regressão Dinâmica usando a série diferenciada (MQO Dinâmico)

lag <- stats::lag
MQODin <- dynlm(dcafe_spot_R ~ lag(dcafe_spot_R, -1))
coef(MQODin)

          (Intercept) lag(dcafe_spot_R, -1) 
          -0.00182699            0.18715316 

summary(MQODin)

Time series regression with "ts" data:
Start = 2014(3), End = 2024(7)

Call:
dynlm(formula = dcafe_spot_R ~ lag(dcafe_spot_R, -1))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.151371 -0.045392 -0.000234  0.040528  0.154448 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)           -0.001827   0.005873  -0.311   0.7562  
lag(dcafe_spot_R, -1)  0.187153   0.084348   2.219   0.0283 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06566 on 123 degrees of freedom
Multiple R-squared:  0.03849,   Adjusted R-squared:  0.03067 
F-statistic: 4.923 on 1 and 123 DF,  p-value: 0.02833

Resíduos da regressão

residuosMQO_SPOT <- residuals(MQODin)
residuosMQO_SPOT
chart.TimeSeries(residuosMQO_SPOT)

Histograma

chart.Histogram(residuosMQO_SPOT)

Gráfico de histograma com curva de distribuição normal sobreposta

par(mfrow = c(1, 1))
hist(residuosMQO_SPOT, main = "", col = "cadetblue", prob = TRUE, 
     xlab = names(residuosMQO_SPOT)[1], breaks = 30)
curve(expr = dnorm(x, mean = mean(residuosMQO_SPOT), sd = sd(residuosMQO_SPOT)), 
      col = "red", add = TRUE, lwd = 2)

Gráfico QQ-Norm (Normal Quantile-Quantile Plot) para verificar a normalidade dos resíduos

qqnorm(residuosMQO_SPOT, col = "blue")
qqline(residuosMQO_SPOT, col = "red")

Testes de Normalidade

## Teste de Jarque-Bera
jarqueberaTest(residuosMQO_SPOT)

Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 0.8121
  P VALUE:
    Asymptotic p Value: 0.6663 

Teste de Shapiro-Wilk

shapiro.test(residuosMQO_SPOT)

    Shapiro-Wilk normality test

data:  residuosMQO_SPOT
W = 0.99413, p-value = 0.8862

Testes de Heterocedasticidade

length(residuosMQO_SPOT)

[1] 125

length(residuosMQO_SPOT) * 0.15

[1] 18.75

Teste de Goldfeld-Quandt

gqtest(MQODin, fraction = 18.75, alternative = "greater")

    Goldfeld-Quandt test

data:  MQODin
GQ = 1.4353, df1 = 52, df2 = 51, p-value = 0.09935
alternative hypothesis: variance increases from segment 1 to 2

Teste de Breusch-Pagan

bptest(MQODin)

    studentized Breusch-Pagan test

data:  MQODin
BP = 0.40405, df = 1, p-value = 0.525

Teste White

white_test(MQODin)

White's test results

Null hypothesis: Homoskedasticity of the residuals
Alternative hypothesis: Heteroskedasticity of the residuals
Test Statistic: 6.21
P-value: 0.044873

Teste ARCH

ArchTest(residuosMQO_SPOT, lags = 4)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  residuosMQO_SPOT
Chi-squared = 1.1526, df = 4, p-value = 0.8858

Testes de autocorrelação dos resíduos

## Teste de Durbin-Watson
dwtest(MQODin)

    Durbin-Watson test

data:  MQODin
DW = 2.0792, p-value = 0.6646
alternative hypothesis: true autocorrelation is greater than 0

Teste de Breusch-Godfrey

bgtest(MQODin, order = 4)

    Breusch-Godfrey test for serial correlation of order up to 4

data:  MQODin
LM test = 9.0328, df = 4, p-value = 0.06029

Teste Box-Pierce

Box.test(residuosMQO_SPOT, lag = 12, type = "Box-Pierce")

    Box-Pierce test

data:  residuosMQO_SPOT
X-squared = 12.95, df = 12, p-value = 0.3727

Função de Autocorrelação dos resíduos com 36 defasagens

par(mfrow = c(1, 1))
acf(residuosMQO_SPOT, lag = 36, col = "red", lwd = 2, 
    main = "ACF com 36 Defasagem", xlab = "Defasagem")

ESTIMAÇÃO GARCH / RETORNO E VARIÂNCIA

Estimação do modelo GARCH(1,1) com uma estrutura ARMA(1,1) na média

Garch11 <- garchFit(~ arma(1,1) + garch(1,1),
                    data = diff(log(cafe_spot_R)), 
                    include.mean = FALSE, 
                    trace = FALSE)


summary(Garch11)

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(1, 1) + garch(1, 1), data = diff(log(cafe_spot_R)), 
    include.mean = FALSE, trace = FALSE) 

Mean and Variance Equation:
 data ~ arma(1, 1) + garch(1, 1)
<environment: 0x0000017f75b12ec0>
 [data = diff(log(cafe_spot_R))]

Conditional Distribution:
 norm 

Coefficient(s):
        ar1          ma1        omega       alpha1        beta1  
 0.50699728  -0.40141863   0.00036306   0.01756702   0.89411389  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value  Pr(>|t|)    
ar1     0.5069973   0.1412834    3.589  0.000333 ***
ma1    -0.4014186   0.1616030   -2.484  0.012992 *  
omega   0.0003631   0.0008688    0.418  0.676042    
alpha1  0.0175670   0.0515071    0.341  0.733058    
beta1   0.8941139   0.2228359    4.012 0.0000601 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
 167.7009    normalized:  1.33096 

Description:
 Tue Aug 20 13:41:55 2024 by user: josue 


Standardised Residuals Tests:
                                 Statistic   p-Value
 Jarque-Bera Test   R    Chi^2   0.5258746 0.7687901
 Shapiro-Wilk Test  R    W       0.9942114 0.8890764
 Ljung-Box Test     R    Q(10)   9.6155678 0.4748428
 Ljung-Box Test     R    Q(15)  18.6576292 0.2296697
 Ljung-Box Test     R    Q(20)  25.6214741 0.1786804
 Ljung-Box Test     R^2  Q(10)   3.6898209 0.9602568
 Ljung-Box Test     R^2  Q(15)   5.7709158 0.9833441
 Ljung-Box Test     R^2  Q(20)  11.4703920 0.9331003
 LM Arch Test       R    TR^2    5.7624978 0.9275777

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-2.582555 -2.470004 -2.585547 -2.536829 

round(Garch11@fit$matcoef, 6)

        Estimate  Std. Error   t value Pr(>|t|)
ar1     0.506997    0.141283  3.588512 0.000333
ma1    -0.401419    0.161603 -2.483980 0.012992
omega   0.000363    0.000869  0.417871 0.676042
alpha1  0.017567    0.051507  0.341060 0.733058
beta1   0.894114    0.222836  4.012431 0.000060

vconds1 = Garch11@h.t

Gráfico da série temporal da variância condicional GARCH(1,1)

plot.ts(vconds1, type = "l",
        main = "Variância Condicional GARCH(1,1)",
        xlab = "",
        ylab = "Variância Condicional")

plot(Garch11, which = 3)

plot.ts(vconds1)

Cálculo do desvio-padrão condicional

DP <- sqrt(vconds1)

Cálculo da volatilidade anualizada assumindo 252 dias úteis no ano

vol <- DP * sqrt(252)

plot.ts(vol)

Cálculo da média da volatilidade anualizada

mean(vol)

[1] 1.015833

Estimação do modelo MQO Dinâmico - CAFÉ ÍNDICE

dcafe_indice_R <- diff(log(cafe_indice_R))
cbind(cafe_indice_R, dcafe_indice_R)

Regressão Dinâmica usando a série diferenciada (MQO Dinâmico)

lag <- stats::lag
MQODin <- dynlm(dcafe_indice_R ~ lag(dcafe_indice_R, -1))
coef(MQODin)

            (Intercept) lag(dcafe_indice_R, -1) 
             -0.0034158               0.1701110 

summary(MQODin)

Time series regression with "ts" data:
Start = 2014(3), End = 2024(7)

Call:
dynlm(formula = dcafe_indice_R ~ lag(dcafe_indice_R, -1))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.12265 -0.01634  0.00230  0.01720  0.07520 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)  
(Intercept)             -0.003416   0.002652  -1.288   0.2002  
lag(dcafe_indice_R, -1)  0.170111   0.088079   1.931   0.0557 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02944 on 123 degrees of freedom
Multiple R-squared:  0.02943,   Adjusted R-squared:  0.02154 
F-statistic:  3.73 on 1 and 123 DF,  p-value: 0.05574

Resíduos da regressão

residuosMQO_IND <- residuals(MQODin)
residuosMQO_IND
chart.TimeSeries(residuosMQO_IND)

Histograma

chart.Histogram(residuosMQO_IND)

Gráfico de histograma com curva de distribuição normal sobreposta

par(mfrow = c(1, 1))
hist(residuosMQO_IND, main = "", col = "cadetblue", prob = TRUE, 
     xlab = names(residuosMQO_IND)[1], breaks = 30)
curve(expr = dnorm(x, mean = mean(residuosMQO_IND), sd = sd(residuosMQO_IND)), 
      col = "red", add = TRUE, lwd = 2)

Gráfico QQ-Norm (Normal Quantile-Quantile Plot) para verificar a normalidade dos resíduos

qqnorm(residuosMQO_IND, col = "blue")
qqline(residuosMQO_IND, col = "red")

Testes de Normalidade

Teste de Jarque-Bera

jarqueberaTest(residuosMQO_IND)

Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 28.2717
  P VALUE:
    Asymptotic p Value: 0.0000007259 

Teste de Shapiro-Wilk

shapiro.test(residuosMQO_IND)

    Shapiro-Wilk normality test

data:  residuosMQO_IND
W = 0.97111, p-value = 0.008818

Testes de Heterocedasticidade

length(residuosMQO_IND)

[1] 125

length(residuosMQO_IND) * 0.15

[1] 18.75

Teste de Goldfeld-Quandt

gqtest(MQODin, fraction = 18.75, alternative = "greater")

    Goldfeld-Quandt test

data:  MQODin
GQ = 3.2616, df1 = 52, df2 = 51, p-value = 0.0000207
alternative hypothesis: variance increases from segment 1 to 2

Teste de Breusch-Pagan

bptest(MQODin)

    studentized Breusch-Pagan test

data:  MQODin
BP = 12.089, df = 1, p-value = 0.0005073

Teste White

white_test(MQODin)

White's test results

Null hypothesis: Homoskedasticity of the residuals
Alternative hypothesis: Heteroskedasticity of the residuals
Test Statistic: 15.72
P-value: 0.000386

Teste ARCH

ArchTest(residuosMQO_IND, lags = 4)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  residuosMQO_IND
Chi-squared = 9.7829, df = 4, p-value = 0.04425

Testes de autocorrelação dos resíduos

Teste de Durbin-Watson

dwtest(MQODin)

    Durbin-Watson test

data:  MQODin
DW = 1.9874, p-value = 0.4651
alternative hypothesis: true autocorrelation is greater than 0

Teste de Breusch-Godfrey

bgtest(MQODin, order = 4)

    Breusch-Godfrey test for serial correlation of order up to 4

data:  MQODin
LM test = 6.6054, df = 4, p-value = 0.1583

Teste Box-Pierce

Box.test(residuosMQO_IND, lag = 12, type = "Box-Pierce")

    Box-Pierce test

data:  residuosMQO_IND
X-squared = 11.703, df = 12, p-value = 0.4698

Função de Autocorrelação dos resíduos com 36 defasagens

par(mfrow = c(1, 1))
acf(residuosMQO_IND, lag = 36, col = "red", lwd = 2, 
    main = "ACF com 36 Defasagem", xlab = "Defasagem")

ESTIMAÇÃO GARCH / RETORNO E VARIÂNCIA

Estimação do modelo GARCH(1,1) com uma estrutura ARMA(1,1) na média

Garch11 <- garchFit(~ arma(1,1) + garch(1,1),
                    data = diff(log(cafe_indice_R)), 
                    include.mean = FALSE, 
                    trace = FALSE)


summary(Garch11)

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(1, 1) + garch(1, 1), data = diff(log(cafe_indice_R)), 
    include.mean = FALSE, trace = FALSE) 

Mean and Variance Equation:
 data ~ arma(1, 1) + garch(1, 1)
<environment: 0x0000017f78c18e18>
 [data = diff(log(cafe_indice_R))]

Conditional Distribution:
 norm 

Coefficient(s):
        ar1          ma1        omega       alpha1        beta1  
 0.37063646  -0.24383438   0.00061679   0.27494167   0.00370377  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)   
ar1     0.3706365   0.2733318    1.356  0.17510   
ma1    -0.2438344   0.2935942   -0.831  0.40625   
omega   0.0006168   0.0002382    2.590  0.00961 **
alpha1  0.2749417   0.1385271    1.985  0.04717 * 
beta1   0.0037038   0.2785586    0.013  0.98939   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
 270.852    normalized:  2.149619 

Description:
 Tue Aug 20 13:41:56 2024 by user: josue 


Standardised Residuals Tests:
                                 Statistic     p-Value
 Jarque-Bera Test   R    Chi^2  11.1949842 0.003707149
 Shapiro-Wilk Test  R    W       0.9830267 0.116082807
 Ljung-Box Test     R    Q(10)  11.0045961 0.357160033
 Ljung-Box Test     R    Q(15)  15.4570635 0.419022312
 Ljung-Box Test     R    Q(20)  17.0139690 0.652066402
 Ljung-Box Test     R^2  Q(10)  16.3517463 0.089996290
 Ljung-Box Test     R^2  Q(15)  18.1087996 0.256988858
 Ljung-Box Test     R^2  Q(20)  21.1070652 0.390855045
 LM Arch Test       R    TR^2   14.1171359 0.293290838

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-4.219873 -4.107322 -4.222865 -4.174147 

round(Garch11@fit$matcoef, 6)

        Estimate  Std. Error   t value Pr(>|t|)
ar1     0.370636    0.273332  1.355995 0.175101
ma1    -0.243834    0.293594 -0.830515 0.406248
omega   0.000617    0.000238  2.589658 0.009607
alpha1  0.274942    0.138527  1.984749 0.047172
beta1   0.003704    0.278559  0.013296 0.989391

vconds1 = Garch11@h.t

Gráfico da série temporal da variância condicional GARCH(1,1)

plot.ts(vconds1, type = "l",
        main = "Variância Condicional GARCH(1,1)",
        xlab = "",
        ylab = "Variância Condicional")

plot(Garch11, which = 3)

plot.ts(vconds1)

Cálculo do desvio-padrão condicional

DP <- sqrt(vconds1)

Cálculo da volatilidade anualizada assumindo 252 dias úteis no ano

vol <- DP * sqrt(252)
plot.ts(vol)

Cálculo da média da volatilidade anualizada

mean(vol)

[1] 0.4549016