Pendugaan Model Data Harga Gula Dunia 2012-2022

Kelompok 6 P2 MPDW

Anggota:

Pendahuluan

Latar Belakang

Gula merupakan salah satu pemegang peranan penting dan merupakan komoditas strategis pada bidang pertanian (Sugiyanto 2017). Gula juga digunakan sebagai bahan baku pembuat alkohol, bahan pengawet makanan dan pencampur obat - obatan (Syafa’at 2005). Konsumsi gula diperkirakan cenderung bertambah tahun ke tahun seiring meningkatnya jumlah penduduk, pertumbuhan ekonomi, dan pertumbuhan industri makanan dan minuman ( Maria 2009).

Produksi gula mengalami fluktuasi karena permintaan yang cenderung meningkat, oleh karena itu untuk memenuhi kebutuhan permintaan gula pemerintah melakukan opsi untuk impor gula (Safrida 2004). Permintaan gula dipengaruhi oleh harga gula yang berkaitan dengan daya beli masyarakat (Kurniawati 2005).

Di pasar produk pertanian, harga gula menjadi salah satu indikator terbaik untuk mengukur dan mengetahui kondisi pasar dunia, regional, dan domestik. Harga gula menjadi sangat penting dalam membantu menguraikan strategi ekonomi serta membuat keputusan yang tepat dan tepat waktu mengenai produksi (untuk petani), pemasaran, serta konsumsi terutama untuk negara-negara dengan rasio produk pertanian terhadap total impor sangat tinggi (Apergis dan Rezitis 2011).

Harga gula dari tahun ke tahun mengalami kenaikkan dan penurunan yang cukup bervariatif sehingga memengaruhi daya jual beli masyarakat terhadap gula. Untuk itu sangat penting untuk mengetahui prediksi dari harga gula yang akan datang pada periode tertentu. Salah satu metode yang digunakan untuk memprediksi harga gula adalah Metode Double Exponential Smoothing (DES).

Double Exponential Smoothing (DES) merupakan suatu metode prediksi dengan memberi nilai pembobot pada beberapa periode atau pengamatan sebelumnya untuk memprediksi nilai pada periode yang akan datang (Hilmy et al. 2021).

Tujuan Penelitian

Tujuan dari penelitian ini adalah sebagai berikut:

  1. Membangun model ARIMA yang menjelaskan hubungan antara harga gula dunia saat ini dengan harga gula dunia pada periode-periode sebelumnya
  2. Melakukan peramalan harga gula dunia untuk periode Agustus 2022 - Juli 2023 berdasarkan model terbaik

Tinjauan Pustaka

Gula

Gula merupakan suatu karbohidrat sederhana karena dapat larut dalam air dan langsung diserap tubuh untuk diubah menjadi energi. Secara umum, gula dapat dibedakan menjadi dua, yaitu monosakarida dan disakarida (Darwin 2013). Gula merupakan komoditas utama perdagangan di Indonesia. Gula merupakan salah satu pemanis yang umum dikonsumsi masyarakat. Gula biasa digunakan sebagai pemanis di makanan maupun minuman, dalam bidang makanan, selain sebagai pemanis, gula juga digunakan sebagai stabilizer dan pengawet. Kebutuhan akan gula setiap tahunnya meningkat, selain karena pertamabahan jumlah penduduk juga karena semakin banyak industry yang memerlukan gula sebagai bahan baku sehingga harga gula pasir terus mengalami perubahan setiap tahunnya.

Analisis Deret Waktu

Analisis deret waktu adalah suatu metode kuantitatif untuk menentukan pola data masa lalu yang telah dikumpulkan secara teratur. Data dapat berupa data harian, mingguan, bulanan, kuartal, dan tahunan. Analisis deret waktu adalah salah satu prosedur statistika yang ditetapkan untuk meramalkan struktur probalistik keadaan yang akan terjadi pada masa akan datang dalam pengambilan keputusan (Aswi dan Sukarna 2006).

Analis deret waktu didasarkan dua domain yaitu domain waktu dan domain frekuensi. Model dalam domain waktu antara lain Autoregressive Integrated Moving Average (ARIMA). Analisis pada domain frekuensi yaitu analisis deret waktu yang dianggap sebagai akibat dari adanya komponen siklus pada frekuensi berbeda yang sulit diperoleh dalam domain waktu (Habinuddin et al. 2019).

Kestasioneran

Asumsi yang paling penting pada analisis deret waktu adalah stasioneritas data. Data yang stasioner adalah data yang nilai rata-rata dan variannya tidak mengalami perubahan secara sistematik sepanjang waktu atau yang memiliki rata- rata dan variannya konstan. Untuk melihat adanya stasioneritas tersebut, dapat dilihat dengan grafik atau plot time series menggunakan Scatter Plot. Grafik tersebut menjelaskan observasi dengan waktu. Jika terlihat memiliki rata-rata dan varians konstan, data tersebut dapat disimpulkan stasioner. Jika kondisi stasioner dalam rata-rata tidak terpenuhi, maka dilakukan proses pembedaan (differencing). Apabila kondisi stasioner dalam variansi tidak terpenuhi, maka dilakukan transformasi pangkat (Power Transformation) (Aswi dan Sukarna, 2006).

Autokorelasi dan Autokorelasi Parsial

Metode grafik dengan Scatter Plot memiliki kelemahan dalam objektivitas peneliti. Setiap peneliti dapat memiliki pandangan yang berbeda-beda sehingga dibutuhkan uji formal yang akan menguatkan keputusan secara ilmiah. Salah satu uji formal tersebut adalah korelogram. Korelogram merupakan teknik identifikasi stasioner data deret waktu melalui fungsi autokorelasi (ACF) dan fungsi autokorelasi parsial (PACF) (Habinuddin et al. 2019).

