Pemodelan Data Saham PTBA Mingguan Periode Tahun 2020 - 2022

Kelompok 4

P1-MPDW 2022

Anggota:

Duwi Sulistianingsih (G1401201033)
Indri Ramdani (G1401201036)
Muhammad Fadhlan (G1401201070)
Naura Tirza Ardhani (G1401201073)
Faisal Arkan (G1401201077)

Pendahuluan

A. Latar Belakang

Perang yang terjadi antara Ukraina dan Rusia membawa dampak pada perekonomian Indonesia terutama pada sektor perdagangan. Kondisi pasar bisa merespon negatif jika harga berbagai komoditas mengalami peningkatan yang signifikan, terutama pada sektor pertambangan. Hal ini terjadi karena negara Rusia merupakan salah satu negara yang memproduksi berbagai bahan baku dan energi seperti bahan baku pertanian, logam, dan lain-lain.

Kenaikan harga komoditas yang terjadi akibat adanya konflik Rusia dan Ukraina ini akan berdampak pada perekonomian dunia, dengan melambatnya pertumbuhan ekonomi global akibat kenaikan harga komoditas, khususnya pada industri perminyakan dan pertambangan. Dampak negatif dirasakan pada sektor minyak. Harga minyak yang meningkat memberikan dampak negatif karena beban subsidi negara Indonesia akan meningkat akibat tekanan fiskal. Hal ini berbanding terbalik ada sektor pertambangan memberi dampak positif karena negara Indonesia merupakan negara eksportir batu bara. Contohnya yang terjadi pada perusahaan batu bara PT Bukit Asam Tbk yang merasakan keuntungan dengan meroketnya harga batu bara. Hal ini disebabkan karena akibat konflik tersebut harga batu bara mencapai nilai USD446 per metrik. Ini berbanding lurus dengan harga saham PTBA yang juga mengalami peningkatan hingga +3,95 persen sejumlah Rp 3.680/unit (Cahyakemala dan Putri, 2022).

Dengan melambungnya harga batu bara menjadi sebuah peluang baik untuk negara Indonesia memulihkan kondisi ekonomi. Konflik yang terjadi antara Ukraina-Rusia yang menyebabkan harga batu bara terus meroket. Manajer diperlukan untuk mengelola proyeksi terkait gencar melakukan ekspor batu bara sebagai bentuk pemanfaat peluang ini. PTBA juga berupaya untuk tetap mempertahankan kebutuhan dalam negeri sebesar 25% atas domestik market obligation (DMO) (Cahyakemala dan Putri, 2022).

Peneliti memilih topik saham PTBA karena harga batu bara sekarang ini didorong oleh meningkatnya permintaan. Permintaan diyakini akan melonjak setelah negara-negara Eropa memutuskan untuk menggunakan kembali batu bara sebagai sumber pembangkit listrik meraka. Keputusan tersebut diambil setelah Rusia (memulai penyerangan terhadap Ukraina) memangkas pasokan gas alam cair ke banyak negara Eropa pada akhir bulan Februari 2022. PTBA adalah Badan Usaha Milik Negara (BUMN) yang fokus pada bisnis pertambangan batu bara di Indonesia. PT. Bukit ASAM (Persero) Tbk merupakan perusahaan industri yang berada di Tanjung Enim, Sumatera Selatan. Saham yang berada di lingkup pertambangan saat ini menjadi komoditas yang paling tidak pernah sepi dari pandangan investor di (BEI) Bursa Efek Indonesia, sebab harga saham batubara saat ini sangat menjanjikan keuntungan. Akan tetapi, di samping itu, berinvestasi pada sektor pertambangan pasti tidak terlepas dari namanya resiko, yang dimana komoditas saham tambang ini sangat ditentukan oleh pasokan atau tersedianya bahan baku tambang, demand dari barang tambang tersebut, dan yang tidak kalah penting peran pemerintah dalam menentukan kebijakan terkait harga barang tambang yang mengakibatkan sahamsektor pertambangan harganya cenderung naik turun (Junaid et al. 2020).

Kinerja saham PTBA menjadi cerminan kepercayaan masyarakat dunia terhadap keberadaan mata uang rupiah. Pada negara berkembang seperti Indonesia, kebijakan pemerintah dinilai memiliki kontrol besar, terutama pada perusahaan BUMN. Oleh sebab itu, Harga Batu Bara Acuan (HBA) dianggap dapat mewakili sisi makroekonomi Indonesia, yaitu penetapan harga tersebut sebagai acuan yang harus dipatuhi oleh seluruh perusahaan pertambangan batu bara di tanah air (Anindita, 2017; Pratama & Yulianto, 2016; Sundari, 2015). Pada saat Harga Batu Bara (HBA) naik, maka saham pada sektor-sektor batu bara secara otomatis akan meningkat. Sangatlah penting untuk mengetahui pergerakan HBA karena dapat memprediksi iklim investasi. Selain itu, HBA juga menentukan potensi investasi suatu negara, terutama Indonesia. Investor menginginkan kepastian HBA pada kisaran harga tertentu untuk melihat apakah kegiatan investasi yang akan mereka lakukan dalam perusahaan pertambangan batu bara tersebut menjanjikan hasil sesuai yang diharapkan.

B. Tinjauan Pustaka

A. Data Stasioner

Suatu data dikatakan memiliki sifat stasioner apabila data tersebut tidak memiliki kecenderungan terhadap trend tertentu. Atau dengan kata lain, apabila fluktuasi data berada disekitar suatu nilai rata-rata yang konstan, tidak tergantung pada waktu dan variansi dari fluktuasi tersebut.

B. Single Moving Average (SMA)

Metode single moving average adalah metode peramalan yang menggunakan sejumlah data aktual permintaan yang baru untuk membangkitkan nilai ramalan untuk permintaan dimasa yang akan datang. Metode moving average cocok digunakan untuk data jangka panjang (Maulina et al. 2019).

C. Double Moving Average (DMA)

Metode Double Moving Average adalah suatu variasi dari prosedur rata-rata bergerak yang diharapkan dapat mengatasi adanya faktor tren secara lebih baik (Iskandar dan Layakana 2020).

D. Single Exponential Smoothing (SES)

Metode Single Exponential Smoothing digunakan untuk peramalan jangka pendek, biasanya hanya satu bulan ke depan. Model ini mengasumsikan bahwa data berfluktuasi dan rata-rata cukup stabil (tidak ada pola pertumbuhan tren atau konsisten). Metode Single Exponential Smoothing lebih cocok digunakan untuk meramalkan hal-hal yang fluktuasinya secara random atau tidak teratur (Jonnius 2016)

E. Double Exponential Smoothing (DES)

Peramalan eksponensial ganda yang berbasis Double Exponential Smoothing (DES) model time series tertentu menggunakan persamaan regresi linier sederhana di mana mereka mencegat intercept b0 dan slope kemiringan b1 yang bervariasi perlahan dari waktu ke waktu. Metode ini digunakan ketika data menunjukkan trend. Pemulusan eksponensial ini lebih banyak dibandingkan smoothing sederhana kecuali bahwa dua komponen harus diperbarui setiap periode level dan trend. Level dimaksud adalah pemulusan estimasi dari nilai data pada akhir setiap periode. Sementara yang dimaksud trend adalah pemulusan estimasi dari pertumbuhan rata-rata pada setiap akhir periode (Jonnius 2016)

F. Model Autoregressive Integrated Moving Average (ARIMA)

