Aktifkan Library

library(readxl)    # Baca file Excel
library(zoo)       # Labeling data time series
library(forecast)  # Model & peramalan
library(ggplot2)   # Visualisasi data

Impor Data

Data yang digunakan dalam penelitian ini adalah data harga komoditas beras di pasar tradisional Provinsi Jawa Tengah periode Januari 2021 hingga Februari 2026. Berdasarkan data tersebut, akan dilakukan peramalan harga hingga Desember 2026.

Data diperoleh dari website Pusat Informasi Harga Pangan Strategis (PIHPS) Nasional yang dikelola oleh Bank Indonesia dan dapat diakses melalui: https://www.bi.go.id/hargapangan/TabelHarga/PasarTradisionalKomoditas

Tabel_Harga_Berdasarkan_Komoditas <- read_excel("D:/KP/ANALISIS DATA/Tabel Harga Berdasarkan Komoditas.xlsx", sheet = "Beras")
Harga_Beras <- ts(
  Tabel_Harga_Berdasarkan_Komoditas[["Harga Beras Provinsi Jawa Tengah"]][1:62],
  start = c(2021, 1),
  frequency = 12)
Harga_Beras
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2021 10950 10950 10900 10800 10850 10850 10800 10700 10750 10750 10750 10850
## 2022 11000 11000 10950 10950 11000 10950 10950 11100 11300 11350 11400 11900
## 2023 12250 12800 12650 12550 12650 12600 12400 12800 13500 13750 13750 13950
## 2024 14050 15400 15400 15000 14500 14450 14550 14550 14650 14650 14500 14550
## 2025 14550 14500 14750 14750 14750 14900 14950 15000 15100 15100 15050 15050
## 2026 15000 15200

Plot Time Series

# Plot utama 
plot.ts(Harga_Beras,
        col = "blue",
        ylab = "Harga (Rp/kg)",
        xlab = "Waktu",
        xaxt = "n")

title(main = "Data Harga Beras Provinsi Jawa Tengah",
      cex.sub = 0.8)

# Titik data
points(Harga_Beras,
       pch = 20,
       col = "red")

# Ambil waktu dalam format yearmon
waktu_full <- as.yearmon(time(Harga_Beras))

idx <- seq(1, length(waktu_full), by = 12)

axis(1,
     at = waktu_full[idx],
     labels = format(waktu_full[idx], "%b %Y"),
     cex.axis = 0.8)

Keterangan: Data harga secara umum memiliki kecenderungan meningkat dari waktu ke waktu. Pola tersebut mengindikasikan adanya komponen trend dalam data. Karena data mengandung unsur tren tanpa pola musiman yang jelas, metode yang sesuai untuk digunakan adalah Double Exponential Smoothing Holt, yang dirancang untuk menangkap pola tren pada data deret waktu.

Pembagian Data Training dan Testing

Data harga beras yang terdiri dari 62 periode observasi dibagi menjadi data training dan data testing dengan proporsi 80:20. Sebanyak 50 periode pertama digunakan sebagai data training untuk membangun model peramalan, sedangkan 12 periode terakhir digunakan sebagai data testing untuk menguji performa model terhadap data yang tidak terlibat dalam proses estimasi parameter. Pendekatan ini digunakan untuk mengidentifikasi potensi overfitting atau underfitting pada model yang dibangun.

# Training 50 data pertama
train_beras <- window(Harga_Beras, end = time(Harga_Beras)[50])

# Testing 12 data terakhir
test_beras <- window(Harga_Beras, start = time(Harga_Beras)[51])

length(train_beras) 
## [1] 50
length(test_beras) 
## [1] 12

Plot Pembagian Data

plot(train_beras,
     type = "o",
     pch = 16,
     cex = 0.6,
     lwd = 1.5,
     xlim = range(time(Harga_Beras)),
     ylim = range(Harga_Beras),
     main = "Pembagian Data Training (80%) dan Testing (20%)",
     ylab = "Harga (Rp/kg)",
     xlab = "Waktu",
     xaxt = "n")