Model Deret Waktu ARIMA

Model deret waktu ARIMA (Autoregressive Integrated Moving Average) merupakan gabungan dari model AutoregressiveI (AR) ordo p, model Moving Average (MA) ordo q, dan diikuti proses differencing ordo d. Model ini dapat digunakan untuk data deret waktu yang bersifat univariat. Model ARIMA (p, q, d) dapat ditulis dalam persamaan (Wei 2006):

Data dan Analisis Eksplorasi

Metodologi

library(dplyr)
library(readxl)
library(forecast)
library(TTR)
library(graphics)
require(smooth)
require(Mcomp)
library(fpp2)
library(tseries)
library(readr)
library(TSA)
library(ggplot2)
library(MLmetrics)
library(cowplot)
library(gridExtra)
library(gtable)
library(grid)
library(MASS)

Data yang kami gunakan pada analisis kali ini adalah data harga gula dunia dalam US cents per bulan sejak Januari 2000 hingga Juni 2022. Data terdiri dari 271 amatan dengan 2 variabel yaitu variabel waktu dan harga gula. Data bersumber dari IMF, yang dapat diakses melalui https://www.imf.org/en/Research/commodity-prices

#Input Data
data<- read_xlsx("C:/Users/Alam/Downloads/bisa mpdw - Copy.xlsx")
head(data)
## # A tibble: 6 x 2
##   t                    Gula
##   <dttm>              <dbl>
## 1 2000-01-01 00:00:00  5.63
## 2 2000-02-01 00:00:00  5.35
## 3 2000-03-01 00:00:00  5.11
## 4 2000-04-01 00:00:00  6.01
## 5 2000-05-01 00:00:00  6.78
## 6 2000-06-01 00:00:00  8.33
str(data)
## tibble [271 x 2] (S3: tbl_df/tbl/data.frame)
##  $ t   : POSIXct[1:271], format: "2000-01-01" "2000-02-01" ...
##  $ Gula: num [1:271] 5.63 5.35 5.11 6.01 6.78 ...
dim(data)
## [1] 271   2
Yt<-data$Gula

Pembagian Data Training dan Data Testing

Data training adalah bagian dari data set yang digunakan untuk membuat prediksi atau dalam tahapan modelling. Data yang digunakan sebagai data training sebanyak 70% dari keseluruhan data set, yaitu sebanyak 191 amatan. Sedangkan data testing adalah bagian dari data set yang digunakan untuk melihat keakuratan dari data training. Data yang digunakan sebagai data testing sebanyak 45% dari keseluruhan data set, yaitu sebanyak 81 amatan.

test <- tail(data,81)
train <- head(data,190)

#data time series
train.ts<-ts(train$Gula)
test.ts<-ts(test$Gula, start = 191)
data.ts<-ts(data$Gula)

Eksplorasi Data

plot(x = data$t,
     y = Yt,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Harga Gula",
     main = "Time Series Plot of Harga Gula Dunia")

plot(x = train$t,
     y = train$Gula,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Harga Gula Dunia",
     main = "Time Series Plot of Harga Gula Dunia (Data Training)")

plot(x = test$t,
     y = test$Gula,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Harga Gula Dunia",
     main = "Time Series Plot of Harga Gula Dunia (Data Testing)")

Terlihat bahwa pola data harga gula membentuk pola tren menurun. Berdasarkan time series plot, harga gula dunia mengalami peningkatan yang signifikan pada tahun 2016 hingga 2017, dan tahun 2020 hingga sekarang, yakni pada masa pandemi.

Pemulusan

Berdasarkan hasil pengecekan data yang tidak stasioner dan memiliki tren negatif, digunakan Double Moving Average dan Double Exponential Smoothing sebagai metode pemulusan terhadap data. Pemulusan akan dilakukan terhadap data training, dan dilakukan perhitungan akurasi forecasting terhadap data testing.

Double Moving Average

Dilakukan proses pemulusan dengan metode double moving average (DMA) dengan optimisasi nilai m.

