1 Import Library

> library(tseries)
> library(TSA)
> library(lmtest)
> library(forecast)
> library(MLmetrics)
> library(xts)
> library(readxl)
> library(tsoutliers)
> library(knitr)
> library(FSA)

2 Import Data

Data yang digunakan merupakan data Indeks Harga Saham Gabungan (IHSG) dari 12 Desember 2023 sampai 9 April 2025. Data ini mencerminkan pergerakan harga saham di pasar modal Indonesia yang dipengaruhi oleh berbagai faktor ekonomi makro, termasuk inflasi, kebijakan moneter, dan sentimen pasar global.

> # Import data
> IHSG <- read_excel("C:/Users/ahmad/Downloads/Ahmad Fathan Muzhaffar - ihsg 2.xlsx")
> 
> # Konversi kolom tanggal
> IHSG$Tanggal <- as.Date(IHSG$Tanggal, format = "%d-%b-%y")
> 
> # Definisi variabel
> tanggal <- IHSG$Tanggal
> H.saham <- IHSG$`Harga Saham Penutup`

3 Eksplorasi Data Deret Waktu

Sebelum melakukan analisis deret waktu, langkah awal yang perlu dilakukan adalah eksplorasi data untuk mengetahui pola umum, tren, serta indikasi adanya outlier atau gangguan pada data.

> plot(tanggal, H.saham,
+      type = "o", col = "blue", lwd = 1,
+      xlab = "Tanggal", ylab = "Harga Saham",
+      main = "Plot Indeks Harga Saham Harian")

Berdasarkan grafik di atas, pada awal tahun 2024 pergerakan harga terlihat relatif stabil di kisaran 7.000–7.300. Harga kemudian meningkat hingga mencapai titik tertinggi sekitar 7.700–7.800 pada pertengahan tahun. Setelah itu, terjadi penurunan secara bertahap, dan memasuki awal tahun 2025 penurunannya menjadi lebih tajam hingga mendekati 6.000. Fluktuasi harga yang semakin besar di akhir periode mengindikasikan peningkatan volatilitas. Secara umum, data ini belum bersifat stasioner karena adanya perubahan tren yang cukup jelas.

> plot(tanggal, H.saham,
+      type = "o", col = "blue", lwd = 2,
+      xlab = "Tanggal", ylab = "Harga Saham",
+      main = "Plot IHSG dengan Titik Intervensi")
> abline(v = tanggal[c(150, 290, 310)],
+        col = "red", lty = 3, lwd = 1.5)
> legend("topright",
+        legend = c("IHSG", "Titik Intervensi"),
+        col    = c("blue", "red"),
+        lty    = c(1, 3), lwd = c(2, 1.5),
+        bty    = "n")

4 Uji Stasioneritas Rata-rata

Stasioneritas merupakan asumsi penting dalam analisis deret waktu, khususnya pada model ARIMA. Data dikatakan stasioner apabila memiliki rata-rata dan varians yang konstan sepanjang waktu. Pengujian stasioneritas rata-rata dilakukan menggunakan uji Augmented Dickey-Fuller (ADF) dengan model sebagai berikut:

\[\Delta Z_t = \alpha + \beta t + \delta Z_{t-1} + \sum_{j=1}^{p} \gamma_j \Delta Z_{t-j} + \varepsilon_t\]

Hipotesis:

\[H_0 : \delta = 0 \quad \text{(data tidak stasioner)}\] \[H_1 : \delta < 0 \quad \text{(data stasioner)}\]

Apabila data tidak stasioner, dilakukan diferensiasi orde pertama:

\[\nabla Z_t = Z_t - Z_{t-1} = (1-B)Z_t\]

4.1 Uji Box-Cox — ARIMA

> BoxCox.lambda(H.saham)
[1] 1.999924

4.2 Uji Box-Cox — Pra-Intervensi

> # Dijalankan setelah objek pra.intervensi terbentuk (lihat Bagian 3.3)
> BoxCox.lambda(pra.intervensi)

5 Identifikasi Model

5.1 Model ARIMA

Identifikasi model bertujuan untuk menentukan orde model ARIMA\((p, d, q)\), yaitu parameter AR (\(p\)), differencing (\(d\)), dan MA (\(q\)). Proses ini dilakukan dengan mengamati pola plot ACF dan PACF.

Persamaan umum model ARIMA\((p, d, q)\):

\[\phi(B)(1-B)^d Z_t = \theta(B)\varepsilon_t\]

dengan komponen:

AR(\(p\)): \[Z_t = \phi_1 Z_{t-1} + \cdots + \phi_p Z_{t-p} + \varepsilon_t\]

