PERSIAPAN DATA

Data yang digunakan adalah data pengeluaran alat kontrasepsi suntik 3 bulan progestin sebayak 50 observasi, mencakup periode Februari 2022 hingga Maret 2026 (data bulanan). Data kemudian dikonversi menjadi objek time series dengan frekuensi 12 (bulanan).

# Library
library(readxl)
library(forecast)

## Warning: package 'forecast' was built under R version 4.5.1

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(tseries)

## Warning: package 'tseries' was built under R version 4.5.1

library(ggplot2)

# Fungsi untuk menghitung MSE langsung dari data aktual dan hasil forecast
# (tidak perlu mengonversi dari RMSE)
hitung_mse <- function(actual, forecast_obj) {
  mean((actual - forecast_obj$mean)^2)
}

# Import data
data <- read_excel("C:/Users/LENOVO/Downloads/alokon.xlsx")
data

# Tampilkan tabel data secara rapi (Tabel 1)
knitr::kable(data, caption = "Tabel 1. Data Pengeluaran Alokon Suntik 3 Bulan Progestin Februari 2022 - Maret 2026")

Tabel 1. Data Pengeluaran Alokon Suntik 3 Bulan Progestin Februari 2022 - Maret 2026
Bulan	Suntik 3 Bulan
2022-02-01	2233
2022-03-01	2597
2022-04-01	2859
2022-05-01	2797
2022-06-01	2545
2022-07-01	2681
2022-08-01	2724
2022-09-01	2773
2022-10-01	2845
2022-11-01	2533
2022-12-01	2575
2023-01-01	3858
2023-02-01	3721
2023-03-01	4580
2023-04-01	4231
2023-05-01	4781
2023-06-01	4726
2023-07-01	4951
2023-08-01	5532
2023-09-01	5504
2023-10-01	5283
2023-11-01	5782
2023-12-01	5227
2024-01-01	6440
2024-02-01	6063
2024-03-01	5995
2024-04-01	6458
2024-05-01	6146
2024-06-01	6149
2024-07-01	6415
2024-08-01	6071
2024-09-01	6307
2024-10-01	6340
2024-11-01	5519
2024-12-01	6761
2025-01-01	6650
2025-02-01	6097
2025-03-01	6327
2025-04-01	6224
2025-05-01	6592
2025-06-01	6427
2025-07-01	5595
2025-08-01	5863
2025-09-01	6612
2025-10-01	6823
2025-11-01	6667
2025-12-01	6148
2026-01-01	6145
2026-02-01	6291
2026-03-01	6525

# Ubah kolom bulan jadi date
data$Bulan <- as.Date(paste0("01-", data$Bulan), format="%d-%b-%y")

# Buat time series (Feb 2022 – Mar 2026 = bulanan)
ts_data <- ts(data$`Suntik 3 Bulan`, start=c(2022,2), frequency=12)

# Statistik Deskriptif
summary(ts_data)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2233    3951    5822    5140    6322    6823

sd(ts_data)

## [1] 1513.968

statdesk <- data.frame(
  Rata_rata = mean(ts_data),
  Median    = median(ts_data),
  Minimum   = min(ts_data),
  Maksimum  = max(ts_data),
  Std_Dev   = sd(ts_data)
)

knitr::kable(statdesk, caption = "Tabel 2. Statistik Deskriptif Data Pengeluaran Alokon Suntik 3 Bulan Progestin", digits = 2)

Tabel 2. Statistik Deskriptif Data Pengeluaran Alokon Suntik 3 Bulan Progestin
Rata_rata	Median	Minimum	Maksimum	Std_Dev
5139.76	5822.5	2233	6823	1513.97

plot(ts_data, main="Data Pengeluaran Suntik 3 Bulan", ylab="Jumlah", xlab="Waktu")

Gambar 1. Plot Data Pengeluaran Alokon Suntik 3 Bulan Progestin Februari 2022 - Maret 2026