akurasi = data.frame()
for(n in 1:95){
  dma.step1<-SMA(train.ts, n=n)
  dma <- SMA(dma.step1, n=n)
  pemulusan.DMA <- c(dma, rep(NA,81))
  At <- 2*dma.step1 - dma
  Bt <- 2/(n-1)*(dma.step1 - dma)
  data.dma<- At+Bt
  f = c()
  for(i in 1:81){
   f[i] = At[length(At)] + Bt[length(Bt)]*(i)
  }
  data.ramal2<- c(NA, data.dma)
  ramalan2 = c(data.ramal2, f[-1])
  
  
  error.dma = test.ts-ramalan2[191:length(ramalan2)]
  SSE.dma = sum(error.dma^2)
  MSE.dma = SSE.dma/length(test.ts)
  RMSE.dma = sqrt(MSE.dma)
  MAPE.dma = mean(abs((error.dma/test.ts)*100))
  
  akurasi.dma <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
  row.names(akurasi.dma)<- c("SSE", "MSE", "MAPE")
  colnames(akurasi.dma) <- c(paste("Akurasi m =", n))
  akurasi = rbind(akurasi, c(n=n, akurasi=akurasi.dma[3]))
}
akurasi[order(akurasi[,2]),]
##    X1      Inf.
## 52 52  18.68041
## 51 51  19.21765
## 53 53  19.50629
## 54 54  20.36333
## 55 55  21.82757
## 50 50  23.63784
## 56 56  23.93290
## 4   4  26.86796
## 57 57  27.09642
## 49 49  29.12170
## 58 58  31.03847
## 48 48  35.24623
## 59 59  35.47231
## 5   5  37.22702
## 60 60  40.11844
## 47 47  40.89457
## 61 61  44.92639
## 46 46  45.98237
## 62 62  48.20328
## 63 63  50.03205
## 45 45  51.47619
## 64 64  51.53315
## 65 65  52.55266
## 66 66  53.17128
## 67 67  54.36658
## 68 68  55.96303
## 44 44  56.59479
## 69 69  58.75480
## 21 21  61.28032
## 70 70  61.49829
## 43 43  61.89787
## 25 25  62.36716
## 20 20  62.65885
## 22 22  62.72118
## 24 24  63.40591
## 23 23  63.64989
## 26 26  64.02930
## 71 71  65.12246
## 42 42  66.31795
## 27 27  66.97993
## 19 19  67.49214
## 72 72  68.07570
## 41 41  69.08096
## 28 28  69.90148
## 73 73  71.12152
## 40 40  71.41621
## 29 29  73.16353
## 6   6  73.24203
## 74 74  74.23434
## 18 18  74.59172
## 39 39  75.25395
## 30 30  76.31989
## 75 75  76.81165
## 38 38  77.23697
## 76 76  78.25992
## 37 37  78.75513
## 31 31  79.27086
## 77 77  79.34220
## 95 95  79.36948
## 94 94  79.88060
## 85 85  79.88493
## 92 92  80.02280
## 91 91  80.03161
## 84 84  80.07475
## 90 90  80.09897
## 78 78  80.10154
## 93 93  80.15563
## 83 83  80.23436
## 86 86  80.24110
## 79 79  80.32033
## 80 80  80.40117
## 36 36  80.46399
## 87 87  80.50170
## 89 89  80.50188
## 82 82  80.52112
## 81 81  80.59034
## 88 88  80.63749
## 32 32  81.13528
## 35 35  81.37653
## 34 34  82.19451
## 33 33  82.33926
## 17 17  83.45730
## 16 16  94.09805
## 7   7 101.41075
## 15 15 108.06681
## 14 14 120.03503
## 8   8 127.55118
## 13 13 128.12117
## 3   3 131.04358
## 12 12 136.99141
## 9   9 137.51868
## 10 10 141.63837
## 11 11 143.74276
## 2   2 479.21584
## 1   1       Inf

Didapatkan M = 43 sebagai nilai M optimum.

val = 52
dma.step1<-SMA(train.ts, n=val)
dma <- SMA(dma.step1, n=val)
pemulusan.DMA <- c(dma, rep(NA,81))
At <- 2*dma.step1 - dma
Bt <- 2/(val-1)*(dma.step1 - dma)
data.dma<- At+Bt
f = c()
for(i in 1:81){
  f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}
data.ramal2<- c(NA, data.dma)
ramalan2 = c(data.ramal2, f[-1])

error.dma = test.ts-ramalan2[191:length(ramalan2)]
  SSE.dma = sum(error.dma^2)
  MSE.dma = SSE.dma/length(test.ts)
  RMSE.dma = sqrt(MSE.dma)
  MAPE.dma = mean(abs((error.dma/test.ts)*100))
accu <- data.frame(SSE.dma, MSE.dma, RMSE.dma, MAPE.dma)
accu
##    SSE.dma  MSE.dma RMSE.dma MAPE.dma
## 1 1268.909 15.66555 3.957973 18.68041
ts.plot(data.ts,xlab="Time Period ", ylab="Harga Gula", main= "DMA M=52 Harga Gula")
lines(pemulusan.DMA,col="green",lwd=2)
lines(ramalan2, col= "red",lwd=2)

Double Exponential Smoothing

Dilakukan proses pemulusan dengan metode double exponential smoothing (DES) dengan optimisasi nilai alpha dan beta.