MA(\(q\)): \[Z_t = \varepsilon_t - \theta_1 \varepsilon_{t-1} - \cdots - \theta_q \varepsilon_{t-q}\]

Panduan identifikasi:

ACF → menentukan orde MA(\(q\)): cut-off pada lag ke-\(q\)
PACF → menentukan orde AR(\(p\)): cut-off pada lag ke-\(p\)

> # Uji ADF data asli
> adf.test(H.saham)

    Augmented Dickey-Fuller Test

data:  H.saham
Dickey-Fuller = -0.59149, Lag order = 6, p-value = 0.9773
alternative hypothesis: stationary
> 
> # Differencing orde 1
> diff.H_saham <- diff(H.saham)
> adf.test(diff.H_saham)

    Augmented Dickey-Fuller Test

data:  diff.H_saham
Dickey-Fuller = -7.703, Lag order = 6, p-value = 0.01
alternative hypothesis: stationary
> 
> # Plot data setelah differencing
> plot(diff.H_saham, type = "l",
+      main = "Plot Data Setelah Differencing",
+      xlab = "Waktu", ylab = expression(nabla * Z[t]))

> 
> # ACF dan PACF
> acf(diff.H_saham,  main = "ACF Data Setelah Differencing")

> pacf(diff.H_saham, main = "PACF Data Setelah Differencing")

5.2 Deteksi Outlier

Outlier dalam deret waktu adalah observasi yang menyimpang secara signifikan dari pola umum data. Chen dan Liu (1993) mengklasifikasikan empat jenis outlier utama sebagai berikut.

Additive Outlier (AO) — gangguan sesaat pada satu titik waktu \(T\): \[Z_t = Z_t^N + \omega_{AO} \cdot I_t^{(T)}\]

Innovational Outlier (IO) — gangguan yang merambat melalui komponen inovasi: \[\phi(B)(1-B)^d Z_t = \theta(B)\varepsilon_t + \omega_{IO} \cdot \theta(B) \cdot I_t^{(T)}\]

Level Shift (LS) — perubahan permanen pada rata-rata: \[Z_t = Z_t^N + \omega_{LS} \cdot S_t^{(T)}\]

Temporary Change (TC) — gangguan yang meluruh secara eksponensial: \[Z_t = Z_t^N + \omega_{TC} \cdot \delta^{(t-T)} \cdot I_{t \geq T}, \quad 0 < \delta < 1\]

> # Residual dari model ARIMA awal (digunakan untuk deteksi outlier)
> prediksi.arima <- H.saham - resid.arima1
Error: object 'resid.arima1' not found
> 
> # Deteksi outlier dengan TSO
> outlier <- tso(prediksi.arima, types = c("AO", "IO", "LS", "TC"))
Error: object 'prediksi.arima' not found
> print(outlier)
Error: object 'outlier' not found
> 
> # Tabel outlier yang terdeteksi
> outlier$outliers
Error: object 'outlier' not found
> 
> # Plot efek outlier
> plot(outlier$effects,
+      main = "Efek Outlier yang Terdeteksi")
Error: object 'outlier' not found

5.3 Model ARIMA Pra-Intervensi

Model pra-intervensi dibangun menggunakan data sebelum terjadinya intervensi (\(t < T\), dengan \(T = 151\)) untuk menangkap struktur stokastik murni dari deret waktu:

\[\phi(B)(1-B)^d Z_t = \theta(B)a_t, \quad t < T\]

> pra.intervensi <- H.saham[1:150]

5.4 Model ARIMA Intervensi

Model intervensi merupakan perluasan ARIMA dengan menambahkan komponen deterministik yang merepresentasikan dampak kejadian eksternal (Box & Tiao, 1975):

\[Z_t = m_t + N_t\]

di mana \(m_t\) adalah fungsi respons intervensi dan \(N_t \sim \text{ARIMA}(p,d,q)\).

6 Pendugaan Parameter

6.1 Model ARIMA

> # Model ARIMA(19,1,0) dengan fixed parameter
> arima.model1 <- arima(H.saham,
+                       order          = c(19, 1, 0),
+                       method         = "ML",
+                       transform.pars = FALSE,
+                       include.mean   = FALSE,
+                       fixed = c(
+                         0, NA, rep(0, 7),   # ar1 (fixed), ar2 (estimasi), ar3-ar9 (fixed)
+                         NA,                  # ar10 (estimasi)
+                         rep(0, 8),           # ar11-ar18 (fixed)
+                         NA                   # ar19 (estimasi)
+                       ))
> 
> # Model ARIMA(0,1,19) sebagai pembanding
> arima.model2 <- arima(H.saham,
+                       order          = c(0, 1, 19),
+                       method         = "ML",
+                       transform.pars = TRUE,
+                       include.mean   = FALSE,
+                       fixed = c(
+                         0, NA, rep(0, 7),
+                         NA,
+                         rep(0, 8),
+                         NA
+                       ))
> 
> arima.model1