Dari plot data keseluruhan, terlihat adanya tren naik yang jelas sepanjang periode pengamatan. Nilai awal berada di kisaran 2.000-3.000 pada awal 2022, kemudian meningkat drastis mulai pertengahan 2023 dan stabil di kisaran 5.500-6.800 pada 2024-026. Pola ini mengindikasikan data tidak stasioner dan mengandung tren positif.

# PARTISI DATA
training <- window(ts_data, end=c(2024,12))
testing  <- window(ts_data, start=c(2025,1))

plot(training, main="Data Training")

Gambar 2. Plot Data Training (Februari 2022 - Desember 2024)

plot(testing, main="Data Testing")

Gambar 3. Plot Data Testing (Januari 2025 - Maret 2026)

Data dibagi menjadi dua bagian: - Data Training: Februari 2022 - Desember 2024 (36 observasi), digunakan untuk membangun model - Data Testing: Januari 2025 - Maret 2026 (15 observasi), digunakan untuk untuk mengevaluasi performa model

Dari plot data training, terlihat tren naik yang konsisten dari ~2.200 hingga ~6.500. Dari plot data testing, nilai berfluktuasi di kisaran 5.600 - 6.900, menunjukkan data masih dalam tren yang relatif stabil di level tinggi.

AUTOREGRESSIVE INTEGRATED MOVING AVERAGE

acf(training, main="ACF Data Training")

Gambar 4. Plot ACF Data Training

Plot ACF menunjukkan autokorelasi yang turun sangat lambat dan tetap signifikan hingga lag jauh. Ini merupakan indikasi kuat bahwa data tidak stasioner dalam mean, sehingga diperlukan proses differencing.

# Uji ADF
adf.test(training)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training
## Dickey-Fuller = -1.5762, Lag order = 3, p-value = 0.7376
## alternative hypothesis: stationary

Hasil uji Augmented Dickey-Fuller menghasilkan nilai p-value = 0,7376 > 0,05, sehingga gagal menolak H₀. Artinya, data training belum stasioner dan harus dilakukan differencing.

# Differencing Ordo 1
diff1 <- diff(training, differences=1)

plot(diff1, main="Differencing Orde 1")

Gambar 5. Plot Differencing Orde 1

acf(diff1) # Cek Kestasioneran

Gambar 6. Plot ACF dan PACF Differencing Orde 1

pacf(diff1)

Gambar 6. Plot ACF dan PACF Differencing Orde 1

