Kestasioneran Data Deret Waktu pada Data Minyak Mentah WTI

Package

Terdapat 6 packages yang digunakan, yaitu:

readxl, TTR, TSA, tseries, graphics, dan forecast

Data

Harga Minyak Mentah WTI (West Texas Index Crude Oil Price)

AlMadi dan Zhang (2011), Pasar minyak mentah dunia memiliki 4 tolak ukur (benchmark) internasional utama yang digunakan untuk penetapan harga minyak mentah lainnya, yaitu 2 di West - West Texas Intermediate (WTI) di AS dan Brent di Eropa, dan 2 di Asia - Dubai dan Oman. Benchmark digunakan untuk menentukan harga ekspor minyak mentah setelah melakukan penyesuaian terkait perbedaan kualitas dan lokasi.

Deskripsi Data dan Time Series Plot

Data yang digunakan dalam pemodelan deret waktu ini adalah sebanyak 128 amatan yaitu Harga Minyak WTI untuk January 2012 hingga Agustus 2022 dengan periode pengamatan bulanan. Dengan Time Series Plot sebagai berikut:

data <- read_excel("D:/MPDW/WTI Crude Oil (2).xlsx", sheet = 2)
head(data)

## # A tibble: 6 x 1
##   Terakhir
##      <dbl>
## 1     98.5
## 2    107. 
## 3    103. 
## 4    105. 
## 5     86.5
## 6     85.0

ts.dt <- ts(data)
ts.plot(ts.dt, xlab = "Time Period", ylab = "Crude Oil Prices",
        main = "Time Series Plot of WTI Crude Oil Prices")

Pergerakan harga minyak mentah cenderung nonstasioner serta dipengaruhi oleh banyak faktor. Berbagai metode peramalan dapat memberikan hasil yang baik jika deretan harga yang terbentuk merepresentasikan fungsi yang cenderung linier. Namun, fluktuasi harga riil minyak mentah cenderung memiliki karakteristik nonlinier dan nonstasioner. Hal ini disebabkan oleh kompleksitas intrinsik dari pasar minyak (Herawati dan Djunaidy 2014).

Eksplorasi Data

par(mfrow = c(2,2))
ts.plot(ts.dt[1:30,], xlab = "Time Period", ylab = "Crude Oil Prices",
        main = "Time Series Plot of WTI Crude Oil Prices ", sub = "(Januari 2012 - June 2014)")

ts.plot(ts.dt[31:96,], xlab = "Time Period", ylab = "Crude Oil Prices",
        main = "Time Series Plot of WTI Crude Oil Prices", sub = "(Juli 2014 - Desember 2019)")
ts.plot(ts.dt[97:128,], xlab = "Time Period", ylab = "Crude Oil Prices",
        main = "Time Series Plot of WTI Crude Oil Prices", sub = "(Januari 2020 - Agustus 2022)")

Berdasarkan plot time series di atas, terlihat bahwa data merupakan data gabungan stasioner dan trend naik di akhir periode.

Terlihat pula dari ketiga plot diatas bahwa harga minyak WTI berfluktuatif. Seperti yang ditunjukkan pada plot untuk periode waktu Juli 2014-Desember 2019, dimana pada awal periode terjadi penurunan yang begitu signifikan. Lonjakan pasokan minyak dan minimnya permintaan merupakan faktor utama di belakang kejatuhan harga minyak semenjak pertengahan 2014. Ditambah lagi bangkitnya minyak shale dari Amerika Serikat dan juga penolakan OPEC untuk menurunkan produksinya. Selain itu pada awal tahun 2016 juga terjadi penurunan, dimana pada Februari 2016 harga minyak WTI hanya sebesar 33.75 USD. Yang mana jika dibandingkan sejak Juli 2014, di mana harga minyak berada di puncaknya, nilai ini sudah jatuh 70 persen (Bareksa.com)

