O agronegócio representa uma boa porcentagem do PIB brasileiro, por ser um mercado competitivo e com valores voláteis, a maioria dos dados não são públicos. Este estudo utiliza as poucas bases públicas referentes ao comércio de carne bovina, por serem escassas, os dados são providos de estados diferentes. Abaixo é apresentado os dados e suas fontes:

Bibliotecas

library(tseries)
library(urca)
library(readxl)
library(vars)
library(tsDyn)
library(MTS)
library(MVN)

Tratamento dos dados

data <- read_excel("D:/Estatística/Bases de cursos/Frigorífico/Preço de compra e venda Boi Gordo mensal 2014-2025.xlsx", 
                    skip = 4, sheet = "Planilha4")

tempo  <- data$Mês_Ano
carne  <- ts(data$carne_ponderada , start = c(2014), end = c(2025, 9), frequency = 12)
boi    <- ts(data$quilo_boi_gordo , start = c(2014), end = c(2025, 9), frequency = 12)
volume <- ts(data$Volume          , start = c(2014), end = c(2025, 9), frequency = 12)
dolar  <- ts(data$Dólar           , start = c(2014), end = c(2025, 9), frequency = 12)
safra  <- ts(data$D_Safra         , start = c(2014), end = c(2025, 9), frequency = 12)

A série do preço da carne no atacado foi calculada usando uma média ponderada do preço da carne no atacado da parte dianteira e traseira do boi. A conta leva em conta os valores não contabilizados da agulha bovina e das vendas de subprodutos do boi. Abaixo está a fórmula improvisada para correção da falta de informações:

\[ carne = 1.15 \ \cdot \ (0.48 \ \cdot traseira + 0.52 \ \cdot dianteira) \]

Este processo foi realizado no excel.

log_carne  <- log(carne)
log_boi    <- log(boi)
log_dolar  <- log(dolar)
log_volume <- log(volume)

Os dados possuem escalas diferentes, como a estimação do modelos é feita através de uma matriz variância-covariância, para não enviesar o modelo é preciso aplicar uma transformação nos dados. Portanto, as variáveis foram transformadas em logaritmos naturais para que os coeficientes estimados representem elasticidades, permitindo a análise de variações percentuais e estabilizando a variância da série de abate.

endogenas <- cbind(log_boi, log_carne, log_volume)
exogenas  <- cbind(log_dolar, safra)

O modelo diferencia as variáveis em endógenas e exógenas. As endógenas possuem influência mútua: o preço do boi, o valor da carne no atacado e o volume de abates interagem dinamicamente entre si. Já as exógenas exercem influência unidirecional sobre o sistema. O dólar, por exemplo, impacta os preços da cadeia, mas não é alterado por eles. Da mesma forma, os períodos de safra e entressafra determinam a oferta e o preço do boi gordo devido às condições das pastagens, porém o comportamento do mercado pecuário não interfere na ocorrência desses ciclos climáticos

O período de safra e entressafra é sazonal, mas não exercem o mesmo impacto todos os anos. As mudanças climáticas afetaram o período de chuvas e o intervalo entre elas, porém ainda a chuva ainda impacta de certa forma o preço dos bois. Os meses de safra e entressafra são:

  • Safra: Novembro - Dezembro - Janeiro - Fevereiro - Março - Abril

  • Entressafra: Maio - Junho - Julho - Agosto - Setembro - Outubro

Visualização das séries

plot.ts(log_carne, main = "Preço médio de venda bovina Atacado",
        ylab = "R$ - Real", xlab = "Tempo")

Nota-se que a série é não estacionária, ela varia muito ao decorrer do tempo e muda de média dependendo do intervalo de tempo utilizado. O valor da carne possui um grande aumento no ano de 2020, este que ocorreu o início da pandemia de Covid-19. Durante dois anos os preços se estabilizaram e houve uma queda, mas seguida de um novo aumento. Os fatores que levaram à dinâmica destes preços, principamente foi a pandemia e o ciclo pecuário de retenção, onde o bezerro ficou super valorizado e as vacas foram retidas para procriação, diminuindo o volume de abates.

plot.ts(log_boi, main = "Preço médio do quilo do Boi Gordo", ylab = "R$ - Real",
        xlab = "Tempo")