lines(test_beras,
      type = "o",
      col = "red",
      pch = 16,
      cex = 0.6,
      lwd = 1.5)

legend("topleft",
       legend = c("Training", "Testing"),
       col = c("black", "red"),
       lty = 1,
       pch = 16,
       pt.cex = 0.6)

# Ambil waktu full 
waktu_full <- as.yearmon(time(Harga_Beras))

idx <- seq(1, length(waktu_full), by = 12)

axis(1,
     at = waktu_full[idx],
     labels = format(waktu_full[idx], "%b %Y"),
     cex.axis = 0.8)

# Legend
legend("topleft",
       legend = c("Training", "Testing"),
       col = c("black", "red"),
       lty = 1,
       pch = 16,
       pt.cex = 0.6)

Analisis dengan Metode Double Exponential Smoothing Holt

Mencari parameter terbaik menggunakan data training

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = train_beras, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0.04011698
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a 14500.00000
## b    61.00513
model_DES$fitted
##              xhat level      trend
## Mar 2021 10950.00 10950   0.000000
## Apr 2021 10897.99 10900  -2.005849
## May 2021 10794.06 10800  -5.937078
## Jun 2021 10846.31 10850  -3.693051
## Jul 2021 10846.46 10850  -3.544897
## Aug 2021 10794.59 10800  -5.408536
## Sep 2021 10690.80 10700  -9.203259
## Oct 2021 10743.17 10750  -6.828203
## Nov 2021 10743.45 10750  -6.554277
## Dec 2021 10743.71 10750  -6.291339
## Jan 2022 10847.97 10850  -2.027252
## Feb 2022 11004.07 11000   4.071622
## Mar 2022 11003.91 11000   3.908281
## Apr 2022 10951.75 10950   1.745644
## May 2022 10951.68 10950   1.675614
## Jun 2022 11003.61 11000   3.614242
## Jul 2022 10951.46 10950   1.463401
## Aug 2022 10951.40 10950   1.404693
## Sep 2022 11107.37 11100   7.365888
## Oct 2022 11315.09 11300  15.093786
## Nov 2022 11366.49 11350  16.494118
## Dec 2022 11417.84 11400  17.838272
## Jan 2023 11937.18 11900  37.181143
## Feb 2023 12299.73 12250  49.730490
## Mar 2023 12869.80 12800  69.799790
## Apr 2023 12710.98 12650  60.982087
## May 2023 12604.52 12550  54.523972
## Jun 2023 12706.35 12650  56.348333
## Jul 2023 12652.08 12600  52.081959
## Aug 2023 12441.97 12400  41.969193
## Sep 2023 12856.33 12800  56.332307
## Oct 2023 13582.15 13500  82.154308
## Nov 2023 13838.89 13750  88.887770
## Dec 2023 13835.32 13750  85.321861
## Jan 2024 14039.92 13950  89.922401
## Feb 2024 14140.33 14050  90.326684
## Mar 2024 15540.86 15400 140.860969
## Apr 2024 15535.21 15400 135.210053
## May 2024 15113.74 15000 113.739044
## Jun 2024 14589.12 14500  89.117689
## Jul 2024 14533.54 14450  83.536708
## Aug 2024 14634.20 14550  84.197165
## Sep 2024 14630.82 14550  80.819430
## Oct 2024 14731.59 14650  81.588896
## Nov 2024 14728.32 14650  78.315796
## Dec 2024 14569.16 14500  69.156457
## Jan 2025 14618.39 14550  68.387958
## Feb 2025 14615.64 14550  65.644440

Plot fitted vs aktual (training DES Holt)

plot(train_beras,
     type = "o",
     pch = 16,
     cex = 0.6,
     col = "black",
     main = "Plot Fitted vs Aktual (Training DES Holt)",
     ylab = "Harga (Rp/kg)",
     xlab = "Waktu",
     xaxt = "n")  