#plot DES
df.master <- data.frame()
for(i in 1:10){
  x = i/10
  for(j in 1:10){
    y = j/10
    optimumdes <- HoltWinters(train.ts, alpha =  y, beta = x,  gamma = FALSE)
    ramal_des <- forecast(optimumdes, h = 81)
    
    error.des1 <- (ramal_des$mean-test.ts)
    ssetest <- sum(error.des1^2)
    msetest <- ssetest/length(error.des1)
    rmsetest <- sqrt(msetest)
    mapetest <- sum(abs(error.des1/test.ts))/length(error.des1)*100
    ak <- data.frame(Alpha=y,Beta=x,SSE=ssetest,MSE=msetest,RMSE=rmsetest, MAPE = mapetest)
    df.master <- rbind(df.master,ak)
  }
}
df.master[order(df.master[,6]),]
##     Alpha Beta         SSE        MSE      RMSE      MAPE
## 10    1.0  0.1   1228.3520   15.16484  3.894206  16.86720
## 16    0.6  0.2    880.6980   10.87281  3.297395  17.02722
## 24    0.4  0.3    966.7278   11.93491  3.454694  18.06583
## 9     0.9  0.1   1895.8140   23.40511  4.837883  20.72893
## 42    0.2  0.5   1953.2913   24.11471  4.910673  22.68673
## 33    0.3  0.4   1411.3732   17.42436  4.174250  24.35295
## 15    0.5  0.2   2454.0640   30.29709  5.504279  25.89437
## 91    0.1  1.0   1825.0587   22.53159  4.746745  26.30653
## 8     0.8  0.1   3064.5945   37.83450  6.150975  29.06497
## 17    0.7  0.2   1765.7098   21.79889  4.668928  29.66708
## 7     0.7  0.1   4889.9452   60.36969  7.769794  40.19100
## 25    0.5  0.3   4308.1850   53.18747  7.292974  46.76203
## 18    0.8  0.2   4294.1980   53.01479  7.281126  46.87910
## 23    0.3  0.3   6070.0305   74.93865  8.656711  48.18877
## 52    0.2  0.6   4819.7053   59.50253  7.713789  49.05287
## 6     0.6  0.1   7538.3320   93.06583  9.647063  52.63900
## 14    0.4  0.2   7919.8190   97.77554  9.888152  55.33913
## 19    0.9  0.2   7789.3547   96.16487  9.806369  61.98918
## 5     0.5  0.1  11195.1524  138.21176 11.756350  66.25707
## 34    0.4  0.4   9990.8141  123.34338 11.106007  69.24324
## 43    0.3  0.5  11618.6611  143.44026 11.976655  74.22599
## 20    1.0  0.2  11579.7772  142.96021 11.956597  74.58938
## 26    0.6  0.3  12007.7457  148.24377 12.175540  75.49600
## 32    0.2  0.4  14984.1214  184.98915 13.601072  80.18097
## 4     0.4  0.1  16002.5491  197.56234 14.055687  80.84636
## 1     0.1  0.1  18431.9234  227.55461 15.084913  87.34807
## 62    0.2  0.7  18498.3575  228.37478 15.112074  92.03466
## 13    0.3  0.2  20373.9200  251.52988 15.859693  92.72027
## 81    0.1  0.9  19196.3715  236.99224 15.394552  93.20881
## 3     0.3  0.1  21492.3261  265.33736 16.289179  94.81128
## 2     0.2  0.1  23818.2006  294.05186 17.147940  99.95991
## 27    0.7  0.3  22713.2220  280.41015 16.745452 101.70797
## 11    0.1  0.2  25773.4582  318.19084 17.837905 102.48873
## 35    0.5  0.4  24978.4482  308.37590 17.560635 106.12205
## 53    0.3  0.6  29396.7473  362.92281 19.050533 114.39739
## 44    0.4  0.5  31142.5876  384.47639 19.608070 117.61519
## 22    0.2  0.3  35062.4962  432.87032 20.805536 123.59680
## 28    0.8  0.3  35178.9865  434.30848 20.840069 125.42367
## 72    0.2  0.8  37107.1050  458.11241 21.403561 127.55971
## 71    0.1  0.8  38836.1077  479.45812 21.896532 131.92168
## 61    0.1  0.7  40583.2504  501.02778 22.383650 134.45916
## 12    0.2  0.2  42289.7355  522.09550 22.849409 135.09985
## 36    0.6  0.4  45497.4075  561.69639 23.700135 141.43478
## 63    0.3  0.7  48148.6568  594.42786 24.380891 144.68725
## 29    0.9  0.3  47921.1365  591.61897 24.323219 145.89831
## 82    0.2  0.9  54242.2561  669.65748 25.877741 153.03911
## 51    0.1  0.6  53760.7303  663.71272 25.762623 154.77766
## 54    0.4  0.6  58560.1950  722.96537 26.888015 159.26800
## 45    0.5  0.5  59517.4787  734.78369 27.106894 160.83823
## 30    1.0  0.3  59277.8005  731.82470 27.052259 162.04613
## 92    0.2  1.0  65281.9283  805.94973 28.389254 167.53486
## 73    0.3  0.8  65711.8297  811.25716 28.482576 168.25289
## 37    0.7  0.4  69949.2502  863.57099 29.386578 174.74738
## 21    0.1  0.3  82610.6767 1019.88490 31.935637 188.05973
## 83    0.3  0.9  83509.0699 1030.97617 32.108818 189.22546
## 41    0.1  0.5  87699.7854 1082.71340 32.904611 197.24629
## 64    0.4  0.7  91464.6333 1129.19300 33.603467 198.21151
## 46    0.6  0.5  96546.7724 1191.93546 34.524418 204.39097
## 38    0.8  0.4  95839.9198 1183.20889 34.397804 204.55044
## 55    0.5  0.6 105025.0371 1296.60540 36.008407 212.67193
## 93    0.3  1.0 106058.3887 1309.36282 36.185119 212.86275
## 31    0.1  0.4 119023.2943 1469.42339 38.333059 228.34938
## 39    0.9  0.4 120029.9684 1481.85146 38.494824 229.23854
## 74    0.4  0.8 132591.7327 1636.93497 40.459053 238.42450
## 47    0.7  0.5 138254.4694 1706.84530 41.313984 245.06034
## 40    1.0  0.4 139081.8348 1717.05969 41.437419 247.03752
## 65    0.5  0.7 161690.1546 1996.17475 44.678571 263.93480
## 56    0.6  0.6 162505.3440 2006.23882 44.791057 265.14533
## 48    0.8  0.5 180000.0566 2222.22292 47.140459 280.03684
## 84    0.4  0.9 184403.1736 2276.58239 47.713545 281.35048
## 49    0.9  0.5 215933.7291 2665.84851 51.631856 307.05347
## 57    0.7  0.6 224319.5338 2769.37696 52.624870 312.04072
## 75    0.5  0.8 229260.1328 2830.37201 53.201241 314.40184
## 66    0.6  0.7 242114.3013 2989.06545 54.672346 323.65334
## 50    1.0  0.5 240000.8560 2962.97353 54.433202 323.97833
## 94    0.4  1.0 246159.7355 3039.00908 55.127208 325.20942
## 58    0.8  0.6 283080.2733 3494.81819 59.116987 350.93426
## 85    0.5  0.9 305745.7278 3774.63861 61.438088 363.15651
## 67    0.7  0.7 326125.7803 4026.24420 63.452693 376.10675
## 59    0.9  0.6 328731.4549 4058.41302 63.705675 378.51299
## 76    0.6  0.8 334077.5360 4124.41402 64.221601 380.15646
## 60    1.0  0.6 351960.5890 4345.19246 65.918074 391.91760
## 95    0.5  1.0 389049.6375 4803.08194 69.304271 409.69724
## 68    0.8  0.7 401453.0173 4956.21009 70.400356 417.71187
## 86    0.6  0.9 437637.7367 5402.93502 73.504660 435.06136
## 77    0.7  0.8 442118.5726 5458.25398 73.879997 437.81225
## 69    0.9  0.7 451904.0633 5579.06251 74.693122 443.51575
## 70    1.0  0.7 465088.4356 5741.83254 75.774881 450.20299
## 78    0.8  0.8 531305.7148 6559.32981 80.989690 480.36634
## 96    0.6  1.0 553098.9191 6828.38172 82.634023 489.06112
## 87    0.7  0.9 570875.3995 7047.84444 83.951441 497.40467
## 80    1.0  0.8 569179.4729 7026.90707 83.826649 497.80459
## 79    0.9  0.8 578281.1013 7139.27286 84.494218 501.48986
## 90    1.0  0.9 654015.7315 8074.26829 89.856932 533.45088
## 88    0.8  0.9 667907.8322 8245.77571 90.806254 538.44384
## 89    0.9  0.9 699812.2741 8639.65771 92.949759 551.50348
## 97    0.7  1.0 710398.3084 8770.34949 93.650144 554.78741
## 100   1.0  1.0 710068.4278 8766.27689 93.628398 555.74325
## 98    0.8  1.0 805211.7881 9940.88627 99.703993 591.08751
## 99    0.9  1.0 808006.5975 9975.39009 99.876875 592.48072

