Input Data

bbca <- read.csv("BBCA.JK.csv")
bbca <- bbca[,-c(2,3,4,6,7)]
View(bbca)
bbca$Date <- as.POSIXct(bbca$Date)
bbca$Close <- as.numeric(bbca$Close)
ts.bbca <- ts(bbca)
View(ts.bbca)
kableExtra::kable(head(bbca) ,caption = 'Subset Data Harian Harga Tutupan Saham')
Subset Data Harian Harga Tutupan Saham
Date Close
2013-03-04 2140
2013-03-05 2100
2013-03-06 2140
2013-03-07 2150
2013-03-08 2160
2013-03-11 2150

Eksplorasi data

str(bbca)
## 'data.frame':    2489 obs. of  2 variables:
##  $ Date : POSIXct, format: "2013-03-04" "2013-03-05" ...
##  $ Close: num  2140 2100 2140 2150 2160 2150 2140 2130 2150 2130 ...
summary(bbca)
##       Date                         Close     
##  Min.   :2013-03-04 00:00:00   Min.   :1710  
##  1st Qu.:2015-09-21 00:00:00   1st Qu.:2650  
##  Median :2018-03-15 00:00:00   Median :4510  
##  Mean   :2018-03-09 12:22:16   Mean   :4603  
##  3rd Qu.:2020-08-25 00:00:00   3rd Qu.:6275  
##  Max.   :2023-03-03 00:00:00   Max.   :9300
bbca$Date <- as.Date(bbca$Date)
View(bbca)
ts.bbca <- ts(bbca)
str(bbca)
## 'data.frame':    2489 obs. of  2 variables:
##  $ Date : Date, format: "2013-03-03" "2013-03-04" ...
##  $ Close: num  2140 2100 2140 2150 2160 2150 2140 2130 2150 2130 ...
# Boxplot
y <- bbca$Date
x <- bbca$Close
boxplot(x, xlab="Harga Tutupan Saham BBCA", col="lightblue")

Data tidak terlihat memiliki pencilan

# plot harian
ggplot(bbca, aes(x=Date, y=Close)) +
  geom_line() +
  labs(title = "Plot Data Harian Harga Tutupan Saham BBCA",
       subtitle = "1 Januari 2012 - 31 Desember 2022") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA))

Berdasarkan plot time series yang ditampilkan di atas, pola data memiliki tren naik

p=(2142/2489)
f=as.integer(p*2489)
jumlah.train <- as.integer(p*nrow(bbca))
ggplot(bbca, aes(x=Date, y=Close)) +
  geom_line() + 
  labs(title = "Plot Data Harian Harga Tutupan Saham BBCA",
       subtitle = "1 Januari 2012 - 31 Desember 2022") +
  geom_vline(aes(xintercept = Date[jumlah.train], 
                 col="Pemisah_data_train"), lty=2, lwd=.7) +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA)) +
  scale_color_manual(name = "", values = c(Pemisah_data_train="blue"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.

Berdasarkan hasil eksplorasi dan pertimbangan terkait pola data,didapatkan titik potong yang dinilai tepat dalam membagi data training dan data testing, yaitu pada nilai proporsi data training sebesar 86.05%. Titik pemisah data training dan data testing tersebut bertepatan pada tanggal 10 November 2021.

bbca[jumlah.train,1]
## [1] "2021-10-10"

Pembagian Data Training dan Data Testing

train=1:f
training.ts <- ts(bbca$Close[train])
testing.ts <- ts(bbca$Close[-train], start = f+1)
#Plot data training
ggplot(bbca[1:2206,], aes(x=Date, y=Close)) +
  geom_line() + 
  labs(title = "Plot Time Series Data Training Harga Tutupan Saham BBCA",
       subtitle = "(Januari 2012 - November 2021)") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA))

#Plot data testing
ggplot(bbca[2143:2489,], aes(x=Date, y=Close)) +
  geom_line() + 
  labs(title = "Plot Time Series Data Testing Harga Tutupan Saham BBCA",
       subtitle = "(November 2021 - Desember 2022)") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom",
        panel.background = element_rect(fill=NA))