Call:
arima(x = H.saham, order = c(19, 1, 0), include.mean = FALSE, transform.pars = FALSE, 
    fixed = c(0, NA, rep(0, 7), NA, rep(0, 8), NA), method = "ML")

Coefficients:
      ar1      ar2  ar3  ar4  ar5  ar6  ar7  ar8  ar9    ar10  ar11  ar12  ar13
        0  -0.1234    0    0    0    0    0    0    0  0.1706     0     0     0
s.e.    0   0.0595    0    0    0    0    0    0    0  0.0624     0     0     0
      ar14  ar15  ar16  ar17  ar18     ar19
         0     0     0     0     0  -0.1914
s.e.     0     0     0     0     0   0.0632

sigma^2 estimated as 5209:  log likelihood = -1766.92,  aic = 3539.84
> arima.model2

Call:
arima(x = H.saham, order = c(0, 1, 19), include.mean = FALSE, transform.pars = TRUE, 
    fixed = c(0, NA, rep(0, 7), NA, rep(0, 8), NA), method = "ML")

Coefficients:
      ma1      ma2  ma3  ma4  ma5  ma6  ma7  ma8  ma9    ma10  ma11  ma12  ma13
        0  -0.0994    0    0    0    0    0    0    0  0.1256     0     0     0
s.e.    0   0.0560    0    0    0    0    0    0    0  0.0636     0     0     0
      ma14  ma15  ma16  ma17  ma18     ma19
         0     0     0     0     0  -0.1730
s.e.     0     0     0     0     0   0.0613

sigma^2 estimated as 5273:  log likelihood = -1768.67,  aic = 3543.34

6.2 Model ARIMA + Outlier

> n <- length(H.saham)
> 
> # --- Variabel Outlier ---
> # AO (Additive Outlier) → pulse pada satu titik
> ao1 <- rep(0, n); ao1[291] <- 1   # AO pada t = 291
> ao2 <- rep(0, n); ao2[311] <- 1   # AO pada t = 311
> 
> # TC (Temporary Change) → decay eksponensial mulai t = 151
> tc1 <- rep(0, n)
> for (i in 151:n) {
+   tc1[i] <- 0.7^(i - 151)
+ }
> 
> # Gabungkan sebagai matriks xreg
> xreg <- cbind(tc1, ao1, ao2)
> 
> # --- Fixed Parameter ---
> # Urutan: ar1..ar19 (19 parameter) + tc1, ao1, ao2 (3 parameter) = 22 total
> fixed_vec        <- rep(0, 22)
> fixed_vec[2]     <- NA    # ar2  → diestimasi
> fixed_vec[10]    <- NA    # ar10 → diestimasi
> fixed_vec[19]    <- NA    # ar19 → diestimasi
> fixed_vec[20:22] <- NA    # tc1, ao1, ao2 → diestimasi
> 
> # Verifikasi
> cat("Jumlah parameter:", length(fixed_vec), "\n")
Jumlah parameter: 22 
> cat("Parameter aktif (NA):", sum(is.na(fixed_vec)), "\n")
Parameter aktif (NA): 6 
> 
> # --- Fitting Model ---
> arima.tc <- arima(H.saham,
+                   order          = c(19, 1, 0),
+                   xreg           = xreg,
+                   method         = "ML",
+                   transform.pars = FALSE,
+                   include.mean   = FALSE,
+                   fixed          = fixed_vec)
> arima.tc

Call:
arima(x = H.saham, order = c(19, 1, 0), xreg = xreg, include.mean = FALSE, transform.pars = FALSE, 
    fixed = fixed_vec, method = "ML")

Coefficients:
      ar1      ar2  ar3  ar4  ar5  ar6  ar7  ar8  ar9    ar10  ar11  ar12  ar13
        0  -0.1523    0    0    0    0    0    0    0  0.1352     0     0     0
s.e.    0   0.0609    0    0    0    0    0    0    0  0.0649     0     0     0
      ar14  ar15  ar16  ar17  ar18     ar19      tc1       ao1       ao2
         0     0     0     0     0  -0.1601  16.6978  153.7717  -14.1055
s.e.     0     0     0     0     0   0.0659  64.5262   51.7529   71.6574

sigma^2 estimated as 5061:  log likelihood = -1762.29,  aic = 3536.59

6.3 Model ARIMA Pra-Intervensi