Model ARIMA memberikan pendekatan lain untuk peramalan deret waktu. Pemulusan eksponensial dan model ARIMA adalah dua pendekatan yang paling banyak digunakan untuk peramalan deret waktu dan memberikan pendekatan pelengkap untuk masalah tersebut. Sementara model pemulusan eksponensial didasarkan pada deskripsi tren dan musiman dalam data, model ARIMA bertujuan untuk menggambarkan autokorelasi dalam data (Hyndman dan Athanasopoulos 2018).

Deret waktu {Yt} dikatakan mengikuti model rata-rata bergerak autoregresif terintegrasi jika perbedaan ke-d Wt = dYt adalah proses ARMA stasioner. Jika {Wt} mengikuti model ARMA(p,q), kita katakan bahwa {Yt} adalah proses ARIMA(\(p\), \(d\), \(q\)). Untungnya, untuk tujuan praktis, kita biasanya dapat mengambil \(d\) = 1 atau paling banyak 2.

G. Teknik Interpolasi

Interpolasi merupakan suatu metode atau fungsi matematika untuk menduga nilai pada lokasi-lokasi yang datanya tidak tersedia. Menurut Burrough and McDonell (1998), interpolasi adalah proses memprediksi nilai pada suatu titik yang bukan merupakan titik sampel, berdasarkan pada nilai-nilai dari titik-titik di sekitarnya yang berkedudukan sebagai sampel.

Interpolasi spasial mempunyai dua asumsi yakni atribut data bersifat kontinu di dalam ruang (space) dan atribut tersebut saling berhubungan (dependence) secara spasial (Anderson, 2001). Kedua asumsi tersebut berimplikasi pada logika bahwa pendugaan atribut data dapat dilakukan berdasarkan data dari lokasi-lokasi di sekitarnya dan nilai pada titik-titik yang berdekatan akan lebih mirip daripada nilai dari titik-titik yang berjauhan (Prasasti, Wijayanto, Christanto, 2005).

G. Teknik Evaluasi Kesesuaian Model

Sum of Square Error (SSE)

Sum of Square Error (SSE) merupakan hasil penjumlahan dari seluruh jarak masing-masing data dengan titik pusat clusternya. Semakin kecil nilai SSE yang didapat, semakin seragam data yang ada didalam masing-masing cluster, semakin baik cluster yang dihasilkan.

Mean Square Error (MSE)

Mean Squared Error (MSE) adalah Rata-rata Kesalahan kuadrat antara nilai aktual dan nilai peramalan. Metode Mean Squared Error secara umum digunakan untuk mengecek estimasi berapa nilai kesalahan pada peramalan. Nilai Mean Squared Error yang rendah atau nilai mean squared error mendekati nol menunjukkan bahwa hasil peramalan sesuai dengan data aktual dan bisa dijadikan untuk perhitungan peramalan di periode mendatang.

Mean Absolute Percentage Error

Mean Absolut Percentage error (MAPE) adalah persentase kesalahan rata-rata secara mutlak (absolut). Pengertian Mean Absolute Percentage Error adalah Pengukuran statistik tentang akurasi perkiraan (prediksi) pada metode peramalan. Metode Mean Absolute Percentage Error (MAPE) memberikan informasi seberapa besar kesalahan peramalan dibandingkan dengan nilai sebenarnya dari series tersebut. Semakin kecil nilai presentasi kesalahan (percentage error) pada MAPE maka semakin akurat hasil peramalan tersebut.

Root Mean Square Error (RMSE)

Root Mean Square Error (RMSE) merupakan besarnya tingkat kesalahan hasil prediksi, dimana semakin kecil (mendekati 0) nilai RMSE maka hasil prediksi akan semakin akurat. Root Mean Squared Error (RMSE) merupakan salah satu cara untuk mengevaluasi model regresi linear dengan mengukur tingkat akurasi hasil perkiraan suatu model.

Data

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)
library(lmtest)
library(GGally)
library(aTSA)
library(bsts)
library(keras)
library(dLagM)
library(dynlm)

Dataset yang digunakan adalah data harga penutupan saham PTBA dengan rentang mingguan (17 Agustus 2020 - 15 Agustus 2022) dengan jumlah 105 amatan. Data ini diperolah dari Yahoo Finance <https://finance.yahoo.com/quote/PTBA.JK>.

#Input Data
ptba<- read_csv("C:/SMT  5/MPDW/1/Tugas Projek/PTBA.JK (1).csv")
head(ptba)

## # A tibble: 6 × 7
##   Date        Open  High   Low Close `Adj Close`    Volume
##   <date>     <dbl> <dbl> <dbl> <dbl>       <dbl>     <dbl>
## 1 2020-08-17  2110  2110  2110  2110       1742.         0
## 2 2020-08-24  2090  2170  2070  2130       1758. 152538300
## 3 2020-08-31  2150  2170  2020  2100       1733. 185101500
## 4 2020-09-07  2100  2130  1810  1970       1626. 178527600
## 5 2020-09-14  2000  2070  2000  2000       1651. 103866200
## 6 2020-09-21  2020  2040  1940  1990       1643.  69160500

str(ptba)

## spec_tbl_df [105 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Date     : Date[1:105], format: "2020-08-17" "2020-08-24" ...
##  $ Open     : num [1:105] 2110 2090 2150 2100 2000 ...
##  $ High     : num [1:105] 2110 2170 2170 2130 2070 2040 2020 2040 2070 2080 ...
##  $ Low      : num [1:105] 2110 2070 2020 1810 2000 ...
##  $ Close    : num [1:105] 2110 2130 2100 1970 2000 ...
##  $ Adj Close: num [1:105] 1742 1758 1733 1626 1651 ...
##  $ Volume   : num [1:105] 0.00 1.53e+08 1.85e+08 1.79e+08 1.04e+08 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Date = col_date(format = ""),
##   ..   Open = col_double(),
##   ..   High = col_double(),
##   ..   Low = col_double(),
##   ..   Close = col_double(),
##   ..   `Adj Close` = col_double(),
##   ..   Volume = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

dim(ptba)

## [1] 105   7

ptba$Date <- as.Date(ptba$Date,"%y/%m/%d")
Yt<-ptba$Close

Analisis Eksplorasi data

A. Pembagian Data Latih dan Data Uji

Dilakukan pembagian data latih ( Data Training ) dan data uji ( Data Test ). Pada data harga saham PTBA ini akan dilakukan pembagian jumlah data latih dan data uji dengan proporsi 85:15, yaitu data latih sebanyak 86 observasi dan 19 observasi sisanya menjadi data uji.

test <- tail(ptba,19)
train <- head(ptba,(length(ptba$Close)-length(test$Close)))

# Data Time Series
train.ts<-ts(train$Close)
test.ts<-ts(test$Close)
data.ts<-ts(ptba$Close)

B. Plot Time Series

# Plot Time Series
plot(x = ptba$Date,
     y = ptba$Close,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Close Stock PTBA",
     main = "Time Series Plot of Stock PTBA")

# Plot TIme Series Data Training
plot(x = train$Date,
     y = train$Close,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Close Stock PTBA",
     main = "Time Series Plot of Stock PTBA (Data Training)")

# Plot Time Series Data Testing
plot(x = test$Date,
     y = test$Close,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Tahun",
     ylab = "Close Stock PTBA",
     main = "Time Series Plot of Stock PTBA (Data Testing)")

Berdasarkan Gambar 1, pada akhir 2020 menuju 2021 terjadi kenaikan harga saham PTBA, tetapi pada 2021 sempat terjadi penurunan karena efek konflik Rusia dan Ukraina. Pada 2022 harga sahah PTBA mengalami peningkatan yang ekstrem karena pemulihan ekonomi global dan konflik Rusia dan Ukraina yang mereda.

