Pendahuluan

Dalam pasar saham yang dinamis, kemampuan untuk meramalkan harga saham menjadi sangat penting bagi para investor, analis, dan manajer keuangan. Dalam upaya untuk mencapai tujuan tersebut, metode peramalan time series telah menjadi alat yang populer dan efektif dalam mengidentifikasi tren dan pola dalam data historis, serta memberikan estimasi yang berguna tentang pergerakan harga saham di masa depan.

Dalam artikel ini, kami akan fokus pada peramalan harga closing saham PT Aneka Tambang Tbk (ANTAM) menggunakan model ARIMA (AutoRegressive Integrated Moving Average). ANTAM merupakan perusahaan pertambangan terkemuka di Indonesia yang terdaftar di Bursa Efek Indonesia (BEI). Harga closing saham ANTAM menggambarkan harga penutupan saham pada akhir setiap sesi perdagangan.

Pada artikel ini, kami akan menjelaskan langkah-langkah yang diambil dalam peramalan harga closing saham ANTAM menggunakan model ARIMA. Kami akan memulai dengan melakukan eksplorasi data dan analisis visual untuk memahami karakteristik time series dan mengidentifikasi tren serta pola yang mungkin ada. Selanjutnya, kami akan memastikan bahwa data time series memenuhi asumsi stasioneritas, dan jika tidak memenuhi, kami akan menerapkan transformasi atau differencing untuk mencapainya. Kami akan membangun model ARIMA dengan memilih parameter yang optimal menggunakan teknik seperti identifikasi model, estimasi parameter, dan pengujian model. Setelah model ARIMA yang tepat ditemukan, kami akan menggunakan model tersebut untuk meramalkan harga closing saham ANTAM untuk periode masa depan.

Tujuan dari artikel ini adalah untuk memberikan pemahaman yang mendalam tentang penerapan model ARIMA dalam peramalan harga closing saham ANTAM. Diharapkan bahwa hasil dari penelitian ini akan memberikan wawasan yang berharga bagi para investor, analis pasar, dan praktisi keuangan dalam pengambilan keputusan investasi yang lebih informasional dan akurat.

Namun, perlu diingat bahwa peramalan harga saham melibatkan ketidakpastian dan risiko inheren, dan hasil peramalan tidak dapat dijamin secara mutlak. Oleh karena itu, sangat penting untuk melihat hasil peramalan sebagai alat bantu pengambilan keputusan, dan tetap mempertimbangkan faktor lain seperti analisis fundamental, kondisi pasar saat ini, dan risiko yang terkait dengan investasi saham.

Data Preparation

Library

library(quantmod) # quantitive Financial Modelling (mendapatkan data saham)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lubridate) # date manipulation
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(plotly) # interactive Plotting
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(dplyr) # data wrangling library
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(forecast) # time series library
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
library(tseries) # adf.test
library(fpp) # usconsumtion
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
library(TSstudio) # timeseries Visualization
library(ggplot2) # data visualization
library(tidyr) # transforming data
library(padr) # complete data frame
library(zoo) # missing value imputation
library(data.table) # data manipulation & analysis
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:xts':
## 
##     first, last
library(glue) # for creating formatted strings into a template string
library(scales) # for formatting numeric and date/time values

Import Dataset

df <- getSymbols("ANTM.JK", auto.assign = FALSE, from=as.Date("2017-01-01"), to = as.Date("2022-12-31"))
## Warning: ANTM.JK contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.
df <- as.data.frame(df)

df <- setDT(df, keep.rownames = TRUE)[]

Data Cleansing

Setelah dataset diunduh, dilakukan pengecekkan Nullitas pada data untuk melihat apakah terdapat data yang kosong.

anyNA(df)
## [1] TRUE
sum(is.na(df))
## [1] 6

Ternyata terdapat 6 data yang kosong atau tidak memilki nilai, selanjutnya dilakukan pengecekkan untuk melihat jumlah baris yang kosong.

sum(!complete.cases(df))
## [1] 1

Ternyata hanya memiliki satu baris yang kosong, hal ini tidak berpengaruh terhadap data. Sehingga dilakukan penghapusan data yang kosong tersebut.

df <- na.omit(df)  # Remove rows with NA values

Data Pre-Processing

Selanjutnya, dilakukan manipulasi tipe data pada kolom tanggal yang semula bertipe data char menjadi tipe data date dan juga membuat kolom baru yaitu tahun, bulan, dan hari.