Dan pada plot untuk periode waktu Januari 2020-Agustus 2022, dimana masa pandemic yang ikut memengaruhi harga minyak WTI. Menurut Widyastuti N L dan Nugroho Hanan (Juni 2020), yaitu Permintaan minyak mentah yang sudah mulai terjadi pada kuartal pertama tahun 2020, akan terus menurun hingga mencapai titik terbawahnya pada akhir kuartal ke-2, kemudian diperkirakan naik lagi dan mencapai titik keseimbangannya dengan pasokan (supply) minyak bumi pada akhir kuartal ke-3 tahun 2020. Dalam masa penurunan permintaan tersebut, minyak yang dikonsumsi akan berkurang hingga lebih 10 juta barel per hari, atau bahkan 15-20 juta barel per hari. Atau sekitar 20%-30% menurut beberapa analisis yang lain. Produksi minyak mentah tidak mungkin bisa diturunkan seketika, dalam masa penurunan permintaan yang tajam tersebut akan terbentuk stock minyak yang cukup besar. Hal ini membuat harga minyak turun.

Namun turunnya harga minyak bukan saja karena dampak Covid-19, tapi juga didorong oleh “pertikaian” (conflict) dalam industri minyak sendiri. Pertikaian, khususnya dalam kelompok produsen, diawali oleh Saudi Arabia dan Rusia, menyangkut kesepakatan mengenai berapa jumlah minyak mentah yang seharusnya diproduksi.

Smoothing (Pemulusan)

Splitting Data

#SPLIT data 80:20
train <- data[1:102,]
test <- data[103:128,]

train.ts <- ts(train)
test.ts <- ts(test)

par(mfrow = c(1,2))
ts.plot(train.ts, xlab = "Time Period", ylab = "Crude Oil Prices",
        main = "Time Series Plot of Data Training")
ts.plot(test.ts, xlab = "Time Period", ylab = "Crude Oil Prices",
        main = "Time Series Plot of Data Testing")

Double Moving Average

Pola data yang cenderung menunjukkan pola trend, sehingga untuk pemulusan digunakan Double Moving Average(DMA).

Penentuan n optimum

Dalam pemulusan dengan menggunakan Double Moving Average perlu ditentukan nilai n optimum, agar model ini menjadi model yang baik, ditunjukkan dengan ukuran keakuratan seperti SSE, MSE, RMSE, MAD, atau MAPE yang minimum

iterasi_dma <- function(x,min,max,metode=c("SSE","MSE","RMSE","MAD","MAPE")){
  df.master <- data.frame()
  z <- max-min+1
  temp <- sqrt(z)
  a <- round(temp) 
  b=a
  if(a*b<z){
    b=b+1
  }
  par(mfrow=c(a,b))
  for(i in min:max) {
    df_ts_sma <- TTR::SMA(x,  n=i)
    df_ts_dma <- TTR::SMA(df_ts_sma, n=i)
    At <- 2*df_ts_sma-df_ts_dma
    Bt <- df_ts_sma-df_ts_dma
    pemulusan_dma <- At+Bt
    ramal_dma <- c(NA, pemulusan_dma)
    df_dma <- cbind(df_aktual=c(x,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
    error.dma <- df_dma[, 1] - df_dma[, 3]
    SSE.dma <- sum(error.dma^2, na.rm = T)
    MSE.dma <- mean(error.dma^2, na.rm = T)
    RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))
    MAD.dma <- mean(abs(error.dma), na.rm = T)
    r.error.dma <- (error.dma/df_dma[, 1])*100 
    MAPE.dma <- mean(abs(r.error.dma), na.rm = T)
    ak <- data.frame(n=i,SSE=SSE.dma,MSE=MSE.dma,RMSE=RMSE.dma,
                     MAD=MAD.dma,MAPE=MAPE.dma)
    df.master <- rbind(df.master,ak)
    ts.plot(x,main=paste("Double Moving Average n =",i), lwd = 1)
    lines(pemulusan_dma,col="green",lwd=0.5)
    lines(ramal_dma,col="red",lwd= 0.5)
  }
  opt <- df.master$n[which.min(df.master[,metode])]
  return(list(Akurasi=df.master,n.optimum=paste("n optimum adalah",opt)))
}

dma_iter<-iterasi_dma(train.ts,2,10, metode= "MAPE")

tibble::as_tibble(dma_iter[[1]])