> # ARIMA(18,1,0) — kandidat 1
> model.pra1 <- arima(pra.intervensi,
+                     order          = c(18, 1, 0),
+                     method         = "ML",
+                     transform.pars = TRUE)
> 
> # ARIMA(0,1,18) — kandidat 2
> model.pra2 <- arima(pra.intervensi,
+                     order  = c(0, 1, 18),
+                     method = "ML")
> 
> model.pra1

Call:
arima(x = pra.intervensi, order = c(18, 1, 0), transform.pars = TRUE, method = "ML")

Coefficients:
          ar1      ar2      ar3     ar4     ar5      ar6     ar7      ar8
      -0.0434  -0.1543  -0.0962  0.1282  0.0670  -0.0411  0.1574  -0.0135
s.e.   0.0849   0.0854   0.0855  0.0866  0.0875   0.0875  0.0865   0.0875
          ar9    ar10    ar11     ar12     ar13     ar14     ar15     ar16
      -0.0611  0.0001  0.1066  -0.1373  -0.0522  -0.0494  -0.0077  -0.1159
s.e.   0.0887  0.0868  0.0870   0.0864   0.0877   0.0871   0.0869   0.0863
        ar17     ar18
      0.0670  -0.2320
s.e.  0.0871   0.0872

sigma^2 estimated as 2686:  log likelihood = -800.72,  aic = 1637.43
> model.pra2

Call:
arima(x = pra.intervensi, order = c(0, 1, 18), method = "ML")

Coefficients:
          ma1      ma2      ma3     ma4     ma5      ma6     ma7      ma8
      -0.0936  -0.2036  -0.1271  0.1458  0.1281  -0.1551  0.0706  -0.0489
s.e.   0.1697   0.1646   0.1411  0.1367  0.1520   0.1607  0.1385   0.1482
          ma9     ma10    ma11     ma12     ma13     ma14    ma15     ma16
      -0.0630  -0.1044  0.2184  -0.1536  -0.1698  -0.0165  0.1626  -0.2337
s.e.   0.1396   0.1441  0.1218   0.1438   0.1270   0.1179  0.1283   0.1167
         ma17     ma18
      -0.1162  -0.2197
s.e.   0.1330   0.1114

sigma^2 estimated as 2462:  log likelihood = -796.91,  aic = 1629.82

6.4 Model ARIMA Intervensi

> # Residual untuk identifikasi orde intervensi
> resid.arima1    <- residuals(arima.model1)
> resid.modelpra2 <- residuals(model.pra2)
> 
> n           <- length(H.saham)
> break_point <- 151
> 
> # Residual gabungan (pra dan pasca intervensi)
> error_idintv                    <- rep(NA, n)
> error_idintv[1:(break_point-1)] <- resid.arima1[1:(break_point-1)]
> error_idintv[break_point:n]     <- resid.modelpra2[break_point:n]
> 
> # --- Variabel Input Intervensi ---
> # AO (Additive Outlier) → pulse
> AO291 <- c(rep(0, 290), 1, rep(0, n - 291))
> AO311 <- c(rep(0, 310), 1, rep(0, n - 311))
> 
> # TC (Temporary Change / Step)
> TC151 <- c(rep(0, 150), rep(1, n - 150))
> 
> # Matriks xreg intervensi
> X <- cbind(TC151, AO291, AO311)
> 
> # --- Fitting Model Intervensi ---
> model.intervensi <- arima(H.saham,
+                           order  = c(0, 1, 18),
+                           xreg   = X,
+                           method = "ML")
> 
> model.intervensi

Call:
arima(x = H.saham, order = c(0, 1, 18), xreg = X, method = "ML")

Coefficients:
         ma1      ma2     ma3     ma4     ma5      ma6      ma7     ma8
      0.0691  -0.1686  0.0346  0.1291  0.0472  -0.1078  -0.1936  0.1047
s.e.  0.0614   0.0605  0.0650  0.0692  0.0632   0.0738   0.0681  0.0660
          ma9    ma10    ma11     ma12    ma13    ma14    ma15     ma16
      -0.0288  0.0395  0.0390  -0.1796  0.0428  0.1628  0.0994  -0.0511
s.e.   0.0779  0.0737  0.0891   0.0636  0.0813  0.0669  0.0771   0.0841
         ma17     ma18    TC151     AO291    AO311
      -0.1369  -0.0214  38.0870  239.4615   1.3891
s.e.   0.0644   0.0775  65.7987   48.7688  77.2433

sigma^2 estimated as 4618:  log likelihood = -1748.68,  aic = 3539.36
> coeftest(model.intervensi)

z test of coefficients:

        Estimate Std. Error z value  Pr(>|z|)    