Pemulusan

Parameter Optimum

# Fungsi Mencari Parameter Optimum
smooth <- function(x,min,max,method=c("SMA","DMA"),
                      metrik=c("SSE","MSE","RMSE","MAD","MAPE")){
  df.master <- data.frame()
  if(method=="SMA"){
    for(i in min:max) {
      sma <- TTR::SMA(x,i)
      forecast <- c(NA,sma)
      df <- cbind(aktual=c(x,NA),mulus=c(sma,NA),forecast)
      error <- df[,1]-df[,3]
      sse <- sum(error^2,na.rm=T)
      mse <- mean(error^2,na.rm=T)
      rmse <- sqrt(mse)
      mad <- mean(abs(error),na.rm=T)
      rerror <- error/df[,1]*100
      mape <- mean(abs(rerror),na.rm=T)
      ak <- data.frame(n=i,SSE=sse,MSE=mse,RMSE=rmse,MAD=mad,MAPE=mape)
      df.master <- rbind(df.master,ak)
       ts.plot(x,col="gray",main=paste("Plot SMA n =",i),xlab="Periode Waktu",ylab="Close PTBA")
      lines(sma,col="blue",lwd=2)
      lines(forecast,col="red",lwd=1)
    }
  }
  else if(method=="DMA") {
    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
      smooth_dma <- At+Bt
      forecast_dma <- c(NA, smooth_dma)
      df_dma <- cbind(df_aktual=c(x,NA), smooth_dma=c(smooth_dma,NA), forecast_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("Plot DMA n =",i),xlab="Periode Waktu",ylab="Close PTBA")
      lines(smooth_dma,col="blue",lwd=2)
      lines(forecast_dma,col="red",lwd= 2)
    }
  }
  opt <- df.master$n[which.min(df.master[,metrik])]
  return(list(Akurasi=df.master,n.optimum=paste("n optimum adalah",opt)))
}

A. Single Moving Average (SMA)

smooth(train.ts,min=2,max=10,method = "SMA",metrik = "MAPE")

## $Akurasi
##    n     SSE      MSE     RMSE      MAD     MAPE
## 1  2 1679388 19992.71 141.3956 103.2738 3.984750
## 2  3 2124200 25592.77 159.9774 118.6345 4.561457
## 3  4 2546875 31059.45 176.2369 128.5671 4.929348
## 4  5 3023834 37331.28 193.2131 140.1975 5.365266
## 5  6 3553860 44423.25 210.7682 154.8542 5.878733
## 6  7 4086804 51731.70 227.4460 166.8535 6.293753
## 7  8 4575364 58658.51 242.1952 178.1571 6.694927
## 8  9 5038742 65438.21 255.8089 190.2453 7.140052
## 9 10 5476861 72063.96 268.4473 201.3816 7.543250
## 
## $n.optimum
## [1] "n optimum adalah 2"

Parameter optimum yang didapat untuk single moving average (SMA) berdasarkan metrik MSE adalah n = 2.

Pemulusan dan Peramalan SMA n=2 Data Training

data.sma<-TTR::SMA(train.ts,  n=2)
data.ramal<-c(NA,data.sma)
data.gab<-cbind(Aktual=c(train.ts, NA), Pemulusan=c(data.sma, NA), Peramalan=data.ramal)

# Plot Time Series SMA n=2
ts.plot(data.gab[,1], xlab="Periode Waktu", ylab="Close PTBA", main= "Plot SMA n=2 (Data Training)")
points(data.gab[,1])
lines(data.gab[,2],col="blue",lwd=2)
lines(data.gab[,3],col="red",lwd=2)
legend("bottomright",c("Data Aktual","Data Pemulusan","Ramalan"), lty=c(3,1,1), col=c("black","blue","red"), cex=0.7)

Berdasarkan plot pemulusan SMA dengan n=2 terlihat bahwa plot pemulusan cukup baik karena plot yang dihasilkan sudah cukup mulus dan optimum.

Menghitung Nilai Keakuratan n=2

error.sma <- data.gab[, 1] - data.gab[, 3] #SSE (Sum Square Error)
SSE.sma <- sum(error.sma^2, na.rm = T) # MSE (Mean Squared Error)
MSE.sma <- mean(error.sma^2, na.rm = T) # RMSE (Root Mean Square Error)
RMSE.sma <- sqrt(mean(error.sma^2, na.rm = T)) # MAD (Mean Absolute Deviation)
MAD.sma <- mean(abs(error.sma), na.rm = T) # MAPE (Mean Absolute Percentage Error)
r.error.sma <- (error.sma/data.gab[, 1])*100 # Relative Error
MAPE.sma <- mean(abs(r.error.sma), na.rm = T)

Pemulusan dan Peramalan SMA n=2 Data Testing

data.sma1<-TTR::SMA(test.ts,  n=2)
data.ramal1<-c(NA,data.sma1)
data.gab1<-cbind(Aktual=c(test.ts, NA), Pemulusan=c(data.sma1, NA), Peramalan=data.ramal1)

#Plot Time Series Data Testing n=2
ts.plot(data.gab1[,1], xlab="Periode Waktu", ylab="Close PTBA", main= "Plot SMA n=2 (Data Testing)")
points(data.gab1[,1])
lines(data.gab1[,2],col="blue",lwd=2)
lines(data.gab1[,3],col="red",lwd=2)
legend("bottomright",c("Data Aktual","Data Pemulusan","Ramalan"), lty=c(3,1,1), col=c("black","blue","red"), cex=0.7)

Berdasarkan plot pemulusan SMA dengan n=2 terlihat bahwa plot pemulusan cukup baik karena plot yang dihasilkan sudah cukup mulus dan optimum.

Menghitung Nilai Keakuratan n=2

error.sma1 <- data.gab1[, 1] - data.gab1[, 3] 
SSE.sma1 <- sum(error.sma1^2, na.rm = T) # SSE (Sum Square Error)
MSE.sma1 <- mean(error.sma1^2, na.rm = T) # MSE (Mean Squared Error)
RMSE.sma1 <- sqrt(mean(error.sma1^2, na.rm = T)) # RMSE (Root Mean Square Error)
MAD.sma1 <- mean(abs(error.sma1), na.rm = T) # MAD (Mean Absolute Deviation)
r.error.sma1 <- (error.sma1/data.gab1[, 1])*100 # Relative Error
MAPE.sma1 <- mean(abs(r.error.sma1), na.rm = T) # MAPE (Mean Absolute Percentage Error)

Perbandingan Training dan Testing

MSESMA <-matrix(c(SSE.sma1, MSE.sma1, MAPE.sma1,SSE.sma, MSE.sma, MAPE.sma ),nrow=3,ncol=2)
row.names(MSESMA)<- c("SSE","MSE","MAPE")
colnames(MSESMA) <- c("Testing","Training")
MSESMA

##           Testing     Training
## SSE  1.373100e+06 1.679388e+06
## MSE  8.077059e+04 1.999271e+04
## MAPE 5.289504e+00 3.984750e+00

B. Double Moving Average (DMA)

smooth(train.ts,min=2,max=10,method = "DMA",metrik = "MAPE")