lines(model_DES$fitted[, "xhat"],
      type = "o",
      pch = 16,
      cex = 0.6,
      col = "blue")

# Ambil waktu training
waktu_train <- as.yearmon(time(train_beras))

# Tampilkan label waktu tiap 8 bulan
idx <- seq(1, length(waktu_train), by = 8)

axis(1,
     at = waktu_train[idx],
     labels = format(waktu_train[idx], "%b %Y"),
     cex.axis = 0.8)

legend("topleft",
       legend = c("Aktual Training", "Fitted DES Holt"),
       col = c("black", "blue"),
       lty = 1,
       pch = 16,
       pt.cex = 0.6)

Hitung akurasi training

# Ambil fitted
fitted_train_DES <- model_DES$fitted[, "xhat"]

# Tidak mengikutkan 2 data awal yang NA
actual_train_valid <- train_beras[-c(1,2)]

error_train_DES <- actual_train_valid - fitted_train_DES

MSE_train_DES  <- mean(error_train_DES^2)
MAD_train_DES  <- mean(abs(error_train_DES))
MAPE_train_DES <- mean(abs(error_train_DES / actual_train_valid)) * 100

round(c(MSE=MSE_train_DES, MAD=MAD_train_DES, MAPE=MAPE_train_DES),2)
##      MSE      MAD     MAPE 
## 79825.56   167.44     1.26

Keterangan: Akurasi training akan dibandingkan dengan testing

Forecast ke data testing

forecast_test_DES <- forecast(model_DES, h = length(test_beras))
forecast_test_DES
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Mar 2025       14561.01 14197.40 14924.61 14004.92 15117.09
## Apr 2025       14622.01 14097.38 15146.64 13819.65 15424.37
## May 2025       14683.02 14027.64 15338.39 13680.71 15685.32
## Jun 2025       14744.02 13972.36 15515.68 13563.87 15924.18
## Jul 2025       14805.03 13925.53 15684.52 13459.96 16150.09
## Aug 2025       14866.03 13884.16 15847.90 13364.39 16367.68
## Sep 2025       14927.04 13846.49 16007.59 13274.48 16579.60
## Oct 2025       14988.04 13811.40 16164.68 13188.53 16787.55
## Nov 2025       15049.05 13778.15 16319.94 13105.38 16992.71
## Dec 2025       15110.05 13746.20 16473.90 13024.23 17195.88
## Jan 2026       15171.06 13715.16 16626.96 12944.45 17397.66
## Feb 2026       15232.06 13684.71 16779.41 12865.60 17598.53

Plot forecast vs aktual testing

plot(forecast_test_DES,
     include = 0,
     xaxt = "n",
     main = "Plot Forecast DES Holt vs Aktual Testing",
     xlim = range(time(test_beras)),
     ylab = "Harga (Rp/kg)", xlab = "Waktu")

lines(test_beras,
      type = "o",
      col = "black",
      pch = 16,
      cex = 0.6)

# Ambil waktu testing
waktu_test <- as.yearmon(time(test_beras))

# Tampilkan tiap 3 bulan
idx <- seq(1, length(waktu_test), by = 3)

axis(1,
     at = waktu_test[idx],
     labels = format(waktu_test[idx], "%b %Y"),
     cex.axis = 0.9)

legend("topleft", 
       legend = c("Forecast DES Holt", "Aktual Testing"), 
       col = c("blue", "black"), 
       lty = 1, 
       pch = 16)

Hitung akurasi testing

actual_test <- test_beras
pred_test_DES   <- forecast_test_DES$mean

error_test_DES <- actual_test - pred_test_DES

MSE_test_DES <- mean(error_test_DES^2)
MAD_test_DES <- mean(abs(error_test_DES))
MAPE_test_DES <- mean(abs(error_test_DES / actual_test)) * 100

round(c(MSE=MSE_test_DES, MAD=MAD_test_DES, MAPE=MAPE_test_DES),2)
##      MSE      MAD     MAPE 
## 16352.41   113.99     0.76