# Tabel 3. Hasil Uji ADF Setelah Differencing Orde 1
adf.test(diff1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1
## Dickey-Fuller = -2.9694, Lag order = 3, p-value = 0.1974
## alternative hypothesis: stationary

Setelah dilakukan differencing orde 1, plot menunjukkan data sudah lebih berfluktuasi di sekitar nol. Namun, uji ADF menghasilkan p-value = 0,1974 > 0,05, yang berarti data masih belum stasioner pada orde 1. Oleh karena itu, dilakukan differencing orde 2.

# Differencing Ordo 2 (karena belum stasioner pada ordo 1)
diff2 <- diff(training, differences=2)

plot(diff2, main="Differencing Orde 2")

Gambar 7. Plot Differencing Orde 2

acf(diff2) # Cek Kestasioneran

Gambar 8. Plot ACF dan PACF Differencing Orde 2

pacf(diff2)

Gambar 8. Plot ACF dan PACF Differencing Orde 2

# Tabel 4. Hasil Uji ADF Setelah Differencing Orde 2
adf.test(diff2)

## Warning in adf.test(diff2): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff2
## Dickey-Fuller = -4.4544, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Setelah differencing orde 2, uji ADF menghasilkan p-value = 0,01 ≤ 0,05, sehingga H₀ ditolak. Data sudah stasioner pada differencing orde 2, yang berarti nilai d = 2 untuk model ARIMA.

acf(diff2, main="ACF diff2")

Gambar 9. Plot ACF Data Differencing Orde 2 (Identifikasi Model)

pacf(diff2, main="PACF diff2")

Gambar 10. Plot PACF Data Differencing Orde 2 (Identifikasi Model)

Karena keterbatasan jumlah data, metode EACF tidak digunakan. Identifikasi dilakukan berdasarkan pola ACF dan PACF setelah differencing orde 2: - ACF diff2: Terdapat 1 lag signifikan yang dominan di awal, lalu cut off -> mengindikasikan q = 1 (komponen MA) - PACF diff2: Terdapat beberapa lag signifikan di awal -> mengindikasikan p = 1 atau 2 (komponen AR)

Berdasarkan hal tersebut, kandidat model tentatif yang diusulkan adalah: ARIMA(1,2,0), ARIMA(1,2,1), ARIMA(0,2,1), dan ARIMA(2,2,0).

# Kandidat Model 
arima120 <- arima(training, order=c(1,2,0))
arima121 <- arima(training, order=c(1,2,1))
arima021 <- arima(training, order=c(0,2,1))
arima220 <- arima(training, order=c(2,2,0))

# Uji Signifikansi Parameter
arima120

## 
## Call:
## arima(x = training, order = c(1, 2, 0))
## 
## Coefficients:
##           ar1
##       -0.7650
## s.e.   0.1248
## 
## sigma^2 estimated as 330724:  log likelihood = -256.96,  aic = 517.93

lmtest::coeftest(arima120)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.76501    0.12483 -6.1285 8.869e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima121

## 
## Call:
## arima(x = training, order = c(1, 2, 1))
## 
## Coefficients:
##           ar1      ma1
##       -0.5481  -1.0000
## s.e.   0.1572   0.2038
## 
## sigma^2 estimated as 173795:  log likelihood = -248.28,  aic = 502.55

lmtest::coeftest(arima121)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.54811    0.15718 -3.4872  0.000488 ***
## ma1 -1.00000    0.20375 -4.9080 9.202e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima021

## 
## Call:
## arima(x = training, order = c(0, 2, 1))
## 
## Coefficients:
##          ma1
##       -1.000
## s.e.   0.081
## 
## sigma^2 estimated as 242720:  log likelihood = -253.18,  aic = 510.37

lmtest::coeftest(arima021)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.999998   0.080976 -12.349 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima220

## 
## Call:
## arima(x = training, order = c(2, 2, 0))
## 
## Coefficients:
##           ar1      ar2
##       -1.2237  -0.6223
## s.e.   0.1499   0.1473
## 
## sigma^2 estimated as 215655:  log likelihood = -250.38,  aic = 506.76

lmtest::coeftest(arima220)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -1.22365    0.14994 -8.1608 3.327e-16 ***
## ar2 -0.62233    0.14734 -4.2238 2.402e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Tabel 5. Ringkasan Uji Signifikansi Parameter Kandidat Model
tabel_signifikansi <- data.frame(
  Model = c("ARIMA(1,2,0)", "ARIMA(1,2,1)", "ARIMA(1,2,1)", "ARIMA(0,2,1)", "ARIMA(2,2,0)", "ARIMA(2,2,0)"),
  Parameter = c("ar1", "ar1", "ma1", "ma1", "ar1", "ar2"),
  Estimate = c(-0.7650, -0.5481, -1.0000, -1.0000, -1.2237, -0.6223),
  Std_Error = c(0.1248, 0.1572, 0.2038, 0.0810, 0.1499, 0.1473),
  p_value = c(8.869e-10, 0.000488, 9.202e-07, 2.2e-16, 3.327e-16, 2.402e-05)
)
knitr::kable(tabel_signifikansi, caption = "Tabel 5. Ringkasan Uji Signifikansi Parameter Kandidat Model ARIMA", digits = 4)

Tabel 5. Ringkasan Uji Signifikansi Parameter Kandidat Model ARIMA
Model	Parameter	Estimate	Std_Error	p_value
ARIMA(1,2,0)	ar1	-0.7650	0.1248	0e+00
ARIMA(1,2,1)	ar1	-0.5481	0.1572	5e-04
ARIMA(1,2,1)	ma1	-1.0000	0.2038	0e+00
ARIMA(0,2,1)	ma1	-1.0000	0.0810	0e+00
ARIMA(2,2,0)	ar1	-1.2237	0.1499	0e+00
ARIMA(2,2,0)	ar2	-0.6223	0.1473	0e+00

Seluruh kandidat model memiliki parameter yang signifikan secara statistik (p-value < 0,05).

# Bandingkan AIC
tabel_aic <- data.frame(
  Model=c("ARIMA(1,2,0)","ARIMA(1,2,1)","ARIMA(0,2,1)","ARIMA(2,2,0)"),
  AIC=c(arima120$aic, arima121$aic, arima021$aic, arima220$aic)
)
knitr::kable(tabel_aic, caption = "Tabel 6. Perbandingan Nilai AIC Kandidat Model ARIMA", digits = 4)

Tabel 6. Perbandingan Nilai AIC Kandidat Model ARIMA
Model	AIC
ARIMA(1,2,0)	517.9282
ARIMA(1,2,1)	502.5525
ARIMA(0,2,1)	510.3651
ARIMA(2,2,0)	506.7584

Model ARIMA(1,2,1) memiliki AIC terkecil (502,55), sehingga dipilih sebagai model terbaik.

# Diagnostik Model
model_terbaik <- forecast::Arima(training, order=c(1,2,1))
res <- residuals(model_terbaik)

par(mfrow=c(2,2))
qqnorm(res); qqline(res)
plot(res, main="Residual")
acf(res)
pacf(res)

Gambar 11. Plot Diagnostik Residual Model ARIMA(1,2,1)

Normal Q-Q Plot: Titik-titik sebagian besar mengikuti garis diagonal, meskipun terdapat sedikit penyimpangan di ekor distribusi. Hal ini dikonfirmasi oleh uji Shapiro-Wilk.
Plot Residual: Residual berfluktuasi acak di sekitar nol tanpa pola tertentu, mengindikasikan model sudah cukup baik menangkap struktur data.
ACF & PACF Residual: Tidak ada lag yang signifikan melewati batas kepercayaan, menunjukkan residual bersifat white noise.

# Uji Non-Autokorelasi dengan Ljung-Box test
ljung <- Box.test(res, type="Ljung")
ljung

## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 0.22795, df = 1, p-value = 0.633

# Uji Normalitas dengan Shapiro-Wilk test
shapiro <- shapiro.test(res)
shapiro

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.94487, p-value = 0.07885

tabel_diagnostik <- data.frame(
  Uji = c("Ljung-Box (non-autokorelasi)", "Shapiro-Wilk (normalitas)"),
  Statistik = c(round(ljung$statistic, 3), round(shapiro$statistic, 3)),
  p_value = c(round(ljung$p.value, 3), round(shapiro$p.value, 3)),
  Keputusan = c("Gagal tolak H0 - tidak ada autokorelasi", "Gagal tolak H0 - residual normal")
)
knitr::kable(tabel_diagnostik, caption = "Tabel 7. Hasil Uji Diagnostik Residual Model ARIMA(1,2,1)")

Tabel 7. Hasil Uji Diagnostik Residual Model ARIMA(1,2,1)
	Uji	Statistik	p_value	Keputusan
X-squared	Ljung-Box (non-autokorelasi)	0.228	0.633	Gagal tolak H0 - tidak ada autokorelasi
W	Shapiro-Wilk (normalitas)	0.945	0.079	Gagal tolak H0 - residual normal

Seluruh asumsi diagnostik terpenuhi, sehingga model ARIMA(1,2,1) valid dan layak digunakan untuk peramalan.

forecast_arima <- forecast(model_terbaik, h=length(testing))

plot(forecast_arima)

Gambar 12. Plot Hasil Peramalan ARIMA(1,2,1) pada Data Testing

akurasi_arima_full <- accuracy(forecast_arima, testing)

mse_arima_training <- mean(residuals(model_terbaik)^2)
mse_arima_testing  <- hitung_mse(testing, forecast_arima)

tabel_akurasi_arima <- data.frame(
  Data = c("Training", "Testing"),
  MAPE = akurasi_arima_full[, "MAPE"],
  MSE  = c(mse_arima_training, mse_arima_testing),
  MAE  = akurasi_arima_full[, "MAE"]
)
knitr::kable(tabel_akurasi_arima, caption = "Tabel 8. Nilai Akurasi Model ARIMA(1,2,1)", digits = 2)

Tabel 8. Nilai Akurasi Model ARIMA(1,2,1)
	Data	MAPE	MSE	MAE
Training set	Training	6.29	163864.3	284.04
Test set	Testing	16.44	1330873.6	1024.45

MAPE data training = 6,29% -> akurasi baik
MAPE data testing = 16,44% -> akurasi cukup, namun terdapat kesenjangan antara training dan testing yang menunjukkan model kurang optimal dalam menggeneralisasi ke data baru

model_final <- forecast::Arima(ts_data, order=c(1,2,1))
forecast_final <- forecast(model_final, h=9)

knitr::kable(as.data.frame(forecast_final), caption = "Tabel 9. Hasil Peramalan ARIMA(1,2,1) April - Desember 2026", digits = 2)

Tabel 9. Hasil Peramalan ARIMA(1,2,1) April - Desember 2026
	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Apr 2026	6535.25	5956.42	7114.09	5650.00	7420.51
May 2026	6623.77	5922.10	7325.45	5550.65	7696.89
Jun 2026	6684.91	5834.46	7535.37	5384.26	7985.57
Jul 2026	6755.63	5783.34	7727.92	5268.64	8242.62
Aug 2026	6823.00	5732.18	7913.83	5154.73	8491.27
Sep 2026	6891.54	5688.94	8094.15	5052.32	8730.77
Oct 2026	6959.67	5648.64	8270.71	4954.62	8964.73
Nov 2026	7027.94	5611.43	8444.46	4861.58	9194.31
Dec 2026	7096.17	5576.22	8616.12	4771.60	9420.73

plot(forecast_final)

Gambar 13. Plot Hasil Peramalan ARIMA(1,2,1) April - Desember 2026

Model ARIMA(1,2,1) dilatih ulang menggunakan seluruh data (50 observasi) dan digunakan untuk meramalkan 9 bulan ke depan. Hasil forecast menunjukkan tren kenaikan yang terus berlanjut dari ~6.535 (April 2026) hingga ~7.096 (Desember 2026). Interval kepercayaan semakin melebar seiring jauhnya horizon peramalan, mencerminkan ketidakpastian yang meningkat untuk prediksi jangka panjang.

DOUBLE EXPONENTIAL SMOOTHING

# DES 1 (alpha=0.2, beta=0.2)
des1 <- HoltWinters(training, beta=0.2, alpha=0.2, gamma=FALSE)
forecast_des1 <- forecast(des1, h=length(testing))

plot(forecast_des1)

Gambar 14. Plot Hasil Peramalan DES 1 (alpha=0.2; beta=0.2) pada Data Testing

DES dengan parameter α=0,2 dan β=0,2 memberikan bobot yang rendah pada data terbaru, sehingga forecast cenderung lebih halus namun lebih lambat beradaptasi terhadap perubahan tren. Dari plot forecast, hasil prediksi relatif stabil dengan tren naik yang landai.

# DES 2 (alpha=0.6, beta=0.3)
des2 <- HoltWinters(training, beta=0.3, alpha=0.6, gamma=FALSE)
forecast_des2 <- forecast(des2, h=length(testing))

plot(forecast_des2)

Gambar 15. Plot Hasil Peramalan DES 2 (alpha=0.6; beta=0.3) pada Data Testing

DES dengan α=0,6 dan β=0,3 memberikan bobot lebih tinggi pada data terbaru, sehingga lebih responsif terhadap perubahan tren. Namun hal ini juga membuat interval kepercayaan jauh lebih lebar dibanding DES 1, menunjukkan ketidakpastian prediksi yang lebih besar.

# Alpha dan Beta Optimum
des_opt <- HoltWinters(training, gamma=FALSE)
des_opt

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = training, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.5936537
##  beta : 0.1661849
##  gamma: FALSE
## 
## Coefficients:
##         [,1]
## a 6390.15753
## b   59.06808

forecast_des_opt <- forecast(des_opt, h=length(testing))

plot(forecast_des_opt)

Gambar 16. Plot Hasil Peramalan DES Optimum pada Data Testing

Parameter optimal yang diperoleh secara otomatis adalah α = 0,5937 dan β = 0,1662. Nilai α yang mendekati 0,6 menunjukkan model memberikan bobot cukup tinggi pada data terkini untuk estimasi level, sementara β yang relatif kecil menunjukkan estimasi tren bersifat lebih stabil dan konservatif. Nilai akhir level (a) = 6.390 dan slope (b) = 59,07 per bulan, artinya setiap bulan diprediksi terjadi kenaikan sekitar 59 unit.

# Akurasi DES (testing) - MSE dihitung langsung, konsisten dengan Tabel 12
akurasi_des1 <- accuracy(forecast_des1, testing)[2,]
akurasi_des2 <- accuracy(forecast_des2, testing)[2,]
akurasi_des_opt <- accuracy(forecast_des_opt, testing)[2,]

tabel_akurasi_des <- data.frame(
  Model = c("DES 1", "DES 2", "DES Optimum"),
  MAPE = c(akurasi_des1["MAPE"], akurasi_des2["MAPE"], akurasi_des_opt["MAPE"]),
  MSE  = c(
    hitung_mse(testing, forecast_des1),
    hitung_mse(testing, forecast_des2),
    hitung_mse(testing, forecast_des_opt)
  ),
  MAE  = c(akurasi_des1["MAE"], akurasi_des2["MAE"], akurasi_des_opt["MAE"])
)
knitr::kable(tabel_akurasi_des, caption = "Tabel 10. Perbandingan Akurasi Tiga Varian Model DES pada Data Testing", digits = 2)

Tabel 10. Perbandingan Akurasi Tiga Varian Model DES pada Data Testing
Model	MAPE	MSE	MAE
DES 1	6.64	256172.2	403.83
DES 2	9.11	450451.9	560.76
DES Optimum	9.05	440765.0	557.07

Berdasarkan Tabel 10, DES 1 menghasilkan nilai MAPE, MSE, dan MAE yang paling kecil di antara ketiga varian model DES pada data testing, sehingga dipilih sebagai model DES terbaik. Hal ini menunjukkan bahwa parameter yang lebih konservatif (α = 0,2; β = 0,2) menghasilkan peramalan yang lebih stabil dan lebih baik dalam menggeneralisasi ke data baru dibandingkan parameter yang lebih responsif (DES 2) maupun parameter optimum yang dipilih otomatis oleh algoritma.

# Forecast DES April - Desember 2026
des_final <- HoltWinters(ts_data, gamma=FALSE)
forecast_des_final <- forecast(des_final, h=9)

knitr::kable(as.data.frame(forecast_des_final), caption = "Tabel 11. Hasil Peramalan DES Optimum April - Desember 2026", digits = 2)

Tabel 11. Hasil Peramalan DES Optimum April - Desember 2026
	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Apr 2026	6435.12	5858.96	7011.27	5553.96	7316.27
May 2026	6451.95	5760.95	7142.96	5395.16	7508.75
Jun 2026	6468.79	5656.63	7280.95	5226.70	7710.89
Jul 2026	6485.63	5546.26	7425.01	5048.98	7922.28
Aug 2026	6502.47	5430.07	7574.87	4862.38	8142.56
Sep 2026	6519.31	5308.30	7730.31	4667.23	8371.38
Oct 2026	6536.15	5181.14	7891.15	4463.85	8608.44
Nov 2026	6552.98	5048.78	8057.18	4252.51	8853.46
Dec 2026	6569.82	4911.39	8228.25	4033.47	9106.17

plot(forecast_des_final)

Gambar 17. Plot Hasil Peramalan DES Optimum April - Desember 2026

Model DES optimal dilatih ulang menggunakan seluruh 50 observasi untuk meramalkan 9 bulan ke depan. Hasil forecast menunjukkan tren kenaikan yang sangat landai dan stabil, dari ~6.435 (April 2026) hingga ~6.570 (Desember 2026), dengan kenaikan sekitar 17 unit per bulan. Dibandingkan ARIMA, interval kepercayaan DES lebih sempit di awal periode namun melebar secara moderat seiring bertambahnya horizon waktu, mencerminkan tingkat keyakinan prediksi yang lebih baik dalam jangka pendek.

PERBANDINGAN ARIMA VS DES

akurasi_arima_testing <- accuracy(forecast_arima, testing)[2,]

tabel_perbandingan <- data.frame(
  Model = c("ARIMA(1,2,1)", "DES 1", "DES 2", "DES Optimum"),
  MAPE = c(
    akurasi_arima_testing["MAPE"],
    accuracy(forecast_des1, testing)[2,"MAPE"],
    accuracy(forecast_des2, testing)[2,"MAPE"],
    accuracy(forecast_des_opt, testing)[2,"MAPE"]
  ),
  MSE = c(
    hitung_mse(testing, forecast_arima),
    hitung_mse(testing, forecast_des1),
    hitung_mse(testing, forecast_des2),
    hitung_mse(testing, forecast_des_opt)
  ),
  MAE = c(
    akurasi_arima_testing["MAE"],
    accuracy(forecast_des1, testing)[2,"MAE"],
    accuracy(forecast_des2, testing)[2,"MAE"],
    accuracy(forecast_des_opt, testing)[2,"MAE"]
  )
)
knitr::kable(tabel_perbandingan, caption = "Tabel 12. Perbandingan Akurasi ARIMA dan DES Berdasarkan Nilai MAPE, MSE, dan MAE pada Data Testing", digits = 2)

Tabel 12. Perbandingan Akurasi ARIMA dan DES Berdasarkan Nilai MAPE, MSE, dan MAE pada Data Testing
Model	MAPE	MSE	MAE
ARIMA(1,2,1)	16.44	1330873.6	1024.45
DES 1	6.64	256172.2	403.83
DES 2	9.11	450451.9	560.76
DES Optimum	9.05	440765.0	557.07

autoplot(ts_data) +
  autolayer(forecast_final$mean, series="ARIMA") +
  autolayer(forecast_des_final$mean, series="DES") +
  ggtitle("Perbandingan Forecast ARIMA vs DES")

Gambar 18. Plot Perbandingan Hasil Peramalan ARIMA vs DES April - Desember 2026

Berdasarkan perbandingan nilai MAPE pada data testing, DES 1 dengan α=0,2 dan β=0,2 merupakan model terbaik dengan MAPE terendah sebesar 6,64%. Seluruh model DES mengungguli ARIMA(1,2,1) dalam hal akurasi prediksi pada data testing. Dari plot perbandingan forecast, ARIMA menghasilkan proyeksi dengan tren kenaikan yang lebih curam hingga ~7.096 di akhir 2026, sementara DES menghasilkan proyeksi yang lebih datar dan konservatif di kisaran ~6.570. Mengingat pola data testing yang berfluktuasi relatif stabil, proyeksi DES tampak lebih realistis dan sesuai dengan kondisi data terkini.

KESIMPULAN

Metode Double Exponential Smoothing DES 1 (α=0,2, β=0,2) direkomendasikan untuk peramalan pengeluaran suntik 3 bulan karena menghasilkan akurasi tertinggi pada data testing dengan MAPE = 6,64%, yang termasuk dalam kategori peramalan baik (MAPE < 10%).

Perbandingan ARIMA dan DES

Patrecia Daniela Butarbutar

2026-05-07

PERSIAPAN DATA

AUTOREGRESSIVE INTEGRATED MOVING AVERAGE

DOUBLE EXPONENTIAL SMOOTHING

PERBANDINGAN ARIMA VS DES

KESIMPULAN