## # A tibble: 9 x 6
##       n    SSE   MSE  RMSE   MAD  MAPE
##   <int>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     2  4855.  49.0  7.00  5.23  9.48
## 2     3  8020.  82.7  9.09  6.60 12.4 
## 3     4 10169. 107.  10.3   7.72 14.9 
## 4     5 11073. 119.  10.9   8.13 15.9 
## 5     6 12288. 135.  11.6   8.68 16.8 
## 6     7 13447. 151.  12.3   8.99 17.7 
## 7     8 14729. 169.  13.0   9.57 19.0 
## 8     9 16004. 188.  13.7  10.3  20.6 
## 9    10 16890. 203.  14.3  11.0  22.3

dma_iter[[2]]

## [1] "n optimum adalah 2"

Dari hasil perbandingan model DMA dengan n dari 2 hingga 10, nilai MAPE minimum didapatkan ketika n bernilai 2

DMA dengan n optimum

Melakukan pemulusan data menggunakan DMA dengan n sebesar 2 sebagai berikut:

# Fungsi Double Moving Average
dma_function <- function(x,n,metode=c("SSE","MSE","RMSE","MAD","MAPE")){
  df.master <- data.frame()
  df_ts_sma <- TTR::SMA(x,  n=n)
  df_ts_dma <- TTR::SMA(df_ts_sma, n=n)
  At <- 2*df_ts_sma-df_ts_dma
  Bt <- df_ts_sma-df_ts_dma
  pemulusan_dma <- At+Bt
  ramal_dma <- c(NA, pemulusan_dma)
  df_dma <- cbind(df_aktual=c(x,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
  error.dma <- df_dma[, 1] - df_dma[, 3]
  SSE.dma <- sum(error.dma^2, na.rm = T)
  MSE.dma <- mean(error.dma^2, na.rm = T)
  RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))
  MAD.dma <- mean(abs(error.dma), na.rm = T)
  r.error.dma <- (error.dma/df_dma[, 1])*100 
  MAPE.dma <- mean(abs(r.error.dma), na.rm = T)
  ak <- data.frame(n=n,SSE=SSE.dma,MSE=MSE.dma,RMSE=RMSE.dma,
                     MAD=MAD.dma,MAPE=MAPE.dma)
  df.master <- rbind(df.master,ak)
  ts.plot(x,main=paste("Double Moving Average n =",n), lwd = 1)
  lines(pemulusan_dma,col="green",lwd=0.5)
  lines(ramal_dma,col="red",lwd= 0.5)
  return(list(df_dma, df.master))
}

smoothing_dma<-dma_function(train.ts, n=2, metode = .)

tibble::as_tibble(tail(smoothing_dma[[1]]))

## # A tibble: 6 x 3
##   df_aktual pemulusan_dma ramal_dma
##       <dbl>         <dbl>     <dbl>
## 1      44.8         40.0      54.5 
## 2      20.5         17.1      40.0 
## 3      18.8          6.70     17.1 
## 4      35.5         34.7       6.70
## 5      39.3         47.6      34.7 
## 6      NA           NA        47.6

tibble::as_tibble(smoothing_dma[[2]])

## # A tibble: 1 x 6
##       n   SSE   MSE  RMSE   MAD  MAPE
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     2 4855.  49.0  7.00  5.23  9.48

Duble Exponential Smoothing

Penentuan parameter optimum

iterasi_des <- function(x,alfa=NULL,beta=NULL,metrik=c("SSE","MSE","RMSE", "MAPE")){
  df.master <- data.frame()
  par(mfrow=c(3,3))
    for (i in alfa) {
      for (j in beta) {
        df_des <- HoltWinters(x, alpha = i, beta=j, gamma=F)
        datades <- data.frame(x,c(NA, NA, df_des$fitted[,1]))
        error<-datades[,1]-datades[,2]
        MAPE<-mean(abs((error[3:length(error)]/datades[3:length(error),1])*100))
        sse <- df_des$SSE 
        mse <- sse/length(x)
        rmse <-sqrt(mse)
        ak <- data.frame(Alpha=i,Beta=j,SSE=sse,MSE=mse,RMSE=rmse, MAPE)
        df.master <- rbind(df.master,ak)
        datades <- data.frame(x, c(NA, NA, df_des$fitted[,1]))
        colnames(datades) = c("y","yhat")
        ts.plot(x,xlab="periode  waktu", ylab="Yt",  col="blue", lty=3,main=paste("DES Alpha=",i,"Beta=",j))
        points(x)
        lines (datades[,2], col="red",lwd=2)
      }
    }
  opt <- df.master$Alpha[which.min(df.master[,metrik])]
    opt2 <- df.master$Beta[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt,
                                                  ",Beta optimum adalah",opt2)))
    
}

