UNIVERSIDADE FEDERAL DA PARAÍBA
remove(list=ls())
ls()character(0)options(scipen=999)
options(max.print = 1000000)
par(mfrow=c(1,1))
setwd("C:/Users/henri/OneDrive/Área de Trabalho/GARCH")
library(tseries)
library(BatchGetSymbols)
library(readxl)
library(tidyverse)
library(ggthemes)
library(deflateBR)
library(knitr)
library(mFilter)
library(forecast)
library(rbcb)preco <- read_excel("BoiMes.xlsx")
# Indice de Commodities Agropecuário do BCB
Indice_AgroPec <- get_series(
27575,
start_date = "2014-01-01",
end_date = "2024-08-01")ultimo_valorf <- tail(preco, 1)
preco |>
ggplot(aes(x = Data, y = Valor_fut)) +
geom_line(size = 1, color = "#132157") + # Linha do gráfico
geom_point(data = ultimo_valorf, aes(x = Data, y = Valor_fut), size = 3, color = "red") + # Ponto no último valor
geom_text(data = ultimo_valorf, aes(x = Data, y = Valor_fut, label = round(Valor_fut, 2)),
vjust = -1, hjust = 0.6, color = "red", size = 3.5) + # Texto do último valor
labs(title = "Boi Gordo Preço Futuro Mensal",
subtitle = "2014/01/01 2024/08/01",
x = "",
y = "Preço",
caption = "Fonte: Investing.com") + # Rótulos e legendas
theme_economist() + # Tema minimalista
theme(plot.title = element_text(hjust = 0.5), # Centralizar título
plot.subtitle = element_text(hjust = 0.5))
retornos <- preco |>
mutate(log_preco = log(Valor_fut),
diff_log_preco = c(NA, diff(log_preco))) |>
select(Data, diff_log_preco) |>
na.omit()
#--------------------------------------#
ultimo_valorr <- tail(retornos, 1)
retornos |>
ggplot(aes(x = Data, y = diff_log_preco)) +
geom_line(size = 1, color = "#132157") + # Linha do gráfico
geom_point(data = ultimo_valorr, aes(x = Data, y = diff_log_preco), size = 3, color = "red") + # Ponto no último valor
geom_text(data = ultimo_valorr, aes(x = Data, y = diff_log_preco, label = round(diff_log_preco, 2)),
vjust = -1, hjust = 0.6, color = "red", size = 3.5) + # Texto do último valor
labs(title = "Boi Gordo Retorno do Preço Futuro Mensal",
subtitle = "2014/01/01 2024/08/01",
x = "",
y = "Preço",
caption = "Fonte: Investing.com") + # Rótulos e legendas
theme_economist() + # Tema minimalista
theme(plot.title = element_text(hjust = 0.5), # Centralizar título
plot.subtitle = element_text(hjust = 0.5))
ultimo_valor <- tail(preco, 1)
preco |>
ggplot(aes(x = Data, y = Valor)) +
geom_line(size = 1, color = "#132157") + # Linha do gráfico
geom_point(data = ultimo_valor, aes(x = Data, y = Valor), size = 3, color = "red") + # Ponto no último valor
geom_text(data = ultimo_valor, aes(x = Data, y = Valor, label = round(Valor, 2)),
vjust = -1, hjust = 0.6, color = "red", size = 3.5) + # Texto do último valor
labs(title = "Boi Gordo Preço Spot Mensal",
subtitle = "2014/01/01 2024/08/01",
x = "",
y = "Preço",
caption = "Fonte: Cepea") + # Rótulos e legendas
theme_economist() + # Tema minimalista
theme(plot.title = element_text(hjust = 0.5), # Centralizar título
plot.subtitle = element_text(hjust = 0.5))
retornosf <- preco |>
mutate(log_preco = log(Valor),
diff_log_preco = c(NA, diff(log_preco))) |>
select(Data, diff_log_preco) |>
na.omit()
#--------------------------------------#
ultimo_valorf <- tail(retornosf, 1)
retornosf |>
ggplot(aes(x = Data, y = diff_log_preco)) +
geom_line(size = 1, color = "#132157") + # Linha do gráfico
geom_point(data = ultimo_valorf, aes(x = Data, y = diff_log_preco), size = 3, color = "red") + # Ponto no último valor
geom_text(data = ultimo_valorf, aes(x = Data, y = diff_log_preco, label = round(diff_log_preco, 2)),
vjust = -1, hjust = 0.6, color = "red", size = 3.5) + # Texto do último valor
labs(title = "Boi Gordo Retorno do Preço Spot Mensal",
subtitle = "2014/01/01 2024/08/01",
x = "",
y = "Preço",
caption = "Fonte: Cepea") + # Rótulos e legendas
theme_economist() + # Tema minimalista
theme(plot.title = element_text(hjust = 0.5), # Centralizar título
plot.subtitle = element_text(hjust = 0.5))
times = seq(as.Date("2014/1/1"), by = "month", length.out = length(preco$Data))
Boif <- ts(preco$Valor_fut, frequency = 12, start = c(2014,01,01)) |>
deflate(nominal_dates = times, real_date = "07/2024", index = "igpm")
Bois <- ts(preco$Valor, frequency = 12, start = c(2014,01,01)) |>
deflate(nominal_dates = times, real_date = "07/2024", index = "igpm")
ts.plot(Boif, Bois, col = c("blue", "red"), lty = 1:2, gpars = list(xlab = "Ano", ylab = "Índice", main = "Spot x Futuro"))
legend("topright", legend = c("Boi Fut", "Boi Spot"), col = c("blue", "red"), lty = 1:2)
Boi_Hp <- hpfilter(Boif, freq = 14400)
plot(Boi_Hp$cycle, main = "Ciclo Boi Fut", col = "blue4");abline(h = 0, col ="grey")
Bois_Hp <- hpfilter(Bois, freq = 14400)
plot(Bois_Hp$cycle, main = "Ciclo Boi Spot", col = "red4");abline(h = 0, col ="grey")
ts.plot(Boi_Hp$cycle, Bois_Hp$cycle, col = c("blue4", "red4"), lty = 1:2, gpars = list(xlab = "Ano", ylab = "Índice", main = "Spot x Futuro"))
legend("topright", legend = c("Boi Fut", "Boi Spot"), col = c("blue4", "red4"), lty = 1:2);abline(h = 0, col = "grey")
ultimo_valor <- tail(Indice_AgroPec, 1)
Indice_AgroPec |>
ggplot(aes(x = date, y = `27575`)) +
geom_line(size = 1, color = "#132157") + # Linha do gráfico
geom_point(data = ultimo_valor, aes(x = date, y = `27575`), size = 3, color = "red") + # Ponto no último valor
geom_text(data = ultimo_valor, aes(x = date, y = `27575`, label = round(`27575`, 2)),
vjust = -0.5, hjust = 0.6, color = "red", size = 3.5) + # Texto do último valor
labs(title = "Índice de Commodities - Brasil - Agropecuária",
subtitle = "2014/01/01 2024/08/01",
x = "",
y = "Pontos",
caption = "Fonte: BCB") + # Rótulos e legendas
theme_economist() + # Tema minimalista
theme(plot.title = element_text(hjust = 0.5), # Centralizar título
plot.subtitle = element_text(hjust = 0.5))
df <- data.frame(Data = preco$Data, Boif, Bois)retornos <- df |>
mutate(
rboif = c(NA, diff(log(Boif))),
rbois = c(NA, diff(log(Bois)))) |>
select(Data, rbois, rboif) |>
na.omit()
# Calcular indicadores de desempenho
estatisticas <- function(retorno) {
media <- mean(retorno)
desvio <- sd(retorno)
list(media = media, desvio = desvio)
}
# Calcular indicadores para todas as empresas
Estatisticas <- lapply(retornos, estatisticas)
df2 <- data.frame(Boi_Spot = Estatisticas$rbois,
Boi_Fut = Estatisticas$rboif)
kable(df2)| Boi_Spot.media | Boi_Spot.desvio | Boi_Fut.media | Boi_Fut.desvio |
|---|---|---|---|
| -0.0003491 | 0.0416452 | -0.0001882 | 0.0544451 |
library(moments)
kurtosis(data.frame(Boi_Spot = retornos$rbois))Boi_Spot
7.444821 kurtosis(data.frame(Boi_Fut = retornos$rboif)) Boi_Fut
12.80386 skewness(data.frame(Boi_Spot = retornos$rbois)) Boi_Spot
0.7096975 skewness(data.frame(Boi_Fut = retornos$rboif)) Boi_Fut
1.676998 library(fBasics)
jarqueberaTest(retornos$rbois)
Title:
Jarque - Bera Normalality Test
Test Results:
STATISTIC:
X-squared: 115.2055
P VALUE:
Asymptotic p Value: < 0.00000000000000022 jarqueberaTest(retornos$rboif)
Title:
Jarque - Bera Normalality Test
Test Results:
STATISTIC:
X-squared: 568.14
P VALUE:
Asymptotic p Value: < 0.00000000000000022 #jarqueberaTest(diff(log(Boif)))acf(retornos$rbois)
acf(retornos$rboif)
pacf(retornos$rbois)
pacf(retornos$rboif)
Box.test((retornos$rbois), lag=1)
Box-Pierce test
data: (retornos$rbois)
X-squared = 1.1163, df = 1, p-value = 0.2907Box.test((retornos$rbois), lag=6)
Box-Pierce test
data: (retornos$rbois)
X-squared = 5.6411, df = 6, p-value = 0.4646Box.test((retornos$rbois), lag=10)
Box-Pierce test
data: (retornos$rbois)
X-squared = 9.7801, df = 10, p-value = 0.46Box.test((retornos$rbois), lag=24)
Box-Pierce test
data: (retornos$rbois)
X-squared = 16.255, df = 24, p-value = 0.8786Box.test((retornos$rboif), lag=1)
Box-Pierce test
data: (retornos$rboif)
X-squared = 1.158, df = 1, p-value = 0.2819Box.test((retornos$rboif), lag=6)
Box-Pierce test
data: (retornos$rboif)
X-squared = 10.101, df = 6, p-value = 0.1205Box.test((retornos$rboif), lag=10)
Box-Pierce test
data: (retornos$rboif)
X-squared = 12.519, df = 10, p-value = 0.2518Box.test((retornos$rboif), lag=24)
Box-Pierce test
data: (retornos$rboif)
X-squared = 26.482, df = 24, p-value = 0.3292Arima21<-arima(Boif,order=c(2,1,0),method=c("CSS"))
Arima21
Call:
arima(x = Boif, order = c(2, 1, 0), method = c("CSS"))
Coefficients:
ar1 ar2
-0.1247 -0.1863
s.e. 0.0869 0.0868
sigma^2 estimated as 233.8: log likelihood = -526.58, aic = NAResiduos<-residuals(Arima21)
par(mfrow=c(1,1))
#chart.TimeSeries(Residuos)
#chart.ACF(Residuos)
library(stats)
tsdiag(Arima21)
library(urca)
ddicky2 <- ur.df(Residuos, lags = 0, type="drift")
summary(ddicky2)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression drift
Call:
lm(formula = z.diff ~ z.lag.1 + 1)
Residuals:
Min 1Q Median 3Q Max
-43.004 -6.826 -1.618 5.246 97.667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.26476 1.35695 -0.195 0.846
z.lag.1 -0.99028 0.08948 -11.067 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.29 on 125 degrees of freedom
Multiple R-squared: 0.4949, Adjusted R-squared: 0.4909
F-statistic: 122.5 on 1 and 125 DF, p-value: < 0.00000000000000022
Value of test-statistic is: -11.0672 61.2421
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81library(fGarch)
Garch21<-garchFit(~ arma(2,1) + garch(2,1),
data=retornos$rbois,
include.mean = F,
trace=F)
summary(Garch21)
Title:
GARCH Modelling
Call:
garchFit(formula = ~arma(2, 1) + garch(2, 1), data = retornos$rbois,
include.mean = F, trace = F)
Mean and Variance Equation:
data ~ arma(2, 1) + garch(2, 1)
<environment: 0x00000195f2bfe2e8>
[data = retornos$rbois]
Conditional Distribution:
norm
Coefficient(s):
ar1 ar2 ma1 omega alpha1 alpha2
-0.15371412 -0.12317226 0.17467483 0.00046728 0.41683260 0.73741687
beta1
0.00000001
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
ar1 -0.15371412 0.20434082 -0.752 0.4519
ar2 -0.12317226 NaN NaN NaN
ma1 0.17467483 0.21198876 0.824 0.4100
omega 0.00046728 NaN NaN NaN
alpha1 0.41683260 0.25796804 1.616 0.1061
alpha2 0.73741687 0.38795916 1.901 0.0573 .
beta1 0.00000001 NaN NaN NaN
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log Likelihood:
235.6627 normalized: 1.855611
Description:
Sat Aug 10 12:53:58 2024 by user: henri
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 135.1792067 0.0000000000
Shapiro-Wilk Test R W 0.9401118 0.0000263034
Ljung-Box Test R Q(10) 14.8138911 0.1389958079
Ljung-Box Test R Q(15) 19.5074086 0.1916540735
Ljung-Box Test R Q(20) 21.5790772 0.3637806752
Ljung-Box Test R^2 Q(10) 1.9563733 0.9966637077
Ljung-Box Test R^2 Q(15) 5.9528857 0.9805283165
Ljung-Box Test R^2 Q(20) 6.8565256 0.9971297336
LM Arch Test R TR^2 1.9389920 0.9994937829
Information Criterion Statistics:
AIC BIC SIC HQIC
-3.600987 -3.444221 -3.606650 -3.537295 round(Garch21@fit$matcoef, 6) Estimate Std. Error t value Pr(>|t|)
ar1 -0.153714 0.204341 -0.752244 0.451904
ar2 -0.123172 NaN NaN NaN
ma1 0.174675 0.211989 0.823982 0.409950
omega 0.000467 NaN NaN NaN
alpha1 0.416833 0.257968 1.615830 0.106131
alpha2 0.737417 0.387959 1.900759 0.057334
beta1 0.000000 NaN NaN NaNvconds1=Garch21@h.t
plot.ts(vconds1,type="l",
main="Variância Condicional GARCH(2,1)",
xlab="",
ylab="Variância Condicional")
plot(Garch21, which = 3)
vconds1 [1] 0.0023751108 0.0023751108 0.0004672786 0.0004872066 0.0007419877
[6] 0.0008937121 0.0005667418 0.0012638324 0.0023397679 0.0026026211
[11] 0.0037246832 0.0037699861 0.0005042042 0.0004928726 0.0005579916
[16] 0.0006730990 0.0008699492 0.0008655283 0.0011525719 0.0015517798
[21] 0.0006718330 0.0006484979 0.0006906585 0.0007889026 0.0009564414
[26] 0.0006614274 0.0006289586 0.0004878901 0.0007272039 0.0008886243
[31] 0.0008107885 0.0014273870 0.0012553758 0.0005149393 0.0005256122
[36] 0.0005636765 0.0005435699 0.0009568424 0.0011814598 0.0016380409
[41] 0.0023254791 0.0016408716 0.0027310195 0.0030882062 0.0061792650
[46] 0.0037320126 0.0006476689 0.0006731994 0.0007658540 0.0005541621
[51] 0.0006243476 0.0006365045 0.0009841302 0.0012513061 0.0011419885
[56] 0.0005444358 0.0008976927 0.0011875993 0.0009845907 0.0011343316
[61] 0.0012583944 0.0008232330 0.0005641780 0.0006250624 0.0011882321
[66] 0.0016730019 0.0008020893 0.0005199840 0.0009396498 0.0016819270
[71] 0.0187269574 0.0325389291 0.0055788930 0.0066301463 0.0007645066
[76] 0.0005309023 0.0005765189 0.0010415436 0.0020632663 0.0015805421
[81] 0.0021405624 0.0033445772 0.0018043497 0.0066122638 0.0107507228
[86] 0.0052302414 0.0005239160 0.0005828707 0.0009518992 0.0014014637
[91] 0.0007372301 0.0006619261 0.0017565008 0.0071355108 0.0125059791
[96] 0.0078582228 0.0045155941 0.0028806765 0.0004984054 0.0014245002
[101] 0.0031881225 0.0025659060 0.0009991980 0.0011545298 0.0017843308
[106] 0.0009668645 0.0012988018 0.0020912919 0.0018017714 0.0012865809
[111] 0.0010897590 0.0013111924 0.0029811075 0.0051133213 0.0015829510
[116] 0.0073718376 0.0125421631 0.0046924692 0.0070945749 0.0025923259
[121] 0.0036628884 0.0012564982 0.0018753589 0.0007256978 0.0007157081
[126] 0.0013162919 0.0018011972plot.ts(vconds1)
DP<-sqrt(vconds1)
vol <- DP*sqrt(252)
vol [1] 0.7736459 0.7736459 0.3431533 0.3503942 0.4324129 0.4745687 0.3779139
[8] 0.5643454 0.7678682 0.8098522 0.9688241 0.9746982 0.3564540 0.3524257
[15] 0.3749852 0.4118506 0.4682171 0.4670258 0.5389324 0.6253387 0.4114631
[22] 0.4042542 0.4171881 0.4458738 0.4909412 0.4082643 0.3981175 0.3506398
[29] 0.4280834 0.4732159 0.4520163 0.5997512 0.5624542 0.3602287 0.3639427
[36] 0.3768905 0.3701076 0.4910441 0.5456445 0.6424845 0.7655199 0.6430394
[43] 0.8295884 0.8821723 1.2478681 0.9697769 0.4039957 0.4118813 0.4393122
[50] 0.3736962 0.3966555 0.4004986 0.4979968 0.5615418 0.5364523 0.3704022
[57] 0.4756244 0.5470604 0.4981133 0.5346509 0.5631300 0.4554720 0.3770582
[64] 0.3968825 0.5472061 0.6493046 0.4495848 0.3619889 0.4866125 0.6510343
[71] 2.1723704 2.8635311 1.1856985 1.2925931 0.4389256 0.3657696 0.3811598
[78] 0.5123173 0.7210708 0.6311074 0.7344533 0.9180596 0.6743116 1.2908487
[85] 1.6459593 1.1480509 0.3633550 0.3832537 0.4897740 0.5942801 0.4310243
[92] 0.4084181 0.6653106 1.3409507 1.7752484 1.4072214 1.0667379 0.8520155
[99] 0.3543983 0.5991444 0.8963297 0.8041196 0.5017947 0.5393899 0.6705605
[106] 0.4936090 0.5720997 0.7259515 0.6738296 0.5694018 0.5240413 0.5748221
[113] 0.8667405 1.1351462 0.6315882 1.3629758 1.7778147 1.0874292 1.3370987
[120] 0.8082488 0.9607538 0.5627055 0.6874521 0.4276399 0.4246863 0.5759389
[127] 0.6737223plot.ts(ts(vol, freq = 12, start = c(2014,01,01)),
ylab = "Vol")
O parâmetro Beta do Modelo GARCH representa a persistência da volatilidade ao longo do tempo, ou seja, o tempo que a série leva para retornar ao seu estado de equilíbrio após um choque. Além disso, pode ser utilizado como um indicador da eficiência do mercado. Quanto maior o Beta, maior o tempo de recuperação pós choque.
Por outro lado, o parâmetro Alfa reflete a reação imediata da volatilidade a um choque, indicando quão acentuada será essa volatilidade logo após o evento. Quanto maior o alfa, mais pontiagudo será.
Quando a soma dos parâmetros Alfa e Beta se aproxima de 1, isso sugere que o tempo necessário para a dissipação do efeito do choque será maior.
Beta = 0.00000001
Alfa = 0.41683260
mean(vol)*100[1] 67.3201seasonplot(Boif,
col = rainbow(12),
year.labels = T,
type = "o",
pch = 16)
library(seasonal)
# Análise sazonal para preços spot
modelo_spot <- seas(Bois)
summary(modelo_spot)
Call:
seas(x = Bois)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
LS2019.Nov 0.17781 0.02486 7.151 0.000000000000859 ***
LS2020.Jan -0.12874 0.02521 -5.106 0.000000328907105 ***
AO2020.Dec -0.09707 0.01560 -6.223 0.000000000488334 ***
AO2021.Oct -0.08478 0.01535 -5.522 0.000000033505226 ***
LS2023.Aug -0.11402 0.02490 -4.580 0.000004660359967 ***
LS2023.Oct 0.12595 0.02547 4.945 0.000000759966103 ***
AR-Seasonal-12 0.22581 0.09081 2.487 0.012892 *
MA-Nonseasonal-01 -0.30133 0.08456 -3.563 0.000366 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SEATS adj. ARIMA: (0 1 1)(1 0 0) Obs.: 128 Transform: log
AICc: 885, BIC: 909 QS (no seasonality in final): 0
Box-Ljung (no autocorr.): 11.1 Shapiro (normality): 0.9902 plot(modelo_spot)
# Decomposição usando STL para preços spot
stl_decomposicao_spot <- stl(Bois, s.window = "periodic")
plot(stl_decomposicao_spot)
seasonplot(Bois,
col = rainbow(12),
year.labels = T,
type = "o",
pch = 16)
# Análise sazonal para preços spot
modelo_fut <- seas(Boif)
summary(modelo_fut)
Call:
seas(x = Boif)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
LS2017.Aug 0.14265 0.02938 4.855 0.00000120526850633 ***
LS2019.Nov 0.32557 0.02939 11.076 < 0.0000000000000002 ***
AO2020.Dec -0.09610 0.01813 -5.302 0.00000011466976606 ***
AO2021.Oct -0.15895 0.01812 -8.772 < 0.0000000000000002 ***
AO2023.Aug -0.14069 0.01818 -7.740 0.00000000000000993 ***
AR-Seasonal-12 0.10099 0.09158 1.103 0.270104
MA-Nonseasonal-01 -0.31546 0.08403 -3.754 0.000174 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SEATS adj. ARIMA: (0 1 1)(1 0 0) Obs.: 128 Transform: log
AICc: 923.2, BIC: 944.7 QS (no seasonality in final): 0
Box-Ljung (no autocorr.): 17.51 Shapiro (normality): 0.9831
Messages generated by X-13:
Warnings:
- Series should not be a candidate for seasonal adjustment
because the spectrum of the prior adjusted series (Table B1)
has no visually significant seasonal peaks.plot(modelo_fut)
# Decomposição usando STL para preços spot
stl_decomposicao_fut <- stl(Boif, s.window = "periodic")
plot(stl_decomposicao_fut)