Cek Kestasioneran Data

par(mfrow=c(1,2))
acf(training.ts, main = "ACF Plot for MA(q)")
pacf(training.ts, main = "PACF Plot for AR(p)")

Hipotesis yang diuji: H0 : data tidak stasioner H1 : data stasioner

adf.test(training.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.ts
## Dickey-Fuller = -3.1837, Lag order = 12, p-value = 0.09071
## alternative hypothesis: stationary

Hasil di atas menunjukkan bahwa \(P-value\) hasil uji di atas lebih besar dari \(alpha=0.05\), sehingga terima H0. Artinya, pada taraf nyata 5% tidak cukup bukti untuk menyatakan bahwa data telah stasioner.

Penanganan Ketidakstasioneran

Differencing

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

Uji Kembali Kestasioneran Data

par(mfrow=c(1,2))
acf(train.diff, main = "ACF Plot for MA(q)")
pacf(train.diff, main = "PACF Plot for AR(p)")

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

Hasil di atas menunjukkan bahwa \(P-value=0.01\) kurang dari \(alpha=0.05\), sehingga tolak H0. Artinya, pada taraf nyata 5% cukup bukti untuk menyatakan bahwa data telah stasioner.

Pemodelan Arima

A. Penentuan Ordo Model Arima

Penentuan ordo model ARIMA dapat dilakukan dengan menampilkan plot ACF,PACF, dan EACF sebagai berikut:

#Identifikasi Model ARIMA dengan ACF, PACF, dan EACF
acf(train.diff, main="ACF data harga tutupan saham")

pacf(train.diff, main="PACF data harga tutupan saham")

eacf(train.diff)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x 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 x x o o o o o o o x  o  o  o 
## 4 x x x x o o o o o o o  o  o  o 
## 5 x x x x x o o o o o o  o  o  o 
## 6 x x x x x x o o o o x  o  o  o 
## 7 x x x x x o x o o o x  o  o  o

Berdasarkan plot di atas, didapatkan kandidat model ARIMA sebagai berikut: ARIMA(1,1,1) ARIMA(0,1,3) ARIMA(2,1,2) ARIMA(3,1,3) ARIMA(1,1,2)

B. Pemilihan Model Terbaik

Model terbaik dipilih berdasarkan AIC terkecil

model1 <- Arima(train.diff,order = c(1,1,1),method ="ML")
summary(model1)
## Series: train.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.0638  -1.0000
## s.e.   0.0216   0.0017
## 
## sigma^2 = 4507:  log likelihood = -12041.86
## AIC=24089.72   AICc=24089.73   BIC=24106.72
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE         ACF1
## Training set 1.165697 67.09058 41.12208 NaN  Inf 0.6417125 -0.002503759
model2 <- Arima(train.diff,order = c(0,1,3),method ="ML")
summary(model2)
## Series: train.diff 
## ARIMA(0,1,3) 
## 
## Coefficients:
##           ma1     ma2     ma3
##       -1.0635  0.0436  0.0199
## s.e.   0.0219  0.0337  0.0223
## 
## sigma^2 = 4507:  log likelihood = -12041.24
## AIC=24090.48   AICc=24090.5   BIC=24113.15
## 
## Training set error measures:
##                  ME     RMSE      MAE MPE MAPE      MASE         ACF1
## Training set 1.1938 67.07031 41.15324 NaN  Inf 0.6421987 -0.001702531
model3 <- Arima(train.diff,order = c(2,1,2),method ="ML")
summary(model3)
## Series: train.diff 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2      ma1      ma2
##       -0.6289  -0.0784  -0.4362  -0.5638
## s.e.   0.1762   0.0223   0.1757   0.1757
## 
## sigma^2 = 4500:  log likelihood = -12039.05
## AIC=24088.11   AICc=24088.14   BIC=24116.45
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE     MASE        ACF1
## Training set 1.192773 67.00182 41.27808 NaN  Inf 0.644147 0.001209986
model4 <- Arima(train.diff,order = c(3,1,3),method = "ML")
summary(model4)
## Series: train.diff 
## ARIMA(3,1,3) 
## 
## Coefficients:
##           ar1      ar2     ar3      ma1      ma2      ma3
##       -0.4445  -0.0999  0.0347  -0.6188  -0.3328  -0.0484
## s.e.   0.7451   0.4694  0.0618   0.7453   0.4063   0.4235
## 
## sigma^2 = 4500:  log likelihood = -12037.98
## AIC=24089.96   AICc=24090.02   BIC=24129.64
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE          ACF1
## Training set 1.155287 66.96921 41.31899 NaN  Inf 0.6447854 -0.0006955678
model5 <- Arima(train.diff,order = c(1,1,2),method = "ML")
summary(model5)
## Series: train.diff 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.1080  -1.1741  0.1741
## s.e.  0.2049   0.2027  0.2027
## 
## sigma^2 = 4508:  log likelihood = -12041.48
## AIC=24090.97   AICc=24090.98   BIC=24113.64
## 
## Training set error measures:
##                    ME    RMSE      MAE MPE MAPE      MASE         ACF1
## Training set 1.182305 67.0784 41.12664 NaN  Inf 0.6417837 0.0002000068
Model <- c("ARIMA (1,1,1)","ARIMA (0,1,3)","ARIMA (2,1,2)","ARIMA (3,1,3)","ARIMA(1,1,2)")
AIC <- c(model1$aic,model2$aic,model3$aic,model4$aic,model5$aic)
BIC <- c(model1$bic,model2$bic,model3$bic,model4$bic,model5$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi)
Model AIC BIC
ARIMA (1,1,1) 24089.72 24106.72
ARIMA (0,1,3) 24090.48 24113.15
ARIMA (2,1,2) 24088.11 24116.45
ARIMA (3,1,3) 24089.96 24129.64
ARIMA(1,1,2) 24090.97 24113.64 
paste("Model yang terbaik berdasarkan akurasi adalah model",Akurasi$Model[which.min(Akurasi[,"AIC"])])
## [1] "Model yang terbaik berdasarkan akurasi adalah model ARIMA (2,1,2)"

C. Overfitting

model6 <- Arima(train.diff, order=c(3,1,2),method="ML") #intersep signifikan
lmtest::coeftest(model6)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ar1 -0.419644   0.286043 -1.4671  0.14236  
## ar2 -0.049953   0.030289 -1.6492  0.09910 .
## ar3  0.038221   0.026881  1.4219  0.15506  
## ma1 -0.643680   0.285640 -2.2535  0.02423 *
## ma2 -0.356317   0.285635 -1.2475  0.21223  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model7 <- Arima(train.diff, order=c(2,1,3),method="ML") #intersep signifikan
lmtest::coeftest(model7)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.456587   0.435160 -1.0492   0.2941
## ar2 -0.124639   0.229814 -0.5423   0.5876
## ma1 -0.610164   0.437727 -1.3939   0.1633
## ma2 -0.331578   0.508030 -0.6527   0.5140
## ma3 -0.058217   0.234351 -0.2484   0.8038

Model terbaik

Model <- c("ARIMA (2,1,2)","ARIMA (3,1,2)","ARIMA(2,1,3)")
AIC <- c(model3$aic, model6$aic, model7$aic)
BIC <- c(model3$bic, model6$bic, model7$bic)
Akurasi <- data.frame(Model,AIC,BIC)
kableExtra::kable(Akurasi)
Model AIC BIC
ARIMA (2,1,2) 24088.11 24116.45
ARIMA (3,1,2) 24088.01 24122.02
ARIMA(2,1,3) 24090.40 24124.42

Model terbaik adalah ARIMA(3,1,2)

D. Analisis Sisaan Model Tentatif Arima(3,1,2)

Eksplorasi

sisaan <- model6$residuals
tsdiag(model6)

par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(sisaan, type="o", 
     ylab = "Sisaan", xlab = "Order", main = "Sisaan Model vs Order")
abline(h = 0, col='red')
acf(sisaan)
pacf(sisaan)

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

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

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

Uji Formal

1) Sisaan Menyebar Normal