iterasi_des(train.ts, alfa = seq(0.1,0.9, 0.1), beta = seq(0.1,0.9, 0.1), metrik = "MAPE")

## $Akurasi
##    Alpha Beta       SSE       MSE      RMSE      MAPE
## 1    0.1  0.1 60420.812 592.36090 24.338465 33.434716
## 2    0.1  0.2 33078.812 324.30208 18.008389 22.082936
## 3    0.1  0.3 29967.096 293.79506 17.140451 22.082751
## 4    0.1  0.4 30409.124 298.12867 17.266403 23.556232
## 5    0.1  0.5 29909.648 293.23184 17.124014 23.747225
## 6    0.1  0.6 28048.092 274.98130 16.582560 23.098588
## 7    0.1  0.7 25525.562 250.25060 15.819311 21.312703
## 8    0.1  0.8 23415.429 229.56303 15.151338 20.726247
## 9    0.1  0.9 22419.044 219.79454 14.825469 20.649836
## 10   0.2  0.1 20926.938 205.16606 14.323619 17.917175
## 11   0.2  0.2 16177.007 158.59811 12.593574 16.284840
## 12   0.2  0.3 15127.378 148.30762 12.178162 15.875399
## 13   0.2  0.4 14582.928 142.96988 11.957001 15.982064
## 14   0.2  0.5 14573.741 142.87981 11.953234 16.240125
## 15   0.2  0.6 14997.938 147.03861 12.125948 16.565347
## 16   0.2  0.7 15552.777 152.47821 12.348207 17.428767
## 17   0.2  0.8 15951.655 156.38878 12.505550 18.029325
## 18   0.2  0.9 16083.510 157.68147 12.557128 18.149239
## 19   0.3  0.1 12950.227 126.96301 11.267786 14.497017
## 20   0.3  0.2 11041.455 108.24956 10.404305 13.406524
## 21   0.3  0.3 10752.123 105.41297 10.267082 13.487621
## 22   0.3  0.4 10855.285 106.42436 10.316219 13.905927
## 23   0.3  0.5 11070.139 108.53078 10.417811 14.563476
## 24   0.3  0.6 11238.682 110.18316 10.496817 15.071510
## 25   0.3  0.7 11345.880 111.23412 10.546759 15.327087
## 26   0.3  0.8 11443.540 112.19157 10.592052 15.406869
## 27   0.3  0.9 11567.802 113.40983 10.649405 15.578168
## 28   0.4  0.1  9635.440  94.46509  9.719316 12.765276
## 29   0.4  0.2  8704.898  85.34213  9.238081 12.173503
## 30   0.4  0.3  8677.171  85.07030  9.223356 12.555884
## 31   0.4  0.4  8814.234  86.41406  9.295916 13.034166
## 32   0.4  0.5  8964.480  87.88706  9.374810 13.296079
## 33   0.4  0.6  9103.907  89.25399  9.447433 13.372129
## 34   0.4  0.7  9242.447  90.61223  9.519046 13.396941
## 35   0.4  0.8  9381.186  91.97241  9.590225 13.569095
## 36   0.4  0.9  9513.008  93.26479  9.657370 13.638916
## 37   0.5  0.1  7828.788  76.75282  8.760869 11.737221
## 38   0.5  0.2  7330.037  71.86311  8.477211 11.619116
## 39   0.5  0.3  7384.401  72.39608  8.508589 11.899783
## 40   0.5  0.4  7521.560  73.74078  8.587245 12.004506
## 41   0.5  0.5  7663.024  75.12769  8.667623 12.106865
## 42   0.5  0.6  7796.641  76.43766  8.742863 12.134801
## 43   0.5  0.7  7918.880  77.63608  8.811134 12.075750
## 44   0.5  0.8  8025.032  78.67678  8.869993 11.949409
## 45   0.5  0.9  8111.501  79.52452  8.917652 11.974163
## 46   0.6  0.1  6684.713  65.53640  8.095456 11.060043
## 47   0.6  0.2  6402.985  62.77437  7.923028 10.972442
## 48   0.6  0.3  6485.970  63.58795  7.974205 11.034780
## 49   0.6  0.4  6613.416  64.83741  8.052168 11.096977
## 50   0.6  0.5  6734.786  66.02732  8.125719 11.045529
## 51   0.6  0.6  6839.505  67.05397  8.188649 11.028261
## 52   0.6  0.7  6923.790  67.88029  8.238950 10.972964
## 53   0.6  0.8  6986.404  68.49416  8.276120 11.080729
## 54   0.6  0.9  7028.616  68.90800  8.301084 11.292594
## 55   0.7  0.1  5891.656  57.76133  7.600087 10.448477
## 56   0.7  0.2  5728.019  56.15705  7.493801 10.318806
## 57   0.7  0.3  5816.930  57.02873  7.551737 10.297504
## 58   0.7  0.4  5927.507  58.11282  7.623176 10.374900
## 59   0.7  0.5  6024.605  59.06475  7.685359 10.464480
## 60   0.7  0.6  6101.760  59.82118  7.734415 10.491950
## 61   0.7  0.7  6159.109  60.38342  7.770677 10.631558
## 62   0.7  0.8  6199.578  60.78017  7.796164 10.830501
## 63   0.7  0.9  6227.496  61.05388  7.813698 10.984395
## 64   0.8  0.1  5315.039  52.10823  7.218603  9.869692
## 65   0.8  0.2  5221.450  51.19069  7.154767  9.753680
## 66   0.8  0.3  5310.310  52.06186  7.215390  9.946628
## 67   0.8  0.4  5409.216  53.03153  7.282275 10.073540
## 68   0.8  0.5  5494.695  53.86956  7.339588 10.099574
## 69   0.8  0.6  5565.530  54.56402  7.386746 10.229840
## 70   0.8  0.7  5626.041  55.15727  7.426794 10.369036
## 71   0.8  0.8  5682.511  55.71090  7.463973 10.509952
## 72   0.8  0.9  5741.932  56.29345  7.502896 10.677383
## 73   0.9  0.1  4892.924  47.96984  6.926026  9.331836
## 74   0.9  0.2  4849.030  47.53951  6.894890  9.489583
## 75   0.9  0.3  4944.648  48.47694  6.962538  9.676914
## 76   0.9  0.4  5048.804  49.49808  7.035487  9.746427
## 77   0.9  0.5  5147.791  50.46854  7.104121  9.914294
## 78   0.9  0.6  5245.522  51.42669  7.171240 10.047833
## 79   0.9  0.7  5350.214  52.45308  7.242450 10.249832
## 80   0.9  0.8  5471.457  53.64173  7.324052 10.522653
## 81   0.9  0.9  5619.447  55.09262  7.422440 10.805432
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.9 ,Beta optimum adalah 0.1"

