Receita da planta
William Yuji Kohatsu Tanigava
Bibliotecas
library(tseries)
library(urca)
library(readxl)
library(vars)
library(forecast)
library(car)
library(ggplot2)Importando as séries
data <- read_excel("D:/Estatística/Bases de cursos/Frigorífico/Preço de compra e venda Boi Gordo mensal 2014-2025.xlsx",
sheet = "RL_planta")names(data)## [1] "ID_semana" "N_semana" "CEPEA_quilo" "IPCA_alimentos"
## [5] "Volume_semanal" "Receita_liquida" "Dólar" "Ano"cepea <- ts(data$CEPEA_quilo, start = c(2023, 14), end = c(2026, 13), frequency = 50)
ipca <- ts(data$IPCA_alimentos, start = c(2023, 14), end = c(2026, 13), frequency = 50)
volume <- ts(data$Volume_semanal, start = c(2023, 14), end = c(2026, 13), frequency = 50)
dolar <- ts(data$Dólar, start = c(2023, 14), end = c(2026, 13), frequency = 50)
rl <- ts(data$Receita_liquida,start = c(2023, 14), end = c(2026, 13), frequency = 50)Aplicando uma diferença caso necessite posteriormente.
cepea_dif <- diff(cepea)
ipca_dif <- diff(ipca)
volume_dif <- diff(volume)
dolar_dif <- diff(dolar)
rl_dif <- diff(rl)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(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(ipca , main = "Variação do IPCA alimentação e bebida" , ylab = "Índice IPCA")
adf.test(ipca)##
## Augmented Dickey-Fuller Test
##
## data: ipca
## Dickey-Fuller = -1.8493, Lag order = 5, p-value = 0.6392
## alternative hypothesis: stationaryadf.test(ipca_dif)##
## Augmented Dickey-Fuller Test
##
## data: ipca_dif
## Dickey-Fuller = -1.745, Lag order = 5, p-value = 0.6826
## 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(rl , main = "Série semanal da receita líquida da planta", ylab = "Real (R$)")
adf.test(rl)##
## Augmented Dickey-Fuller Test
##
## data: rl
## Dickey-Fuller = -3.383, Lag order = 5, p-value = 0.06008
## alternative hypothesis: stationaryadf.test(rl_dif)##
## Augmented Dickey-Fuller Test
##
## data: rl_dif
## Dickey-Fuller = -9.436, 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(rl ~ cepea + ipca + volume + dolar)
summary(model)##
## Call:
## lm(formula = rl ~ cepea + ipca + volume + dolar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3277226 -1020208 -26532 1122649 2779362
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.452e+07 2.077e+06 -11.802 < 2e-16 ***
## cepea 8.234e+05 4.872e+04 16.900 < 2e-16 ***
## ipca 3.116e+05 3.508e+05 0.888 0.375832
## volume 5.472e+03 3.437e+02 15.922 < 2e-16 ***
## dolar 1.508e+06 4.094e+05 3.684 0.000324 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1359000 on 145 degrees of freedom
## Multiple R-squared: 0.8479, Adjusted R-squared: 0.8437
## F-statistic: 202 on 4 and 145 DF, p-value: < 2.2e-16Agora vamos avaliar os resíduos e a multicolinearidade das variáveis.
vif(model)## cepea ipca volume dolar
## 1.564260 1.034749 1.010452 1.602036As 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* = 28.962, df = 10, p-value = 0.001264
##
## 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(rl_dif ~ cepea_dif + volume_dif + dolar_dif)
summary(model_dif)##
## Call:
## lm(formula = rl_dif ~ cepea_dif + volume_dif + dolar_dif)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4687548 -1509673 -97867 1060702 4512840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33313.3 155820.7 0.214 0.831
## cepea_dif 537154.0 441364.5 1.217 0.226
## volume_dif 5272.3 331.7 15.897 <2e-16 ***
## dolar_dif -1858622.1 2441646.7 -0.761 0.448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1897000 on 145 degrees of freedom
## Multiple R-squared: 0.642, Adjusted R-squared: 0.6346
## F-statistic: 86.66 on 3 and 145 DF, p-value: < 2.2e-16vif(model_dif)## cepea_dif volume_dif dolar_dif
## 1.007201 1.011925 1.019135res_dif <- model_dif$residuals
checkresiduals(res_dif)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 76.749, df = 10, p-value = 2.17e-12
##
## Model df: 0. Total lags used: 10adf.test(res_dif)##
## Augmented Dickey-Fuller Test
##
## data: res_dif
## Dickey-Fuller = -7.8891, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarymix <- cepea_dif * (1 + ipca_dif)
#md <- lm(rec_dif ~ mix + volume_dif + dolar_dif)
md <- lm(rl_dif ~ volume_dif)
summary(md)##
## Call:
## lm(formula = rl_dif ~ volume_dif)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4614960 -1509663 -116616 1134159 4563614
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45699.2 155332.1 0.294 0.769
## volume_dif 5296.3 329.6 16.069 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1896000 on 147 degrees of freedom
## Multiple R-squared: 0.6372, Adjusted R-squared: 0.6348
## F-statistic: 258.2 on 1 and 147 DF, p-value: < 2.2e-16res_md <- md$residuals
checkresiduals(res_md)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 76.66, df = 10, p-value = 2.259e-12
##
## Model df: 0. Total lags used: 10adf.test(res_md)##
## Augmented Dickey-Fuller Test
##
## data: res_md
## Dickey-Fuller = -6.4316, 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.
exogenas <- cbind(volume)
endogena <- rln <- length(rl)
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 : 3889.948
## Regression with ARIMA(0,1,0) errors : 3950.86
## Regression with ARIMA(1,1,0) errors : 3912.691
## Regression with ARIMA(0,1,1) errors : 3893.456
## Regression with ARIMA(0,1,0) errors : 3948.799
## Regression with ARIMA(1,1,2) errors : 3892.793
## Regression with ARIMA(2,1,1) errors : 3892.573
## Regression with ARIMA(3,1,2) errors : 3892.188
## Regression with ARIMA(2,1,3) errors : 3888.04
## Regression with ARIMA(1,1,3) errors : 3890.174
## Regression with ARIMA(3,1,3) errors : Inf
## Regression with ARIMA(2,1,4) errors : 3889.419
## Regression with ARIMA(1,1,4) errors : 3892.454
## Regression with ARIMA(3,1,4) errors : Inf
## Regression with ARIMA(2,1,3) errors : 3886.016
## Regression with ARIMA(1,1,3) errors : 3888.108
## Regression with ARIMA(2,1,2) errors : 3887.929
## Regression with ARIMA(3,1,3) errors : 3892.017
## Regression with ARIMA(2,1,4) errors : 3887.313
## Regression with ARIMA(1,1,2) errors : 3890.75
## Regression with ARIMA(1,1,4) errors : 3890.353
## Regression with ARIMA(3,1,2) errors : 3890.135
## Regression with ARIMA(3,1,4) errors : Inf
##
## Best model: Regression with ARIMA(2,1,3) 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 70849.34 1441174 1209092 -0.2372608 7.271453 0.2914418
## Test set 1553879.55 2162998 1790463 6.8555023 8.383905 0.4315765
## ACF1 Theil's U
## Training set -0.02929369 NA
## Test set 0.19024234 0.4893716Testando a estabilidade das previsões
f_arimax <- function(x, h, xreg, newxreg) {
forecast(Arima(x, order=c(2,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: 1847091 cat("MAE: " , MAE_h5 , "\t \n")## MAE: 1425768 cat("MAPE: ", MAPE_h5, "% \t \n")## MAPE: 8.6328 % 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: 1601131cat("MAE: ", MAE_h1, "\n")## MAE: 1265801cat("MAPE:", MAPE_h1, "%\n")## MAPE: 7.606996 %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: 1615295cat("MAE: ", MAE_h3, "\n")## MAE: 1276670cat("MAPE:", MAPE_h3, "%\n")## MAPE: 7.594499 %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.
arimax_treino <- Arima(y_treino, xreg = x_treino, order = c(2,1,2))summary(arimax_treino)## Series: y_treino
## Regression with ARIMA(2,1,2) errors
##
## Coefficients:
## ar1 ar2 ma1 ma2 xreg
## -1.0629 -0.3646 0.3004 -0.4046 5049.2526
## s.e. 0.1366 0.1137 0.1346 0.1185 365.3134
##
## sigma^2 = 2.263e+12: log likelihood = -1937.61
## AIC=3887.21 AICc=3887.93 BIC=3904.13
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 69342.99 1467699 1206764 -0.2416699 7.255559 0.2908806
## ACF1
## Training set 3.967226e-05coeftest(arimax_treino)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.06287 0.13656 -7.7834 7.06e-15 ***
## ar2 -0.36465 0.11373 -3.2063 0.0013446 **
## ma1 0.30044 0.13459 2.2323 0.0255983 *
## ma2 -0.40455 0.11852 -3.4134 0.0006416 ***
## xreg 5049.25257 365.31338 13.8217 < 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(2,1,2) errors
## Q* = 17.26, df = 21, p-value = 0.6952
##
## Model df: 4. 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 = -4.2077, 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.98345, p-value = 0.1307qqnorm(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("[Planta] 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 = "[Planta] 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,2))checkresiduals(arimax_final)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(2,1,2) errors
## Q* = 21.884, df = 26, p-value = 0.695
##
## Model df: 4. Total lags used: 30summary(arimax_final)## Series: endogena
## Regression with ARIMA(2,1,2) errors
##
## Coefficients:
## ar1 ar2 ma1 ma2 xreg
## -1.0504 -0.3320 0.3090 -0.4304 5271.026
## s.e. 0.1227 0.1077 0.1171 0.1056 328.546
##
## sigma^2 = 2.141e+12: log likelihood = -2324.55
## AIC=4661.09 AICc=4661.68 BIC=4679.11
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 113103.5 1433602 1170826 0.002940608 6.898036 0.3079363
## ACF1
## Training set -0.002933945summary(arimax_treino)## Series: y_treino
## Regression with ARIMA(2,1,2) errors
##
## Coefficients:
## ar1 ar2 ma1 ma2 xreg
## -1.0629 -0.3646 0.3004 -0.4046 5049.2526
## s.e. 0.1366 0.1137 0.1346 0.1185 365.3134
##
## sigma^2 = 2.263e+12: log likelihood = -1937.61
## AIC=3887.21 AICc=3887.93 BIC=3904.13
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 69342.99 1467699 1206764 -0.2416699 7.255559 0.2908806
## ACF1
## Training set 3.967226e-05As 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.
Testar
Testar fazer a modelagem da receita líquida da planta apenas com um modelo univariado como Box-Jenkins ou ETS