Keterangan: Berdasarkan hasil evaluasi model DES Holt, diperoleh nilai MAPE pada data training sebesar 1,26% dan pada data testing sebesar 0,76%. Perbedaan nilai yang relatif kecil serta nilai MAPE testing yang lebih rendah menunjukkan bahwa model tidak mengalami overfitting dan memiliki kemampuan generalisasi yang baik terhadap data yang belum diamati.

Aplikasi model DES Holt pada seluruh data

Peramalan hingga Desember 2026 dilakukan menggunakan parameter terbaik yang diperoleh dari model training sebelumnya.
Semua nilai smoothing dan trend yang optimal dari tahap training diterapkan pada seluruh dataset untuk menghasilkan prediksi masa depan.

alpha_opt <- model_DES$alpha
beta_opt  <- model_DES$beta

model_DES_full <- HoltWinters(Harga_Beras,
                       alpha = alpha_opt,
                       beta  = beta_opt,
                       gamma = FALSE)
model_DES_full
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = Harga_Beras, alpha = alpha_opt, beta = beta_opt,     gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0.04011698
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a 15200.00000
## b    58.72108

Perhitungan manual model Double Exponential Smoothing Holt

Diketahui parameter hasil estimasi: \[ \alpha = 1 \]

\[ \beta = 0{,}04011698 \]

Inisialisasi

Nilai awal level ditentukan sebagai: \[ S_1 = X_1 = 10950 \]

Nilai awal trend dihitung sebagai: \[ T_1 = X_2 - X_1 = 10950 - 10950 = 0 \]

Periode ke-2 (Februari 2021)

Nilai pemulusan eksponensial tunggal pada periode ke-2: \[ S_2 = \alpha X_2 + (1 - \alpha)(S_1 + T_1) \]

\[ S_2 = 1(10950) + (1 - 1)(10950 + 0) \]

\[ S_2 = 10950 \]

Nilai trend periode ke-2: \[ T_2 = \beta (S_2 - S_1) + (1 - \beta)T_1 \]

\[ T_2 = 0{,}04011698(10950 - 10950) + (1 - 0{,}04011698)(0) \]

\[ T_2 = 0 \]

Hasil peramalan untuk periode ke-3: \[ F_3 = S_2 + T_2 = 10950 + 0 = 10950 \]

Periode ke-3 (Maret 2021)

Nilai pemulusan eksponensial tunggal pada periode ke-3: \[ S_3 = 1(10900) + (1 - 1)(10950 + 0) \]

\[ S_3 = 10900 \]

Nilai trend periode ke-3: \[ T_3 = 0{,}04011698(10900 - 10950) + (1 - 0{,}04011698)(0) \]

\[ T_3 = 0{,}04011698(-50) \]

\[ T_3 = -2{,}005849 \]

Hasil peramalan untuk periode ke-4: \[ F_4 = S_3 + T_3 \]

\[ F_4 = 10900 - 2{,}005849 = 10897{,}994151 \]

Untuk perhitungan periode berikutnya adalah sebagai berikut:

model_DES_full$fitted
##              xhat level      trend
## Mar 2021 10950.00 10950   0.000000
## Apr 2021 10897.99 10900  -2.005849
## May 2021 10794.06 10800  -5.937078
## Jun 2021 10846.31 10850  -3.693051
## Jul 2021 10846.46 10850  -3.544897
## Aug 2021 10794.59 10800  -5.408536
## Sep 2021 10690.80 10700  -9.203259
## Oct 2021 10743.17 10750  -6.828203
## Nov 2021 10743.45 10750  -6.554277
## Dec 2021 10743.71 10750  -6.291339
## Jan 2022 10847.97 10850  -2.027252
## Feb 2022 11004.07 11000   4.071622
## Mar 2022 11003.91 11000   3.908281
## Apr 2022 10951.75 10950   1.745644
## May 2022 10951.68 10950   1.675614
## Jun 2022 11003.61 11000   3.614242
## Jul 2022 10951.46 10950   1.463401
## Aug 2022 10951.40 10950   1.404693
## Sep 2022 11107.37 11100   7.365888
## Oct 2022 11315.09 11300  15.093786
## Nov 2022 11366.49 11350  16.494118
## Dec 2022 11417.84 11400  17.838272
## Jan 2023 11937.18 11900  37.181143
## Feb 2023 12299.73 12250  49.730490
## Mar 2023 12869.80 12800  69.799790
## Apr 2023 12710.98 12650  60.982087
## May 2023 12604.52 12550  54.523972
## Jun 2023 12706.35 12650  56.348333
## Jul 2023 12652.08 12600  52.081959
## Aug 2023 12441.97 12400  41.969193
## Sep 2023 12856.33 12800  56.332307
## Oct 2023 13582.15 13500  82.154308
## Nov 2023 13838.89 13750  88.887770
## Dec 2023 13835.32 13750  85.321861
## Jan 2024 14039.92 13950  89.922401
## Feb 2024 14140.33 14050  90.326684
## Mar 2024 15540.86 15400 140.860969
## Apr 2024 15535.21 15400 135.210053
## May 2024 15113.74 15000 113.739044
## Jun 2024 14589.12 14500  89.117689
## Jul 2024 14533.54 14450  83.536708
## Aug 2024 14634.20 14550  84.197165
## Sep 2024 14630.82 14550  80.819430
## Oct 2024 14731.59 14650  81.588896
## Nov 2024 14728.32 14650  78.315796
## Dec 2024 14569.16 14500  69.156457
## Jan 2025 14618.39 14550  68.387958
## Feb 2025 14615.64 14550  65.644440
## Mar 2025 14561.01 14500  61.005135
## Apr 2025 14818.59 14750  68.587037
## May 2025 14815.84 14750  65.835533
## Jun 2025 14813.19 14750  63.194410
## Jul 2025 14966.68 14900  66.676788
## Aug 2025 15016.01 14950  66.007766
## Sep 2025 15065.37 15000  65.365582
## Oct 2025 15166.76 15100  66.755011
## Nov 2025 15164.08 15100  64.077001
## Dec 2025 15109.50 15050  59.500577
## Jan 2026 15107.11 15050  57.113594
## Feb 2026 15052.82 15000  52.816520

Keterangan: xhat = hasil peramalan, level = S, trend = T

Perhitungan manual peramalan harga beras

Berdasarkan output model diperoleh koefisien:

  • \(a = 15200\)
  • \(b = 58{,}72108\)

dengan keterangan:

  • \(a\) merupakan estimasi level terakhir atau \(S_{62}\)
  • \(b\) merupakan estimasi trend terakhir atau \(T_{62}\)

Sehingga dapat dituliskan:

\[ S_{62} = 15200 \]

\[ T_{62} = 58{,}72108 \]

Peramalan Maret 2026 (\(m = 1\))

\[ F_{63} = 15200 + (58{,}72108 \times 1) \]

\[ F_{63} = 15258{,}72 \]

Peramalan April 2026 (\(m = 2\))

\[ F_{64} = 15200 + (58{,}72108 \times 2) \]

\[ F_{64} = 15317{,}44 \]

Untuk perhitungan periode berikutnya adalah sebagai berikut:

Proyeksi 10 bulan kedepan

forecast_DES <- forecast(model_DES_full, h = 10)
forecast_DES
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Mar 2026       15258.72 14929.06 15588.38 14754.55 15762.90
## Apr 2026       15317.44 14841.78 15793.10 14589.99 16044.90
## May 2026       15376.16 14781.97 15970.36 14467.42 16284.90
## Jun 2026       15434.88 14735.26 16134.51 14364.90 16504.87
## Jul 2026       15493.61 14696.21 16291.00 14274.10 16713.11
## Aug 2026       15552.33 14662.11 16442.54 14190.86 16913.79
## Sep 2026       15611.05 14631.37 16590.73 14112.76 17109.34
## Oct 2026       15669.77 14602.97 16736.57 14038.24 17301.29
## Nov 2026       15728.49 14576.24 16880.74 13966.27 17490.71
## Dec 2026       15787.21 14550.68 17023.74 13896.10 17678.32