ma1     0.069093   0.061397  1.1253  0.260441    
ma2    -0.168595   0.060478 -2.7877  0.005308 ** 
ma3     0.034562   0.065023  0.5315  0.595046    
ma4     0.129118   0.069201  1.8659  0.062061 .  
ma5     0.047219   0.063245  0.7466  0.455297    
ma6    -0.107841   0.073776 -1.4617  0.143811    
ma7    -0.193553   0.068058 -2.8439  0.004456 ** 
ma8     0.104733   0.066047  1.5857  0.112797    
ma9    -0.028845   0.077899 -0.3703  0.711168    
ma10    0.039515   0.073714  0.5361  0.591913    
ma11    0.039011   0.089078  0.4379  0.661427    
ma12   -0.179563   0.063602 -2.8232  0.004754 ** 
ma13    0.042771   0.081305  0.5261  0.598851    
ma14    0.162844   0.066864  2.4355  0.014873 *  
ma15    0.099429   0.077073  1.2901  0.197032    
ma16   -0.051076   0.084076 -0.6075  0.543521    
ma17   -0.136910   0.064357 -2.1273  0.033392 *  
ma18   -0.021359   0.077464 -0.2757  0.782755    
TC151  38.087037  65.798746  0.5788  0.562696    
AO291 239.461463  48.768839  4.9101 9.101e-07 ***
AO311   1.389134  77.243314  0.0180  0.985652    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7 Pengujian Asumsi

Diagnostik model dilakukan untuk memverifikasi bahwa residual memenuhi asumsi white noise: tidak berkorelasi, berdistribusi normal, dan memiliki varians homogen.

Uji Autokorelasi — Ljung-Box:

\[Q^* = n(n+2)\sum_{k=1}^{m} \frac{r_k^2(\hat{\varepsilon})}{n-k} \sim \chi^2(m-p-q)\]

\[H_0 : \rho_1 = \rho_2 = \cdots = \rho_m = 0 \quad \text{(residual white noise)}\] \[H_1 : \exists\, k \text{ dengan } \rho_k \neq 0 \quad \text{(residual tidak independen)}\]

Uji Normalitas — Jarque-Bera:

\[JB = \frac{n}{6}\left(S^2 + \frac{(K-3)^2}{4}\right) \sim \chi^2(2)\]

\[H_0 : \text{residual berdistribusi normal}\] \[H_1 : \text{residual tidak berdistribusi normal}\]

7.1 Asumsi Model ARIMA

> resid.arima1 <- residuals(arima.model1)
> resid.arima2 <- residuals(arima.model2)
> 
> # Uji autokorelasi
> Box.test(resid.arima1, lag = 20, type = "Ljung-Box")

    Box-Ljung test

data:  resid.arima1
X-squared = 19.039, df = 20, p-value = 0.5193
> Box.test(resid.arima2, lag = 20, type = "Ljung-Box")

    Box-Ljung test

data:  resid.arima2
X-squared = 19.755, df = 20, p-value = 0.4734
> 
> # Uji normalitas
> jarque.bera.test(resid.arima1)

    Jarque Bera Test

data:  resid.arima1
X-squared = 223.39, df = 2, p-value < 2.2e-16
> jarque.bera.test(resid.arima2)

    Jarque Bera Test

data:  resid.arima2
X-squared = 349.47, df = 2, p-value < 2.2e-16
> 
> # Perbandingan AIC
> cat("AIC Model 1:", arima.model1$aic, "\n")
AIC Model 1: 3539.839 
> cat("AIC Model 2:", arima.model2$aic, "\n")
AIC Model 2: 3543.337

7.2 Asumsi Model ARIMA + Outlier

> resid.arimatc <- residuals(arima.tc)
> 
> # Uji autokorelasi
> Box.test(resid.arimatc, lag = 20, type = "Ljung-Box")

    Box-Ljung test

data:  resid.arimatc
X-squared = 24.867, df = 20, p-value = 0.2066
> 
> # Uji normalitas
> jarque.bera.test(resid.arimatc)

    Jarque Bera Test

data:  resid.arimatc
X-squared = 411.91, df = 2, p-value < 2.2e-16

7.3 Asumsi Model ARIMA Pra-Intervensi

> resid.modelpra1 <- residuals(model.pra1)
> resid.modelpra2 <- residuals(model.pra2)
> 
> # Uji autokorelasi
> Box.test(resid.modelpra1, lag = 20, type = "Ljung-Box")

    Box-Ljung test

data:  resid.modelpra1
X-squared = 4.3475, df = 20, p-value = 0.9999
> Box.test(resid.modelpra2, lag = 20, type = "Ljung-Box")

    Box-Ljung test

data:  resid.modelpra2
X-squared = 0.83876, df = 20, p-value = 1
> 
> # Uji normalitas
> jarque.bera.test(resid.modelpra1)

    Jarque Bera Test

