library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
library(ggplot2)
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'stringr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(nlme)
## Warning: package 'nlme' was built under R version 4.3.3
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:forecast':
##
## getResponse
##
## The following object is masked from 'package:dplyr':
##
## collapse
data <- read_excel("C:/Users/Lenovo/Documents/Sinar Matahari per Bulan di Kota Bandung (Persen).xlsx")
View(data)
summary(data)
## Bulan Tahun 2019 tahun 2020 Tahun 2021
## Length:12 Min. :48.00 Min. :33.66 Min. :38.50
## Class :character 1st Qu.:58.75 1st Qu.:49.92 1st Qu.:49.30
## Mode :character Median :68.00 Median :58.14 Median :54.95
## Mean :68.92 Mean :60.81 Mean :55.97
## 3rd Qu.:83.25 3rd Qu.:77.27 3rd Qu.:63.40
## Max. :86.00 Max. :81.50 Max. :77.90
## Tahun 2022 tahun 2023 Tahun 2024
## Min. :38.00 Min. :34.00 Min. :39.00
## 1st Qu.:45.58 1st Qu.:48.50 1st Qu.:48.50
## Median :54.35 Median :51.50 Median :63.00
## Mean :53.48 Mean :58.92 Mean :60.58
## 3rd Qu.:61.15 3rd Qu.:74.75 3rd Qu.:71.50
## Max. :67.70 Max. :90.00 Max. :78.00
sum(duplicated(data))
## [1] 0
colSums(is.na(data))
## Bulan Tahun 2019 tahun 2020 Tahun 2021 Tahun 2022 tahun 2023 Tahun 2024
## 0 0 0 0 0 0 0
data_long <- data %>%
pivot_longer(
cols = -Bulan,
names_to = "Tahun",
values_to = "Produksi"
) %>%
mutate(
# Ambil angka tahun dari nama kolom (misal "Tahun 2019" → 2019)
Tahun = str_extract(Tahun, "\\d+") |> as.numeric(),
# Ubah nama bulan ke angka (1 = Januari, dst)
Bulan = match(Bulan, month.name)
) %>%
arrange(Tahun, Bulan)
#Buat boxplot per tahun
ggplot(data_long, aes(x = as.factor(Tahun), y = Produksi, group = Tahun)) +
geom_boxplot(fill = "skyblue", color = "darkblue") +
labs(
title = "Boxplot Produksi per Tahun",
x = "Tahun",
y = "Produksi"
) +
theme_minimal()
#Boxplot per bulan
ggplot(data_long, aes(x = factor(Bulan, labels = month.name), y = Produksi)) +
geom_boxplot(fill = "lightgreen", color = "darkgreen") +
labs(
title = "Boxplot Produksi per Bulan (Gabungan Semua Tahun)",
x = "Bulan",
y = "Produksi"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#Boxplot keseluruhan
ggplot(data_long, aes(y = Produksi)) +
geom_boxplot(fill = "orange", color = "darkred") +
labs(
title = "Boxplot Produksi Keseluruhan (2019–2024)",
x = "",
y = "Produksi"
) +
theme_minimal()
# objek time series (2019 Jan - 2024 Des, frekuensi bulanan = 12)
ts_data <- ts(data_long$Produksi, start = c(2019, 1), frequency = 12)
plot(decompose(ts_data))
plot(ts_data, main = "plot awal")
# Uji Augmented Dickey-Fuller (ADF)
train <- window(ts_data, end = c(2023, 12))
test <- window(ts_data, start = c(2024, 1))
adf_test <- adf.test(train)
## Warning in adf.test(train): p-value smaller than printed p-value
adf_test
##
## Augmented Dickey-Fuller Test
##
## data: train
## Dickey-Fuller = -5.3257, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
diff_train <- diff(train)
# ACF dan PACF
par(mfrow = c(1, 2))
acf(diff_train, main = "ACF Plot")
pacf(diff_train, main = "PACF Plot")
#HOLT winter (data trend dan musiman)
hw_add <- hw(train, seasonal = "additive", h = length(test))
hw_mul <- hw(train, seasonal = "multiplicative", h = length(test))
accuracy(hw_add, test)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.5432959 9.703822 7.575665 -3.869459 13.82597 0.7090553
## Test set 10.4339870 13.682548 11.144196 15.359650 17.18070 1.0430571
## ACF1 Theil's U
## Training set 0.3519338 NA
## Test set -0.2713243 1.108512
accuracy(hw_mul, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2400055 9.637644 7.457130 -2.231937 13.34665 0.6979609
## Test set 7.3236643 11.453074 9.286543 9.739277 14.63761 0.8691874
## ACF1 Theil's U
## Training set 0.293569 NA
## Test set -0.264863 0.9176958
#holt winter multiplicative karena RMSE dan MAPE lebih kecil dibandingkan dengan additive
hw_model_mult <- HoltWinters(train, seasonal = "multiplicative")
hw_model_mult
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = train, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.00823339
## beta : 1
## gamma: 0.3385724
##
## Coefficients:
## [,1]
## a 53.4102778
## b 0.4554378
## s1 0.9476000
## s2 0.7856314
## s3 0.8845751
## s4 0.9697713
## s5 0.9938133
## s6 0.9439605
## s7 1.3092273
## s8 1.4054740
## s9 1.3126250
## s10 1.1345236
## s11 0.9220484
## s12 0.8523128
# Plot fitting model terhadap data training
plot(hw_model_mult, main = "Fitting Holt-Winters Multiplicative") #garis pink hasil peramalan
#melaukan peramalan untuk data testing
forecast_hw_mult <- forecast(hw_model_mult, h = length(test))
# Plot hasil forecast dan bandingkan dengan data aktual
autoplot(forecast_hw_mult) +
autolayer(test, series = "Actual", color = "green") +
ggtitle("Peramalan Holt–Winters Multiplicative") +
ylab("Produksi") +
xlab("Tahun") +
theme_minimal()
# Hitung RMSE
rmse_mult <- sqrt(mean((test - forecast_hw_mult$mean)^2))
# Hitung MAPE
mape_mult <- mean(abs((test - forecast_hw_mult$mean) / test)) * 100
# Tampilkan hasil evaluasi
cat("Evaluasi Model Holt–Winters (Multiplicative):\n")
## Evaluasi Model Holt–Winters (Multiplicative):
cat("RMSE =", round(rmse_mult, 2), "\n")
## RMSE = 10.5
cat("MAPE =", round(mape_mult, 2), "%\n")
## MAPE = 16.22 %
#menggunakan regresi time series karena data trend dan musiman
data_long$t <- 1:nrow(data_long) # variabel trend
data_long$bulan_factor <- as.factor(data_long$Bulan) # variabel musiman (bulan)
# Model regresi dengan trend dan musiman
model_ts <- lm(Produksi ~ t + bulan_factor, data = data_long)
summary(model_ts)
##
## Call:
## lm(formula = Produksi ~ t + bulan_factor, data = data_long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.5906 -7.5412 0.4657 5.7436 19.7132
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.27885 4.70318 11.754 < 2e-16 ***
## t -0.11867 0.06071 -1.955 0.055367 .
## bulan_factor2 -1.73799 6.09581 -0.285 0.776556
## bulan_factor3 -1.16265 6.09671 -0.191 0.849414
## bulan_factor4 5.42935 6.09823 0.890 0.376912
## bulan_factor5 10.53469 6.10034 1.727 0.089418 .
## bulan_factor6 10.72003 6.10306 1.757 0.084192 .
## bulan_factor7 25.20870 6.10638 4.128 0.000117 ***
## bulan_factor8 25.87904 6.11030 4.235 8.11e-05 ***
## bulan_factor9 22.14938 6.11483 3.622 0.000609 ***
## bulan_factor10 13.25139 6.11995 2.165 0.034422 *
## bulan_factor11 -0.29994 6.12567 -0.049 0.961113
## bulan_factor12 -3.97793 6.13198 -0.649 0.519037
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.56 on 59 degrees of freedom
## Multiple R-squared: 0.5534, Adjusted R-squared: 0.4626
## F-statistic: 6.093 on 12 and 59 DF, p-value: 8.405e-07
res <- resid(model_ts)
# Plot residual untuk melihat pola
plot(res, type = "l", col = "darkblue",
main = "Residual Regresi Time Series",
ylab = "Residuals", xlab = "Waktu")
# Uji autokorelasi residual (Durbin-Watson)
dwtest(model_ts)
##
## Durbin-Watson test
##
## data: model_ts
## DW = 1.3218, p-value = 0.003005
## alternative hypothesis: true autocorrelation is greater than 0
# Plot ACF residual untuk cek pola autokorelasi
acf(res, main = "ACF Residual Regresi Time Series")
#karena terdapat autokorelasi maka harus di atasi
# Kita fit model gls dengan korelasi AR(1) terhadap urutan waktu 't'
model_gls_ar1 <- gls(Produksi ~ t, correlation = corAR1(form = ~ t), data = data_long)
summary(model_gls_ar1)
## Generalized least squares fit by REML
## Model: Produksi ~ t
## Data: data_long
## AIC BIC logLik
## 563.1964 572.1904 -277.5982
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.6176659
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 63.55220 6.938891 9.158841 0.0000
## t -0.11544 0.163382 -0.706567 0.4822
##
## Correlation:
## (Intr)
## t -0.859
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.8901968 -0.7164912 -0.1072045 0.8726974 2.2001326
##
## Residual standard error: 14.95931
## Degrees of freedom: 72 total; 70 residual
# Residual GLS
res_gls <- residuals(model_gls_ar1, type = "normalized")
# Plot & ACF residual GLS
plot(res_gls, type = "l", col = "darkgreen",
main = "Residual GLS AR(1)",
ylab = "Residuals", xlab = "Waktu (t)")
acf(res_gls, main = "ACF Residual GLS AR(1)")
# Ljung-Box untuk residual GLS
Box.test(res_gls, lag = 12, type = "Ljung-Box", fitdf = length(coef(model_gls_ar1)))
##
## Box-Ljung test
##
## data: res_gls
## X-squared = 24.939, df = 10, p-value = 0.005463
# Untuk GLS: fitted values (predict)
data_long$Fitted_GLS <- fitted(model_gls_ar1)
# Hitung RMSE dan MAPE (hanya pada baris yg tidak NA pada fitted masing2)
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2, na.rm = TRUE))
mape <- function(obs, pred) mean(abs((obs - pred) / obs) * 100, na.rm = TRUE)
cat("Evaluasi Model (2019-2024):\n")
## Evaluasi Model (2019-2024):
cat("GLS - RMSE =", round(rmse(data_long$Produksi, data_long$Fitted_GLS), 3),
" | MAPE =", round(mape(data_long$Produksi, data_long$Fitted_GLS), 3), "%\n")
## GLS - RMSE = 14.167 | MAPE = 21.09 %
# Visualisasi: Actual vs Fitted (GLS)
plot_df <- data_long %>%
select(t, Produksi, Fitted_GLS) %>%
pivot_longer(cols = starts_with("Fitted"), names_to = "Model", values_to = "Fitted")
ggplot(plot_df, aes(x = t)) +
geom_line(aes(y = Produksi, color = "Aktual"), size = 1) +
geom_line(aes(y = Fitted, color = Model), size = 0.9, linetype = "dashed") +
scale_color_manual(values = c("Aktual" = "yellow",
"Fitted_GLS" = "darkgreen")) +
labs(title = "Actual vs Fitted (OLS, GLS-AR1)",
x = "Waktu (t)", y = "Produksi") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# MODEL ARIMA & SARIMA
# Buat data training dan testing kembali
train <- window(ts_data, end = c(2023, 12))
test <- window(ts_data, start = c(2024, 1))
# Plot awal untuk cek pola musiman
plot(train, main = "Data Training (2019–2023)", ylab = "Produksi", xlab = "Tahun")
# Differencing musiman (lag = 12)
diff_seasonal <- diff(train, lag = 12)
# Plot hasil differencing musiman
plot(diff_seasonal, main = "Differencing Musiman (lag = 12)",
ylab = "Differenced Produksi", xlab = "Tahun")
# Uji stasioneritas pada data yang sudah di-differencing
adf.test(na.omit(diff_seasonal))
## Warning in adf.test(na.omit(diff_seasonal)): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: na.omit(diff_seasonal)
## Dickey-Fuller = -4.6305, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
# ACF dan PACF untuk data yang sudah differencing
par(mfrow = c(1, 2))
acf(na.omit(diff_seasonal), main = "ACF Setelah Differencing (12)")
pacf(na.omit(diff_seasonal), main = "PACF Setelah Differencing (12)")
par(mfrow = c(1, 1))
# Tentukan model ARIMA
model_auto <- auto.arima(train, D = 1, seasonal = TRUE,
stepwise = FALSE, approximation = FALSE,
trace = TRUE)
##
## ARIMA(0,1,0)(0,1,0)[12] : 399.4989
## ARIMA(0,1,0)(0,1,1)[12] : 389.3651
## ARIMA(0,1,0)(1,1,0)[12] : 384.0676
## ARIMA(0,1,0)(1,1,1)[12] : 384.6412
## ARIMA(0,1,1)(0,1,0)[12] : 383.5969
## ARIMA(0,1,1)(0,1,1)[12] : 381.4137
## ARIMA(0,1,1)(1,1,0)[12] : 381.3311
## ARIMA(0,1,1)(1,1,1)[12] : 383.4443
## ARIMA(0,1,2)(0,1,0)[12] : 385.2921
## ARIMA(0,1,2)(0,1,1)[12] : 381.036
## ARIMA(0,1,2)(1,1,0)[12] : 379.9841
## ARIMA(0,1,2)(1,1,1)[12] : 382.2802
## ARIMA(0,1,3)(0,1,0)[12] : 385.0766
## ARIMA(0,1,3)(0,1,1)[12] : 381.6265
## ARIMA(0,1,3)(1,1,0)[12] : 380.669
## ARIMA(0,1,3)(1,1,1)[12] : 383.0933
## ARIMA(0,1,4)(0,1,0)[12] : 384.616
## ARIMA(0,1,4)(0,1,1)[12] : 382.1461
## ARIMA(0,1,4)(1,1,0)[12] : 380.905
## ARIMA(0,1,5)(0,1,0)[12] : 384.3008
## ARIMA(1,1,0)(0,1,0)[12] : 389.5613
## ARIMA(1,1,0)(0,1,1)[12] : 384.7863
## ARIMA(1,1,0)(1,1,0)[12] : 382.8023
## ARIMA(1,1,0)(1,1,1)[12] : 384.3958
## ARIMA(1,1,1)(0,1,0)[12] : 385.0839
## ARIMA(1,1,1)(0,1,1)[12] : 380.2037
## ARIMA(1,1,1)(1,1,0)[12] : 378.5873
## ARIMA(1,1,1)(1,1,1)[12] : 380.6633
## ARIMA(1,1,2)(0,1,0)[12] : 386.3088
## ARIMA(1,1,2)(0,1,1)[12] : 383.7197
## ARIMA(1,1,2)(1,1,0)[12] : 381.0827
## ARIMA(1,1,2)(1,1,1)[12] : 383.2974
## ARIMA(1,1,3)(0,1,0)[12] : 388.8198
## ARIMA(1,1,3)(0,1,1)[12] : 385.936
## ARIMA(1,1,3)(1,1,0)[12] : 382.6101
## ARIMA(1,1,4)(0,1,0)[12] : 385.3885
## ARIMA(2,1,0)(0,1,0)[12] : 387.0342
## ARIMA(2,1,0)(0,1,1)[12] : 385.5068
## ARIMA(2,1,0)(1,1,0)[12] : 384.4291
## ARIMA(2,1,0)(1,1,1)[12] : 386.3231
## ARIMA(2,1,1)(0,1,0)[12] : 386.6892
## ARIMA(2,1,1)(0,1,1)[12] : 382.4589
## ARIMA(2,1,1)(1,1,0)[12] : 381.0724
## ARIMA(2,1,1)(1,1,1)[12] : 383.2959
## ARIMA(2,1,2)(0,1,0)[12] : 388.8195
## ARIMA(2,1,2)(0,1,1)[12] : 385.0786
## ARIMA(2,1,2)(1,1,0)[12] : 383.6928
## ARIMA(2,1,3)(0,1,0)[12] : Inf
## ARIMA(3,1,0)(0,1,0)[12] : 387.2514
## ARIMA(3,1,0)(0,1,1)[12] : 387.4847
## ARIMA(3,1,0)(1,1,0)[12] : 386.8387
## ARIMA(3,1,0)(1,1,1)[12] : 388.9414
## ARIMA(3,1,1)(0,1,0)[12] : 386.1482
## ARIMA(3,1,1)(0,1,1)[12] : 387.441
## ARIMA(3,1,1)(1,1,0)[12] : 382.7638
## ARIMA(3,1,2)(0,1,0)[12] : 386.3502
## ARIMA(4,1,0)(0,1,0)[12] : 378.4864
## ARIMA(4,1,0)(0,1,1)[12] : 379.7115
## ARIMA(4,1,0)(1,1,0)[12] : 379.7309
## ARIMA(4,1,1)(0,1,0)[12] : 378.9049
## ARIMA(5,1,0)(0,1,0)[12] : 380.4918
##
##
##
## Best model: ARIMA(4,1,0)(0,1,0)[12]
summary(model_auto)
## Series: train
## ARIMA(4,1,0)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4
## -0.4469 -0.2491 -0.0494 -0.5044
## s.e. 0.1287 0.1615 0.1574 0.1363
##
## sigma^2 = 151.6: log likelihood = -183.51
## AIC=377.02 AICc=378.49 BIC=386.27
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.130434 10.4231 7.251782 0.7133535 12.67464 0.678741 -0.07150646
# Diagnostik residual
checkresiduals(model_auto)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(4,1,0)(0,1,0)[12]
## Q* = 14.038, df = 8, p-value = 0.08079
##
## Model df: 4. Total lags used: 12
# Forecast ke periode testing (2024)
forecast_arima <- forecast(model_auto, h = length(test))
# Evaluasi performa model
rmse_arima <- sqrt(mean((forecast_arima$mean - test)^2))
mape_arima <- mean(abs((forecast_arima$mean - test) / test)) * 100
cat("\nEvaluasi Model ARIMA:\n")
##
## Evaluasi Model ARIMA:
cat("RMSE =", round(rmse_arima, 3), "\n")
## RMSE = 20.125
cat("MAPE =", round(mape_arima, 3), "%\n")
## MAPE = 32.254 %
# Visualisasi hasil forecast ARIMA
autoplot(forecast_arima) +
autolayer(test, series = "Data Aktual", color = "green") +
ggtitle("Peramalan Menggunakan ARIMA (dengan differencing musiman 12)") +
ylab("Produksi") +
xlab("Tahun") +
theme_minimal()
# SARIMA (Seasonal ARIMA)
# Misalnya kita tentukan SARIMA(1,1,1)(1,1,1)[12] sebagai contoh manual
model_sarima <- Arima(train,
order = c(1,1,1), # komponen ARIMA non-musiman (p,d,q)
seasonal = list(order = c(1,1,1), period = 12)) # komponen musiman (P,D,Q)[s]
summary(model_sarima)
## Series: train
## ARIMA(1,1,1)(1,1,1)[12]
##
## Coefficients:
## ar1 ma1 sar1 sma1
## 0.4870 -0.9322 -0.7514 0.3092
## s.e. 0.1679 0.0877 0.2802 0.4277
##
## sigma^2 = 144.6: log likelihood = -184.6
## AIC=379.2 AICc=380.66 BIC=388.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.597271 10.17823 7.190001 1.151485 12.74855 0.6729585 -0.03723596
checkresiduals(model_sarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,1,1)[12]
## Q* = 17.826, df = 8, p-value = 0.02257
##
## Model df: 4. Total lags used: 12
# Forecast SARIMA
forecast_sarima <- forecast(model_sarima, h = length(test))
# Evaluasi SARIMA
rmse_sarima <- sqrt(mean((forecast_sarima$mean - test)^2))
mape_sarima <- mean(abs((forecast_sarima$mean - test) / test)) * 100
cat("\nEvaluasi Model SARIMA:\n")
##
## Evaluasi Model SARIMA:
cat("RMSE =", round(rmse_sarima, 3), "\n")
## RMSE = 13.965
cat("MAPE =", round(mape_sarima, 3), "%\n")
## MAPE = 20.825 %
# Plot hasil SARIMA
autoplot(forecast_sarima) +
autolayer(test, series = "Data Aktual", color = "green") +
ggtitle("Peramalan Menggunakan SARIMA(1,1,1)(1,1,1)[12]") +
ylab("Produksi") +
xlab("Tahun") +
theme_minimal()
# Bandingkan performa model
cat("\nPerbandingan Model:\n")
##
## Perbandingan Model:
cat("Holt–Winters (Multiplicative): RMSE =", round(rmse_mult, 3),
"| MAPE =", round(mape_mult, 3), "%\n")
## Holt–Winters (Multiplicative): RMSE = 10.502 | MAPE = 16.22 %
cat("ARIMA :", round(rmse_arima, 3),
"| MAPE =", round(mape_arima, 3), "%\n")
## ARIMA : 20.125 | MAPE = 32.254 %
cat("SARIMA :", round(rmse_sarima, 3),
"| MAPE =", round(mape_sarima, 3), "%\n")
## SARIMA : 13.965 | MAPE = 20.825 %