Observa-se que o comportamento do preço do boi gordo é muito semelhante ao preço da carne, pois a mudança de preço pago pelo frigorífico é repassado diretamente ao valor de venda da carne e seus subprodutos. Isto nada mais é que a dinâmica do comércio, logo as mesmas questões abordadas anteriormente explica também esta série.

plot.ts(log_volume, main = "Capacidade do volume de bovinos abatidos", ylab = "Quantidade de abates",
        xlab = "Tempo")

A série do volume de abates é bem comportada, porém no ano de 2022 ela sofre um grande crescimento afetando a média e variância neste período. O fator principal deste comportamento, foi a fase de retenção em 2020, deste modo, as vacas que foram retidas naquele momento começaram a ser abatidas, consequentemente aumentando o volume de abate.

plot.ts(log_dolar, main = "Cotação do dólar comercial", ylab = "R$ - Real",
        xlab = "Tempo")

O dólar é a principal moeda do comércio internacional, sendo afetada por diversos fatores, por isso que ela varia bastante ao decorrer dos anos. Em 2018 houve uma alta do dólar, devido a uma desvalorização do real (greve dos caminhoneiros e política conturbada acabou afetando a “credibilidade” da moeda) e a alta dos juros no Banco Central dos EUA. Em 2020 houve outra alta, mas essa foi ocasionada principalmente pela pandemia, afetando globalmente.

Divisão da base

Diferente dos modelos convencionais de machine learning, a validação em séries temporais deve respeitar a ordem cronológica dos dados para evitar o vazamento de informações do futuro para o passado. Neste estudo, a divisão entre treino e teste foi definida com base na volatilidade recente das séries, reservando-se os últimos 18 meses para a base de teste (avaliação out-of-sample). Para a estimação, foram consideradas duas abordagens de validação cruzada:

A técnica de janela expansiva foi a escolhida para este projeto. Essa abordagem é superior para lidar com a alta volatilidade e as mudanças de regime observadas nos últimos anos (como a transição para a fase de liquidação no ciclo pecuário), pois permite que o modelo VECM reajuste seus coeficientes de correção de erros à medida que novos dados reais são disponibilizados.

n <- nrow(endogenas)
n_teste <- 18  # Separando os últimos 18 meses
n_treino <- n - n_teste

# Treino
endog_treino <- endogenas[1:n_treino, ]
exog_treino  <- exogenas[1:n_treino, ]

# Teste 
endog_teste <- endogenas[(n_treino + 1):n, ]
exog_teste  <- exogenas[(n_treino + 1):n, ]

Verificação dos pressupostos do VEC

Os principais pressupostos do VEC é a não estacionaridade das séries e a relação de longo prazo entre as séries (cointegração). Deste modo, a primeira análise a ser realizada é o teste de estacionaridade das séries, se alguma delas for estacionária, ela precisa ser retirada do modelo. Após a verificação da estacionaridade das séries é necessário identificar o lag correto, ou seja, é necessário encontrar quantos tempos passados o modelo necessita para explicar o presente e prever melhor o futuro. O lag será abordado melhor mais para frente.

Teste de estacionaridade

O teste clássico de estacionaridade é o Augmented Dick-Fuller, onde ele testa se a variação de um \(y_t\) depende do valor de \(y_{t-1}\), ou seja, se o valor da série hoje não depende do valor de ontem a série irá caminhar sem um rumo, se o valor da série hoje depende do valor de ontem a série está caminhando sobre uma média.

Como em qualquer teste estatístico, existe uma hipótese nula e alternativa, sendo neste caso:

\[ \left\{\begin{matrix} H_0: A \ série\ não \ é \ estacionária.\\ H_a: A \ série \ é \ estacionária. \end{matrix}\right. \] , ou seja, se \(p_{valor} \ \leq\) 0.05 rejeita-se \(H_0\) em prol da \(H_a\), assumindo com nível de confiança de 95% que a série é estacionária. Caso \(p_{valor} \ \geq\) 0.05, de modo analógo assume-se que a série é não estacionária.

boi_treino <- endog_treino[,1]
adf.test(boi_treino)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  boi_treino
## Dickey-Fuller = -1.1232, Lag order = 4, p-value = 0.9154
## alternative hypothesis: stationary
carne_treino <- endog_treino[,2]
adf.test(carne_treino)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  carne_treino
## Dickey-Fuller = -1.2377, Lag order = 4, p-value = 0.8932
## alternative hypothesis: stationary
volume_treino <- endog_treino[,3]
adf.test(volume_treino)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  volume_treino
## Dickey-Fuller = -1.6495, Lag order = 4, p-value = 0.7222
## alternative hypothesis: stationary
dolar_treino <- exog_treino[,1]
adf.test(dolar_treino)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dolar_treino
## Dickey-Fuller = -2.156, Lag order = 4, p-value = 0.5118
## alternative hypothesis: stationary

Como mostra os testes, as séries utilizadas na construção do modelo são todas não estacionárias.

Seleção de lags

Com as séries validadas, o próximo passo é encontrar valores de lags onde exista cointegração. Uma função para encontrar um chute inicial é a VARselect, ela calcula modelos com diferentes lags, depois retorna os valores de lags com as melhores métricas calculadas, geralmente se leva mais em conta o AIC.

lag_select <- VARselect(endog_treino, lag.max = 12, type = "const")
print(lag_select$selection)
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      2      2      3

A função retornou o lag 3, sendo assim testaremos esse valor para verificar se existe ou não cointegração nas séries.

Cointegração

O teste multivariado de cointegração clássico é o teste de Johansen, ele funciona como uma espécie de “PCA do equilíbrio”, utilizando a decomposição de autovalores (eigenvalues) para identificar quantas relações estáveis de longo prazo existem entre variáveis não estacionárias, como o preço do boi e da carne. Em vez de focar na variância, como o PCA tradicional, ele extrai vetores que representam combinações lineares estacionárias, testando sucessivamente se a força dessas conexões (os autovalores) é estatisticamente significativa em relação a valores críticos de 1%, 5% ou 10%.

Dessa forma, o teste determina o “posto” (rank) da matriz de impacto do modelo, definindo se as séries vagam de forma independente ou se estão presas por uma “coleira” econômica invisível. Se o teste rejeitar a hipótese de zero vetores (\(r=0\)), ele valida matematicamente que o spread de abate não é uma coincidência estatística, mas sim um fundamento sólido que permite a construção segura de um modelo VEC para previsões de margem.

teste_johansen <- ca.jo(endog_treino, ecdet = "const", type = "trace",
                        K = 3, spec = "longrun",
                        dumvar = exog_treino) 

summary(teste_johansen)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 1.946509e-01 8.748095e-02 2.850361e-02 4.168827e-18
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 2 |  3.47  7.52  9.24 12.97
## r <= 1 | 14.46 17.85 19.96 24.60
## r = 0  | 40.43 32.00 34.91 41.07
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               log_boi.l3 log_carne.l3 log_volume.l3    constant
## log_boi.l3     1.0000000   1.00000000      1.000000  1.00000000
## log_carne.l3  -1.1592901  -0.15282972     -1.278463 -0.98650961
## log_volume.l3  0.2256127   0.05162147     -1.360908  0.03905656
## constant      -1.4081741  -1.07829528     11.420980 -1.13438080
## 
## Weights W:
## (This is the loading matrix)
## 
##              log_boi.l3 log_carne.l3 log_volume.l3      constant
## log_boi.d     0.1249808  -0.07869571    0.01077478 -5.645598e-15
## log_carne.d   0.3151658  -0.03217141    0.01041388  6.469768e-16
## log_volume.d -0.3992854   0.01548744    0.09182137 -8.419831e-14

Como a estatística do teste é maior que o valor crítico de 5%, pode-se assumir com um nível de confiança de 95% que existe ao menos um vetor de cointegração, neste caso existe apenas um pois o teste falha para mais de um vetor.

Portanto, os pressupostos iniciais para construção do modelo VEC foram atendidas e pode-se seguir com a contrução usando lag = 2 (K-1) e r = 1.

Construção do modelo

modelo_treino <- VECM(endog_treino, 
                      lag = 2, r = 1, 
                      include = "const", 
                      exogen = exog_treino, 
                      estim = "ML")

summary(modelo_treino)
## #############
## ###Model VECM 
## #############
## Full sample size: 123    End sample size: 120
## Number of variables: 3   Number of estimated slope parameters 30
## AIC -2368.983    BIC -2279.783   SSR 0.7503165
## Cointegrating vector (estimated by ML):
##    log_boi log_carne log_volume
## r1       1 -1.161495  0.2263903
## 
## 
##                     ECT                Intercept           log_boi -1        
## Equation log_boi    0.1347(0.1280)     -0.1990(0.1844)     0.3126(0.1828).   
## Equation log_carne  0.3183(0.0823)***  -0.4527(0.1186)***  0.3187(0.1175)**  
## Equation log_volume -0.4027(0.2058).   0.5734(0.2964).     0.3469(0.2938)    
##                     log_carne -1        log_volume -1       log_boi -2         
## Equation log_boi    -0.2560(0.2413)     0.0526(0.0570)      -0.3308(0.1752).   
## Equation log_carne  -0.1096(0.1551)     -0.0584(0.0367)     -0.3216(0.1126)**  
## Equation log_volume -0.1657(0.3878)     -0.6142(0.0917)***  -0.1202(0.2816)    
##                     log_carne -2        log_volume -2       log_dolar          
## Equation log_boi    0.4493(0.1811)*     -0.0196(0.0607)     0.0265(0.0259)     
## Equation log_carne  0.2534(0.1164)*     -0.0148(0.0391)     0.0538(0.0167)**   
## Equation log_volume -0.2370(0.2910)     -0.2692(0.0976)**   -0.0518(0.0416)    
##                     safra             
## Equation log_boi    0.0105(0.0088)    
## Equation log_carne  0.0051(0.0056)    
## Equation log_volume -0.0321(0.0141)*

Com o modelo construído será analisado os pressupostos necessários para que as predições do modelo sejam válidas e consistentes. Para as predições do VEC ter um aval estatístico, os resíduso do modelo necessitam ser ruído branco, ou seja, precisam ser independentes (sem autocorrelação), ter distribuição normal e ser homocedástica (variância controlada). A seguir será avaliado todos esses tópicos, se algum teste falhar é necessário escolher outro valor de K e identificar o melhor r até que o modelo passe por todos os testes.

res_VECM <- residuals(modelo_treino)

modelo_cajo <- ca.jo(endog_treino, ecdet = "const", type = "trace",
                        K = 3, spec = "longrun",
                        dumvar = exog_treino) 

modelo_var <- vec2var(teste_johansen, r = 1)

Para a construção dos testes é preciso colocar as variáveis em formatos que as funções implementadas nas bibliotecas aceitem.

Autocorrelação dos resíduos

O teste de autocorrelação dos resíduos é fundamental para garantir que o modelo produza erros aleatórios. Quando os resíduos não são aleatórios, significa que o modelo foi mal especificado e deixou de capturar informações ou comportamentos relevantes das séries temporais. Para verificar a autocorrelação multivariada, utiliza-se o teste de Portmanteau, que avalia se os resíduos no instante t são influenciados por seus valores passados (t - n) por meio de suas autocovariâncias. As hipóteses do teste são:

\[ \left\{\begin{matrix} H_0: Independência\ dos\ resíduos\ (ausência\ de\ autocorrelação).\\ H_a: Dependência\ dos\ resíduos\ (presença\ de\ autocorrelação). \end{matrix}\right. \]

serial.test(modelo_var, lags.pt = 7, type = "PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 39.79, df = 39, p-value = 0.4347
serial.test(modelo_var, lags.pt = 12, type = "PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 100.35, df = 84, p-value = 0.1078
serial.test(modelo_var, lags.pt = 15, type = "PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 127.09, df = 111, p-value = 0.141
serial.test(modelo_var, lags.pt = 23, type = "PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 201.42, df = 183, p-value = 0.1668

O teste foi aplicado para diferentes defasagens e, em todas, o \(p_{valor}\) foi superior a 0,05. Desse modo, não se rejeita \(H_0\), assumindo-se que os resíduos são independentes com um nível de significância de 5% (ou nível de confiança de 95%).

Normalidade dos resíduos

O teste de normalidade usualmente utilizado para dados que não são séries temporais é o de Shapiro-Wilk, enquanto em séries univariadas utiliza-se o de Jarque-Bera. No entanto, quando trabalhamos com múltiplas séries, é necessário verificar a normalidade multivariada. Isso ocorre porque séries que são individualmente normais podem não apresentar normalidade quando analisadas em conjunto.

O teste clássico para essa finalidade é o Jarque-Bera multivariado. Ele segue os princípios da versão univariada, avaliando a assimetria ( skewness ) e a curtose ( kurtosis ) dos resíduos do modelo para compará-las com os parâmetros de uma distribuição normal. As hipóteses do teste são:

\[ \left\{\begin{matrix} H_0: Os\ resíduos\ seguem\ uma\ distribuição\ normal\ (assimetria\ e\ curtose\ condizentes\ com\ a\ normalidade).\\ H_a: Os\ resíduos\ não\ seguem\ uma\ distribuição\ normal. \end{matrix}\right. \]

# Testa assimetria e curtose multivariada
norm_test <- normality.test(modelo_var)
print(norm_test)
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 110.61, df = 6, p-value < 2.2e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 25.116, df = 3, p-value = 1.46e-05
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 85.495, df = 3, p-value < 2.2e-16
## $jb.mul
## $jb.mul$JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 110.61, df = 6, p-value < 2.2e-16
## 
## 
## $jb.mul$Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 25.116, df = 3, p-value = 1.46e-05
## 
## 
## $jb.mul$Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 85.495, df = 3, p-value < 2.2e-16

Nos vários testes de assimetria e curtose realizados, em todos o \(p_{valor}\) foi inferior a 0,05. Desse modo, rejeita-se \(H_0\) com um nível de confiança de 95% e assume-se que os resíduos não possuem uma distribuição normal. Ao final desta seção, as implicações desse resultado serão discutidas com maior detalhamento.

Homocesticidade dos resíduos

Para verificar a estabilidade da variância dos erros ao longo do tempo, utiliza-se o teste de efeito ARCH (Autoregressive Conditional Heteroskedasticity). Em modelos multivariados, esse teste avalia se a volatilidade (variância) dos resíduos em um período é influenciada pela volatilidade de períodos anteriores. A presença de efeitos ARCH indicaria que os erros aparecem em ‘agrupamentos’ (períodos de alta e baixa variação), o que invalidaria os erros-padrão das estimativas. As hipóteses do teste são:

\[ \left\{\begin{matrix} H_0: Os \ resíduos \ são \ homocedásticos \ (a \ variância \ é \ constante \ e \ não \ há \ efeitos \ ARCH).\\ H_a: Os\ resíduos\ são\ heterocedásticos\ (existe\ o\ efeito\ ARCH,\ ou \ seja,\ a \ variância \ é \ dependente \ do \ tempo). \end{matrix}\right. \]

arch_test <- arch.test(modelo_var, lags.multi = 5)
arch_test
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 170.42, df = 180, p-value = 0.6839

A interpretação segue a regra usual: um \(p_{valor}\) indica que não há evidências de heterocedasticidade, confirmando que o modelo capturou adequadamente a variabilidade dos dados.

Estabilidade da série

Para garantir que os parâmetros estimados pelo modelo permaneçam constantes durante todo o período analisado, utiliza-se o teste de estabilidade CUSUM (Soma Cumulativa). Esse teste é fundamental para identificar se houve alguma quebra estrutural nos dados — como uma crise econômica ou mudança de política — que tenha alterado a relação entre as variáveis. O teste monitora a soma acumulada dos resíduos recursivos e verifica se essa soma ultrapassa limites críticos preestabelecidos. As hipóteses do teste são:

\[ \left\{\begin{matrix} H_0: Os \ parâmetros \ do \ modelo \ são \ estáveis \ ao \ longo \ do \ tempo \ (ausência \ de \ quebra \ estrutural.\\ H_a: Os \ parâmetros \ do \ modelo \ são \ instáveis \ (presença \ de \ quebra \ estrutural \ em \ algum \ ponto \ da \ amostra). \end{matrix}\right. \]

library(strucchange)

# Extraindo os resíduos do modelo estável que você já tem
# Vamos testar a estabilidade da primeira equação (ex: Preço do Boi)
res_eq1 <- residuals(modelo_var)[, 1]

# Calculando o processo CUSUM
cusum_teste <- efp(res_eq1 ~ 1, type = "OLS-CUSUM")

# Plotando o resultado
plot(cusum_teste)

res_eq2 <- residuals(modelo_var)[, 2]

# Calculando o processo CUSUM
cusum_teste2 <- efp(res_eq2 ~ 1, type = "OLS-CUSUM")

# Plotando o resultado
plot(cusum_teste2)

res_eq3 <- residuals(modelo_var)[, 3]

# Calculando o processo CUSUM
cusum_teste3 <- efp(res_eq3 ~ 1, type = "OLS-CUSUM")

# Plotando o resultado
plot(cusum_teste3)

Visualmente, o teste é interpretado por meio de um gráfico onde a soma acumulada deve permanecer dentro de duas faixas (limites de confiança de 95%). Se a linha do CUSUM ultrapassar esses limites, o modelo é considerado instável naquele período específico.

Portanto, como em nenhum gráfico das séries foi ultrapassados os limites, assume-se que as séries das variáveis endógenas são estáveis e não possuem quebra significativa da estrutura.

Resíduos sem Normalidade

No contexto econômico, é raro que séries de commodities gerem resíduos perfeitamente normais, dada a frequência de mudanças drásticas de preços em curtos intervalos de tempo. Conforme demonstrado nos gráficos de resíduos a seguir, embora a distribuição central se assemelhe a uma normal, a presença de outliers severos eleva a assimetria e a curtose, levando à rejeição da hipótese nula nos testes formais.

Contudo, a literatura especializada — destacando-se Kilian & Lütkepohl (2017) — aponta que a falha na normalidade não invalida o modelo para previsões de curto e longo prazo, desde que os estimadores sejam consistentes. Os autores enfatizam que a independência dos resíduos (ausência de autocorrelação) é um pressuposto muito mais crítico para a robustez do modelo do que a normalidade. Portanto, embora as bandas de confiança devam ser interpretadas com cautela, a validade estrutural do modelo permanece sustentada pela qualidade dos demais diagnósticos.

# Instale se necessário: install.packages("GGally")
library(GGally)
## Carregando pacotes exigidos: ggplot2
ggpairs(res_VECM, 
        upper = list(continuous = wrap("density", alpha = 0.5)),
        lower = list(continuous = wrap("points", alpha = 0.3, size=0.5)),
        diag = list(continuous = wrap("densityDiag"))) +
  theme_minimal() +
  labs(title = "Distribuição e Densidade dos Resíduos")

Portanto, o modelo construído foi validado estatisticamente. Agora é preciso analisar suas predições e IRF.

Predições

O modelo VEC estimado cumpre os principais pressupostos econométricos, com exceção da normalidade dos resíduos, fator que exige cautela na interpretação das variâncias, mas não inviabiliza a capacidade preditiva das séries. Para testar a robustez das projeções, aplicamos duas metodologias de validação distintas: o Hold-out Simples e a Janela Expansiva. Enquanto o primeiro método utiliza um corte estático entre as bases de treino e teste, a abordagem de janelas expansivas simula o comportamento real do mercado ao incorporar novos dados a cada passo temporal. Esta última técnica é a referência utilizada no mercado para setores dinâmicos, como a pecuária, uma vez que permite ao modelo capturar informações e variações estruturais mais recentes.

Inicialmente, apresentamos os resultados do método Hold-out acompanhados de suas métricas de erro, como o MAPE e o MAE, além de um gráfico de estresse que plota os valores reais e estimados sob uma banda de confiança. Na sequência, aplicamos a mesma análise para a validação por janela expansiva. Observa-se que a técnica de janelas apresenta resultados superiores, com uma aderência muito maior entre os dados reais e as estimativas, o que valida a eficácia do modelo em cenários de atualização periódica.

Validação Hold-out

previsao <- predict(modelo_treino, n.ahead = n_teste, exoPred = exog_teste)

# Extraindo os valores previstos (eles vêm em log)
pred_log_boi    <- previsao[, "log_boi"]
pred_log_carne  <- previsao[, "log_carne"]
pred_log_volume <- previsao[, "log_volume"]

# Exemplo para a Carne
erro_carne <- endog_teste[, "log_carne"] - pred_log_carne
rmse_carne <- sqrt(mean(erro_carne^2))
mape_carne <- mean(abs(erro_carne / endog_teste[, "log_carne"])) * 100

erro_boi <- endog_teste[, "log_boi"] - pred_log_boi
rmse_boi <- sqrt(mean(erro_boi^2))
mape_boi <- mean(abs(erro_boi / endog_teste[, "log_boi"])) * 100

erro_volume <- endog_teste[, "log_volume"] - pred_log_volume
rmse_volume <- sqrt(mean(erro_volume^2))
mape_volume <- mean(abs(erro_volume / endog_teste[, "log_volume"])) * 100

cat("O erro médio (MAPE) da carne   foi de:", mape_carne, "%")
## O erro médio (MAPE) da carne   foi de: 4.100327 %
cat("O erro médio (MAPE) do boi     foi de:", mape_boi, "%")
## O erro médio (MAPE) do boi     foi de: 5.273435 %
cat("O erro médio (MAPE) do volume  foi de:", mape_volume, "%")
## O erro médio (MAPE) do volume  foi de: 0.572873 %
# Reestimando o modelo (mantendo seus parâmetros de Rank e Lag)
modelo_validacao <- VECM(endog_treino, lag = 2, r = 1, 
                         include = "const", estim = "ML", 
                         exogen = exog_treino)


# Previsão de 6 meses à frente
previsoes <- predict(modelo_validacao, n.ahead = n_teste, exoPred = exog_teste)

# Transformando em data frame para facilitar a análise
df_validacao <- as.data.frame(previsoes)
colnames(df_validacao) <- c("prev_boi", "prev_carne", "prev_vol")

# Margem Real vs Margem Prevista (em Log)
margem_teste <- endog_teste[, "log_carne"] - endog_teste[, "log_boi"]
margem_prev <- df_validacao$prev_carne - df_validacao$prev_boi

# Calculando o Erro Médio Absoluto (MAE) da Margem
mae_margem <- mean(abs(margem_teste - margem_prev))

# Calculando o RMSE (Penaliza erros grandes)
rmse_margem <- sqrt(mean((margem_teste - margem_prev)^2))

cat("Erro Médio da Margem (MAE):", round(mae_margem * 100, 2), "%\n")
## Erro Médio da Margem (MAE): 3.26 %
previsao_banda <- predict(modelo_validacao, n.ahead = n_teste, 
                          exoPred = exog_teste)
# O output 'previsao_banda' agora contém as colunas 'fcst', 'lower' e 'upper'
# Vamos focar na Carne e no Boi para reconstruir a Margem

# Extraindo os valores previstos e seus intervalos
p_carne <- previsao_banda[, "log_carne"]
p_boi   <- previsao_banda[, "log_boi"]

# Erro padrão aproximado da margem (simplificado)
# Margem Prevista
margem_f <- p_carne - p_boi

# Calculando bandas de 95% para a margem
# Nota: Usamos o erro padrão do modelo para criar a "sombra"
se_margem <- sd(resid(modelo_validacao)[,2] - resid(modelo_validacao)[,1])
margem_up <- margem_f + (1.96 * se_margem)
margem_lw <- margem_f - (1.96 * se_margem)

par(mfrow = c(1,1))
# Gráfico de Validação
plot(1:n_teste, margem_teste, type="n", ylim=range(c(margem_lw, margem_up, margem_teste)),
     main="Teste de Stress: Margem Real vs Banda de Confiança",
     xlab="Meses", ylab="Log Margem")

# Adiciona a sombra da confiança (Área cinza)
polygon(c(1:n_teste, rev(1:n_teste)), c(margem_up, rev(margem_lw)), 
        col = rgb(0.8, 0.8, 0.8, 0.5), border = NA)

lines(1:n_teste, margem_teste, col="blue", lwd=2, type="b", pch=19) # REAL
lines(1:n_teste, margem_f, col="red", lwd=2, lty=2)               # PREVISTO
legend("bottomleft", legend=c("Real", "Previsto", "Confiança 95%"), 
       col=c("blue", "red", "gray"), lty=c(1,2,1), pch=c(19,NA,NA))

Validação em Janela Expansiva

# 1. Configurações Iniciais baseadas no seu código
n <- nrow(endogenas)
n_teste <- 18
n_treino_inicial <- n - n_teste

# Dataframe para guardar os resultados de cada "passo" da janela
resultados_janela <- data.frame()

# 2. Loop da Janela Expansiva
# O loop vai percorrer os 18 meses do período de teste
for (i in 0:(n_teste - 1)) {
  
  # A cada volta, o treino aumenta em 1 observação (i)
  fim_treino <- n_treino_inicial + i
  
  # Define treino expansivo
  endog_treino_exp <- endogenas[1:fim_treino, ]
  exog_treino_exp  <- exogenas[1:fim_treino, ]
  
  # Define o ponto exato do teste (apenas o mês n+1)
  # Precisamos da exógena do mês que queremos prever
  exog_teste_ponto <- exogenas[(fim_treino + 1), , drop = FALSE]
  real_ponto       <- endogenas[(fim_treino + 1), , drop = FALSE]
  
  # 3. Reestimar o modelo VECM com os parâmetros que você definiu
  modelo_exp <- VECM(endog_treino_exp, 
                     lag = 2, r = 1, 
                     include = "const", 
                     exogen = exog_treino_exp, 
                     estim = "ML")
  
  # 4. Previsão para 1 passo à frente (o mês seguinte imediato)
  prev_ponto <- predict(modelo_exp, n.ahead = 1, exoPred = exog_teste_ponto)
  
  # 5. Organizar os resultados
  res_temp <- data.frame(
    mes = fim_treino + 1,
    real_boi   = real_ponto[, "log_boi"],
    prev_boi   = prev_ponto[, "log_boi"],
    real_carne = real_ponto[, "log_carne"],
    prev_carne = prev_ponto[, "log_carne"],
    real_vol   = real_ponto[, "log_volume"],
    prev_vol   = prev_ponto[, "log_volume"]
  )
  
  resultados_janela <- rbind(resultados_janela, res_temp)
}

# --- 6. Cálculos de Precisão da Margem (Sua métrica final) ---

# Margem Real vs Margem Prevista (em Log)
resultados_janela$margem_real <- resultados_janela$real_carne - resultados_janela$real_boi
resultados_janela$margem_prev <- resultados_janela$prev_carne - resultados_janela$prev_boi

# Erro da Margem
erro_margem_exp <- resultados_janela$margem_real - resultados_janela$margem_prev

# Métricas Finais
mae_margem_final  <- mean(abs(erro_margem_exp))
rmse_margem_final <- sqrt(mean(erro_margem_exp^2))

# Exibição dos Resultados
cat("\n--- RESULTADOS JANELA EXPANSIVA (18 MESES) ---\n")
## 
## --- RESULTADOS JANELA EXPANSIVA (18 MESES) ---
cat("MAE da Margem: ", round(mae_margem_final * 100, 2), "%\n")
## MAE da Margem:  2.64 %
cat("RMSE da Margem:", round(rmse_margem_final * 100, 2), "%\n")
## RMSE da Margem: 3.32 %
# 1. Calcular o Erro e o Desvio Padrão para as Bandas
# Vamos usar os resultados da janela expansiva (resultados_janela)
resultados_janela$erro_margem <- resultados_janela$margem_real - resultados_janela$margem_prev

# Desvio padrão dos erros de previsão (medida de incerteza)
sd_erro <- sd(resultados_janela$erro_margem)

# 2. Definir os limites das bandas (ex: 95% de confiança = 1.96 * SD)
# Você pode ajustar o multiplicador para simular cenários de maior estresse
multiplicador <- 1.96 

resultados_janela$banda_sup <- resultados_janela$margem_prev + (multiplicador * sd_erro)
resultados_janela$banda_inf <- resultados_janela$margem_prev - (multiplicador * sd_erro)

ggplot(resultados_janela, aes(x = mes)) +
  # Área das bandas de estresse
  geom_ribbon(aes(ymin = banda_inf, ymax = banda_sup), fill = "gray80", alpha = 0.5) +
  # Linha da Margem Real
  geom_line(aes(y = margem_real, color = "Margem Real"), size = 1) +
  # Linha da Margem Prevista (Modelo)
  geom_line(aes(y = margem_prev, color = "Previsto (VECM)"), linetype = "dashed", size = 1) +
  # Estilização
  scale_color_manual(values = c("Margem Real" = "black", "Previsto (VECM)" = "red")) +
  labs(title = "Validação VECM: Margem de Abate vs. Bandas de Estresse",
       subtitle = "Janela Expansiva (18 meses de teste)",
       x = "Meses de Teste",
       y = "Margem (Log)",
       color = "Legenda") +
  theme_minimal() +
  theme(legend.position = "bottom")