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)
# Import data
data <- read_excel("C:/Users/LENOVO/Downloads/alokon.xlsx")
data
## # A tibble: 50 × 2
##    Bulan               `Suntik 3 Bulan`
##    <dttm>                         <dbl>
##  1 2022-02-01 00:00:00             2233
##  2 2022-03-01 00:00:00             2597
##  3 2022-04-01 00:00:00             2859
##  4 2022-05-01 00:00:00             2797
##  5 2022-06-01 00:00:00             2545
##  6 2022-07-01 00:00:00             2681
##  7 2022-08-01 00:00:00             2724
##  8 2022-09-01 00:00:00             2773
##  9 2022-10-01 00:00:00             2845
## 10 2022-11-01 00:00:00             2533
## # ℹ 40 more rows
# 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)
# Plot data
plot(ts_data, main="Data Pengeluaran Suntik 3 Bulan", ylab="Jumlah", xlab="Waktu")

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")

plot(testing, main="Data Testing")

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

# Cek Stasioneritas
# Plot ACF
acf(training, main="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 
# Differencing Ordo 1
diff1 <- diff(training, differences=1)
plot(diff1, main="Differencing Orde 1")

acf(diff1) # Cek Kestasioneran

pacf(diff1)

adf.test(diff1) # Augmented Dickey-Fuller
## 
##  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")

acf(diff2) # Cek Kestasioneran

pacf(diff2)

adf.test(diff2) # Augmented Dickey-Fuller
## 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.

# Identifikasi Model Tentatif
# Plot ACF
acf(diff2, main="ACF diff2")

# Plot PACF
pacf(diff2, main="PACF diff2")

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

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

# Bandingkan AIC
data.frame(
  Model=c("120","121","021","220"),
  AIC=c(arima120$aic, arima121$aic, arima021$aic, arima220$aic)
)
##   Model      AIC
## 1   120 517.9282
## 2   121 502.5525
## 3   021 510.3651
## 4   220 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)
# Eksplorasi
par(mfrow=c(2,2))
qqnorm(res); qqline(res)
plot(res, main="Residual")
acf(res)
pacf(res)

# Uji Non-Autokorelasi dengan Ljung-Box test
Box.test(res, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 0.22795, df = 1, p-value = 0.633

Uji Ljung-Box: p-value = 0,633 > 0,05 -> gagal menolak H0, artinya residual tidak terdapat autokorelasi (asumsi terpenuhi).

# Uji Normalitas dengan Shapiro-Wilk test
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.94487, p-value = 0.07885

Uji Shapiro-Wilk: p-value = 0,07885 > 0,05 -> gagal menolak H0, artinya residual berdistribusi normal (asumsi terpenuhi).

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

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

plot(forecast_arima)

accuracy(forecast_arima, testing)
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set  -27.60887  404.8015  284.0442  -1.007558  6.288565 0.1549210
## Test set     -973.02696 1153.6350 1024.4452 -15.667452 16.440659 0.5587441
##                     ACF1 Theil's U
## Training set -0.07736216        NA
## Test set      0.52396743   2.96728
# Forecast ARIMA April - Desember 2026
model_final <- forecast::Arima(ts_data, order=c(1,2,1))
forecast_final <- forecast(model_final, h=9)
forecast_final
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Apr 2026       6535.254 5956.417 7114.091 5650.000 7420.509
## May 2026       6623.772 5922.098 7325.447 5550.653 7696.891
## Jun 2026       6684.914 5834.463 7535.365 5384.262 7985.567
## Jul 2026       6755.632 5783.341 7727.924 5268.641 8242.624
## Aug 2026       6823.001 5732.176 7913.826 5154.727 8491.274
## Sep 2026       6891.541 5688.936 8094.146 5052.315 8730.767
## Oct 2026       6959.671 5648.637 8270.706 4954.617 8964.726
## Nov 2026       7027.945 5611.432 8444.457 4861.576 9194.314
## Dec 2026       7096.168 5576.217 8616.119 4771.604 9420.733
plot(forecast_final)

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 
forecast_des1 <- forecast(des1, h=length(testing))

plot(forecast_des1)

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)

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.

# Lamda 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)

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
accuracy(forecast_des1, testing)
##                     ME     RMSE      MAE       MPE      MAPE      MASE
## Training set -258.9249 632.8842 540.5035 -8.640747 14.333485 0.2947968
## Test set     -383.5258 506.1346 403.8250 -6.335381  6.637448 0.2202508
##                   ACF1 Theil's U
## Training set 0.5125035        NA
## Test set     0.2332825  1.313956
accuracy(forecast_des2, testing)
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  -50.47299 473.8278 362.1408 -2.077451 8.168828 0.1975157
## Test set     -529.03593 671.1571 560.7556 -8.631129 9.108116 0.3058425
##                    ACF1 Theil's U
## Training set -0.2110817        NA
## Test set      0.4139070  1.734288
accuracy(forecast_des_opt, testing)
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  -93.66219 470.1652 375.5136 -3.363125 8.941217 0.2048094
## Test set     -530.30215 663.9013 557.0721 -8.652185 9.054741 0.3038335
##                    ACF1 Theil's U
## Training set -0.1685509        NA
## Test set      0.3903089  1.717643

DES 1 menghasilkan MAPE testing terkecil (6,64%), meskipun MAPE training-nya paling besar. Ini menunjukkan DES 1 lebih baik dalam generalisasi ke data baru, kemungkinan karena parameter yang lebih konservatif menghasilkan peramalan yang lebih stabil.

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

forecast_des_final <- forecast(des_final, h=9)

forecast_des_final
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Apr 2026       6435.116 5858.962 7011.270 5553.965 7316.267
## May 2026       6451.954 5760.952 7142.956 5395.158 7508.750
## Jun 2026       6468.792 5656.630 7280.955 5226.698 7710.887
## Jul 2026       6485.631 5546.256 7425.005 5048.981 7922.280
## Aug 2026       6502.469 5430.072 7574.865 4862.380 8142.558
## Sep 2026       6519.307 5308.300 7730.313 4667.232 8371.382
## Oct 2026       6536.145 5181.142 7891.148 4463.847 8608.443
## Nov 2026       6552.983 5048.783 8057.184 4252.507 8853.459
## Dec 2026       6569.821 4911.389 8228.254 4033.468 9106.175
plot(forecast_des_final)

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

data.frame(
  Model=c("ARIMA","DES 1","DES 2","DES Optimum"),
  MAPE=c(
    16.44,
    accuracy(forecast_des1, testing)[2,"MAPE"],
    accuracy(forecast_des2, testing)[2,"MAPE"],
    accuracy(forecast_des_opt, testing)[2,"MAPE"]
  )
)
##         Model      MAPE
## 1       ARIMA 16.440000
## 2       DES 1  6.637448
## 3       DES 2  9.108116
## 4 DES Optimum  9.054741
autoplot(ts_data) +
  autolayer(forecast_final$mean, series="ARIMA") +
  autolayer(forecast_des_final$mean, series="DES") +
  ggtitle("Perbandingan Forecast ARIMA vs DES")

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%).