## $Akurasi
##    n     SSE      MSE     RMSE      MAD     MAPE
## 1  2 1552438 18704.07 136.7628 103.9157 4.029404
## 2  3 2133605 26340.80 162.2985 121.6187 4.658076
## 3  4 2634138 33343.51 182.6021 139.9842 5.355202
## 4  5 3314083 43040.04 207.4609 160.7065 6.085056
## 5  6 4137062 55160.83 234.8634 185.7741 7.046425
## 6  7 4906215 67208.42 259.2459 205.7898 7.824720
## 7  8 5310166 74791.07 273.4796 210.9595 8.032777
## 8  9 5418149 78523.91 280.2212 211.9216 8.135828
## 9 10 5559545 82978.28 288.0595 216.8507 8.402904
## 
## $n.optimum
## [1] "n optimum adalah 2"

Parameter optimum yang didapat untuk double moving average (DMA) berdasarkan metrik MSE adalah n = 2.

Pemulusan dan Peramalan Data Training

#SMA
data.sma3<-TTR::SMA(train.ts,  n=2)

#DMA N=2
data.dma3 <- TTR::SMA(data.sma3, n=2)
At3 <- 2*data.sma3-data.dma3
Bt3 <- data.sma3-data.dma3
pemulusan_dma3 <- At3+Bt3
ramal_dma3 <- c(NA, pemulusan_dma3)
df_dma3 <- cbind(df_aktual=c(train.ts,NA), pemulusan_dma=c(pemulusan_dma3,NA), ramal_dma3)

Plot DMA N=2

ts.plot(df_dma3[,1], xlab="Periode Waktu", ylab="Close PTBA",  main= "Plot DMA n=2 (Data Training)")
points(df_dma3[,1])
lines(df_dma3[,2],col="blue",lwd=2)
lines(df_dma3[,3],col="red",lwd= 2)
legend("bottomright", c("Data Aktual","Pemulusan","Ramalan"),lty=c(3,1,1),col=c ("black","blue","red"), cex=0.8)

Ukuran Keakuratan

error.dma3 <- df_dma3[, 1] - df_dma3[, 3] 
SSE.dma3 <- sum(error.dma3^2, na.rm = T) # SSE (Sum Square Error)
MSE.dma3 <- mean(error.dma3^2, na.rm = T)# MSE (Mean Squared Error)
RMSE.dma3 <- sqrt(mean(error.dma3^2, na.rm = T))# RMSE (Root Mean Square Error) 
MAD.dma3 <- mean(abs(error.dma3), na.rm = T)# MAD (Mean Absolute Deviation) 
r.error.dma3 <- (error.dma3/df_dma3[, 1])*100 # Relative Error
MAPE.dma3 <- mean(abs(r.error.dma3), na.rm = T) # MAPE (Mean Absolute Percentage Error)

Pemulusan dan Peramalan Data Testing

#SMA
data.sma4<-TTR::SMA(test.ts,  n=2)

#Double Moving  Average dengan N=2
data.dma4 <- TTR::SMA(data.sma4, n=2)
At4 <- 2*data.sma4-data.dma4
Bt4 <- data.sma4-data.dma4
pemulusan_dma4 <- At4+Bt4
ramal_dma4 <- c(NA, pemulusan_dma4)
df_dma4<- cbind(df_aktual=c(test.ts,NA), pemulusan_dma=c(pemulusan_dma4,NA), ramal_dma4)

ts.plot(df_dma4[,1], xlab="Periode Waktu", ylab="Close PTBA",  main= "Plot DMA n=2 (Data Testing)")
points(df_dma4[,1])
lines(df_dma4[,2],col="blue",lwd=2)
lines(df_dma4[,3],col="red",lwd= 2)
legend("bottomright", c("Data Aktual","Pemulusan","Ramalan"),lty=8,col=c ("black","blue","red"), cex=0.7)

Ukuran Keakuratan

error.dma4 <- df_dma4[, 1] - df_dma4[, 3] 
SSE.dma4 <- sum(error.dma4^2, na.rm = T) # SSE (Sum Square Error)
MSE.dma4 <- mean(error.dma4^2, na.rm = T) # MSE (Mean Squared Error)
RMSE.dma4 <- sqrt(mean(error.dma4^2, na.rm = T)) # RMSE (Root Mean Square Error)
MAD.dma4 <- mean(abs(error.dma4), na.rm = T) # MAD (Mean Absolute Deviation)
r.error.dma4 <- (error.dma4/df_dma4[, 1])*100 # Relative Error
MAPE.dma4 <- mean(abs(r.error.dma4), na.rm = T) # MAPE (Mean Absolute Percentage Error)

Perbandingan Training dan Testing

MSEDMA <-matrix(c(SSE.dma4, MSE.dma4,MAPE.dma4, SSE.dma3, MSE.dma3, MAPE.dma3),nrow=3,ncol=2)
row.names(MSEDMA)<- c("SSE","MSE","MAPE")
colnames(MSEDMA) <- c("Testing","Training")
MSEDMA

##           Testing     Training
## SSE  1.883575e+06 1.552438e+06
## MSE  1.177234e+05 1.870407e+04
## MAPE 6.162142e+00 4.029404e+00

C. Single Exponential Smoothing

pemulusan2 <- function(x,alfa=NULL,beta=NULL,method=c("SES","DES"),
                       metrik=c("SSE","MSE","RMSE")){
  df.master <- data.frame()
  if(method=="SES"){
    for(i in alfa){
      df_ses <- HoltWinters(x, alpha = i, beta=F, gamma=F)
      sse <- df_ses$SSE 
      mse <- sse/length(x)
      rmse <-sqrt(mse)
      ak <- data.frame(Alpha=i,SSE=sse,MSE=mse,RMSE=rmse)
      df.master <- rbind(df.master,ak)
      datases <- data.frame(x, c(NA, df_ses$fitted[,1]))
      colnames(datases) = c("y","yhat")
      ts.plot(x, xlab="Periode  waktu", ylab="Close PTBA", col="blue", lty=3,main=paste("Plot SES alpha=",i))
      points(x)
      lines(datases[,2], col="red",lwd=2) #nilai dugaan
    }
  }
  else if(method=="DES"){
    for (i in alfa) {
      for (j in beta) {
        df_des <- HoltWinters(x, alpha = i, beta=j, gamma=F)
        sse <- df_des$SSE 
        mse <- sse/length(x)
        rmse <-sqrt(mse)
        ak <- data.frame(Alpha=i,Beta=j,SSE=sse,MSE=mse,RMSE=rmse)
        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="Close PTBA",  col="blue", lty=3,main=paste("Plot DES alpha=",i,"beta=",j))
        points(x)
        lines (datades[,2], col="red",lwd=2)
      }
    }
  }
  if(method=="SES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt)))
  }
  else if(method=="DES"){
    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)))
  }
}

pemulusan2(x=train.ts,alfa=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
           method="SES",metrik="MSE")

## $Akurasi
##   Alpha     SSE      MSE     RMSE
## 1   0.1 6707042 77988.86 279.2649
## 2   0.2 3911282 45480.02 213.2605
## 3   0.3 2751679 31996.26 178.8750
## 4   0.4 2152952 25034.33 158.2224
## 5   0.5 1799164 20920.52 144.6393
## 6   0.6 1571777 18276.48 135.1905
## 7   0.7 1417351 16480.83 128.3777
## 8   0.8 1308612 15216.42 123.3548
## 9   0.9 1230692 14310.38 119.6260
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.9"

Dengan melakukan looping mulai dari 0.1 hingga 0.9 diperoleh nilai \(\alpha\) optimum = 0.9.

#SES dengan alpha=0.9
df_ses <- HoltWinters(train.ts, alpha = 0.9, beta=F, gamma=F)

#nilai keakuratan
sse2 <- df_ses$SSE 
mse2 <- sse2/length(train.ts)
rmse2 <-sqrt(mse2)
akurasi <- c("SSE"=sse2, "MSE"=mse2, "RMSE"=rmse2)
akurasi