Plot peramalan metode Double Exponential Smoothing Holt

# Data Aktual
df_actual <- data.frame(
  Waktu = as.yearmon(time(Harga_Beras)),
  Harga = as.numeric(Harga_Beras)
)

# Fitted training
fitted_train <- model_DES$fitted[, "xhat"]

# Samakan panjang waktu dengan fitted
time_fitted_train <- tail(time(train_beras), length(fitted_train))

df_fitted_train <- data.frame(
  Waktu = as.yearmon(time_fitted_train),
  Fitted = as.numeric(fitted_train)
)

# Forecast testing
n_test <- length(test_beras)

df_forecast_test <- data.frame(
  Waktu = as.yearmon(time(test_beras)),
  Fitted = as.numeric(forecast_test_DES$mean[1:n_test])
)

# Forecast
df_forecast_full <- data.frame(
  Waktu = as.yearmon(time(forecast_DES$mean)),
  Forecast = as.numeric(forecast_DES$mean),
  Lo80 = forecast_DES$lower[,"80%"],
  Hi80 = forecast_DES$upper[,"80%"],
  Lo95 = forecast_DES$lower[,"95%"],
  Hi95 = forecast_DES$upper[,"95%"]
)

# Gabungkan plot
ggplot() +
  
  # Confidence Interval
  geom_ribbon(data = df_forecast_full,
              aes(x = Waktu, ymin = Lo95, ymax = Hi95),
              fill = "blue", alpha = 0.2) +
  
  geom_ribbon(data = df_forecast_full,
              aes(x = Waktu, ymin = Lo80, ymax = Hi80),
              fill = "lightblue", alpha = 0.4) +
  
  # Aktual Full
  geom_line(data = df_actual,
            aes(x = Waktu, y = Harga, color = "Aktual"),
            linewidth = 1) +
  
  # Fitted Training
  geom_line(data = df_fitted_train,
            aes(x = Waktu, y = Fitted, color = "Fitted Training"),
            linetype = "dashed",
            linewidth = 1) +
  
  # Fitted Testing
  geom_line(data = df_forecast_test,
            aes(x = Waktu, y = Fitted, color = "Forecast Testing"),
            linetype = "dotted",
            linewidth = 1) +
  
  # Forecast Full
  geom_line(data = df_forecast_full,
            aes(x = Waktu, y = Forecast, color = "Forecast"),
            linewidth = 1) +
  
  labs(
    title = "Peramalan Harga Beras Provinsi Jawa Tengah dengan Metode DES Holt",
    x = "Waktu",
    y = "Harga (Rp/kg)",
    color = "Keterangan"
  ) +
  
  scale_color_manual(values = c(
    "Aktual" = "black",
    "Fitted Training" = "blue",
    "Forecast Testing" = "red",
    "Forecast" = "purple"
  )) +
  
  theme_minimal() +
  theme(legend.position = "right")

Kesimpulan

Berdasarkan evaluasi model Double Exponential Smoothing (DES) Holt, diperoleh nilai MAPE pada data training sebesar 1,26% dan pada data testing sebesar 0,76%. Nilai kesalahan yang rendah serta perbedaan yang kecil antara keduanya menunjukkan bahwa model memiliki kemampuan prediksi yang baik dan tidak mengalami overfitting, sehingga layak digunakan untuk peramalan.

Hasil peramalan menunjukkan bahwa harga beras pada Desember 2026 diperkirakan mencapai Rp15.787,21 per kilogram, yang mengindikasikan adanya kecenderungan peningkatan harga. Jika kenaikan ini terjadi, maka dapat meningkatkan beban pengeluaran rumah tangga dan berpotensi memengaruhi stabilitas inflasi, sehingga perlu diantisipasi melalui kebijakan pengendalian harga dan pasokan.