didapatkan alpha = 1 dan beta = 0.1 sebagai alpha dan beta optimum

optimumdes <- HoltWinters(train.ts, alpha = 1, beta = 0.1, gamma = FALSE)
ramal_des <- forecast(optimumdes, h = 81)

error.des1 <- (ramal_des$mean-test.ts)
sse.des <- sum(error.des1^2)
mse.des <- sse.des/length(test.ts)
rmse.des <- sqrt(mse.des)
mape.des <- sum(abs(error.des1/test.ts))/length(error.des1)*100
aku <- data.frame(SSE=sse.des,MSE=mse.des,RMSE=rmse.des, MAPE = mape.des)
aku
##        SSE      MSE     RMSE    MAPE
## 1 1228.352 15.16484 3.894206 16.8672
ts.plot(data.ts,xlab="Time Period ", ylab="Harga Gula", main= "DES Harga Gula", sub = "beta, alpha optimum")
lines(optimumdes$fitted[,1], col ="blue", lwd = 2)
lines(ramal_des$mean, col="red", lwd =2)

plot(ramal_des, xlab="Bulanan", ylab="US Cents", main="Plot Pemulusan Metode DES Ramalan 1")
lines((length(train.ts)+1):(length(train.ts)+81), test.ts,col="black")

Perbandingan Akurasi Smoothing

plot(1:length(train.ts), train.ts, type="b", col="black", xlim=c(1,length(train.ts)+81),
     ylab="Data Nilai Tukar",xlab="Bulanan", main="Perbandingan Metode Ramalan vs Aktual\nDilihat setelah data training")
lines((length(train.ts)+1):(length(train.ts)+81), test.ts,col="black")
lines((length(train.ts)+1):(length(train.ts)+81), ramalan2[191:271],col="blue", lwd = 2)
lines((length(train.ts)+1):(length(train.ts)+81), ramal_des$mean,col="red", lwd =2)
legend("topleft",c("data aktual", "Ramalan DMA", "Ramalan DES"), lty=1, 
       col=c("black", "blue", "red"))

metric<- matrix(c(SSE.dma,MSE.dma, RMSE.dma, MAPE.dma,
         sse.des,mse.des,rmse.des,mape.des),nrow=4,ncol=2)
colnames(metric) <-  c("DMA","DES")
rownames(metric) <- c("SSE", "MSE", "RMSE", "MAPE")
metric
##              DMA         DES
## SSE  1268.909313 1228.351980
## MSE    15.665547   15.164839
## RMSE    3.957973    3.894206
## MAPE   18.680412   16.867202

Berdasarkan perbandingan tersebut, didapatkan bahwa DES memiliki nilai SSE, MSE, RMSE, dan MAPE lebih kecil, sehingga digunakan metode DES sebagai metode pemulusan.

Pengecekan Kestasioneran Data

Plot ACF

acf(train.ts, lag.max = 50)

Berdasarkan plot ACF di atas dapat diketahui bahwa plot tersebut membentuk Tails off (slowlu)y yang artinya data tidak stasioner dan perlu proses pembedaan (differencing) terlebih dahulu sampai data menjadi stasioner.

Uji Formal ADF

Kestasioneran data diuji menggunakan ADF-Test dengan hipotesis sebagai berikut: H0 : Data tidak stasioner H1 : Data stasioner

