Load library yang akan digunakan.
library(tidyverse) #data manipulation
library(lubridate) # date manipulation
library(forecast) # time series library
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
library(tseries) # adf.test
library(fpp)
library(padr)
library(zoo)
library(dplyr)
library(ggplot2)
library(gridExtra)
library(imputeTS)
library(quantmod)
library(TSA)
library(nortest)
library(prettydoc)Data Proyek
IHSG adalah indeks harga saham gabungan secara umum dan berfungsi sebagai pengukuran kinerja suatu saham gabungan di Bursa Efek Indonesia (BEI). IHSG dipengaruhi oleh faktor eksternal yaitu tingkat suku bunga luar negeri dan indeks harga saham luar negeri serta faktor internal yaitu nilai tukar rupiah terhadap USD, inflasi, suku bunga deposito, suku bunga Sertifikat Bank Indonesia, dan jumlah uang yang ada di masyarakat (Purba 2014). Data yang digunakan adalah data IHSG sektoral keuangan (JKFINA) periode Januari 2012 sampai dengan Maret 2022 yang diperoleh dari www.yahoo.finance.com.
Load data yang akan digunakan
start <- as.POSIXct("2012-01-01")
end <- as.POSIXct("2022-03-31")
getSymbols(Symbols = "^JKFINA",src = "yahoo", from = start, to = end)[1] “^JKFINA”
class(JKFINA)[1] “xts” “zoo”
head(JKFINA) JKFINA.Open JKFINA.High JKFINA.Low JKFINA.Close JKFINA.Volume2012-01-02 492.05 493.50 488.52 491.10 84100000 2012-01-03 491.82 496.68 491.77 496.62 269100000 2012-01-04 495.90 501.61 495.85 499.21 254800000 2012-01-05 498.94 499.97 494.75 496.29 545200000 2012-01-06 496.01 496.85 490.24 492.98 204600000 2012-01-09 492.26 498.93 488.91 498.93 539600000 JKFINA.Adjusted 2012-01-02 491.10 2012-01-03 496.62 2012-01-04 499.21 2012-01-05 496.29 2012-01-06 492.98 2012-01-09 498.93
colnames(JKFINA)[1] “JKFINA.Open” “JKFINA.High” “JKFINA.Low” “JKFINA.Close”
[5] “JKFINA.Volume” “JKFINA.Adjusted”
df <- data.frame(JKFINA)
df <- cbind(Date=rownames(df),df)
rownames(df) <- 1:nrow(df)Tabel 1 memperlihatkan seluruh variabel yang akan digunakan. Data yang akan diguakan hanya harga saat closing, sehingga hanya kolom tanggal dan kolom harga closing yang akan diambil. Tabel 2 memperlihatkan sebagian dari data yang akan digunakan.
df1 <- df %>%
select(c(Date,JKFINA.Close))
df1$Date <- ymd(df1$Date)
head(df1) %>% knitr::kable(
align = "c",
caption = "Lima data pertama *subset* data yang digunakan"
)| Date | JKFINA.Close |
|---|---|
| 2012-01-02 | 491.10 |
| 2012-01-03 | 496.62 |
| 2012-01-04 | 499.21 |
| 2012-01-05 | 496.29 |
| 2012-01-06 | 492.98 |
| 2012-01-09 | 498.93 |
Metode ARIMA
Salah satu metode yang banyak digunakan dalam peramalan deret waktu adalah metode ARIMA. Metode ARIMA merupakan gabungan dari proses AR (Autoregressive) dan MA (Moving Average) (Rakhmawati et al. 2020). Metode ini menghiraukan variabel bebas karena metode ini hanya bergantung pada data historis (Fauziyah et al. 2021). Salah satu keunggulan dari metode ARIMA adalah bahwa hanya data yang akan diprediksi yang dibutuhkan sehingga akan memudahkan dalam memprediksi data dalam jangka panjang (Meyler et al. 1998).
cariin rumus sama parameternya apa aja dong. Sama kalo ada yang punya ide lagi boleh banget nambahin di sini.
Eksplorasi Data
Pada data tersebut, terlihat bahwa terdapat beberapa tanggal yang terlewati, seperti pada tanggal 1 Januari 2012, 7 Januari 2012, dan 8 Januari 2012 dan masih banyak lagi. Setelah dilakukan investigasi, ternyata tanggal-tanggal ini terletak pada hari Sabtu dan Minggu, dan hari libur nasional, sehingga tanggal-tanggal yang dilewati ini harus dimasukkan dalam dataset tersebut seperti pada Tabel 3.
date.seq <- data.frame("Date" = seq(ymd(start), ymd(end), by = "days"))
df1 <- merge(date.seq, df1, by = "Date", all.x = T)
knitr::kable(head(df1),
caption = "Data dengan tanggal yang benar", align = "c",
label = "data_tanggal_bener")| Date | JKFINA.Close |
|---|---|
| 2012-01-01 | NA |
| 2012-01-02 | 491.10 |
| 2012-01-03 | 496.62 |
| 2012-01-04 | 499.21 |
| 2012-01-05 | 496.29 |
| 2012-01-06 | 492.98 |
Pada analisis data deret waktu, tidak boleh ada data yang kosong dan periode deret waktu harus seragam untuk seluruh data. Untuk itu, pengecekan data kosong atau missing-values harus dilakukan terlebih dahulu.
sum(is.na(df1))## [1] 1248Terdapat 1248 data kosong yang berada pada dataset tersebut, sehingga harus dilakukan interpolasi terlebih dahulu. Interpolasi yang digunakan adalah spline interpolation dimana data dibagi menjadi beberapa bagian lalu data kosong disubstitusikan berdasarkan data pada bagian tersebut. Metode ini relatif mudah digunakan dan pola dan karakteristik data tetap dipertahankan (Rappai 2015).
Menurut Wang (2013), model matematika untuk spline interpolation dapat digambarkan sebagai berikut:
Untuk i = 1,2,…,n titik data, interpolasi untuk semua pasangan knots \((x_{i-1},y_{i-1})\) dan \((x_i,y_i)\) dengan polinomial :
\[y=q_i(x) ;i = 1,2,…,n\]
kelengkungan kurva fungsi \(y=f(x)\) adalah
\[k = \frac{y"}{(1+y'^2)^\frac{2}{3}}\]
Data yang sudah dilakukan imputasi dapat dilihat pada Tabel 4.
dataimp = na_interpolation(df1, option = "spline")
head(dataimp) %>% knitr::kable(caption = "Data yang sudah diimputasi")| Date | JKFINA.Close |
|---|---|
| 2012-01-01 | 491.10 |
| 2012-01-02 | 491.10 |
| 2012-01-03 | 496.62 |
| 2012-01-04 | 499.21 |
| 2012-01-05 | 496.29 |
| 2012-01-06 | 492.98 |
sum(is.na(dataimp))## [1] 0Pada dataset tersebut sudah tidak ada lagi missing values sehingga analisis deret waktu dapat dilanjutkan. Setelah itu, data akan diagregasi menjadi data bulanan.
Plot data sebelum interpolasi dan setelah interpolasi dapat dilihat pada Gambar 1.
ggplot_na_imputations(df1$JKFINA.Close, dataimp$JKFINA.Close,
title = NULL,
subtitle = NULL,
size_points = 1.5,
size_imputations = 0.75,
xlab = "\nTahun",
ylab = "Harga saat closing\n",
x_axis_labels = dataimp$Date,
label_known = "Sebelum imputasi",
label_imputations = "Setelah imputasi",
shape_points = 20,
shape_imputations = 20,
legend_size = 5,
theme = theme_classic())
Plot data interpolasi dan data asli.
Terlihat bahwa imputasi data dengan metode spline interpolation dapat memberikan nilai-nilai yang memiliki pola yang sama dengan data aslinya. Setelah melakukan imputasi, maka data akan diagregasikan menjadi data rataan per harga closing IHSG per bulannya. Data setelah diagregasi dapat dilihat pada Tabel 5.
dataimp$Date <- floor_date(ymd(dataimp$Date),"month")
df2 <- dataimp %>%
group_by(Date) %>%
summarize(Close = mean(JKFINA.Close))
head(df2) %>% knitr::kable(caption = "Data agregasi per bulan.")| Date | Close |
|---|---|
| 2012-01-01 | 499.9487 |
| 2012-02-01 | 487.7797 |
| 2012-03-01 | 496.7415 |
| 2012-04-01 | 522.1882 |
| 2012-05-01 | 501.7686 |
| 2012-06-01 | 484.4287 |
Statistika deskriptif singkat dari data dapat dilihat pada Tabel 6.
summary(df2)## Date Close
## Min. :2012-01-01 Min. : 484.4
## 1st Qu.:2014-07-16 1st Qu.: 662.9
## Median :2017-02-01 Median : 824.5
## Mean :2017-01-30 Mean : 926.3
## 3rd Qu.:2019-08-16 3rd Qu.:1203.4
## Max. :2022-03-01 Max. :1615.8Didapatkan nilai rata-rata sebesar 926.3 dan median sebesar 824.5. Untuk nilai minimum data tersebut yaitu sebesar 484.429 terjadi pada bulan Juni tahun 2012 sedangkan untuk nilai maksimum data tersebut sebesar 1615.753 terjadi pada bulan Maret tahun 2022.
Pada bulan Desember 2019 IHSG bergerak menguat hingga 4,79% hal tersebut sudah biasa terjadi karena sejak tahun 2011 sudah 7 kali IHSG selalu meningkat pada penutupan perdagangan (akhir tahun). Selain peristiwa akhir tahun, peristiwa Idul Adha di Indonesia pun bisa memberikan efek positif terhadap IHSG (Nugroho. 2017).
Salah satu metode eksplorasi yang paling cepat untuk dilakukan adalah dengan melakukan plot time-series. Plot ini memperlihatkan secara cepat bagaimana pola dan trend dari data itu sendiri.
ggplot(df2, aes(x = Date, y = Close)) +
geom_line(width = 1.5, col = "steelblue") +
geom_point(size = 1.5, col = "steelblue") +
scale_x_date(date_breaks = "year", date_labels = "%Y") +
scale_y_continuous(breaks = seq(200,
1500,
by = 200)
) +
xlab("\nPeriode Waktu") +
ylab("Harga saat closing\n") +
theme_classic() 
Plot time series data IHSG.
Pada Grafik 2, seiring dengan berjalannya waktu, harga closing IHSG semakin meningkat. Tentunya hal ini bukan tanpa sebab, karena IHSG sangat terikat dengan fenomena ekonomi makro, baik internal, seperti inflasi, suku bunga, nilai tukar, dan lainnya, serta eksternal, seperti indeks bursa harga luar negeri, harga minyak mentah, harga emas, dan kestabilan politik suatu negara (Wijaya 2015). Penelitian oleh Tambunan dan Aminda (2021) menyatakan bahwa trend naiknya IHSG dapat dijelaskan secara kuat oleh kurs rupiah. Dari plot tersebut tidak terlihat adanya plot musiman. Hal ini ditandai dengan tidak adanya pola siklik pada data tersebut.
Dalam pemodelan dengan ARIMA, perlu diketahui apakah data tersebut stasioner atau tidak. Stasionaritas data dibutuhkan dalam forecasting dengan metode ARIMA, sehingga data yang belum stasioner harus ditransformasikan terlebih dahulu menggunakan proses pembedaan atau differencing (Fattah et al. 2018). Dengan melihat plot time-series data seperti tersebut terlihat bahwa data memiliki trend, sehingga data tersebut tidak stasioner. Ketidakstasioneran data dapat juga dilihat dari plot ACF-nya seperti di bawah ini.
acf(df2$Close, lag.max = 100, plot = F) %>%
autoplot() +
ggtitle("Plot ACF Harga Closing IHSG") +
theme_classic()
Plot ACF memperlihatkan bahwa autokorelasi antar lag menurun secara perlahan, sehingga membuktikan bahwa data tersebut tidak stasioner. Uji formal Augmented Dicky-Fuller dapat membantu untuk memastikan hal tersebut. Mohamed (2008) menyebutkan bahwa terdapat dua hipotesis yang digunakan dalam uji ini yaitu:
- \(H_0\): Data tidak stasioner
- \(H_1\): Data stasioner.
Hasil perhitungan ADF dapat dilihat sebagai berikut.
adf.test(df2$Close)##
## Augmented Dickey-Fuller Test
##
## data: df2$Close
## Dickey-Fuller = -2.6089, Lag order = 4, p-value = 0.3236
## alternative hypothesis: stationaryKarena nilai p-value > 0.05, maka belum cukup bukti untuk menolah \(H_0\), sehingga dapat dikatakan bahwa data memang tidak stasioner.
Seperti yang telah dijelaskan sebelumnya, data non-stasioner dapat distasionerkan dengan melakukan proses differencing. Hasil differencing pertama data harga closing dapat dilihat di bawah ini.
df.diff.close <- cbind.data.frame(date = df2$Date[2:nrow(df2)], diff = diff(df2$Close, 1))
head(df.diff.close) %>% knitr::kable()| date | diff |
|---|---|
| 2012-02-01 | -12.169024 |
| 2012-03-01 | 8.961815 |
| 2012-04-01 | 25.446714 |
| 2012-05-01 | -20.419620 |
| 2012-06-01 | -17.339884 |
| 2012-07-01 | 23.997179 |
Kestasioneran dapat dilihat seperti sebelumnya, yaitu dengan menggunakan plot time-series, plot ACF, dan uji formal ADF.
ts.plot.diff <-
df.diff.close %>%
ggplot(aes(x = date, y = diff))+
geom_line(color = "steelblue") +
geom_point(color = "steelblue") +
geom_hline(yintercept = mean(df.diff.close$diff),
color = "red", linetype = "dashed") +
xlab("Tanggal") +
ylab("Differencing pertama closing IHSG") +
theme_classic()
acf.diff <-
acf(df.diff.close$diff, lag.max = 100, plot = F) %>%
autoplot() +
ggtitle(NULL) +
theme_classic()
grid.arrange(ts.plot.diff, acf.diff, ncol = 2)
adf.test(df.diff.close$diff)##
## Augmented Dickey-Fuller Test
##
## data: df.diff.close$diff
## Dickey-Fuller = -5.1687, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryPada plot time-series, terlihat bahwa data yang sudah dilakukan pembedaan sudah menyebar di sekitar suatu konstanta, sehingga dapat dikatakan bahwa data sudah terlihat stasioner. Akan tetapi, dari plot tersebut masih diragukan apabila ragam pada data tersebut sudah homogen. Pada plot ACF juga masih meragukan, karena terlihat pola siklus pada data tersebut. Akan tetapi hasil pengujian ADF memperlihatkan bahwa sudah stasioner.
Pemodelan data IHSG dengan Metode ARIMA
Penentuan parameter model ARIMA dapat dilihat dari beberapa sudut pandang. Untuk parameter \(p\) dan \(q\) dapat dilihat dari plot ACF, plot PACF, dan plot EACF. Untuk parameter \(q\) dapat dilihat dari apakah data tersebut stasioner atau tidak. Apabila dilakukan differencing ordo 1, maka nilai \(q = 1\), dan seterusnya. Sebelum memulai pemodelan tersebut, data dibagi menjadi data training dan data testing. Data training dan testing dibagi beradasarkan informasi trend dan pola data karena diusahakan pola data pada data training dan data testing sama. Partisi yang akan digunakan sebagai berikut.
ggplot(df2, aes(x = Date, y = Close)) +
geom_line(width = 1.5, col = "steelblue") +
geom_point(size = 1.5, col = "steelblue") +
scale_x_date(date_breaks = "year", date_labels = "%Y") +
scale_y_continuous(breaks = seq(200,
1500,
by = 200)
) +
geom_vline(xintercept = as.numeric(ymd("2019-02-01")), linetype = "dashed", color = "red") +
xlab("\nPeriode Waktu") +
ylab("Harga saat closing\n") +
theme_classic()
Partisi yang ditentukan untuk data training memiliki rentang dari
Januari 2012 sampai dengan Desember 2018. Sedangkan data testing berawal
dari Januari 2019 sampai dengan Maret 2022. Indeks data Desember 2018
dapat dilihat sebagai berikut.
nov_2019 <- which(df2$Date == "2019-11-01")
train_ts <- ts(df2$Close[1:nov_2019])
test_ts <- ts(df2$Close[(nov_2019+1):nrow(df2)])
head(train_ts) %>% knitr::kable(caption = "Lima data pertama training")| x |
|---|
| 499.9487 |
| 487.7797 |
| 496.7415 |
| 522.1882 |
| 501.7686 |
| 484.4287 |
head(test_ts) %>% knitr::kable(caption = "Lima data pertama testing")| x |
|---|
| 1327.0432 |
| 1361.0328 |
| 1326.9846 |
| 1067.3271 |
| 963.6142 |
| 914.9505 |
Pemodelan akan didasarkan pada data training, kemudian validasi model akan didasarkan pada testing. Karena data dalam keadaan tidak stasioner, maka data training dan testing akan dikenakan proses differencing. Kemudian, plot ACF, PACF, dan EACF dari data training dapat dilihat sebagai berikut.
train_ts_diff <- diff(train_ts, 1)
test_ts_diff <- diff(test_ts, 1)
acf.diff <-
acf(train_ts_diff, lag.max = 20, plot = F) %>%
autoplot() +
ggtitle(NULL) +
theme_classic()
pacf.diff <-
pacf(train_ts_diff, lag.max = 20, plot = F) %>%
autoplot() +
ggtitle(NULL) +
theme_classic()
grid.arrange(acf.diff, pacf.diff, ncol = 2)
eacf(train_ts_diff)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 x o o o o o o o o o o o o o
## 3 o o x o o o o o o o o o o o
## 4 x x x x o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 o x o o x o o o o o o o o o
## 7 o x o o x o o o o o o o o oPada plot ACF dan PACF, terlihat bahwa terdapat pola siklus sehingga susah untuk menentukan model tentatif yang dapat digunakan. Dari plot EACF terlihat bahwa model tentatif yang dapat terbentuk adalah model ARIMA(0,1,1), ARIMA(0,1,2), ARIMA(2,1,2), dan ARIMA(1,1,2). Untuk melihat model tentatif mana yang paling efektif untuk melakukan peramalan, dapat dilihat terlebih dahulu model mana yang memiliki seluruh koefisien signifikan.
mod1.ihsg <- Arima(train_ts, order = c(0,1,1), method = "ML")
mod2.ihsg <- Arima(train_ts, order = c(0,1,2), method = "ML")
mod3.ihsg <- Arima(train_ts, order = c(2,1,2), method = "ML")
mod4.ihsg <- Arima(train_ts, order = c(1,1,2), method = "ML")
lmtest::coeftest(mod1.ihsg)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.331246 0.099553 3.3273 0.0008768 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod2.ihsg)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.329745 0.102321 3.2227 0.00127 **
## ma2 0.062118 0.124090 0.5006 0.61666
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(mod3.ihsg)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.20710 0.44438 -0.4660 0.6412
## ar2 0.28122 0.20249 1.3888 0.1649
## ma1 0.56709 0.44632 1.2706 0.2039
## ma2 -0.15294 0.28332 -0.5398 0.5893lmtest::coeftest(mod4.ihsg)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.391738 0.638047 0.6140 0.5392
## ma1 -0.051600 0.634517 -0.0813 0.9352
## ma2 -0.041438 0.266804 -0.1553 0.8766Hanya model ARIMA(0,1,1) yang memiliki seluruh parameter signifikan. Sehingga model ini yang akan digunakan untuk melakukan peramalan.
Diagnostik Model
punten cariin jurnal kenapa diagnostik model itu penting dong
resid.mod1 <- mod1.ihsg$residuals
qqplot_resid <-
resid.mod1 %>%
ggplot(aes(sample = resid.mod1)) +
stat_qq() +
stat_qq_line() +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
theme_classic()
resid_plot <-
resid.mod1 %>%
ggplot(aes(x = df2$Date[1:length(resid.mod1)],
y = resid.mod1)) +
geom_line(color = "steelblue") +
geom_point(color = "steelblue") +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
xlab("Tanggal") +
ylab("Residual") +
theme_classic()
acf_resid <-
acf(resid.mod1, plot = F) %>%
autoplot() +
ggtitle(NULL) +
theme_classic()
pacf_resid <-
pacf(resid.mod1, plot = F) %>%
autoplot() +
ggtitle(NULL)+
theme_classic()
grid.arrange(qqplot_resid, resid_plot, acf_resid, pacf_resid)
Dari hasil diagnostik, terlihat bahwa sisaan model memiliki deviasi dari garis 45 derajat pada kedua ujungnya, sehingga harus dicek menggunakan uji formal untuk melihat apakah sisaan model sudah menyebar normal. Sisaan juga sudah terlihat menyebar di sekitar 0 dan pada plot ACF dan PACF sudah terlihat stasioner, sehingga dapat dikatakan sisaan tersebut sudah bebas satu sama lain. Akan tetapi, untuk lebih memastikan dapat menggunakan uji normal Anderson Darling untuk mengecek kenormalan, uji-t untuk melihat apakah benar nilai sisaan menyebar di sekitar 0, dan apakah sisaan sudah saling bebas dengan menggunakan uji Ljung-Box.
cari kelebihan uji anderson darling dibanding sama ks test sama hipotesisnya.
norm.resid <- ad.test(resid.mod1)
autocor.resid <- Box.test(resid.mod1, type = "Ljung-Box")
ttest.resid <- t.test(resid.mod1, mu = 0, conf.level = 0.95)
results <- data.frame("Statistics" = c(norm.resid$statistic,
autocor.resid$statistic,
ttest.resid$statistic),
"p-value" = c(norm.resid$p.value,
autocor.resid$p.value,
ttest.resid$p.value))
rownames(results) <- c("Anderson-Darling",
"Ljung-Box",
"t-Student")
results %>% knitr::kable(digits = 3)| Statistics | p.value | |
|---|---|---|
| Anderson-Darling | 0.692 | 0.068 |
| Ljung-Box | 0.074 | 0.786 |
| t-Student | 2.091 | 0.039 |
Ternyata dari hasil uji formal, hanya asumsi kebabasan sisaan dan normalitas yang terpenuhi, sehingga harus dilakukan penanganan terhadap hal tersebut. Salah satunya adalah dengan menambahkan drift ke dalam model ARIMA tersebut. cariin drift itu apaan di ARIMA.
model.drift <- Arima(train_ts, order = c(0, 1, 1), include.drift = T, method = "ML")
lmtest::coeftest(model.drift)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.29872 0.10222 2.9223 0.003475 **
## drift 8.21416 3.84807 2.1346 0.032792 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid.drift <- model.drift$residuals
norm.resid.drift <- ad.test(resid.drift)
autocor.resid.drift <- Box.test(resid.drift, type = "Ljung-Box")
ttest.resid.drift <- t.test(resid.drift, mu = 0, conf.level = 0.95)
results <- data.frame("Statistics" = c(norm.resid.drift$statistic,
autocor.resid.drift$statistic,
ttest.resid.drift$statistic),
"p-value" = c(norm.resid.drift$p.value,
autocor.resid.drift$p.value,
ttest.resid.drift$p.value))
rownames(results) <- c("Anderson-Darling",
"Ljung-Box",
"t-Student")
results %>% knitr::kable(digits = 3)| Statistics | p.value | |
|---|---|---|
| Anderson-Darling | 0.659 | 0.083 |
| Ljung-Box | 0.001 | 0.971 |
| t-Student | 0.019 | 0.985 |
Dengan menambahka drift, maka seluruh asumsi sisaan dapat terpenuhi, sehingga model ARIMA(0,1,1) dengan penambahan faktor drift yang akan digunakan.
Overfitting model
Overfitting adalah proses meningkatkan salah satu parameter dari model tentatif untuk melihat apakah terdapat model yang lebih baik di sekitar dugaan model tentatif yang telah ditentukan. Karena model yang digunakan adalah model ARIMA(0,1,1) dengan drift, maka model overfit yang digunakan adalah ARIMA(1,1,1), ARIMA(0,2,1), dan ARIMA(0,1,2).
overfit1 <- Arima(train_ts, order = c(1, 1, 1), include.drift = T, method = "ML")
overfit2 <- Arima(train_ts, order = c(0, 2, 1), include.drift = T, method = "ML")
overfit3 <- Arima(train_ts, order = c(1, 1, 2), include.drift = T, method = "ML")
lmtest::coeftest(overfit1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.094677 0.519592 0.1822 0.85541
## ma1 0.205064 0.525063 0.3906 0.69613
## drift 8.206711 3.941982 2.0819 0.03735 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.000000 0.087327 -11.451 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(overfit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.198024 1.036550 0.1910 0.84849
## ma1 0.102994 1.032386 0.0998 0.92053
## ma2 -0.034583 0.334546 -0.1034 0.91767
## drift 8.205802 3.944825 2.0801 0.03751 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hanya model ARIMA(0,2,1) dengan drift yang signifikan, sehingga akan dilihat AIC dan BICnya. Nilai AIC dan BIC yang lebih kecil menandakan bahwa model tersebut lebih baik untuk digunakan. Cariin juga jurnal disni sabi.
crit <- data.frame(AIC = c(model.drift$aic, overfit2$aic),
BIC = c(model.drift$bic, overfit2$bic))
rownames(crit) <- c("ARIMA(0,1,1)", "ARIMA(0,2,1)")
crit %>% knitr::kable()| AIC | BIC | |
|---|---|---|
| ARIMA(0,1,1) | 904.5861 | 912.2160 |
| ARIMA(0,2,1) | 906.6179 | 911.6831 |
Karena nilai AIC lebih kecil pada model ARIMA(0,1,1) dan BIC lebih kecil pada ARIMA(0,2,1), maka model terbaik tidak dapat ditentukan dari sana. Diagnostik untuk model overfit dapat dilihat sebagai berikut.
resid.mod1 <- overfit2$residuals
qqplot_resid <-
resid.mod1 %>%
ggplot(aes(sample = resid.mod1)) +
stat_qq() +
stat_qq_line() +
xlab("Theoretical Quantiles") +
ylab("Sample Quantiles") +
theme_classic()
resid_plot <-
resid.mod1 %>%
ggplot(aes(x = df2$Date[1:length(resid.mod1)],
y = resid.mod1)) +
geom_line(color = "steelblue") +
geom_point(color = "steelblue") +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
xlab("Tanggal") +
ylab("Residual") +
theme_classic()
acf_resid <-
acf(resid.mod1, plot = F) %>%
autoplot() +
ggtitle(NULL) +
theme_classic()
pacf_resid <-
pacf(resid.mod1, plot = F) %>%
autoplot() +
ggtitle(NULL)+
theme_classic()
grid.arrange(qqplot_resid, resid_plot, acf_resid, pacf_resid)
Sama seperti model tentatif sebelumnya, model overfit ARIMA(0,2,1) terlihat memiliki sisaan normal, menyebar di sekitar 0, dan saling bebas. Hasil uji formalnya dapat dilihat sebagai berikut.
norm.resid <- ad.test(resid.mod1)
autocor.resid <- Box.test(resid.mod1, type = "Ljung-Box")
ttest.resid <- t.test(resid.mod1, mu = 0, conf.level = 0.95)
results <- data.frame("Statistics" = c(norm.resid$statistic,
autocor.resid$statistic,
ttest.resid$statistic),
"p-value" = c(norm.resid$p.value,
autocor.resid$p.value,
ttest.resid$p.value))
rownames(results) <- c("Anderson-Darling",
"Ljung-Box",
"t-Student")
results %>% knitr::kable(digits = 3)| Statistics | p.value | |
|---|---|---|
| Anderson-Darling | 0.759 | 0.047 |
| Ljung-Box | 7.939 | 0.005 |
| t-Student | 0.843 | 0.401 |
Ternyata meskipun terlihat sudah memenuhi semua asumsi, hanya asumsi bahwa sisaan menyebar dengan nilai harapan 0 yang terpenuhi. Sehingga model yang akan digunakan adalah model ARIMA(0,1,1) dengan drift.
Peramalan data
Data peramalan dapat dilihat sebagai berikut.
forecast_arima <- forecast(model.drift, h = length(test_ts))
forecast_arima %>% knitr::kable(digits = 3,
caption = "Hasil peramalan model ARIMA(0,1,1) dengan drift")| Point Forecast | Lo 80 | Hi 80 | Lo 95 | Hi 95 | |
|---|---|---|---|---|---|
| 96 | 1285.835 | 1248.531 | 1323.139 | 1228.783 | 1342.887 |
| 97 | 1294.049 | 1232.904 | 1355.195 | 1200.535 | 1387.563 |
| 98 | 1302.263 | 1224.251 | 1380.276 | 1182.954 | 1421.573 |
| 99 | 1310.477 | 1218.646 | 1402.309 | 1170.033 | 1450.922 |
| 100 | 1318.692 | 1214.864 | 1422.520 | 1159.900 | 1477.483 |
| 101 | 1326.906 | 1212.331 | 1441.481 | 1151.678 | 1502.133 |
| 102 | 1335.120 | 1210.723 | 1459.517 | 1144.871 | 1525.369 |
| 103 | 1343.334 | 1209.836 | 1476.832 | 1139.166 | 1547.502 |
| 104 | 1351.548 | 1209.531 | 1493.566 | 1134.351 | 1568.745 |
| 105 | 1359.762 | 1209.709 | 1509.816 | 1130.275 | 1589.250 |
| 106 | 1367.977 | 1210.296 | 1525.657 | 1126.824 | 1609.129 |
| 107 | 1376.191 | 1211.235 | 1541.147 | 1123.912 | 1628.469 |
| 108 | 1384.405 | 1212.482 | 1556.328 | 1121.471 | 1647.339 |
| 109 | 1392.619 | 1214.000 | 1571.238 | 1119.445 | 1665.793 |
| 110 | 1400.833 | 1215.760 | 1585.906 | 1117.789 | 1683.878 |
| 111 | 1409.047 | 1217.738 | 1600.356 | 1116.466 | 1701.629 |
| 112 | 1417.261 | 1219.913 | 1614.610 | 1115.444 | 1719.079 |
| 113 | 1425.476 | 1222.268 | 1628.684 | 1114.696 | 1736.255 |
| 114 | 1433.690 | 1224.786 | 1642.593 | 1114.200 | 1753.180 |
| 115 | 1441.904 | 1227.456 | 1656.352 | 1113.935 | 1769.873 |
| 116 | 1450.118 | 1230.266 | 1669.970 | 1113.883 | 1786.353 |
| 117 | 1458.332 | 1233.205 | 1683.459 | 1114.030 | 1802.634 |
| 118 | 1466.546 | 1236.266 | 1696.827 | 1114.362 | 1818.731 |
| 119 | 1474.761 | 1239.439 | 1710.083 | 1114.867 | 1834.655 |
| 120 | 1482.975 | 1242.717 | 1723.232 | 1115.533 | 1850.417 |
| 121 | 1491.189 | 1246.095 | 1736.282 | 1116.351 | 1866.027 |
| 122 | 1499.403 | 1249.567 | 1749.239 | 1117.312 | 1881.494 |
| 123 | 1507.617 | 1253.127 | 1762.107 | 1118.408 | 1896.826 |
Plot dari peramalan tersebut dapat dilihat sebagai berikut.
forecast_arima %>% autoplot() +
xlab("Periode Waktu") +
ylab("Harga Closing IHSG") +
ggtitle(NULL) +
theme_classic()
Terlihat bahwa harga closing IHSG semakin meningkat seiring
berjalannya waktu. Peramalan yang berbentuk garis lurus yang naik ini
akibat tidak adanya faktor autoregressive pada model tersebut,
sehingga hanya proses differencing dan moving average
yang mempengaruhi peramalan harga closing IHSG tersebut.
Akurasi dari peramalan tersebut dapat dilihat sebagai berikut.
accuracy(ts(forecast_arima$mean, start = 1), test_ts)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -77.98828 166.1014 127.9881 -7.800949 11.11353 0.8154713 2.372467Terlihat bahwa error dari proses peramalan tersebut relatif kecil, sehingga model tersebut adalah model yang baik untuk digunakan dalam peramalan harga closing IHSG.
Daftar Pustaka
Fattah J, Ezzine L, Aman Z, El Moussami H, Lachhab A. 2018. Forecasting of demand using ARIMA model. Int J Eng Bus Manag. 10:1–9. doi:10.1177/1847979018808673.
Meyler, Aidan dan Kenny, Geoff dan Quinn, Terry. 1998. Munich Personal RePEc Archive Forecasting irish inflation using ARIMA models. Munich Pers RePEc Arch.(11359):1–8.
Rakhmawati D, Wahyudi R, Yuliawan CG. 2020. Pemodelan harga saham IHSG selama pandemi Covid-19 menggunakan ARIMA non-musiman. Jurnal Pro Bisnis. 13(2):39-48.