##         SSE         MSE        RMSE 
## 1230692.499   14310.378     119.626

datases <- data.frame(train.ts, c(NA, df_ses$fitted[,1]))
colnames(datases) = c("y","yhat")

#Plot 
ts.plot(train.ts, xlab="Periode  Waktu", ylab="Close PTBA", col="blue", lty=3)
points(train.ts)
lines(datases[,2], col="red",lwd=2) #nilai dugaan
title("Plot SES alpha=0.9 (Data Training)",cex.main=1, font.main=2 ,col.main="black")
legend("bottomright", c("Data aktual", "Data Pemulusan"), lty=c(3,1),col=c ("blue","red"), cex=0.7)

Berikut merupakan hasil peramalan dengan alpha dan beta optimum untuk 19 periode ke depan.

ramal_ses <- forecast::forecast(df_ses,h=19) #periode ramalan sebanyak 19
df_ramal_ses <- data.frame(ramal_ses)

plot(ramal_ses,xlab="Periode  Waktu", ylab="Close PTBA", main="Plot Peramalan SES alpha=0.9 (Data Training)")
lines(data.ts,lty=1,col="red")
points(data.ts,cex=0.5)

#Gabungan Data Aktual, Pemulusan, dan Ramalan
data.ses <- cbind(aktual=c(train.ts, rep(NA,19)),
            pemulusan=c(df_ses$fitted[,1], as.numeric(df_ses$coefficients), rep(NA,19)),
            ramalan = c(NA, df_ses$fitted[,1], df_ramal_ses$Point.Forecast))
data.ses <- ts(data.ses)

#Plot
plot(data.ses[, 1], 
     xlab = "Periode Waktu", 
     ylab = "Close PTBA", 
     main="Plot SES alpha=0.9 (Data Training)", 
     lty = c(3,1,1))
points(data.ses[, 1])
lines(data.ses[, 2],col="red")
lines(data.ses[, 3],col="blue")
legend("bottomright", legend=c("Data Aktual", "Data Pemulusan", "Data Peramalan"), lty = c(3,1,1), col = c("black", "red", "blue"),cex = 0.7)

Data Testing

df_sestest <- HoltWinters(test.ts, alpha = 0.9, beta=F, gamma=F)

ramal_sestest <- forecast::forecast(df_sestest,h=19) #periode ramalan sebanyak 19
df_ramal_sestest <- data.frame(ramal_sestest)
plot(ramal_sestest,xlab="Periode Waktu",ylab="Close PTBA",main = "Plot Peramalan SES alpha=0.9 (Data Testing)")

#nilai keakuratan
sse5 <- df_ses$SSE 
mse5 <- sse5/length(test.ts)
rmse5 <-sqrt(mse5)

bandingses <-matrix(c(sse2, mse2,rmse2, sse5,mse5,rmse5),nrow=3,ncol=2)
row.names(bandingses)<- c("SSE","MSE","RMSE")
colnames(bandingses) <- c("Training","Testing")
bandingses

##         Training     Testing
## SSE  1230692.499 1230692.499
## MSE    14310.378   64773.289
## RMSE     119.626     254.506

D. Double Exponential Smoothing

pemulusan2(x=train.ts,alfa=c(0.9),
           beta=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),method="DES",metrik="MSE")

## $Akurasi
##   Alpha Beta     SSE      MSE     RMSE
## 1   0.9  0.1 1245296 14480.19 120.3336
## 2   0.9  0.2 1284736 14938.79 122.2244
## 3   0.9  0.3 1329843 15463.29 124.3515
## 4   0.9  0.4 1378288 16026.60 126.5962
## 5   0.9  0.5 1429680 16624.19 128.9348
## 6   0.9  0.6 1484167 17257.76 131.3688
## 7   0.9  0.7 1542011 17930.35 133.9043
## 8   0.9  0.8 1603758 18648.35 136.5590
## 9   0.9  0.9 1670333 19422.48 139.3646
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.9 ,Beta optimum adalah 0.1"

Dari grafik dapat dilihat bahwa dengan menggunakan metode DES pada \(\alpha\) = 0.9 dan \(\beta\) = 0.1 (optimum) telah mengikuti trend dari data sebenarnya.

# Double Exponential Smooting Alpha dan Beta Optimum
df_des <- HoltWinters(train.ts, alpha = 0.9, beta=0.1, gamma=F)

# Nilai Keakuratan
sse7 <- df_des$SSE 
mse7 <- sse7/length(train.ts)
rmse7 <-sqrt(mse7)
akurasi <- c("SSE"=sse7, "MSE"=mse7, "RMSE"=rmse7)
akurasi

##          SSE          MSE         RMSE 
## 1245295.9295   14480.1852     120.3336

datades <- data.frame(train.ts, c(NA, NA, df_des$fitted[,1]))
colnames(datades) = c("y","yhat")
head(datades)

##      y     yhat
## 1 2110       NA
## 2 2130       NA
## 3 2100 2150.000
## 4 1970 2120.500
## 5 2000 1987.005
## 6 1990 2001.825

# Plot Alpha dan Beta Optimum
par(mfrow=c(1,1))
ts.plot(train.ts,xlab="Periode  Waktu", ylab="Yt",  col="blue", lty=3)
points(train.ts)
lines (datades[,2], col="red",lwd=2)
title("Double Exponential Smoothing (Alpha=0.9 dan Beta 0.1)", cex.main=1, font.main=4 ,col.main="black")
legend("topright", c("Data aktual","Fitted DES"), lty=1:3,col=c ("blue","red"),cex=0.7)

Dari grafik dapat dilihat bahwa dengan menggunakan metode DES pada \(\alpha\) = 0.9 dan \(\beta\) = 0.1 (optimum) telah mengikuti trend dari data sebenarnya.

Berikut merupakan hasil peramalan dengan \(\alpha\) dan \(\beta\) optimum untuk 19 periode ke depan.

ramal_des <- forecast::forecast(df_des,h=19) #periode ramalan sebanyak 19
(df_ramal_des <- data.frame(ramal_des))

##     Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
##  87       3633.036 3476.182 3789.890 3393.148 3872.923
##  88       3689.325 3468.606 3910.044 3351.765 4026.885
##  89       3745.614 3467.381 4023.848 3320.093 4171.136
##  90       3801.904 3468.597 4135.210 3292.155 4311.652
##  91       3858.193 3470.697 4245.689 3265.569 4450.817
##  92       3914.482 3472.902 4356.061 3239.144 4589.819
##  93       3970.771 3474.778 4466.765 3212.214 4729.328
##  94       4027.060 3476.058 4578.062 3184.375 4869.745
##  95       4083.349 3476.578 4690.121 3155.372 5011.326
##  96       4139.639 3476.228 4803.049 3125.039 5154.238
##  97       4195.928 3474.937 4916.919 3093.267 5298.589
##  98       4252.217 3472.657 5031.777 3059.983 5444.450
##  99       4308.506 3469.359 5147.653 3025.142 5591.871
## 100       4364.795 3465.023 5264.567 2988.713 5740.878
## 101       4421.085 3459.639 5382.530 2950.680 5891.489
## 102       4477.374 3453.200 5501.547 2911.035 6043.712
## 103       4533.663 3445.707 5621.619 2869.777 6197.548
## 104       4589.952 3437.161 5742.744 2826.910 6352.994
## 105       4646.241 3427.567 5864.916 2782.439 6510.043