Hipotesis yang digunakan adalah sebagai berikut: H0: Sisaan mengikuti sebaran normal H1: Sisaan tidak mengikuti sebaran normal

jarque.bera.test(sisaan)
## 
##  Jarque Bera Test
## 
## data:  sisaan
## X-squared = 19076, df = 2, p-value < 2.2e-16

Berdasarkan Jarque Bera Test, diperoleh p-value < α (0.05), maka Tolak H0. Artinya, cukup bukti untuk menyatakan bahwa sisaan tidak menyebar normal pada taraf nyata 5%.

2) Sisaan saling bebas/tidak ada autokorelasi Hipotesis yang digunakan adalah sebagai berikut: H0: Tidak ada autokorelasi H1: Ada autokorelasi

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

Berdasarkan hasil uji di atas, didapat nilai \(P−Value=0.973\) yang berarti TERIMA H0. Artinya, pada taraf nyata 5%, ada bukti untuk menyatakan bahwa tidak ada autokorelasi pada data.

3) Nilai tengah sisaan sama dengan nol Hipotesis yang digunakan adalah sebagai berikut: H0: Nilai tengah sisaan bernilai nol H1: Nilai tengah sisaan tak bernilai nol

t.test(sisaan, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.79752, df = 2140, p-value = 0.4252
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.684205  3.992970
## sample estimates:
## mean of x 
##  1.154383

Berdasarkan hasil uji di atas, didapat nilai \(P−Value=0.4252\) yang berarti TERIMA H0. Artinya, pada taraf nyata 5%, ada bukti untuk menyatakan bahwa nilai tengah sisaan bernilai nol.

E. Identifikasi Efek ARCH

Identifikasi adanya proses ARCH pada sisaan model ARIMA(3,1,2) dapat dilakukan dengan uji Lagrange Multiplier (LM). Hipotesis yang digunakan adalah sebagai berikut: H0 : Tidak terjadi heteroskedastisitas H1 : Terjadi Heteroskedastisitas

for (i in 1:15) {
  ArchTest <- ArchTest(model5$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 1.464725e-09 
## P Value LM Test lag ke 2 adalah 0 
## P Value LM Test lag ke 3 adalah 0 
## P Value LM Test lag ke 4 adalah 0 
## P Value LM Test lag ke 5 adalah 0 
## P Value LM Test lag ke 6 adalah 0 
## P Value LM Test lag ke 7 adalah 0 
## P Value LM Test lag ke 8 adalah 0 
## P Value LM Test lag ke 9 adalah 0 
## P Value LM Test lag ke 10 adalah 0 
## P Value LM Test lag ke 11 adalah 0 
## P Value LM Test lag ke 12 adalah 0 
## P Value LM Test lag ke 13 adalah 0 
## P Value LM Test lag ke 14 adalah 0 
## P Value LM Test lag ke 15 adalah 0

Berdasarkan hasil uji diperoleh nilai \(P-Value\) sampai lag ke-15 yang didapat lebih kecil dari α = 5%, maka keputusan yang diambil adalah Tolak H0. Dengan kata lain, dapat disimpulkan bahwa terdapat unsur heteroskedastisitas pada sisaan model ARIMA(3,1,2). Kondisi yang menolak H0 ini terjadi sampai dengan lag ke-15 sehingga dapat dikatakan bahwa indikasi pemodelan sisaan ini lebih cocok menggunakan model GARCH dibandingkan model ARCH. Hal tersebut bersesuaian dengan pernyataan Desvina & Meijer (2018) yang menyebutkan bahwa ketika pemeriksaan sisaan yang dilakukan dari lag 1 sampai 12 signifikan, maka model ARCH lebih cocok digunakan. Akan tetapi, jika pemeriksaan sisaan dilakukan hingga lebih dari lag 12, maka model GARCH lebih cocok digunakan. Oleh karena itu, selanjutnya akan dilakukan pendugaan model dengan menggunakan GARCH.

F. Pemodelan GARCH

Pendugaan model GARCH untuk mendapat model terbaik ditentukan dengan cmempertimbangkan nilai \(p-value\) ARCH-LM, nilai AIC terkecil, dan signifikansi parameternya.

Hasil trial and error menghasilkan model ARIMA(3,1,2)-GARCH(3,7) dengan \(p-value\) uji ARCH-LM sebesar 0,9774, artinya efek ARCH sudah teratasi. Selanjutnya akan dilakukan overfitting pada model GARCH.

Dari overfitting didapatkan ringkasan tabel sebagai berikut.

Model Signifikansi Parameter AIC
ARIMA(3,1,2)-GARCH(3,7) Signifikan 10.76075
ARIMA(3,1,2)-GARCH(3,8) Signifikan 10.76196
ARIMA(3,1,2)-GARCH(4,7) Signifikan 10.76162

Berdasarkan tabel di atas, didapat bahwa model ARIMA(3,1,2)-GARCH(3,7) telah menghasilkan nilai AIC paling minimum sehingga dipilih menjadi model terbaik. Selanjutnya, sebelum melakukan peramalan, dilakukan uji diagnostik kembali pada model ARIMA(3,1,2)-GARCH(3,7) dengan uji Shapiro.Wilk dan uji LM ARCH dan hasilnya sebagai berikut:

Jenis Pengujian P-value
Normalitas 0.00
Efek ARCH 0.97

Tabel di atas menunjukkan hasil uji diagnostik model ARIMA(3,1,2) GARCH(3,7) untuk normalitas dan heterokesdatisitas. Berdasarkan nilai \(P-value\) yang didapatkan dari kedua uji tersebut, model tidak memenuhi asumsi kenormalan dan telah memenuhi asumsi kehomogenan ragam karena memiliki p-value yang lebih besar daripada \(α = 5%\). Hasil tersebut menunjukkan bahwa model GARCH sudah mengatasi heteroskesdastisitas tetapi memiliki ketidaknormalan sisaan hasil uji formal.

G. Peramalan

Peramalan dilakukan dengan menggunakan model terbaik yang didapatka dari pendugaan model GARCH. Hasil ramalan harga tutupan saham BBCA menggunakan model ARIMA (3,1,2)-GARCH (3,7) selama 347 periode ke depan ditunjukkan pada plot di bawah:

# Peramalan
garchSpec <- ugarchspec(
  variance.model=list(model="sGARCH",
                      garchOrder=c(3,7)),
  mean.model=list(armaOrder=c(3,2)),
  distribution.model="std")
garchFitt <- ugarchfit(spec=garchSpec, data=train.diff)
forc<- ugarchforecast(fitORspec = garchFitt, data = tlkm$Close, n.ahead = 347, n.roll = 0)
plot(forc, which= 1)

Hasil forecasting di atas masih perlu ditransformasi balik karena adanya prosess differencing. Hasil transformasi balik serta perbandingan data ramalan dengan data aktual dapat dilihat di pada tabel berikut.

pt_1 <- bbca$Close[2142] #nilai akhir data latih
hasil.forc.Diff <- forc@forecast$seriesFor[,1]
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1
perbandingan <- data.frame("Aktual"= testing.ts,
                           "Ramalan" = hasil[-1])
head(perbandingan,10)
##    Aktual  Ramalan
## 1    7320 7245.837
## 2    7525 7243.032
## 3    7750 7238.920
## 4    7650 7238.737
## 5    7525 7237.771
## 6    7500 7238.728
## 7    7400 7239.257
## 8    7525 7240.727
## 9    7525 7241.967
## 10   7525 7243.667
Date<-c(bbca$Date[2143:2489]) 
dataframe <- data.frame(Date, perbandingan) 
T <- nrow(dataframe) 
MAPE <- 1/T*sum(abs((dataframe$Aktual-dataframe$Ramalan)/dataframe$Aktual)*100)
MAPE
## [1] 5.891575

Menurut Marikar (2019), nilai MAPE yang kurang dari 10% dikategorikan sebagai hasil peramalan yang sangat baik. Hasil akurasi peramalan data nilai tukar rupiah menggunakan nilai MAPE menunjukkan angka 5.9%. Hal ini berarti peramalan data nilai tukar rupiah terhadap dolar AS dengan menggunakan model ARIMA(3,1,2)-GARCH(3,7) sudah tepat digunakan karena cukup mendekati data aktual.