DES dengan parameter optimum

des_function<- function(x,alfa=NULL,beta=NULL,metrik=c("SSE","MSE","RMSE")){
  df.master <- data.frame()
        df_des <- HoltWinters(x, alpha = alfa, beta=beta, gamma=F)
        datades <- data.frame(x,c(NA, NA, df_des$fitted[,1]))
        error<-datades[,1]-datades[,2]
        mape<-mean(abs((error[3:length(error)]/datades[3:length(error),1])*100))
        sse <- df_des$SSE 
        mse <- sse/length(x)
        rmse <-sqrt(mse)
        ak <- data.frame(SSE=sse,MSE=mse,RMSE=rmse, MAPE=mape)
        df.master <- rbind(df.master,ak)
        datades <- data.frame(x, c(NA, NA, df_des$fitted[,1]))
        colnames(datades) = c("y","yhat")
        ts.plot(x,xlab="periode  waktu", ylab="Yt",  col="blue", lty=3,main=paste("DES Alpha=",alfa,"Beta=",beta))
        points(x)
        lines (datades[,2], col="red",lwd=2)
         return(df.master)
      }

smoothing_des<-des_function(train.ts, alfa = 0.9, beta = 0.1, metrik = .)

smoothing_des

##        SSE      MSE     RMSE     MAPE
## 1 4892.924 47.96984 6.926026 9.331836