adf.test(train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -2.0546, Lag order = 5, p-value = 0.5533
## alternative hypothesis: stationary

Terlihat p-value dari uji ADF tersebut lebih besar dari alpha = 0.05. Artinya, keputusan uji hipotesis ini adalah tidak tolak H0 atau dapat dikatakan bahwa data tidak stasioner dan perlu dilakukan penanganan ketidakstasioneran data atau differencing.

Penanganan Ketidakstasioneran Data

Differencing Pertama

Gula.dif1 <- diff(train.ts, differences = 1)

Plot Hasil Differencing Pertama

plot(Gula.dif1,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Time",
     ylab = "Data",
     main = "Gula")

Secara eksploratif dapat dilihat bahwa data setelah mengalami differencing dengan ordo 1 sudah lebih stasioner jika dibandingkan dengan data awal. Untuk mengetahuinya secara lebih pasti, perlu dilakukan uji formal menggunakan ADF Test.

Uji ADF Test Hasil Differencing Pertama

adf.test(Gula.dif1)
## Warning in adf.test(Gula.dif1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Gula.dif1
## Dickey-Fuller = -5.6449, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Terlihat p-value dari uji ADF lebih kecil dari alpha = 0.05. Artinya, keputusan uji hipotesis ini adalah tolak H0 atau dapat dikatakan bahwa data sudah stasioner dan tidak diperlukan differencing kedua.

Kandidat Model

Plot ACF

acf(Gula.dif1, lag.max = 20)

## Plot PACF

pacf(Gula.dif1, lag.max = 20)

## EACF

eacf(Gula.dif1)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o x o o o  x  o  o 
## 1 o o o o o o o o o o o  x  o  o 
## 2 o o o o o o o x o o o  x  o  o 
## 3 x x o o o o o o o o o  x  o  o 
## 4 x o x o o o o o o o o  o  o  o 
## 5 x x x x o o o o o o o  o  o  o 
## 6 o x x x o o o o o o o  o  x  o 
## 7 x x x x x x o o o o o  o  x  o

Berdasarkan plot ACF, PACF, dan EACF, didapati kandidat model dapat dibentuk yaitu : 1. ARIMA (0,1,1) 2. ARIMA (1,1,0) 3. ARIMA(3,1,2)
4. ARIMA(3,1,3) 5. ARIMA (4,1,3)

Pengecekan Auto Arima

auto.arima(train.ts)
## Series: train.ts 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.2785
## s.e.  0.0702
## 
## sigma^2 = 1.371:  log likelihood = -297.52
## AIC=599.04   AICc=599.1   BIC=605.52

Menggunakan fungsi auto.arima juga diperoleh kandidat model untuk data ini ialah ARIMA(1,1,0) yang sama seperti model 2.

Pendugaan Parameter dan Penentuan Model Terbaik

Model 1

model1 <- Arima(train.ts, order = c(0,1,1), method = "ML")
summary(model1)
## Series: train.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.2556
## s.e.  0.0647
## 
## sigma^2 = 1.38:  log likelihood = -298.14
## AIC=600.27   AICc=600.34   BIC=606.76
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.03740414 1.168457 0.8326766 0.1714894 6.127103 0.9708736
##                    ACF1
## Training set 0.02234317
lmtest::coeftest(model1)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 0.255599   0.064748  3.9476 7.893e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model pertama ARIMA (0,1,1) memiliki nilai AIC sebesar 600.27 dan parameter yang signifikan seluruhnya.

Model 2

model2 <- Arima(train.ts, order = c(1,1,0), method = "ML")
summary(model2)
## Series: train.ts 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.2785
## s.e.  0.0702
## 
## sigma^2 = 1.371:  log likelihood = -297.52
## AIC=599.04   AICc=599.1   BIC=605.52
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.03535191 1.164612 0.8300669 0.1842793 6.076027 0.9678308
##                     ACF1
## Training set 0.000872468
lmtest::coeftest(model2)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)    
## ar1 0.278457   0.070153  3.9693 7.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model kedua ARIMA (1,1,0) memiliki nilai AIC sebesar 599.1 dan parameter yang signifikan seluruhnya.

Model 3

model3 <- Arima(train.ts, order = c(3,1,2), method = "ML")
summary(model3)
## Series: train.ts 
## ARIMA(3,1,2) 
## 
## Coefficients:
##          ar1     ar2      ar3      ma1      ma2
##       0.3671  0.6672  -0.2724  -0.0915  -0.6975
## s.e.  0.3005  0.3284   0.0934   0.3021   0.2815
## 
## sigma^2 = 1.378:  log likelihood = -296.01
## AIC=604.03   AICc=604.49   BIC=623.48
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.05308206 1.155086 0.8280263 0.2157467 6.048277 0.9654514
##                      ACF1
## Training set -0.008051724
lmtest::coeftest(model3)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ar1  0.367081   0.300513  1.2215 0.221892   
## ar2  0.667226   0.328387  2.0318 0.042171 * 
## ar3 -0.272371   0.093432 -2.9152 0.003555 **
## ma1 -0.091529   0.302115 -0.3030 0.761919   
## ma2 -0.697550   0.281522 -2.4778 0.013220 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model pertama ARIMA (3,2,1) memiliki nilai AIC sebesar 604.03 dan parameter yang tidak signifikan di AR1, dan MA1

Model 4

