Receita Miúdos
William Yuji Kohatsu Tanigava
Bibliotecas
library(tseries)
library(urca)
library(readxl)
library(vars)
library(forecast)
library(car)
library(ggplot2)Importando as séries
dado <- read_excel("D:/Estatística/Bases de cursos/Frigorífico/Preço de compra e venda Boi Gordo mensal 2014-2025.xlsx",
sheet = "RL_miudos (2)")
couro <- ts(dado$Preço_couro, start = c(2023, 14), end = c(2026, 13), frequency = 50)
sebo <- ts(dado$Preço_sebo, start = c(2023, 14), end = c(2026, 13), frequency = 50)
rend <- ts(dado$Peso_miudos, start = c(2023, 14), end = c(2026, 13), frequency = 50)
volume <- ts(dado$Volume_semanal, start = c(2023, 14), end = c(2026, 13), frequency = 50)
dolar <- ts(dado$Dólar, start = c(2023, 14), end = c(2026, 13), frequency = 50)
rec <- ts(dado$Receita_liquida,start = c(2023, 14), end = c(2026, 13), frequency = 50)Aplicando uma diferença caso necessite posteriormente.
couro_dif <- diff(couro)
sebo_dif <- diff(sebo)
rend_dif <- diff(rend)
volume_dif <- diff(volume)
dolar_dif <- diff(dolar)
rec_dif <- diff(rec)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(couro , main = "Série semanal do preço do Couro" , ylab = "Real (R$)")
adf.test(couro)##
## Augmented Dickey-Fuller Test
##
## data: couro
## Dickey-Fuller = -2.6253, Lag order = 5, p-value = 0.3159
## alternative hypothesis: stationaryadf.test(couro_dif)##
## Augmented Dickey-Fuller Test
##
## data: couro_dif
## Dickey-Fuller = -4.9293, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(sebo , main = "Série semanal do preço do Sebo" , ylab = "Real (R$)")
adf.test(sebo)##
## Augmented Dickey-Fuller Test
##
## data: sebo
## Dickey-Fuller = -2.2939, Lag order = 5, p-value = 0.4539
## alternative hypothesis: stationaryadf.test(sebo_dif)##
## Augmented Dickey-Fuller Test
##
## data: sebo_dif
## Dickey-Fuller = -5.1024, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryplot.ts(rend , main = "Série semanal do peso rendido para miúdos" , ylab = "Peso (Kg)")
adf.test(rend)##
## Augmented Dickey-Fuller Test
##
## data: rend
## Dickey-Fuller = -3.6986, Lag order = 5, p-value = 0.02668
## alternative hypothesis: stationaryadf.test(rend_dif)##
## Augmented Dickey-Fuller Test
##
## data: rend_dif
## Dickey-Fuller = -9.4849, 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(rec , main = "Série semanal da receita líquida dos miúdos", ylab = "Real (R$)")
adf.test(rec)##
## Augmented Dickey-Fuller Test
##
## data: rec
## Dickey-Fuller = -2.0686, Lag order = 5, p-value = 0.5478
## alternative hypothesis: stationaryadf.test(rec_dif)##
## Augmented Dickey-Fuller Test
##
## data: rec_dif
## Dickey-Fuller = -8.7344, 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 e a de rendimento dos miúdos 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
model <- lm(rec ~ couro + sebo + rend + volume + dolar)
summary(model)##
## Call:
## lm(formula = rec ~ couro + sebo + rend + volume + dolar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -235861 -47512 1110 51997 214885
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.212e+06 1.764e+05 -29.554 <2e-16 ***
## couro -3.436e+05 2.720e+04 -12.632 <2e-16 ***
## sebo -3.253e+04 1.299e+04 -2.504 0.0134 *
## rend 2.569e+04 1.070e+03 24.019 <2e-16 ***
## volume 5.932e+02 1.982e+01 29.928 <2e-16 ***
## dolar 7.622e+05 2.408e+04 31.655 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 78330 on 144 degrees of freedom
## Multiple R-squared: 0.9534, Adjusted R-squared: 0.9518
## F-statistic: 589.4 on 5 and 144 DF, p-value: < 2.2e-16Agora vamos avaliar os resíduos e a multicolinearidade das variáveis.
vif(model)## couro sebo rend volume dolar
## 2.080823 1.394552 1.015320 1.012081 1.668600As variáveis apresentam VIF próximos de 1, mostrando que elas não são correlacionadas.
residuo <- model$residuals
checkresiduals(residuo)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 7.0771, df = 10, p-value = 0.7181
##
## 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.
model_dif <- lm(rec_dif ~ couro_dif + rend_dif + volume_dif)
summary(model_dif)##
## Call:
## lm(formula = rec_dif ~ couro_dif + rend_dif + volume_dif)
##
## Residuals:
## Min 1Q Median 3Q Max
## -268132 -50144 3311 53273 269701
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3116.08 7340.98 -0.424 0.672
## couro_dif 930325.35 159637.68 5.828 3.5e-08 ***
## rend_dif 25982.20 866.74 29.977 < 2e-16 ***
## volume_dif 602.02 15.62 38.532 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 89340 on 145 degrees of freedom
## Multiple R-squared: 0.9487, Adjusted R-squared: 0.9476
## F-statistic: 893.1 on 3 and 145 DF, p-value: < 2.2e-16vif(model_dif)## couro_dif rend_dif volume_dif
## 1.000784 1.011733 1.012145res_dif <- model_dif$residuals
checkresiduals(res_dif)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 30.468, df = 10, p-value = 0.0007178
##
## Model df: 0. Total lags used: 10adf.test(res_dif)##
## Augmented Dickey-Fuller Test
##
## data: res_dif
## Dickey-Fuller = -5.7894, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarysoma_miudos <- couro_dif + sebo_dif
md <- lm(rec_dif ~ soma_miudos + rend_dif + volume_dif)
summary(md)##
## Call:
## lm(formula = rec_dif ~ soma_miudos + rend_dif + volume_dif)
##
## Residuals:
## Min 1Q Median 3Q Max
## -280886 -57490 4869 50425 255312
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1776.34 7260.12 0.245 0.807
## soma_miudos 231292.52 37958.50 6.093 9.48e-09 ***
## rend_dif 26041.78 859.34 30.304 < 2e-16 ***
## volume_dif 601.84 15.49 38.862 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 88560 on 145 degrees of freedom
## Multiple R-squared: 0.9496, Adjusted R-squared: 0.9485
## F-statistic: 909.8 on 3 and 145 DF, p-value: < 2.2e-16vif(md)## soma_miudos rend_dif volume_dif
## 1.001094 1.012127 1.012018res_md <- md$residuals
checkresiduals(res_md)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 33.529, df = 10, p-value = 0.0002219
##
## Model df: 0. Total lags used: 10adf.test(res_md)##
## Augmented Dickey-Fuller Test
##
## data: res_md
## Dickey-Fuller = -5.6569, 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(res_md, 60)
pacf(res_md, 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_miudos <- couro + sebo
exogenas <- cbind(mix_miudos, rend, volume)
endogena <- recn <- length(rec)
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, ]Agora iremos rodar uma função para buscar o modelo com o menor AIC, este que será considerado como o melhor modelo.
arimax_treino <- auto.arima(y_treino, xreg = x_treino,
stepwise = TRUE, approximation = FALSE,
trace = TRUE, seasonal = TRUE)##
## Regression with ARIMA(2,1,2) errors : 3168.987
## Regression with ARIMA(0,1,0) errors : 3191.775
## Regression with ARIMA(1,1,0) errors : 3169.207
## Regression with ARIMA(0,1,1) errors : 3164.514
## Regression with ARIMA(0,1,0) errors : 3189.659
## Regression with ARIMA(1,1,1) errors : 3165.815
## Regression with ARIMA(0,1,2) errors : 3165.361
## Regression with ARIMA(1,1,2) errors : 3154.627
## Regression with ARIMA(1,1,3) errors : 3156.892
## Regression with ARIMA(0,1,3) errors : 3166.099
## Regression with ARIMA(2,1,1) errors : 3166.828
## Regression with ARIMA(2,1,3) errors : Inf
## Regression with ARIMA(1,1,2) errors : 3153.861
## Regression with ARIMA(0,1,2) errors : 3163.497
## Regression with ARIMA(1,1,1) errors : 3163.962
## Regression with ARIMA(2,1,2) errors : 3167.042
## Regression with ARIMA(1,1,3) errors : 3154.701
## Regression with ARIMA(0,1,1) errors : 3162.714
## Regression with ARIMA(0,1,3) errors : 3164.207
## Regression with ARIMA(2,1,1) errors : 3164.912
## Regression with ARIMA(2,1,3) errors : 3156.832
##
## Best model: Regression with ARIMA(1,1,2) errorsprevisao <- forecast(arimax_treino, h = length(y_teste), xreg = x_teste)
metricas <- accuracy(previsao, y_teste)
print(metricas)## ME RMSE MAE MPE MAPE MASE
## Training set -1528.277 75320.05 59219.72 -0.04224223 3.086637 0.11541793
## Test set 20771.353 55603.58 41126.41 1.36323547 2.438880 0.08015448
## ACF1 Theil's U
## Training set -0.025251122 NA
## Test set 0.002903374 0.1413551Testando a estabilidade das previsões
f_arimax <- function(x, h, xreg, newxreg) {
forecast(Arima(x, order=c(1,1,2), 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: 109632.3 cat("MAE: " , MAE_h5 , "\t \n")## MAE: 86395.77 cat("MAPE: ", MAPE_h5, "% \t \n")## MAPE: 4.414142 % 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: 83504.67cat("MAE: ", MAE_h1, "\n")## MAE: 64873.01cat("MAPE:", MAPE_h1, "%\n")## MAPE: 3.286933 %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: 94155.38cat("MAE: ", MAE_h3, "\n")## MAE: 74325.49cat("MAPE:", MAPE_h3, "%\n")## MAPE: 3.748622 %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)## Series: y_treino
## Regression with ARIMA(1,1,2) errors
##
## Coefficients:
## ar1 ma1 ma2 mix_miudos rend volume
## 0.9769 -1.581 0.6249 235410.71 26115.5703 601.4873
## s.e. 0.0220 0.083 0.0860 37106.12 959.8568 18.1877
##
## sigma^2 = 6.01e+09: log likelihood = -1569.45
## AIC=3152.9 AICc=3153.86 BIC=3172.64
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -1528.277 75320.05 59219.72 -0.04224223 3.086637 0.1154179
## ACF1
## Training set -0.02525112coeftest(arimax_treino)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 9.7687e-01 2.2007e-02 44.3889 < 2.2e-16 ***
## ma1 -1.5810e+00 8.2973e-02 -19.0544 < 2.2e-16 ***
## ma2 6.2488e-01 8.5999e-02 7.2661 3.700e-13 ***
## mix_miudos 2.3541e+05 3.7106e+04 6.3443 2.235e-10 ***
## rend 2.6116e+04 9.5986e+02 27.2078 < 2.2e-16 ***
## volume 6.0149e+02 1.8188e+01 33.0711 < 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$residuals
checkresiduals(arimax_treino)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(1,1,2) errors
## Q* = 15.409, df = 22, p-value = 0.8441
##
## Model df: 3. 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.064, 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.97878, p-value = 0.04639qqnorm(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("[Miúdos] 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)
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 = "[Miúdos] 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(1,1,2))checkresiduals(arimax_final)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(1,1,2) errors
## Q* = 22.438, df = 27, p-value = 0.7148
##
## Model df: 3. Total lags used: 30summary(arimax_final)## Series: endogena
## Regression with ARIMA(1,1,2) errors
##
## Coefficients:
## ar1 ma1 ma2 mix_miudos rend volume
## 0.9635 -1.6341 0.6915 231423.47 25831.2847 586.7171
## s.e. 0.0292 0.0648 0.0590 31785.25 882.5591 16.1114
##
## sigma^2 = 5.448e+09: log likelihood = -1879.16
## AIC=3772.33 AICc=3773.12 BIC=3793.36
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 133.9889 72068.05 56794.79 0.03491869 2.988553 0.1080259
## ACF1
## Training set 0.007864387summary(arimax_treino)## Series: y_treino
## Regression with ARIMA(1,1,2) errors
##
## Coefficients:
## ar1 ma1 ma2 mix_miudos rend volume
## 0.9769 -1.581 0.6249 235410.71 26115.5703 601.4873
## s.e. 0.0220 0.083 0.0860 37106.12 959.8568 18.1877
##
## sigma^2 = 6.01e+09: log likelihood = -1569.45
## AIC=3152.9 AICc=3153.86 BIC=3172.64
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -1528.277 75320.05 59219.72 -0.04224223 3.086637 0.1154179
## ACF1
## Training set -0.02525112As 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.