Forecasting

Double Moving Average

df_ts_sma <- TTR::SMA(train.ts,  n=2)
df_ts_dma <- TTR::SMA(df_ts_sma, n=2)
At <- 2*df_ts_sma-df_ts_dma
Bt <- df_ts_sma-df_ts_dma
pemulusan_dma <- At+Bt
ramal_dma <- c(NA, pemulusan_dma)
df_dma <- cbind(df_aktual=c(train.ts,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
head (df_dma)

##      df_aktual pemulusan_dma ramal_dma
## [1,]     98.48            NA        NA
## [2,]    107.07            NA        NA
## [3,]    103.02       107.315        NA
## [4,]    104.87       102.845   107.315
## [5,]     86.53        87.455   102.845
## [6,]     84.96        75.790    87.455

Eksploratif: DMA Plot

ts.plot(data, xlab="periode waktu", ylab="Yt", col="blue", lty=1,
        xlim=c(0,150),ylim=c(0,160),gpars=list(col.main="black",col.axis="white",col.sub="black"))
points(data,col="black")
par(col = "black")
lines(pemulusan_dma,col="red",lwd=2,lty=1)
par(col = "black")
lines(ramal_dma,col="green",lwd= 2,lty=1)
title("Rataan Bergerak  Berganda N=2", cex.main=1, font.main=4 ,col.main="black")
legend("topright", c("Data Aktual","Pemulusan","Ramalan"),lty=1,col=c ("blue","red","green"))
box(col="black",lwd=2)

Berdasarkan times series plot metode DMA, n yang optimum adalah 2, karena data training DMA dengan n=2 lebih mendekati sebaran pola dari data aktual dibandingkan hasil peramalannya

Double Exponential Smoothing

df_des <- HoltWinters(train.ts, alpha = 0.9, beta=0.1, gamma=F)
datades <- data.frame(train.ts,c(NA, NA, df_des$fitted[,1]))
colnames(datades)=c("y","yhat")
head(datades)

##        y      yhat
## 1  98.48        NA
## 2 107.07        NA
## 3 103.02 115.66000
## 4 104.87 111.73640
## 5  86.53 112.39106
## 6  84.96  93.62303

ramal_des <- forecast::forecast(df_des,h=3)
(df_ramal_des <- data.frame(ramal_des))

##     Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## 103       37.60995 28.70389 46.51602 23.98930 51.23061
## 104       36.68618 24.15391 49.21844 17.51973 55.85263
## 105       35.76240 19.96448 51.56032 11.60157 59.92324

ts.plot(data, xlab="periode  waktu", ylab="Yt", col="blue", lty=1,
        xlim=c(0,150),ylim=c(0,150))
points(data,col="black")
par(col="black")
lines(datades[,2], col="red",lwd=2) #nilai dugaan
par(col="black")
lines(ts(df_ramal_des$Point.Forecast,start = 108), col="green",lwd=2) #nilai dugaan
title("Double Exponential Smoothing Alpha=0.9 Beta = 0.1",cex.main=1, font.main=4 ,col.main="black")
legend("bottomright", c("Data aktual", "Fitted DES","Ramalan"), lty=1,col=c ("blue","red","green"))
box(col="black",lwd=2)

Berdasarkan time series plot DES, plot data training hasil fitted menggunakan metode pemulusan DES dengan α = 0.9 dan β = 0.1 menampilkan plot yang berpola tidak terlalu jauh dari plot data aktual dibandingkan dengan data peramalannya.