plot(ramal_des, xlab="periode  waktu", ylab="Yt", main="Forecast from DES")
lines(data.ts,lty=1,col="red")
points(data.ts,cex=0.5)
legend("topleft", legend = c("Data Peramalan", "Data Aktual"), lty = 1, col=c("blue","red"),cex=0.8)

#Gabungan Data Aktual, Pemulusan, dan Ramalan
data.des <- cbind(aktual=c(train.ts, rep(NA,19)),
                  pemulusan=c(NA, df_des$fitted[,2], as.numeric(df_des$coefficients[1]+df_des$coefficients[2]), rep(NA,19)),
                  ramalan = c(NA, NA, df_des$fitted[,1], df_ramal_des$Point.Forecast))
data.des <- ts(data.des)

#Plot
plot(data.ts, 
     xlab = "Time Period", 
     ylab = "Close PTBA", 
     main="Forecasting Using Double Exponential Smoothing", 
     lty = c(3,1,1),
     col="black")
points(data.ts)
lines(data.des[, 2], col="red")
lines(data.des[, 3], col="blue")

legend("bottomleft", legend=c("Data Aktual", "Data Pemulusan", "Data Peramalan"), lty = c(3,1,1), col = c("black", "red", "blue"),cex = 0.8)

Nilai Akurasi dari hasil peramalan data testing adalah sebagai berikut:

df_des <- HoltWinters(test.ts, alpha = 0.9, beta=0.1, gamma=F)
df_des

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = test.ts, alpha = 0.9, beta = 0.1, gamma = F)
## 
## Smoothing parameters:
##  alpha: 0.9
##  beta : 0.1
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a 4204.258695
## b    8.773859

df_des$fitted

## Time Series:
## Start = 3 
## End = 19 
## Frequency = 1 
##        xhat    level      trend
##  3 3610.000 3690.000 -80.000000
##  4 3737.900 3799.000 -61.100000
##  5 3758.079 3811.790 -53.711000
##  6 3745.870 3795.808 -49.938110
##  7 3997.021 4019.587 -22.566391
##  8 4393.204 4377.702  15.501756
##  9 4643.134 4606.320  36.813413
## 10 3964.045 3992.313 -28.268629
## 11 3704.072 3753.404 -49.332656
## 12 4036.808 4051.407 -14.599120
## 13 3708.569 3751.681 -43.111846
## 14 3884.674 3907.857 -23.183052
## 15 3916.264 3934.467 -18.203698
## 16 4030.459 4036.626  -6.167430
## 17 4291.137 4273.046  18.091266
## 18 4139.803 4137.114   2.688921
## 19 4152.587 4148.980   3.606684

#nilai keakuratan
sse6 <- df_des$SSE 
mse6 <- sse6/length(test.ts)
rmse6 <-sqrt(mse6)
akurasi <- c("SSE"=sse6, "MSE"=mse6, "RMSE"=rmse6)
akurasi

##         SSE         MSE        RMSE 
## 1382402.966   72758.051     269.737

Perbandingan Pemulusan

Akurasi data Testing

MSEfull <-matrix(c(SSE.sma1, MSE.sma1, RMSE.sma1,SSE.dma4, MSE.dma4, RMSE.dma4, sse5, mse5,rmse5, sse6, mse6, rmse6  ),nrow=3,ncol=4)
row.names(MSEfull)<- c("SSE","MSE","MAPE")
colnames(MSEfull) <- c("SMA","DMA", "SES", "DES")
MSEfull

##               SMA          DMA         SES         DES
## SSE  1373100.0000 1883575.0000 1230692.499 1382402.966
## MSE    80770.5882  117723.4375   64773.289   72758.051
## MAPE     284.2017     343.1085     254.506     269.737

Akurasi Data Training

MSEfull <-matrix(c(SSE.sma, MSE.sma, RMSE.sma, SSE.dma3, MSE.dma3, RMSE.dma3, sse2, mse2, rmse2,sse7, mse7, rmse7 ),nrow=3,ncol=4)
row.names(MSEfull)<- c("SSE","MSE","RMSE")
colnames(MSEfull) <- c("SMA","DMA", "SES", "DES")
MSEfull

##               SMA          DMA         SES          DES
## SSE  1679387.5000 1552437.5000 1230692.499 1245295.9295
## MSE    19992.7083   18704.0663   14310.378   14480.1852
## RMSE     141.3956     136.7628     119.626     120.3336

Pemodelan ARIMA

A. Pemeriksaan Kestasioneran Data

Harga penutupan saham PTBA berfluktuatif yang mengindikasikan bahwa data tidak stasioner dalam rataan. Selain tidak stasioner dalam rataan, data saham PTBA juga tidak stasioner dalam ragam karena melebar dan menyebar secara tidak seimbang.

Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\): \(\rho = 0\) (terdapat akar unit dalam data, data tidak stasioner)

\(H_1\): \(\rho \neq 0\) (tidak terdapat akar unit dalam data, data stasioner)

Pada pengujian ini, nilai \(\alpha\) yang akan digunakan sebesar 0.05 / 5%.

Kestasioneran dalam Rataan

#plotacf
par(mfrow=c(1,2))
acf(train.ts, lag.max = 100)
acf(train.ts, lag.max = 20)

#Uji Formal ADF
adf.test(train.ts)

## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag  ADF p.value
## [1,]   0 1.34   0.953
## [2,]   1 1.05   0.919
## [3,]   2 1.27   0.946
## [4,]   3 1.27   0.947
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -0.203   0.930
## [2,]   1 -0.595   0.837
## [3,]   2 -0.314   0.914
## [4,]   3 -0.613   0.831
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -1.06   0.923
## [2,]   1 -1.47   0.790
## [3,]   2 -1.13   0.913
## [4,]   3 -1.31   0.860
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Berdasarkan Gambar 2A, plot ACF dengan lag max 100 terlihat bahwa plotnya berbentuk Tails off (slowly) yang artinya model tidak stasioner dan perlu proses pembedaan atau differencing terlebih dahulu hingga data menjadi stasioner. Hal ini didukung dengan hasil uji formal yang menghasilkan p-value = 0.913 > \(\alpha\) = 0.05 Tak Tolak H0, yang artinya pola data tidak stasioner. Selanjutnya akan dilakukan pengecekan kestasioneran dalam ragam yang dilakukan menggunakan plot Box-Cox sebagai berikut.

Kestasioneran dalam Ragam

#Plot Box-Cox
b <- boxcox(train$Close~train$Date, lambda = seq(2,-5))

#Nilai Rounded Lambda
lambda <- b$x[which.max(b$y)] # -2.383838
lambda

## [1] -0.4747475

Menurut Trianto (2015) dalam Haris dan Arum (2020) dan Wahyuni et al. (2021), data dikatakan stasioner terhadap ragam jika nilai \(\lambda\) mendekati nilai satu atau pada selang kepercayaan 95% mencakup \(\lambda\) yang memuat nilai satu. Berdasarkan Gambar 3 Plot Box-Cox harga saham PTBA, terlihat nilai rounded value (\(\lambda\)) adalah sebesar -0.4747475 dan untuk selang kepercayaan 95% nilai \(\lambda\)-nya jauh dari nilai satu sehingga dapat dikatakan bahwa data harga saham PTBA penutupan tidak stasioner dalam ragam sehingga perlu dilakukan penanganan transformasi.

Penanganan Ketidakstasioneran

Wahyuni et al. (2021) melakukan penanganan ketidakstasioneran dalam ragam terlebih dahulu sebelum melakukan penanganan terhadap ketidakstasioneran dalam rataan. Berdasarkan hal tersebut, akan dilakukan transformasi terlebih dahulu. Transformasi ini akan menggunakan nilai \(\lambda\)=-1 yang masih berada pada selang 95%. Kemudian akan dicek kembali melalui plot Box-Cox.