model4 <- Arima(train.ts, order = c(3,1,3), method = "ML")
summary(model4)
## Series: train.ts 
## ARIMA(3,1,3) 
## 
## Coefficients:
##          ar1     ar2      ar3      ma1      ma2     ma3
##       0.5275  0.6738  -0.4795  -0.2752  -0.7070  0.2348
## s.e.  0.3811  0.3970   0.3216   0.4162   0.3412  0.3467
## 
## sigma^2 = 1.382:  log likelihood = -295.77
## AIC=605.54   AICc=606.16   BIC=628.23
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.05066084 1.153526 0.8252545 0.1920238 6.052164 0.9622196
##                     ACF1
## Training set 0.009429015
lmtest::coeftest(model4)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ar1  0.52746    0.38112  1.3840  0.16636  
## ar2  0.67383    0.39700  1.6973  0.08963 .
## ar3 -0.47952    0.32156 -1.4912  0.13590  
## ma1 -0.27518    0.41624 -0.6611  0.50854  
## ma2 -0.70705    0.34121 -2.0722  0.03825 *
## ma3  0.23484    0.34667  0.6774  0.49813  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model keempat ARIMA (3,1,3) memiliki nilai AIC sebesar 605.54 dan parameter hanya signifikan di MA2.

Model 5

model5 <- Arima(train.ts, order = c(4,1,3), method = "ML")
summary(model5)
## Series: train.ts 
## ARIMA(4,1,3) 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1     ma2      ma3
##       0.4324  -0.2876  0.8054  -0.3183  -0.1486  0.2128  -0.7333
## s.e.  0.2712   0.2091  0.1692   0.0705   0.2851  0.2051   0.1904
## 
## sigma^2 = 1.371:  log likelihood = -294.6
## AIC=605.19   AICc=605.99   BIC=631.13
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.05191168 1.146068 0.8238667 0.2079054 6.140046 0.9606014
##                     ACF1
## Training set -0.01559986
lmtest::coeftest(model5)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.432389   0.271171  1.5945 0.1108180    
## ar2 -0.287554   0.209065 -1.3754 0.1689979    
## ar3  0.805387   0.169170  4.7608 1.928e-06 ***
## ar4 -0.318273   0.070536 -4.5122 6.415e-06 ***
## ma1 -0.148583   0.285106 -0.5211 0.6022627    
## ma2  0.212805   0.205068  1.0377 0.2993963    
## ma3 -0.733346   0.190355 -3.8525 0.0001169 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model kelima ARIMA (4,1,3) memiliki nilai AIC sebesar 605.19 dan parameter tidak signifikan di AR1, AR2, MA1, dan MA2.

Pemilihan Model

Model <- c("ARIMA(0,1,1)", "ARIMA(1,1,0)", "ARIMA(3,1,2)", "ARIMA(3,1,3)", "ARIMA(4,1,3)")
AIC <- c(model1$aic, model2$aic, model3$aic, model4$aic, model5$aic)
Sig.Coef <- c("1 dari 1", "1 dari 1", "3 dari 5", "1 dari 5", "3 dari 7")
tabel <- data.frame(Model, AIC, Sig.Coef)
kableExtra::kable(tabel, align = c(rep('c', times = 4)), 
                  col.names = c("Model", "AIC", "Parameter Signifikan"), digits = 3)
Model AIC Parameter Signifikan
ARIMA(0,1,1) 600.272 1 dari 1
ARIMA(1,1,0) 599.039 1 dari 1
ARIMA(3,1,2) 604.025 3 dari 5
ARIMA(3,1,3) 605.537 1 dari 5
ARIMA(4,1,3) 605.194 3 dari 7

Berdasarkan hasil perbandingan di atas, maka dapat disimpulkan model terbaik dengan nilai AIC terkecil dan seluruh parameter signifikan adalah Model 2 ARIMA (1,1,0).

Diagnostik Model

Diagnostik model dilakukan untuk mengetahui apakah sisaan dari model yang dibangun telah memenuhi asumsi pemodelan. Asumsi sisaan yang harus terpenuhi pada model ARIMA(p,d,q) adalah sebagai berikut:

  1. Sisaan saling bebas (tidak ada autokorelasi) dan identik
  2. Sisaan mengikuti sebaran Normal (0, σ2)

Eksploratif

sisaan <- model2$residuals
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "red", lwd = 2)
plot(sisaan, type="o", 
     ylab = "Sisaan", xlab = "Order", main = "Sisaan M1 vs Order")
abline(h = 0, col='red')
acf(sisaan)
pacf(sisaan)

Normal Q-Q Plot Berdasarkan hasil eksplorasi di atas, terlihat bahwa masih banyak amatan yang cenderung menjauhi garis qq-plot distribusi normal. Sehingga secara eksploratif, dapat disimpulkan bahwa sisaan belum cukup menyebar normal.

Residual vs Order Berdasarkan hasil eksplorasi di atas, titik pada plot kebebasan sisaan mayoritas bergerak di sekitar titik nol. Namun, terdapat beberapa titik amatan yang terletak cukup jauh dari titik nol. Sehingga, belum dapat disimpulkan apakah terdapat autokorelasi atau tidak.

Plot ACF dan PACF Berdasarkan hasil eksplorasi di atas, baik dari plot ACF maupun plot PACF, pada keduanya terdapat garis vertikal di lag tertentu yang melebihi tinggi garis biru horizontal. Artinya, menurut kedua plot ini, terdapat autokorelasi pada model.

Uji Formal

Sisaan Menyebar Normal

Uji formal ini dilakukan dengan Jarque Bera test dan Shapiro-Wilk test.

Hipotesis yang diuji: H0: Sisaan Menyebar Normal H1: Sisaan Tidak Menyebar Normal

jarque.bera.test(sisaan)
## 
##  Jarque Bera Test
## 
## data:  sisaan
## X-squared = 39.636, df = 2, p-value = 2.473e-09
shapiro.test(sisaan)
## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan
## W = 0.95309, p-value = 6.42e-06