glimpse(df)
## Rows: 1,507
## Columns: 7
## $ rn               <chr> "2017-01-02", "2017-01-03", "2017-01-04", "2017-01-05…
## $ ANTM.JK.Open     <dbl> 895, 895, 885, 875, 880, 870, 855, 875, 870, 885, 920…
## $ ANTM.JK.High     <dbl> 895, 895, 905, 885, 885, 875, 875, 880, 880, 920, 930…
## $ ANTM.JK.Low      <dbl> 895, 875, 865, 865, 865, 850, 855, 865, 865, 885, 900…
## $ ANTM.JK.Close    <dbl> 895, 885, 865, 870, 865, 855, 865, 875, 865, 920, 905…
## $ ANTM.JK.Volume   <dbl> 0, 3577900, 34615900, 22442900, 31596400, 28355600, 2…
## $ ANTM.JK.Adjusted <dbl> 820.9537, 811.7811, 793.4357, 798.0221, 793.4357, 784…
df <- df %>% 
  mutate(rn = ymd(rn),
         Date = ymd(rn),
         Year = as.factor(year(Date)),
         Month = as.factor(month(Date)),
         WeekDay = lubridate::wday(Date,label = TRUE))
head(df)
##            rn ANTM.JK.Open ANTM.JK.High ANTM.JK.Low ANTM.JK.Close
## 1: 2017-01-02          895          895         895           895
## 2: 2017-01-03          895          895         875           885
## 3: 2017-01-04          885          905         865           865
## 4: 2017-01-05          875          885         865           870
## 5: 2017-01-06          880          885         865           865
## 6: 2017-01-09          870          875         850           855
##    ANTM.JK.Volume ANTM.JK.Adjusted       Date Year Month WeekDay
## 1:              0         820.9537 2017-01-02 2017     1     Mon
## 2:        3577900         811.7811 2017-01-03 2017     1     Tue
## 3:       34615900         793.4357 2017-01-04 2017     1     Wed
## 4:       22442900         798.0221 2017-01-05 2017     1     Thu
## 5:       31596400         793.4357 2017-01-06 2017     1     Fri
## 6:       28355600         784.2631 2017-01-09 2017     1     Mon

Exploratory Data Analysis

Melihat pergerakan Harga saham

Untuk melihat pergerakan harga penjualan saham Bukalapak cenderung naik atau turun, maka akan disajikan line chart dari Opening Price saham tahun 2021 - 2022

 OpenPlot <- ggplot(df, aes(x = rn, y = ANTM.JK.Open)) +
      geom_line(col="red")+
        labs(title = "Opening Price Saham ANTAM")+
        theme_light()

ggplotly(OpenPlot)

Dan berikut line chart untuk Close Price

 OpenPlot <- ggplot(df, aes(x = rn, y = ANTM.JK.Close)) +
      geom_line(col="blue")+
        labs(title = "Closing Price Saham ANTM")+
        theme_light()

ggplotly(OpenPlot)
 OpenPlot <- ggplot(df, aes(x = rn, y = ANTM.JK.Volume)) +
      geom_line(col="green")+
        labs(title = "Volume Perdagangan Saham ANTM")+
        theme_light()

ggplotly(OpenPlot)

Dari visualisasi line chart data Open Price, Close Price, dan Volume Price diatas, terlihat bahwa harga saham ANTM cenderung fluktuatif dan mengalami kenaikan . Dan dapat dilihat bahwa dari tahun 2020-2021 saham ANTM mengalami peningkatan harga dan volume yang signifikan kemudian sempat turun sedikit hingga kembali mengalami fluktuatif sampai tahun 2022

Rata-rata nilai harga saham perhari

average_open <- df %>%
  group_by(WeekDay) %>% 
  summarise(Average.Perday = mean(ANTM.JK.Open)) %>% 
  ungroup()

average_open
## # A tibble: 5 × 2
##   WeekDay Average.Perday
##   <ord>            <dbl>
## 1 Mon              1299.
## 2 Tue              1287.
## 3 Wed              1277.
## 4 Thu              1290.
## 5 Fri              1281.
average_open <- average_open %>%
  mutate(label = glue("Day: {WeekDay}\n Average Opening Price: {format(Average.Perday, big.mark = ',')}"))

opnPlot <- ggplot(average_open, aes(x = WeekDay, y = Average.Perday, text = label)) +
  geom_col(aes(fill = Average.Perday)) +
  labs(title = "Average Opening Price",
       subtitle = "contoh",
       x = "Day",
       y = NULL) +
  scale_fill_gradient(high= "red",low = "orange") +
  theme_minimal() +
  theme(legend.position = "none") 