trainTrans <- (train$Close)^(-1)
boxcox(trainTrans~train$Date)

Gambar 4 menunjukkan plot Box-Cox hasil data yang telah ditransformasi. Pada plot terlihat nilai \(\lambda\)=1 sudah berapa pada selang 95%. Akan tetapi, transformasi data ini tidak akan menangani ketidakstasioneran dalam rataan. Selanjutnya akan dilakukan pembedaan dengan \(d\)=1 dilanjutkan dengan melihat plot ACF dan hasil uji ADF.

acf(trainTrans)

train.tr.dif<- diff(trainTrans,differences = 1)
acf(train.tr.dif)

adf.test(train.tr.dif)

## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag   ADF p.value
## [1,]   0 -7.32    0.01
## [2,]   1 -6.34    0.01
## [3,]   2 -4.27    0.01
## [4,]   3 -3.92    0.01
## Type 2: with drift no trend 
##      lag   ADF p.value
## [1,]   0 -7.38    0.01
## [2,]   1 -6.45    0.01
## [3,]   2 -4.43    0.01
## [4,]   3 -4.09    0.01
## Type 3: with drift and trend 
##      lag   ADF p.value
## [1,]   0 -7.37  0.0100
## [2,]   1 -6.44  0.0100
## [3,]   2 -4.40  0.0100
## [4,]   3 -4.08  0.0104
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

Berdasarkan Gambar 5, plot ACF terlihat sudah tidak adanya penurunan nilai lag secara perlahan dan hasil uji ADF menunjukkan p-value = 0.03865 < \(\alpha\)=0.05, artinya Tolak Tolak \(H_0\), cukup bukti untuk menyatakan bahwa data stasioner.

plot(train.tr.dif,
     col = "black",
     lwd = 1,
     type = "o",
     xlab = "Time",
     ylab = "Data",
     main = "PTBA")

Berdasarkan plot data hasil differencing, pola data menjadi stasioner dalam rataan dan lebar pita relatif sama sehingga pola data menjadi stasioner dalam ragam.

B. Seleksi Model

acf(train.tr.dif, lag.max = 20)

pacf(train.tr.dif, lag.max = 20)

eacf(train.tr.dif)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o o o o o o o o o o  o  o  o 
## 4 x o o x o o o o o o o  o  o  o 
## 5 x o o o o o o o o o o  o  o  o 
## 6 o o x x x o o o o o o  o  o  o 
## 7 x x o o x o o o o o o  o  o  o

Berdasarkan plot ACF dan PACF tidak ada cut-off maupun tail-off sehingga sulit ditentukan model tentatifnya. Maka digunakan extended ACF(EACF) untuk melihat model.
Pada tabel EACF, kemungkinan model ARIMA yang dapat dibentuk berdasarkan pola segitiga-nol pada tabel EACF adalah : ARIMA(0,1,0), ARIMA(0,1,1), ARIMA(0,1,2), ARIMA(1,1,1), ARIMA(1,1,2).

C. Pendugaan Parameter

arima1 <- Arima(ts(trainTrans), order=c(0,1,0),method = "ML")
arima2 <- Arima(ts(trainTrans), order=c(0,1,1),method = "ML")
arima3 <- Arima(ts(trainTrans), order=c(0,1,2),method = "ML")
arima4 <- Arima(ts(trainTrans), order=c(1,1,1),method = "ML")
arima5 <- Arima(ts(trainTrans), order=c(1,1,2),method = "ML")

D. Signifikansi Tiap Model

coeftest(arima2) # Model Arima 2

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ma1  0.29498    0.12935  2.2805  0.02258 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(arima3) # Model Arima 3

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ma1  0.278394   0.106281  2.6194 0.008808 **
## ma2 -0.121932   0.097214 -1.2543 0.209747   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(arima4) # model Arima 4

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.77277    0.10079 -7.6673 1.757e-14 ***
## ma1  0.99048    0.12108  8.1802 2.833e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(arima5) # model Arima 5

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.274809   0.471312 -0.5831   0.5598
## ma1  0.541002   0.466302  1.1602   0.2460
## ma2 -0.043903   0.180965 -0.2426   0.8083

Berdasarkan hasil di atas, terdapat satu model dengan seluruh parameter berpengaruh. Sementara, model lainnya hanya menghasilkan beberapa parameter yang signifikan, bahkan tidak ada satupun yang signifikan. Akan tetapi, saat menggunakan fungsi auto.ARIMA, didapatkan bahwa model optimum adalah model ARIMA(0,1,1).

E. Perbandingan AIC

aicmodel <- data.frame("Model"=c("Arima 1","Arima 2","Arima 3", "Arima 4","Arima 5"), "p,d,q"=c("0,1,0","0,1,1", "0,1,2","1,1,1", "1,1,2") ,"AIC"=c(AIC(arima1),AIC(arima2),AIC(arima3),AIC(arima4),AIC(arima5)))
aicmodel

##     Model p.d.q       AIC
## 1 Arima 1 0,1,0 -1617.581
## 2 Arima 2 0,1,1 -1620.575
## 3 Arima 3 0,1,2 -1620.107
## 4 Arima 4 1,1,1 -1620.248
## 5 Arima 5 1,1,2 -1618.392

Berdasarkan perbandingan nilai AIC tersebut, didapatkan hasil bahwa model ARIMA (0,1,1) memiliki nilai AIC terkecil. Dikarenakan model ARIMA(0,1,1) merupakan model yang menghasilkan seluruh peubah signifikan dengan nilai AIC terkecil, serta didukung dengan hasil fungsi auto.ARIMA, maka model ARIMA(0,1,1) dijadikan sebagai model tentatif sehingga akan dilakukan analisis lebih lanjut.

F. Overfitting

overfit.1 <- arima(ts(trainTrans),order=c(0,1,2), method = "ML") # Model Overfitting 1
overfit.1

## 
## Call:
## arima(x = ts(trainTrans), order = c(0, 1, 2), method = "ML")
## 
## Coefficients:
##          ma1      ma2
##       0.2784  -0.1219
## s.e.  0.1063   0.0972
## 
## sigma^2 estimated as 2.874e-10:  log likelihood = 813.05,  aic = -1622.11

overfit.2 <- arima(ts(trainTrans), order=c(1,1,1), method="ML") # Model Overfitting 2
overfit.2

## 
## Call:
## arima(x = ts(trainTrans), order = c(1, 1, 1), method = "ML")
## 
## Coefficients:
##           ar1     ma1
##       -0.7728  0.9905
## s.e.   0.1008  0.1211
## 
## sigma^2 estimated as 2.814e-10:  log likelihood = 813.12,  aic = -1622.25

# Perbandingan AIC Model
Perbadingan.Overfit <- data.frame("Model"=c("Model Terbaik Sebelumnya", "Model Overfitting 1", "Model Overfitting 2"), "p,d,q"=c("1,1,1","1,1,2","2,1,1"),
                         "AIC"=c(arima2$aic,overfit.1$aic,overfit.2$aic))
Perbadingan.Overfit

##                      Model p.d.q       AIC
## 1 Model Terbaik Sebelumnya 1,1,1 -1620.575
## 2      Model Overfitting 1 1,1,2 -1622.107
## 3      Model Overfitting 2 2,1,1 -1622.248

G. Analisis Sisaan