Berdasarkan Jarque-Bera test, diperoleh p-value (2.473e-09) < α (0.05). Selain itu, hasil Shapiro-Wilk test, diperoleh p-value (6.42e-06) < α (0.05) maka tolak H0. Artinya, tidak cukup bukti untuk menyatakan bahwa sisaan menyebar normal pada taraf nyata 5%.

Sisaan Saling Bebas atau Tidak ada Autokorelasi

Uji formal ini dilakukan dengan LJung-Box test.

Hipotesis yang diuji: H0: Sisaan antara lag saling bebas H1: Sisaan antara lag tidak saling bebas

Box.test(sisaan, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 0.00014692, df = 1, p-value = 0.9903

Berdasarkan LJung-Box test, diperoleh p-value (0.9903) > α (0.05), maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan antara lag saling bebas atau dapat dikatakan tidak ada autokorelasi antara sisaan lag pada taraf nyata 5%.

Nilai Tengah Sisaan Sama dengan Nol

Uji formal ini dilakukan dengan t-test.

Hipotesis yang diuji: H0: Nilai tengah sisaan sama dengan nol H1: Nilai tengah sisaan tidak sama dengan nol

t.test(sisaan, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.41751, df = 189, p-value = 0.6768
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1316756  0.2023794
## sample estimates:
##  mean of x 
## 0.03535191

Berdasarkan t-test, diperoleh p-value (0.6768) > α (0.05), maka tak tolak H0. Artinya, cukup bukti untuk menyatakan bahwa nilai tengah sisaan sama dengan nol pada taraf nyata 5%.

Berdasarkan eksplorasi dan uji formal yang telah dilakukan, hasil dari analisis model ARIMA (1,1,0) yaitu sisaan tidak menyebar normal, tidak ada autokorelasi, dan nilai tengah sisaan bernilai nol.

Overfitting

Overfitting model dilakukan dengan menggunakan:

  1. ARIMA (1,1,1)

  2. ARIMA (2,1,0)

model1a <- Arima(train.ts, order=c(1,1,1), method="ML")
summary(model1a)
## Series: train.ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1     ma1
##       0.2766  0.0020
## s.e.  0.2269  0.2332
## 
## sigma^2 = 1.378:  log likelihood = -297.52
## AIC=601.04   AICc=601.17   BIC=610.76
## 
## Training set error measures:
##                      ME     RMSE       MAE      MPE     MAPE     MASE
## Training set 0.03536171 1.164611 0.8300886 0.184171 6.076408 0.967856
##                      ACF1
## Training set 0.0007048348
lmtest::coeftest(model1a)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 0.2765734  0.2268816  1.2190   0.2228
## ma1 0.0020274  0.2332476  0.0087   0.9931
model1b <- Arima(train.ts, order=c(2,1,0), method="ML")
summary(model1b)
## Series: train.ts 
## ARIMA(2,1,0) 
## 
## Coefficients:
##          ar1      ar2
##       0.2786  -0.0007
## s.e.  0.0729   0.0730
## 
## sigma^2 = 1.378:  log likelihood = -297.52
## AIC=601.04   AICc=601.17   BIC=610.76
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.03536295 1.164611 0.8300921 0.1841541 6.076469 0.9678601
##                      ACF1
## Training set 0.0006668884
lmtest::coeftest(model1b)
## 
## z test of coefficients:
## 
##        Estimate  Std. Error z value  Pr(>|z|)    
## ar1  0.27863496  0.07288263  3.8231 0.0001318 ***
## ar2 -0.00065545  0.07301091 -0.0090 0.9928371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Karena koefisien penduga parameter pada model ARIMA (1,1,1) dan ARIMA (2,1,0) tidak nyata/tidak signifikan dan nilai AICnya lebih besar, maka model yang dipilih tetap model ARIMA (1,1,0).

Forecasting

ramalan <- forecast(model2, h = 81)
data.ramalan <- ramalan$mean

plot(ramalan)

Perbandingan data aktual dan hasil forecast

accu_mod<-accuracy(data.ramalan, test.ts)
accu_mod
##                ME    RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set 0.314603 3.16632 2.661461 -2.085646 17.53248 0.9361954  2.835419

Perbandingan Hasil Smoothing dan Modeling

perbandingan <- matrix(c(RMSE.dma, MAPE.dma, rmse.des,mape.des, accu_mod[,2], accu_mod[,5]),nrow=2,ncol=3)
colnames(perbandingan) <-  c("DMA","DES", "Arima (1,1,0)")
rownames(perbandingan) <- c("RMSE", "MAPE")
perbandingan
##            DMA       DES Arima (1,1,0)
## RMSE  3.957973  3.894206       3.16632
## MAPE 18.680412 16.867202      17.53248
plot(1:length(train.ts), train.ts, type="b", col="black", xlim=c(1,length(train.ts)+81),
     ylab="Data Nilai Tukar",xlab="Bulanan", main="Perbandingan Metode Ramalan vs Aktual\nDilihat setelah data training")
lines((length(train.ts)+1):(length(train.ts)+81), test.ts,col="black")
lines((length(train.ts)+1):(length(train.ts)+81), ramalan2[191:271],col="blue", lwd = 2)
lines((length(train.ts)+1):(length(train.ts)+81), ramal_des$mean,col="red", lwd =2)
lines((length(train.ts)+1):(length(train.ts)+81), data.ramalan,col="green", lwd =2)
legend("topleft",c("data aktual", "Ramalan DMA", "Ramalan DES", "ramalan ARIMA (1,1,0)"), lty=1, 
       col=c("black", "blue", "red", "green"))