data:  resid.modelpra1
X-squared = 13.372, df = 2, p-value = 0.001248
> jarque.bera.test(resid.modelpra2)

    Jarque Bera Test

data:  resid.modelpra2
X-squared = 9.9704, df = 2, p-value = 0.006839
> 
> # Perbandingan AIC
> cat("AIC Pra Model 1:", model.pra1$aic, "\n")
AIC Pra Model 1: 1637.434 
> cat("AIC Pra Model 2:", model.pra2$aic, "\n")
AIC Pra Model 2: 1629.818

7.4 Asumsi Model ARIMA Intervensi

> resid.model.inter <- residuals(model.intervensi)
> 
> # Uji autokorelasi
> Box.test(resid.model.inter, lag = 10, type = "Ljung-Box")

    Box-Ljung test

data:  resid.model.inter
X-squared = 0.44243, df = 10, p-value = 1
> 
> # Uji normalitas
> jarque.bera.test(resid.model.inter)

    Jarque Bera Test

data:  resid.model.inter
X-squared = 144.38, df = 2, p-value < 2.2e-16

8 Pemilihan Model Terbaik

Model terbaik dipilih berdasarkan nilai AIC dan BIC terkecil, dengan tetap mempertimbangkan prinsip parsimoni.

\[AIC = -2\ln L + 2k\] \[BIC = -2\ln L + k\ln n\]

8.1 Tabel Perbandingan AIC dan BIC

> evaluasi.model <- data.frame(
+   Model = c("ARIMA", "ARIMA + Outlier", "ARIMA Intervensi"),
+   AIC   = round(c(AIC(arima.model1), AIC(arima.tc), AIC(model.intervensi)), 3),
+   BIC   = round(c(BIC(arima.model1), BIC(arima.tc), BIC(model.intervensi)), 3)
+ )
> 
> kable(evaluasi.model,
+       caption = "Perbandingan Nilai AIC dan BIC Antar Model",
+       align   = c("l", "r", "r"))

Perbandingan Nilai AIC dan BIC Antar Model
Model	AIC	BIC
ARIMA	3541.839	3556.786
ARIMA + Outlier	3538.590	3564.746
ARIMA Intervensi	3541.359	3623.563

8.2 Plot Perbandingan Fitted Values

> prediksi.arima       <- H.saham - resid.arima1
> prediksi.arimatc     <- H.saham - resid.arimatc
> prediksi.arima.inter <- H.saham - resid.model.inter
> 
> plot(tanggal, H.saham,
+      type = "l", col = "blue", lwd = 2,
+      xlab = "Tanggal", ylab = "Harga Saham",
+      main = "Perbandingan Fitted Values Antar Model")
> lines(tanggal, prediksi.arima,       col = "red",    lwd = 1.5)
> lines(tanggal, prediksi.arimatc,     col = "green",  lwd = 1.5)
> lines(tanggal, prediksi.arima.inter, col = "orange", lwd = 1.5)
> legend("topright",
+        legend = c("Aktual", "ARIMA", "ARIMA + Outlier", "ARIMA Intervensi"),
+        col    = c("blue", "red", "green", "orange"),
+        lty    = 1, lwd = 2, bty = "n")

9 Evaluasi Prediksi

Evaluasi prediksi menggunakan ukuran akurasi in-sample untuk membandingkan performa ketiga model.

\[RMSE = \sqrt{\frac{1}{n}\sum_{t=1}^{n}(Z_t - \hat{Z}_t)^2}\]

\[MAE = \frac{1}{n}\sum_{t=1}^{n}|Z_t - \hat{Z}_t|\]

\[MAPE = \frac{100\%}{n}\sum_{t=1}^{n}\left|\frac{Z_t - \hat{Z}_t}{Z_t}\right|\]

> evaluasi.prediksi <- data.frame(
+   Model = c("ARIMA", "ARIMA + Outlier", "ARIMA Intervensi"),
+   RMSE  = round(c(
+     sqrt(mean(residuals(arima.model1)^2)),
+     sqrt(mean(residuals(arima.tc)^2)),
+     sqrt(mean(residuals(model.intervensi)^2))
+   ), 4),
+   MAE   = round(c(
+     mean(abs(residuals(arima.model1))),
+     mean(abs(residuals(arima.tc))),
+     mean(abs(residuals(model.intervensi)))
+   ), 4),
+   MAPE  = round(c(
+     mean(abs(residuals(arima.model1) / H.saham)) * 100,
+     mean(abs(residuals(arima.tc)     / H.saham)) * 100,
+     mean(abs(residuals(model.intervensi) / H.saham)) * 100
+   ), 4)
+ )
> 
> kable(evaluasi.prediksi,
+       caption = "Evaluasi Akurasi Prediksi In-Sample",
+       align   = c("l", "r", "r", "r"))