ggplotly(opnPlot, tooltip = "text")
average_close <- df %>%
  group_by(WeekDay) %>% 
  summarise(Average.Perday = mean(ANTM.JK.Close)) %>% 
  ungroup()

average_close
## # A tibble: 5 × 2
##   WeekDay Average.Perday
##   <ord>            <dbl>
## 1 Mon              1296.
## 2 Tue              1280.
## 3 Wed              1279.
## 4 Thu              1287.
## 5 Fri              1279.
average_close <- average_close %>%
  mutate(label = glue("Day: {WeekDay}\n Average Closing Price: {format(Average.Perday, big.mark = ',')}"))

clsPlot <- ggplot(average_close, aes(x = WeekDay, y = Average.Perday, text = label)) +
  geom_col(aes(fill = Average.Perday)) +
  labs(title = "Average Closing Price",
       subtitle = "contoh",
       x = "Day",
       y = NULL) +
  scale_fill_gradient(high= "red",low = "orange") +
  theme_minimal() +
  theme(legend.position = "none") 

ggplotly(clsPlot, tooltip = "text")

Berdasarkan dari visualisasi diatas hari senin menempati posisi teratas dengan rata-rata Opening price dan Closing price terbesar dibandingkan hari yang lain, namun secara umum rata-rata opening price dan closing price perharinya tidak memiliki perbedaan yang signifikan.

Data Processing

Pada kasus kali ini saya akan melakukan forecasting pada pada Closing price, hal ini dikarenakan Closing price merupakan hal yang paling penting untuk menimbang untuk menanam saham bagi investor. Alasan lain karena harga closing price merupakan acuan harga opening price, pada hari kerja selanjutnya.

close <- df %>% 
  select(c(rn,ANTM.JK.Close))

Padding Data

Pada Data time-series dilakukan padding untuk memastikan interval dari sequences data yang digunakan sama

# padding data
c_pad <- pad(
  close,
  interval = NULL,
  start_val = ymd("2017-01-01"),
  end_val = ymd("2022-12-31"),
  by = NULL,
  group = NULL,
  break_above = 1
)
## pad applied on the interval: day
is.null(df)
## [1] FALSE
# mengisi data NA
antm_pad<- c_pad %>% 
  mutate(ANTM.JK.Close = na.fill(ANTM.JK.Close, fill = "extend"))
head(antm_pad)
##            rn ANTM.JK.Close
## 1: 2017-01-01           895
## 2: 2017-01-02           895
## 3: 2017-01-03           885
## 4: 2017-01-04           865
## 5: 2017-01-05           870
## 6: 2017-01-06           865

Decomposition

Sebelum melakukan pemodelan pada permalan forecasting kita perlu mengamati objek timeseries dari hasil decompose. Alasan utama dari decompose adalah untuk menguraikan ketiga komponen dari objek ts (trend, seasonal, residual).

antm_ts <- ts(data = antm_pad$ANTM.JK.Close,
             start = range(antm_pad$rn)[[1]],
             frequency = 365) # yearly seasonality

# melihat plot time series
antm_ts %>%
  decompose()%>%
  autoplot()

Berdasarkan plot dekomposisi dapat diketahui pola data yang digunakan bersifat additive.

Feature Engineering

Splitting Data

Selanjutnya dilakukan pembagian antara data train dan data test.Data train yang digunakan adalah data perhari selama 4 tahun, yaitu data dari tahun 2017-2021. Sedangkan data testing yang digunakan adalah data 365 hari terakhir.

test <- tail(antm_ts, 365)
train <- head(antm_ts, length(antm_ts)-length(test))

Stasioneritas Data

Apabila suatu data deret waktu memiliki kecenderungan untuk bergerak diantara mean-nya1. Sebelum dilakukan pemodelan menggunakan Arima perlu dilakukan pengujian untuk melihat kestasioneran pada data.

adf.test(train)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train
## Dickey-Fuller = -1.9906, Lag order = 12, p-value = 0.5823
## alternative hypothesis: stationary

Adapun Hipotesis pada data: H0 : Data tidak stasioner H1 : Data stasioner

Daerah kritis

p-value > \(\alpha\)

📌 Keputusan: Hasil menunjukkan bahwasannya p-value > \(\alpha\)(0,05) sehingga dapat dikatakan data belom stasioner

Differencing/Integreted (I)

