Receita das carnes
William Yuji Kohatsu Tanigava
Bibliotecas
library(tseries)
library(urca)
library(readxl)
library(vars)
library(forecast)
library(car)
library(ggplot2)Importando as séries
dados <- read_excel("D:/Estatística/Bases de cursos/Frigorífico/Preço de compra e venda Boi Gordo mensal 2014-2025.xlsx",
sheet = "RL_carne (2)")names(dados)## [1] "ID_semana" "N_semana" "Valor_carne" "CEPEA_quilo"
## [5] "Peso_carcaça" "Volume_semanal" "Receita_liquida" "Dólar"
## [9] "Ano"carne <- ts(dados$Valor_carne, start = c(2023, 14), end = c(2026, 13), frequency = 50)
cepea <- ts(dados$CEPEA_quilo, start = c(2023, 14), end = c(2026, 13), frequency = 50)
rendim <- ts(dados$Peso_carcaça, start = c(2023, 14), end = c(2026, 13), frequency = 50)
volume <- ts(dados$Volume_semanal, start = c(2023, 14), end = c(2026, 13), frequency = 50)
dolar <- ts(dados$Dólar, start = c(2023, 14), end = c(2026, 13), frequency = 50)
receita <- ts(dados$Receita_liquida,start = c(2023, 14), end = c(2026, 13), frequency = 50)Aplicando uma diferença caso necessite posteriormente.
carne_dif <- diff(carne)
cepea_dif <- diff(cepea)
rendim_dif <- diff(rendim)
volume_dif <- diff(volume)
dolar_dif <- diff(dolar)
receita_dif <- diff(receita)Visualização das séries e estacionariedade
Lembrando que, no teste de Dicky-Fuller a hipótese alternativa é a de estacionariedade da série. Logo, quando o \(p_{valor}\) for menor que 0.05 isso significa que o teste está apontando estacionariedade da série, caso seja maior ele aponta que a série é não estacionária.
plot.ts(carne , main = "Série semanal do preço do Carne no atacado" , ylab = "Real (R$)")
adf.test(carne)##
## Augmented Dickey-Fuller Test
##
## data: carne
## Dickey-Fuller = -2.7753, Lag order = 5, p-value = 0.2534
## alternative hypothesis: stationaryadf.test(carne_dif)##
## Augmented Dickey-Fuller Test
##
## data: carne_dif
## Dickey-Fuller = -4.0746, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(cepea , main = "Série semanal do Índice do Boi Gordo/kg" , ylab = "Real (R$)")
adf.test(cepea)##
## Augmented Dickey-Fuller Test
##
## data: cepea
## Dickey-Fuller = -2.2754, Lag order = 5, p-value = 0.4617
## alternative hypothesis: stationaryadf.test(cepea_dif)##
## Augmented Dickey-Fuller Test
##
## data: cepea_dif
## Dickey-Fuller = -4.2173, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(rendim , main = "Série semanal do peso rendido em carcaça" , ylab = "Peso (Kg)")
adf.test(rendim)##
## Augmented Dickey-Fuller Test
##
## data: rendim
## Dickey-Fuller = -4.2392, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryadf.test(rendim_dif)##
## Augmented Dickey-Fuller Test
##
## data: rendim_dif
## Dickey-Fuller = -10.812, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(dolar , main = "Série semanal da cotação do Dólar comercial", ylab = "Real (R$)")
adf.test(dolar)##
## Augmented Dickey-Fuller Test
##
## data: dolar
## Dickey-Fuller = -0.85128, Lag order = 5, p-value = 0.9549
## alternative hypothesis: stationaryadf.test(dolar_dif)##
## Augmented Dickey-Fuller Test
##
## data: dolar_dif
## Dickey-Fuller = -4.8312, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(volume, main = "Série semanal da quantiade de bovinos abatidos ", ylab = "Unidade")
adf.test(volume)##
## Augmented Dickey-Fuller Test
##
## data: volume
## Dickey-Fuller = -6.0531, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryadf.test(volume_dif)##
## Augmented Dickey-Fuller Test
##
## data: volume_dif
## Dickey-Fuller = -9.8487, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(receita , main = "Série semanal da receita líquida das carnes", ylab = "Real (R$)")
adf.test(receita)##
## Augmented Dickey-Fuller Test
##
## data: receita
## Dickey-Fuller = -3.464, Lag order = 5, p-value = 0.04812
## alternative hypothesis: stationaryadf.test(receita_dif)##
## Augmented Dickey-Fuller Test
##
## data: receita_dif
## Dickey-Fuller = -9.2089, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryBasicamente, todas as séries diferenciadas em uma ordem são estacionárias, sendo que a de volume, rendimento da carcaça e receita já apresentavam estacionariedade em nível. Logo, posteriormente se os resíduos do modelo de regressão não apresentar os pressupostos do ARIMA com as séries em nível pode ser utilizado as séries diferenciadas.
Modelo de regressão e escolha das variáveis exógenas
modelo <- lm(receita ~ carne + cepea + rendim + volume + dolar)
summary(modelo)##
## Call:
## lm(formula = receita ~ carne + cepea + rendim + volume + dolar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1423630 -407647 -36979 362280 2327884
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.973e+07 1.701e+06 -17.480 <2e-16 ***
## carne 9.524e+05 7.502e+04 12.696 <2e-16 ***
## cepea -1.307e+05 7.758e+04 -1.684 0.0943 .
## rendim 5.195e+04 5.056e+03 10.275 <2e-16 ***
## volume 4.822e+03 1.574e+02 30.624 <2e-16 ***
## dolar -1.899e+05 2.010e+05 -0.945 0.3465
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 622300 on 144 degrees of freedom
## Multiple R-squared: 0.9612, Adjusted R-squared: 0.9599
## F-statistic: 714.4 on 5 and 144 DF, p-value: < 2.2e-16Agora vamos avaliar os resíduos e a multicolinearidade das variáveis.
vif(modelo)## carne cepea rendim volume dolar
## 21.337840 18.925513 1.010640 1.011934 1.843163As variáveis apresentam VIF próximos de 1, mostrando que elas não são correlacionadas.
res <- modelo$residuals
checkresiduals(res)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 12.272, df = 10, p-value = 0.2673
##
## Model df: 0. Total lags used: 10Pelo teste de Ljung-Box, os resíduos não apresentam autocorrelação, logo eles são independentes. Portanto, não tem necessidade da modelagem ARIMA em cima deles. Agora iremos testar com as variáveis diferenciadas, caso algum modelo apresente os pressupostos pedidos, a análise será prosseguida com ele.
modelo_dif <- lm(receita_dif ~ carne_dif + rendim_dif + volume_dif)
summary(modelo_dif)##
## Call:
## lm(formula = receita_dif ~ carne_dif + rendim_dif + volume_dif)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2492777 -604054 25322 652280 2342774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7944.1 77466.4 0.103 0.91846
## carne_dif 584136.6 207038.6 2.821 0.00545 **
## rendim_dif 56442.4 5413.6 10.426 < 2e-16 ***
## volume_dif 4708.2 163.8 28.742 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 938900 on 145 degrees of freedom
## Multiple R-squared: 0.8741, Adjusted R-squared: 0.8715
## F-statistic: 335.7 on 3 and 145 DF, p-value: < 2.2e-16vif(modelo_dif)## carne_dif rendim_dif volume_dif
## 1.003486 1.010597 1.007354resid_dif <- modelo_dif$residuals
checkresiduals(resid_dif)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 56.562, df = 10, p-value = 1.611e-08
##
## Model df: 0. Total lags used: 10adf.test(resid_dif)##
## Augmented Dickey-Fuller Test
##
## data: resid_dif
## Dickey-Fuller = -8.0962, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryAmbos os modelos apresentam resíduos estacionários e autocorrelacionados, deste modo é possível construir um modelo SARIMA. O segundo modelo contendo a soma do couro e sebo carrega um pouco mais de informação, portanto as variáveis exógenas escolhidas para o modelo SARIMAX são: Mix de couro e sebo, volume de abate e rendimento do boi para miúdos. O dólar não apresentou significância para o modelo de regressão, por este motivo que ele ficou de fora.
acf(resid_dif, 60)
pacf(resid_dif, 60)
Modelo SARIMAX
Para construirmos o melhor modelo SARIMAX é necessário seguir alguns procedimentos padrões de qualquer modelo. A divisão da base em treino e teste verifica a capacidade de generalização que o modelo construído possui e também possibilita observar a precisão das predições.
mix_carne <- carne
exogenas <- cbind(mix_carne, rendim, volume)
endogena <- receitan <- length(receita)
teste <- 125
y_treino <- subset(endogena, end = teste)
x_treino <- exogenas[1:teste, ]
y_teste <- subset(endogena, start = teste + 1)
x_teste <- exogenas[(teste + 1):n, ]arimax_treino_carne <- Arima(y_treino, xreg = x_treino, order = c(2,1,4))previsao <- forecast(arimax_treino_carne, h = length(y_teste), xreg = x_teste)
metricas <- accuracy(previsao, y_teste)
print(metricas)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 65848.9 546196.5 425085.9 0.3587394 2.925035 0.1140183 -0.0176694
## Test set 165776.0 737015.8 602707.9 0.5526748 3.311010 0.1616607 -0.1003889
## Theil's U
## Training set NA
## Test set 0.1809108Testando a estabilidade das previsões
f_arimax <- function(x, h, xreg, newxreg) {
forecast(Arima(x, order=c(2,1,4), xreg=xreg), h=h, xreg=newxreg)
}
erros_cv <- tsCV(endogena, f_arimax, h=5, xreg=exogenas)RMSE_h5 <- sqrt(mean(erros_cv[,5]^2, na.rm = TRUE))
MAE_h5 <- mean(abs(erros_cv [,5]), na.rm = TRUE)
MAPE_h5 <- mean(abs(erros_cv [,5] / endogena), na.rm = TRUE) * 100
cat("Métricas calculadas para h=5 \n")## Métricas calculadas para h=5cat("RMSE: ", RMSE_h5, "\t \n")## RMSE: 768281.7 cat("MAE: " , MAE_h5 , "\t \n")## MAE: 579630.6 cat("MAPE: ", MAPE_h5, "% \t \n")## MAPE: 3.970935 % plot(erros_cv[,5], main = "Erros na validação de janela com horizonte de 5 tempos",
ylab = "Resíduo", xlab = "Tempo")
abline(h=0, lwd = 0.5, lty = 2)
abline(h=mean(abs(erros_cv[,5]), na.rm = TRUE), col = "red")
abline(h=mean(erros_cv[,5], na.rm = TRUE), col = "green")
plot(erros_cv[,1], main = "Erros na validação de janela com horizonte de 1 tempo",
ylab = "Resíduo", xlab = "Tempo")
abline(h=0, lwd = 0.5, lty = 2)
abline(h=mean(abs(erros_cv[,1]), na.rm = TRUE), col = "red")
abline(h=mean(erros_cv[,1], na.rm = TRUE), col = "green")
RMSE_h1 <- sqrt(mean(erros_cv[,1]^2, na.rm = TRUE))
MAE_h1 <- mean(abs(erros_cv[,1]), na.rm = TRUE)
MAPE_h1 <- mean(abs(erros_cv[,1] / endogena), na.rm = TRUE) * 100
cat("Resultados para h=1:\n")## Resultados para h=1:cat("RMSE:", RMSE_h1, "\n")## RMSE: 720921.1cat("MAE: ", MAE_h1, "\n")## MAE: 562222.8cat("MAPE:", MAPE_h1, "%\n")## MAPE: 3.750888 %plot(erros_cv[,3], main = "Erros na validação de janela com horizonte de 3 tempos",
ylab = "Resíduo", xlab = "Tempo")
abline(h=0, lwd = 0.5, lty = 2)
abline(h=mean(abs(erros_cv[,3]), na.rm = TRUE), col = "red")
abline(h=mean(erros_cv[,3], na.rm = TRUE), col = "green")
RMSE_h3 <- sqrt(mean(erros_cv[,3]^2, na.rm = TRUE))
MAE_h3 <- mean(abs(erros_cv[,3]), na.rm = TRUE)
MAPE_h3 <- mean(abs(erros_cv[,3] / endogena), na.rm = TRUE) * 100
cat("Resultados para h=3:\n")## Resultados para h=3:cat("RMSE:", RMSE_h3, "\n")## RMSE: 715983.2cat("MAE: ", MAE_h3, "\n")## MAE: 557425.1cat("MAPE:", MAPE_h3, "%\n")## MAPE: 3.736278 %boxplot(erros_cv[,1], erros_cv[,3], erros_cv[,5],
names = c("h=1", "h=3", "h=5"),
main = "Distribuição do Erro por Horizonte",
col = c("lightblue", "orange", "red"),
ylab = "Erro Bruto")
abline(h = 0, lty = 2)
Analisando a validação do modelo para predição da receita em diferentes horizontes de tempo é visto que o modelo para curto prazo se sai melhor do que a longo prazo, porém ele ainda continua confiável.
summary(arimax_treino_carne)## Series: y_treino
## Regression with ARIMA(2,1,4) errors
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4 mix_carne rendim
## -1.6070 -0.9273 0.5569 -0.8265 -0.7859 0.2276 778871.18 47304.497
## s.e. 0.0428 0.0572 0.1070 0.1041 0.0968 0.1049 29157.41 5145.874
## volume
## 4565.1232
## s.e. 142.2514
##
## sigma^2 = 3.243e+11: log likelihood = -1818.23
## AIC=3656.45 AICc=3658.4 BIC=3684.66
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 65848.9 546196.5 425085.9 0.3587394 2.925035 0.1140183 -0.0176694coeftest(arimax_treino_carne)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.6070e+00 4.2819e-02 -37.5305 < 2.2e-16 ***
## ar2 -9.2726e-01 5.7231e-02 -16.2019 < 2.2e-16 ***
## ma1 5.5690e-01 1.0702e-01 5.2037 1.954e-07 ***
## ma2 -8.2647e-01 1.0411e-01 -7.9383 2.049e-15 ***
## ma3 -7.8592e-01 9.6837e-02 -8.1159 4.823e-16 ***
## ma4 2.2755e-01 1.0490e-01 2.1692 0.03007 *
## mix_carne 7.7887e+05 2.9157e+04 26.7126 < 2.2e-16 ***
## rendim 4.7304e+04 5.1459e+03 9.1927 < 2.2e-16 ***
## volume 4.5651e+03 1.4225e+02 32.0919 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Os coeficientes do modelo são todos estatisticamente significativos.
Análise dos pressupostos dos resíduos do modelo
residuo_ARIMA <- arimax_treino_carne$residuals
checkresiduals(arimax_treino_carne)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(2,1,4) errors
## Q* = 27.981, df = 19, p-value = 0.08379
##
## Model df: 6. Total lags used: 25Segundo o teste de Ljung-Box os resíduos do modelo construído não são autocorrelacionados, agora é preciso verificar se são estacionários.
adf.test(residuo_ARIMA)##
## Augmented Dickey-Fuller Test
##
## data: residuo_ARIMA
## Dickey-Fuller = -5.3394, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryDe acordo com o teste de Dickey-Fuller, os resíduos são estacionários. Deste modo, conclui-se que o modelo ARIMAX(1,1,2) possui resíduos ruído branco, sendo assim válido estatisticamente. Por fim, será testado a normalidade dos resíduos, caso ele seja o intervalo de confiança construído terá uma precisão garantido, caso contrário, a banda de confiança não é totalmente seguro.
shapiro.test(residuo_ARIMA)##
## Shapiro-Wilk normality test
##
## data: residuo_ARIMA
## W = 0.98228, p-value = 0.101qqnorm(residuo_ARIMA)
qqline(residuo_ARIMA, col = "red")
Pelo teste de Shapiro-Wilk e observando o QQplot, confirmamos que os resíduos do modelo não são normais. Porém, ele ainda é um modelo válido visto que a série da receita líquida tinha alguns pontos outliers, e muito provalvemente o modelo está errando na previsão destes pontos (isso é um bom sinal, pois mostra que o modelo não superajustou aos dados overfitting).
Como dito anteriormente, as bandas de confiança podem não ser precisas, porém a predição da receita tem seu valor, visto que o modelo acerta bem a média.
Visualização da predição do modelo
autoplot(previsao) +
ggtitle("[Carne] Receita Líquida estimada com bandas de confiança") +
xlab("Ano") +
ylab("Real (R$)") +
theme_minimal() +
scale_fill_manual(values = c("gray80", "gray60")) # Cores das bandas
ajustados <- fitted(arimax_treino_carne)
autoplot(endogena) +
# 1. Definimos a série real com um nome para a legenda
autolayer(endogena, series = "Real") +
# 2. Definimos a série estimada com um nome para a legenda
autolayer(ajustados, series = "Estimado") +
ggtitle("Real vs. Estimado - Receita Líquida na base treino") +
xlab("Ano") + ylab("Real (R$)") +
# 3. Agora o scale_color_manual consegue encontrar os nomes
scale_color_manual(values = c("Estimado" = "red", "Real" = "black")) +
guides(color = guide_legend(title = "", override.aes = list(fill = NA))) +
theme_minimal()
df_plot <- data.frame(
Tempo = 1:length(y_teste),
Real = as.numeric(y_teste),
Previsto = as.numeric(previsao$mean),
Inferior95 = as.numeric(previsao$lower[,2]),
Superior95 = as.numeric(previsao$upper[,2])
)ggplot(df_plot, aes(x = Tempo)) +
# Banda de Confiança
geom_ribbon(aes(ymin = Inferior95, ymax = Superior95), fill = "lightblue", alpha = 0.3) +
# Linhas
geom_line(aes(y = Previsto, color = "Estimado"), size = 1) +
geom_line(aes(y = Real, color = "Real"), size = 1) +
# Cores manuais
scale_color_manual(values = c("Estimado" = "blue", "Real" = "black")) +
guides(color = guide_legend(title = "", override.aes = list(linetype = c(1, 1)))) +
labs(title = "[Carne] Receita Líquida Real x Estimada - Base de teste",
subtitle = "Sombra azul representa o intervalo de confiança de 95%",
x = "Semanas", y = "Real (R$)") +
theme_light() +
theme(legend.position = "bottom", legend.key = element_blank())
Modelo base completa
Agora será construido o modelo com as ordens encontradas na base de treino usando todas as observações, para verificar se o modelo consegue generalizar bem.
arimax_final <- Arima(endogena, xreg = exogenas, order = c(2,1,4))checkresiduals(arimax_final)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(2,1,4) errors
## Q* = 29.585, df = 24, p-value = 0.1989
##
## Model df: 6. Total lags used: 30summary(arimax_final)## Series: endogena
## Regression with ARIMA(2,1,4) errors
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4 mix_carne rendim
## -0.8441 0.0612 -0.2506 -0.9111 0.2884 -0.0510 789984.97 48226.932
## s.e. 1.4564 1.3361 1.4576 0.2711 1.3956 0.3293 31557.49 4978.128
## volume
## 4826.1281
## s.e. 158.0107
##
## sigma^2 = 3.864e+11: log likelihood = -2196.82
## AIC=4413.63 AICc=4415.23 BIC=4443.67
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 63101.18 600537.2 469659.7 0.3112193 3.136457 0.1377485
## ACF1
## Training set -0.01495523summary(arimax_treino_carne)## Series: y_treino
## Regression with ARIMA(2,1,4) errors
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4 mix_carne rendim
## -1.6070 -0.9273 0.5569 -0.8265 -0.7859 0.2276 778871.18 47304.497
## s.e. 0.0428 0.0572 0.1070 0.1041 0.0968 0.1049 29157.41 5145.874
## volume
## 4565.1232
## s.e. 142.2514
##
## sigma^2 = 3.243e+11: log likelihood = -1818.23
## AIC=3656.45 AICc=3658.4 BIC=3684.66
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 65848.9 546196.5 425085.9 0.3587394 2.925035 0.1140183 -0.0176694As métricas de erro do modelo são bem parecidas e com bons resultados.
O que ainda falta?
Um teste de stress e cenário, simular variáveis exógenas com cenários positivos, negativos e neutros.