Evaluasi Akurasi Prediksi In-Sample
Model	RMSE	MAE	MAPE
ARIMA	72.0563	53.3465	0.7561
ARIMA + Outlier	71.0273	52.1974	0.7385
ARIMA Intervensi	67.8496	50.8631	0.7181

10 Peramalan

Peramalan dilakukan menggunakan model ARIMA + Outlier terpilih dengan horizon 30 periode ke depan. Model di-refit menggunakan fungsi Arima() dari package forecast agar kompatibel dengan fungsi forecast().

Ramalan titik \(h\)-langkah ke depan:

\[\hat{Z}_t(h) = E(Z_{t+h} \mid Z_t, Z_{t-1}, \ldots) = \sum_{k=1}^{K} \frac{\hat{\omega}_k(B)}{\hat{\delta}_k(B)} I_{k,t+h} + \hat{Z}_t^N(h)\]

Selang kepercayaan \((1-\alpha) \times 100\%\):

\[\hat{Z}_t(h) \pm z_{\alpha/2} \cdot \hat{\sigma}_a \cdot \sqrt{\sum_{j=0}^{h-1} \psi_j^2}\]

> # ==============================================
> # STEP 1: Refit model dengan Arima() dari forecast
> # ==============================================
> p <- 19; d <- 1; q <- 0
> 
> # Fixed parameter: NA = diestimasi, 0 = dikunci
> # Urutan: ar1..ar19 (19), tc1, ao1, ao2 (3) = 22 total
> fixed_arima <- c(
+   0, NA, rep(0, 7),   # ar1 (fixed), ar2 (estimasi), ar3-ar9 (fixed)
+   NA,                  # ar10 (estimasi)
+   rep(0, 8),           # ar11-ar18 (fixed)
+   NA,                  # ar19 (estimasi)
+   NA, NA, NA           # tc1, ao1, ao2 (estimasi)
+ )
> 
> xreg_train <- cbind(tc1 = tc1, ao1 = ao1, ao2 = ao2)
> 
> arima.tc2 <- Arima(H.saham,
+                    order          = c(p, d, q),
+                    xreg           = xreg_train,
+                    include.mean   = FALSE,
+                    fixed          = fixed_arima,
+                    transform.pars = FALSE,
+                    method         = "ML")
> 
> # Verifikasi koefisien
> summary(arima.tc2)
Series: H.saham 
Regression with ARIMA(19,1,0) errors 

Coefficients:
      ar1      ar2  ar3  ar4  ar5  ar6  ar7  ar8  ar9    ar10  ar11  ar12  ar13
        0  -0.1523    0    0    0    0    0    0    0  0.1352     0     0     0
s.e.    0   0.0609    0    0    0    0    0    0    0  0.0649     0     0     0
      ar14  ar15  ar16  ar17  ar18     ar19      tc1       ao1       ao2
         0     0     0     0     0  -0.1601  16.6978  153.7717  -14.1055
s.e.     0     0     0     0     0   0.0659  64.5262   51.7529   71.6574

sigma^2 = 5161:  log likelihood = -1762.29
AIC=3538.59   AICc=3538.96   BIC=3564.75

Training set error measures:
                    ME     RMSE      MAE         MPE     MAPE      MASE
Training set -4.055064 71.02731 52.19736 -0.06805568 0.738469 0.9635303
                   ACF1
Training set 0.06373397
> 
> # ==============================================
> # STEP 2: Buat xreg untuk periode ramalan
> # ==============================================
> h          <- 30
> last_index <- length(tc1)
> 
> # TC: lanjutkan decay eksponensial dari indeks terakhir
> tc1_future <- 0.7^((last_index + 1):(last_index + h) - 151)
> 
> # AO: bernilai 0 di masa depan (outlier hanya di masa lalu)
> ao1_future <- rep(0, h)
> ao2_future <- rep(0, h)
> 
> xreg_future <- cbind(tc1 = tc1_future,
+                      ao1 = ao1_future,
+                      ao2 = ao2_future)
> 
> # Verifikasi nama kolom harus identik
> cat("Kolom xreg_train  :", colnames(xreg_train),  "\n")
Kolom xreg_train  : tc1 ao1 ao2 
> cat("Kolom xreg_future :", colnames(xreg_future), "\n")
Kolom xreg_future : tc1 ao1 ao2 
> cat("Dimensi xreg_future:", nrow(xreg_future), "x", ncol(xreg_future), "\n")
Dimensi xreg_future: 30 x 3 
> 
> # ==============================================
> # STEP 3: Forecast
> # ==============================================
> forecast_arima_tc <- forecast(arima.tc2,
+                               h    = h,
+                               xreg = xreg_future)
> 
> # ==============================================
> # STEP 4: Visualisasi
> # ==============================================
> plot(forecast_arima_tc,
+      main = "Peramalan IHSG — Model ARIMA + Outlier",
+      xlab = "Waktu",
+      ylab = "Harga Saham",
+      col  = "blue",
+      lwd  = 2,
+      fcol = "red",
+      flwd = 2)
> legend("topright",
+        legend = c("Data Aktual", "Ramalan", "CI 80%", "CI 95%"),
+        col    = c("blue", "red", "gray70", "gray90"),
+        lty    = c(1, 1, 2, 2),
+        lwd    = c(2, 2, 1, 1),
+        bty    = "n")