residual <- arima2$residuals
par(mfrow=c(2,2))
# QQ-Plot
qqnorm(residual)
qqline(residual, col="red", lwd=2, lty=1)
plot(residual, type="o",ylab="Sisaan", xlab="Order",main="Plot Sisaan vs Order")
abline(h=0,col="red")
acf(residual, main="Plot ACF Residual")
pacf(residual, main="Plot PACF Residual")

Normal QQ-Plot
Terlihat bahwa sisaan berada di sepanjang garis normal, sehingga mengindikasikan bahwa sisaan menyebar normal.
Residual vs Order
Terlihat tebaran amatan tidak memiliki pola atau acak, serta sisaan mayoritas berada disekitar titik nol, sehingga dapat diasumsikan bahwa sisaan saling bebas atau tidak terdapat autokorelasi
Plot ACF dan PACF
Berdasarkan plot ACF dan PACF terlihat tidak ada garis vertikal di lag tertentu yang melewati garis biru horizontal, sehingga dapat disimpulkan tidak terdapat autokorelasi pada.

Uji Normalitas Sisaan

\(H_0\): Sisaan mengikuti sebaran normal

\(H_1\): Sisaan tidak mengikuti sebaran normal

# Jarque Bera Test
jarque.bera.test(residual)

## 
##  Jarque Bera Test
## 
## data:  residual
## X-squared = 2.8663, df = 2, p-value = 0.2386

# Shapiro-Wilk Test
shapiro.test(residual)

## 
##  Shapiro-Wilk normality test
## 
## data:  residual
## W = 0.98244, p-value = 0.2941

Berdasarkan uji Jarque Bera dan Shapiro-Wilk, keduanya menghasilkan p-value > \(\alpha\) = 0.05 Tak Tolak H0, yang artinya sisaan mengikuti sebaran normal dalam taraf nyata 5%.

Uji Autokorelasi

\(H_0\): Tidak terdapat autokorelasi

\(H_1\): Terdapat autokorelasi

Box.test(residual,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  residual
## X-squared = 0.22538, df = 1, p-value = 0.635

Berdasarkan uji Box-Ljung yang menghasilkan p-value = 0.635 > \(\alpha\) = 0.05 Tak Tolak H0, yang artinya tidak terdapat autokorelasi dalam taraf nyata 5%.

Uji Nilai Tengah

\(H_0\): Nilai tengah sisaan bernilai 0

\(H_1\): Nilai tengah sisaan tidak bernilai 0

t.test(residual, mu = 0, conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  residual
## t = -0.99106, df = 85, p-value = 0.3245
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -5.465742e-06  1.829431e-06
## sample estimates:
##     mean of x 
## -1.818156e-06

Berdasarkan uji t yang menghasilkan p-value = 0.3245 > \(\alpha\) = 0.05 Tak Tolak H0, yang artinya nilai tengah sisaan bernilai 0 dalam taraf nyata 5%.

H. Identifikasi Efek ARCH

Efek ARCH pada sisaan model ARIMA (0,1,1) diidentifikasi dengan melakukan uji Lagrange Multiplier (LM). Hipotesis yang digunakan adalah sebagai berikut: \(H_0\): Tidak ada heteroskedastisitas pada sisaan

\(H_1\): Ada heteroskedastisitas pada sisaan

model2d <- arima(trainTrans, order=c(0,1,1), method="ML")
aTSA::arch.test(model2d)

## ARCH heteroscedasticity test for residuals 
## alternative: heteroscedastic 
## 
## Portmanteau-Q test: 
##      order    PQ p.value
## [1,]     4  3.03   0.553
## [2,]     8  6.04   0.643
## [3,]    12 10.68   0.556
## [4,]    16 15.59   0.482
## [5,]    20 18.44   0.558
## [6,]    24 20.50   0.668
## Lagrange-Multiplier test: 
##      order     LM  p.value
## [1,]     4 27.832 3.94e-06
## [2,]     8 12.088 9.77e-02
## [3,]    12  5.618 8.98e-01
## [4,]    16  3.265 9.99e-01
## [5,]    20  1.752 1.00e+00
## [6,]    24  0.933 1.00e+00

for (i in 1:15) {
  ArchTest <- ArchTest(model2d$residuals, lags=i, demean=TRUE)
  cat("P Value LM Test lag ke", i,"adalah" , ArchTest$p.value, "\n") }

## P Value LM Test lag ke 1 adalah 0.5545835 
## P Value LM Test lag ke 2 adalah 0.2900247 
## P Value LM Test lag ke 3 adalah 0.483532 
## P Value LM Test lag ke 4 adalah 0.6182282 
## P Value LM Test lag ke 5 adalah 0.7609083 
## P Value LM Test lag ke 6 adalah 0.7465258 
## P Value LM Test lag ke 7 adalah 0.845187 
## P Value LM Test lag ke 8 adalah 0.6064018 
## P Value LM Test lag ke 9 adalah 0.4083108 
## P Value LM Test lag ke 10 adalah 0.4976978 
## P Value LM Test lag ke 11 adalah 0.5497174 
## P Value LM Test lag ke 12 adalah 0.5666264 
## P Value LM Test lag ke 13 adalah 0.5911461 
## P Value LM Test lag ke 14 adalah 0.6126341 
## P Value LM Test lag ke 15 adalah 0.562215

Berdasarkan hasil uji LM diperoleh bahwa semua lag memiliki p-value yang lebih besar dari \(\alpha\) = 0.05 yang artinya Tak Tolak \(H_0\), dapat disimpulkan nahwa tidak terdapat heteroskedastisitas pada sisaan model ARIMA (0,1,1).

I. Forecasting

Tranformasi Balik

trainTrans.balik<-1/trainTrans

fitramal <- forecast::forecast(ts(trainTrans.balik),model=arima2,h=19)
meanramal <- fitramal$mean
plot(fitramal,ylab="Close PTBA Trans",xlab="Periode Waktu",main="Plot Peramalan ARIMA(0,1,1) (Data Training)")
lines((length(trainTrans.balik)+1):(length(trainTrans.balik)+19),test.ts,col="black")

Saran dan Simpulan

Simpulan

Pemulusan terbaik adalah dengan metode Double Exponential Smoothing.

Saran

Daftar Pustaka

Cahyakemala SR, Putri KPR. 2022. Konflik Ukraina Rusia memengaruhi harga saham pertambangan yang disebabkan perilaku investor. Sociability: Social & Humaniora Journal. 1(1): 1-9.

Iskandar S, Layakana M. 2020. Penerapan metode double moving average dan double eksponential smoothing dalam meramalkan jumlah produksi crude palm oil (CPO) pada PT. Perkebunan Nusantara IV unit dolok sinumbah. Karismatika. 6(1): 44-53.

Junaid MT, Juliana A, dan Sabrina H. 2020. Studi perbandingan model arima dan garch untuk memprediksi harga saham pada perusahaan tambang di Indonesia. Jurnal Ilmu Keuangan dan Perbankan (JIKA). 10(1): 83-98.

Maulina D, Hidayat T, Novianti B, Astuti Y. 2019. Penerapan metode single moving average untuk peramalan penjualan mainan anak. SENSITIf: Seminar Nasional Sistem Informasi dan Teknik Informatika. 3(1): 253-261.

Wahyuni T, Indahwati I, Sadik K. 2021. Perbandingan arima dan artificial neural networks dalam peramalan jumlah positif covid-19 Di DKI Jakarta. Xplore: Journal of Statistics. 10(3):289–302. doi:10.29244/xplore.v10i3.846.

Wulandari M. 2016. Analisis kualitas interpolasi terhadap fitur statistik pada citra. Ultimatics: Jurnal Teknik Informatika. 8(2): 139-146.