Setelah diketahui bahwa data belum stasioner, dilakukan differencing untuk membuat data menjadi stasioner dengan cara mengurangi nilai data saat ini dengan nilai data sebelumnya2

I(d): \(Y't = Y_t - Y{t-1}\)

dimana, - \(Y_t\) adalah nilai observasi time series pada waktu ke-t. - \(Y'_{t}\) adalah data yang telah mengalami differencing sebanyak d kali.

# differencing 1 kali
antm_ts1 <- train %>% 
  diff()

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

Adapun Hipotesisnya:

H0 : Data tidak stasioner H1 : Data stasioner

Daerah kritis p-value < \(\alpha\)

📌 Keputusan: Hasil menunjukkan bahwasannya p-value < \(\alpha\) (0,05) sehingga dapat dikatakan data telah stasioner

Pemodelan

Identifikasi Model

Sebelum Melakukan pemodelan, perlu dilakukan Identifikasi model AR(p) dan MA(q) yang akan digunakan dengan melihat plot Partial Autocorrelation Function (PACF) dan auto correlation Function (ACF)

Auto Regressive (AR) Model

Model linear regression berdasarkan lag dari datanya sebagai prediktor3

\[Y_t = \alpha + \beta_1Y_{t-1}+\beta_2Y_{t-2}+...+\beta_pY_{t-p}+ \epsilon_1\]

Dimana,

  • \(Y_{t-1}\) adalah lag1.
  • \(\beta_1\) adalah coefficient dari lag1.
  • \(\alpha\) adalah intercept.

MA(q) : Moving Average

Rata-rata berjalan terhadap data time series itu sendiri.4

\[ Y_t = \alpha + \epsilon_t + \phi_1\epsilon_{t-1}+\phi_2\epsilon_{t-2}+...+\phi_q\epsilon_{t-q} \] Dimana,

  • \(Y_t\) adalah nilai observasi time series pada waktu ke-t.
  • \(\alpha\) adalah konstanta atau intercept.
  • \(\epsilon_t\) adalah error atau residual pada waktu ke-t.
  • \(\phi_1, \phi_2, ..., \phi_q\) adalah koefisien moving average yang menggambarkan pengaruh error pada nilai saat ini dan sebelumnya.
antm_ts1 %>% 
  tsdisplay(lag.max = 20)

Berdasarkan plot ACF (cut off pada lag ke-1 dan ke-2), PACF (cut off pada lag ke-1 dan ke-2), dan ACF, diperoleh kandidat model:

  • P = 0,1,2
  • Q = 0,1,2

Sehingga, dapat dikatakan model yang dapat dibangun adalah5:

  • ARIMA (0,1,1)
  • ARIMA (1,1,0)
  • ARIMA (1,1,1)
  • ARIMA (0,1,2)
  • ARIMA (1,1,2)
  • ARIMA (2,1,0)
  • ARIMA (2,1,2)

Selain menentukan ordo ARIMA secara manual, R menyediakan fitur ‘auto.arima’ yang dapat digunakan untuk menentukan ordo model secara otomatis yang akan kita gunakan dalam pemodelan kali ini.

Pemodelan ARIMA

ARIMA(p, d, q): \[Y't = c + \phi_1Y'{t-1} + \phi_2Y'{t-2} + ... + \phi_pY'{t-p} + \theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2} + ... + \theta_q\varepsilon_{t-q} + \varepsilon_t\]

Dimana,

  • \(Y_t\) adalah nilai observasi time series pada waktu ke-t.
  • \(Y'_t\) adalah data yang telah mengalami differencing (penyesuaian integrasi) sebanyak d kali.
  • \(c\) adalah konstanta atau intercept.
  • \(\phi_1, \phi_2, ..., \phi_p\) adalah koefisien autoregressive yang menggambarkan pengaruh nilai-nilai sebelumnya pada nilai saat ini.
  • \(\theta_1, \theta_2, ..., \theta_q\) adalah koefisien moving average yang menggambarkan pengaruh error pada nilai saat ini dan sebelumnya.
  • \(\varepsilon_t\) adalah error atau residual pada waktu ke-t.

Pemodelan dengan menggunakan STLM (Seasonal Trend with Loess Model)

Pada metode decompose biasa, untuk mendapatkan komponen trend dengan cara central moving average (CMA) diberikan bobot yang sama sesuai ordo yang ditetapkan. Dikarenakan dengan merata-ratakan data tengahnya hasilnya data awal dan akhir akan hilang, sehingga terdapat beberapa informasi yang hilang.

Untuk itu, terdapat salah satu cara untuk mendapatkan decompose data namun tetap mempertahankan informasi dari seluruh data yaitu dengan menggunakan STL(Seasonal Trend with Loess)6.

Secara Konsep, STL akan melakukan metode smoothing terhadap data tetangga dari setiap observasi dengan memberikan bobot yang lebih berat terhadap data yang dekat dengan data yang diobservasi.

Dalam memodelkan Hasil STL, digunakan STLM (Seasonal Trend With Loess Model) dalam menerapkan metode ARIMA. STLM digunakan dalam menangkap seasonal yang belum bisa ditangkap oleh metode ARIMA biasa.

# ARIMA
close_arima <- stlm(train, method = "arima", lambda = 0, s.window = 365)

Parameter stlm() yang harus di atur:

  • y : object time series
  • s.window : pola seasonal yang ingin ditangkap (dalam metode kali ini kita menggunakan 365 hari)
  • method : metode forecast yang akan digunakan, tersedia ets dan arima

Peramalan (Forecasting)

Dalam Peramalan kali ini dilakukan peramalan untuk 365 hari mendatang

# forecasting model arima

forecast_arima <- forecast(close_arima, h = 365)

Evaluasi Model

# result from stlm method
MAPE(forecast_arima$mean, y_true = test)*100
## [1] 15.95365
accuracy(forecast_arima$mean, test)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -84.67965 397.5379 332.5358 -6.328607 15.95365 0.9903767  8.524294

Dalam mengevaluasi model didasarkan pada akurasi. Akurasi yang digunakan adalah MAPE (Mean Absolut Percentage Error) dan RMSE (Root Mean Squared Error). RMSE dan MAPE sama-sama bagus. Alasan digunakan RMSE adalah karena RMSE lebih sensitif pada data pencilan (amatan yang nilainya jauh dari rata-rata) Sementara itu, MAPE memberikan nilai presentase sehingga lebih mudah untuk diinterpretasikan7

Dalam Hal ini model yang dibangun dapat dikatakan cukup baik8 dengan nilai MAPE berada di angka 15.953

Plot Peramalan

arima_p <- autoplot(forecast_arima, series = "ARIMA", fcol = "red") +
  autolayer(antm_ts, series = "Actual", color = "black") + 
  labs(subtitle = "ANTM Stock Price",
       y = "Closing Price") +
  theme_minimal()

arima_p

Dapat dilihat berdasarkan data diatas bahwasannya model yang dibangun dapat dikatakan cukup baik dengan peramalan/prediksi digambarkan dengan garis dan plot penyebaran berwarna merah.

Uji Asumsi

shapiro.test(forecast_arima$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  forecast_arima$residuals
## W = 0.9255, p-value < 2.2e-16

Hipotesis H0 : Data menyebar pada distribusi normal H1 : Data tidak menyebar pada distribusi normal

Daerah kritis p-value < \(\alpha\)

📌 Keputusan: Data tidak menyebar pada distribusi normal

plot(forecast_arima$residuals)

hist(forecast_arima$residuals, breaks = 30)

Box.test(forecast_arima$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  forecast_arima$residuals
## X-squared = 0.00055183, df = 1, p-value = 0.9813

Hipotesis H0 : Tidak terdapat residual autokorelasi H1 : Terdapat residual autokorelasi

Daerah kritis p-value > \(\alpha\)

📌 Keputusan: Data tidak terdapat autokorelasi

Kesimpulan dan Saran

Berdasarkan pengujian asumsi, tidak terdapat autokorelasi pada residu peramalan (p-value > 0,05). Namun, residual peramalan dapat dikatakan tidak terdistribusi secara normal, oleh karena itu residual mungkin tidak muncul di sekitar rata-ratanya seperti yang terlihat pada histogram. Tapi, jika dilakukan pemeriksaan terhadap distribusi residu melalui plot garis, itu sebenarnya mirip dengan plot kesalahan dari dekomposisi objek data deret waktu. Sehingga dapat dikatakan Model Arima ini sudah memprediksi dengan cukup baik.

Model yang dibangun bersifat Analisis Univariat sehingga dapat dikatakan perlu adanya pertimbangan mendalam terhadap peubah-peubah lain dalam melakukan pertimbangan untuk pengambilan keputusan terhadap pembelian/penjualan saham. Dikarenakan fluktuasi harga saham sendiri berasal dari banyak faktor yang terdiri oleh peubah-peubah lain yang mempengaruhi pergerakkan harga saham.

Referensi