> 
> # ==============================================
> # STEP 5: Tabel hasil ramalan
> # ==============================================
> hasil_ramalan <- data.frame(
+   Periode = 1:h,
+   Ramalan = round(as.numeric(forecast_arima_tc$mean),      2),
+   CI80_Lo = round(as.numeric(forecast_arima_tc$lower[, 1]), 2),
+   CI80_Hi = round(as.numeric(forecast_arima_tc$upper[, 1]), 2),
+   CI95_Lo = round(as.numeric(forecast_arima_tc$lower[, 2]), 2),
+   CI95_Hi = round(as.numeric(forecast_arima_tc$upper[, 2]), 2)
+ )
> 
> kable(hasil_ramalan,
+       caption = "Hasil Peramalan IHSG 30 Periode ke Depan",
+       align   = c("c", "r", "r", "r", "r", "r"))

Hasil Peramalan IHSG 30 Periode ke Depan
Periode	Ramalan	CI80_Lo	CI80_Hi	CI95_Lo	CI95_Hi
1	6002.67	5910.60	6094.73	5861.86	6143.47
2	6002.90	5872.70	6133.11	5803.78	6202.03
3	6006.33	5854.53	6158.14	5774.17	6238.50
4	5995.65	5824.96	6166.34	5734.60	6256.70
5	5990.39	5801.81	6178.98	5701.98	6278.81
6	5983.00	5778.08	6187.93	5669.60	6296.41
7	6018.64	5798.71	6238.57	5682.28	6355.00
8	6046.03	5812.05	6280.01	5688.19	6403.88
9	5978.03	5730.78	6225.27	5599.89	6356.16
10	6011.74	5751.90	6271.58	5614.35	6409.13
11	6010.75	5734.98	6286.51	5589.00	6432.49
12	5994.44	5703.63	6285.26	5549.68	6439.20
13	6014.82	5710.82	6318.82	5549.89	6479.75
14	6031.38	5714.74	6348.01	5547.13	6515.63
15	6015.66	5686.63	6344.68	5512.46	6518.85
16	5974.24	5633.29	6315.20	5452.80	6495.69
17	5975.33	5622.89	6327.78	5436.32	6514.35
18	6067.69	5704.12	6431.26	5511.66	6623.72
19	6060.58	5686.21	6434.94	5488.03	6633.12
20	6047.78	5666.06	6429.50	5463.99	6631.57
21	6048.69	5659.43	6437.94	5453.37	6644.00
22	6047.88	5650.35	6445.42	5439.91	6655.86
23	6052.21	5646.72	6457.70	5432.06	6672.36
24	6055.41	5642.31	6468.51	5423.63	6687.20
25	6053.81	5633.20	6474.43	5410.53	6697.09
26	6042.02	5613.98	6470.06	5387.39	6696.64
27	6038.03	5602.71	6473.35	5372.26	6703.79
28	6063.20	5620.72	6505.67	5386.49	6739.91
29	6057.45	5607.92	6506.97	5369.96	6744.93
30	6052.04	5596.26	6507.83	5354.98	6749.10

Berdasarkan grafik hasil peramalan menggunakan model ARIMA dengan outlier, nilai IHSG pada periode mendatang diprediksi cenderung bergerak stabil di kisaran 6.000-an tanpa menunjukkan tren kenaikan maupun penurunan yang signifikan. Setelah mengalami penurunan tajam pada akhir periode historis, model memperkirakan kondisi pasar akan memasuki fase yang relatif tenang. Namun demikian, interval kepercayaan yang semakin melebar seiring bertambahnya horizon peramalan menunjukkan bahwa ketidakpastian prediksi meningkat, sehingga hasil ramalan jangka panjang perlu diinterpretasikan dengan hati-hati.

Perbadingan Model ARIMA, ARIMA Outlier dan ARIMA Intervensi pada Peramalan Data IHSG

Ahmad Fathan Muzhaffar

2026-04-29