LOAD PACKAGE
library(forecast)## Warning: package 'forecast' was built under R version 4.5.3library(tseries)## Warning: package 'tseries' was built under R version 4.5.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(readxl)## Warning: package 'readxl' was built under R version 4.5.3library(TSA)## Warning: package 'TSA' was built under R version 4.5.3## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarIMPORT DATA EXCEL
data <- read_excel(
"~/PROJECT SEMESTER 6/PROJECT STATISTIKA 1/DATA PROJECT UAS JST.xlsx"
)
head(data)## # A tibble: 6 × 2
## Tanggal Terakhir
## <chr> <chr>
## 1 30/12/2025 8.646,94
## 2 29/12/2025 8.644,26
## 3 24/12/2025 8.537,91
## 4 23/12/2025 8.584,78
## 5 22/12/2025 8.645,84
## 6 19/12/2025 8.609,55Kode tersebut digunakan untuk membaca data historis Jakarta Stock Exchange Composite (JST) yang tersimpan dalam format Microsoft Excel menggunakan fungsi read_excel() dari package readxl. Data yang berhasil dibaca kemudian disimpan dalam objek data.
MELIHAT STRUKTUR DATA
str(data)## tibble [474 × 2] (S3: tbl_df/tbl/data.frame)
## $ Tanggal : chr [1:474] "30/12/2025" "29/12/2025" "24/12/2025" "23/12/2025" ...
## $ Terakhir: chr [1:474] "8.646,94" "8.644,26" "8.537,91" "8.584,78" ...summary(data)## Tanggal Terakhir
## Length:474 Length:474
## Class :character Class :character
## Mode :character Mode :characterDATA CLEANING
Data harga penutupan pada file Excel masih menggunakan format numerik Indonesia, yaitu titik sebagai pemisah ribuan dan koma sebagai pemisah desimal. Oleh karena itu dilakukan transformasi agar dapat diproses oleh perangkat lunak R.
Misalnya:
8.710,70
diubah menjadi:
8710.70
sehingga dapat dikenali sebagai bilangan real.
Langkah ini penting karena seluruh perhitungan statistik dan model ARIMA membutuhkan data numerik.
# Mengubah format angka Indonesia menjadi numerik
data$Terakhir <- as.numeric(
gsub(",", ".",
gsub("\\.", "", data$Terakhir))
)
# Mengubah tanggal menjadi format Date
data$Tanggal <- as.Date(
data$Tanggal,
format = "%d/%m/%Y"
)
# Mengurutkan data berdasarkan tanggal
data <- data[order(data$Tanggal), ]
# Mengecek missing value
sum(is.na(data$Terakhir))## [1] 0head(data)## # A tibble: 6 × 2
## Tanggal Terakhir
## <date> <dbl>
## 1 2024-01-15 7224
## 2 2024-01-16 7243.
## 3 2024-01-17 7201.
## 4 2024-01-18 7253.
## 5 2024-01-19 7227.
## 6 2024-01-22 7248.MEMBUAT TIME SERIES
Kode ini mengubah data harga penutupan menjadi objek deret waktu (time series object),dengan urutan observasi berdasarkan waktu.
Transformasi ini memungkinkan fungsi-fungsi time series seperti ACF, PACF, ARIMA, dan SARIMA diterapkan pada data.
ts_data <- ts(data$Terakhir)
plot(
ts_data,
main = "Plot Data IHSG",
ylab = "Harga Penutupan",
xlab = "Waktu",
col = "blue"
)
# PEMBAGIAN DATA TRAINING DAN TESTING Sebanyak 80% data digunakan
sebagai data pelatihan (training data) dan 20% sisanya digunakan sebagai
data pengujian (testing data).
Jika jumlah observasi:
n = 474
maka:
training = 0.8 x 474 = 379
testing = 474 - 379 = 95
Data training digunakan untuk membangun model, sedangkan data testing digunakan untuk mengukur akurasi forecasting.
n <- length(ts_data)
train_size <- round(0.8*n)
train <- ts_data[1:train_size]
test <- ts_data[(train_size+1):n]UJI STASIONERITAS AWAL
Uji Augmented Dickey-Fuller (ADF) digunakan untuk mengetahui apakah data memiliki akar unit (unit root).
Hipotesis yang digunakan:
H0 = Data tidak stasioner
H1 = Data stasioner
adf.test(train)##
## Augmented Dickey-Fuller Test
##
## data: train
## Dickey-Fuller = -1.5623, Lag order = 7, p-value = 0.7619
## alternative hypothesis: stationaryDIFFERENCING
Differencing dilakukan untuk menghilangkan tren sehingga rata-rata data menjadi konstan.
d <- ndiffs(train)
d ## [1] 1train_diff <- diff(train, differences = d)
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 = -6.1835, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationaryJika data masih belum stasioner maka differencing dapat dilakukan lebih dari satu kali.
Pada penelitian ini diperoleh:
d = 1
sehingga data menjadi stasioner setelah differencing pertama.
PLOT ACF DAN PACF
Tahap identifikasi model dilakukan menggunakan plot Autocorrelation Function (ACF), Partial Autocorrelation Function (PACF), dan Extended Autocorrelation Function (EACF) pada data yang telah stasioner. Ketiga metode ini digunakan untuk menentukan kandidat orde model ARIMA yang sesuai.
acf(
train_diff,
main = "Plot ACF"
)
pacf(
train_diff,
main = "Plot PACF"
)
TSA::eacf(train_diff)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 x x o o o o o o o o o o o o
## 3 x x x o o o o o o o o o o o
## 4 x x x o o o o o o o o o o o
## 5 x o x o x o o o o o o o o o
## 6 x o x x x o o o o o o o o o
## 7 x x x o x o o o o o o o o oPEMODELAN ARIMA
Model ARIMA memiliki bentuk umum:
\[ϕ_p(B)(1−B)^dY_t=θ_q(B)ε_t\]
# ARIMA(1,1,0)
m1 <- Arima(train, order = c(1,1,0))
# ARIMA(0,1,1)
m2 <- Arima(train, order = c(0,1,1))
# ARIMA(1,1,1)
m3 <- Arima(train, order = c(1,1,1))
# ==============================
# 2. Tabel Perbandingan Model
# ==============================
perbandingan_arima <- data.frame(
Model = c(
"ARIMA(1,1,0)",
"ARIMA(0,1,1)",
"ARIMA(1,1,1)"
),
AIC = c(AIC(m1), AIC(m2), AIC(m3)),
BIC = c(BIC(m1), BIC(m2), BIC(m3)),
# AICc dari model (aman, tidak error)
AICc = c(m1$aicc, m2$aicc, m3$aicc)
)
# tampilkan tabel
perbandingan_arima## Model AIC BIC AICc
## 1 ARIMA(1,1,0) 4573.825 4581.694 4573.857
## 2 ARIMA(0,1,1) 4573.809 4581.679 4573.841
## 3 ARIMA(1,1,1) 4575.769 4587.574 4575.833# ==============================
# 3. Pilih Model Terbaik (AIC terkecil)
# ==============================
best_index <- which.min(perbandingan_arima$AIC)
best_arima <- list(m1, m2, m3)[[best_index]]
summary(best_arima)## Series: train
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.0291
## s.e. 0.0525
##
## sigma^2 = 10449: log likelihood = -2284.9
## AIC=4573.81 AICc=4573.84 BIC=4581.68
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.3707987 101.9484 67.32092 -0.01615881 0.9330212 1.00034
## ACF1
## Training set 0.0006284642# ==============================
# 4. Urutkan Model (opsional)
# ==============================
perbandingan_arima[order(perbandingan_arima$AIC), ]## Model AIC BIC AICc
## 2 ARIMA(0,1,1) 4573.809 4581.679 4573.841
## 1 ARIMA(1,1,0) 4573.825 4581.694 4573.857
## 3 ARIMA(1,1,1) 4575.769 4587.574 4575.833Hasil perbandingan menunjukkan bahwa model ARIMA(0,1,1) memiliki nilai AIC, BIC, dan AICc terkecil dibandingkan model lainnya, sehingga dipilih sebagai model terbaik.
maka:
\[p = 0, d = 1, q = 1\]
Sehingga model menjadi:
\[(1-B)Y_t=(1+θ_1 B)ε_t\] \[(1-B)Yt=Yt-Y_(t-1)\] \[Y_t-Y_{t-1}=ε_t+θ_1 ε_{t-1}\] Pindahkan Y_(t-1) ke ruas kanan:
\[Y_t=Y_{t-1}+ε_t+θ_1 ε_{t-1}\] Dari hasil estimasi diperoleh:
\[𝜃_1=−0,0161\]
Sehingga model menjadi:
\[𝑌_𝑡=𝑌_{𝑡−1}+𝜀_𝑡−0,0161𝜀_{𝑡−1}\]
PEMODELAN SARIMA
Model SARIMA memiliki bentuk umum:
\[ϕ_p (B)Φ_P (B^s)(1-B)^d (1-B^s )^D Y_t=θ_q (B)Θ_Q (B^s)ε_t\] Karena nilai 𝜃_1 sangat kecil (-0.0161), maka kontribusi komponen MA(1) relatif lemah, namun tetap signifikan dalam membentuk model terbaik berdasarkan kriteria AIC, BIC, dan AICc.
# ==============================
# 1. SARIMA(1,1,0)(1,1,0)[12]
# ==============================
m1 <- Arima(
train,
order = c(1,1,0),
seasonal = list(order = c(1,1,0), period = 12)
)
# ==============================
# 2. SARIMA(0,1,1)(0,1,1)[12]
# ==============================
m2 <- Arima(
train,
order = c(0,1,1),
seasonal = list(order = c(0,1,1), period = 12)
)
# ==============================
# 3. SARIMA(1,1,1)(1,1,0)[12]
# ==============================
m3 <- Arima(
train,
order = c(1,1,1),
seasonal = list(order = c(1,1,0), period = 12)
)
# ==============================
# 4. SARIMA(0,1,1)(1,1,0)[12]
# ==============================
m4 <- Arima(
train,
order = c(0,1,1),
seasonal = list(order = c(1,1,0), period = 12)
)
# ==============================
# 5. SARIMA(1,1,0)(0,1,1)[12]
# ==============================
m5 <- Arima(
train,
order = c(1,1,0),
seasonal = list(order = c(0,1,1), period = 12)
)
perbandingan_sarima <- data.frame(
Model = c(
"SARIMA(1,1,0)(1,1,0)[12]",
"SARIMA(0,1,1)(0,1,1)[12]",
"SARIMA(1,1,1)(1,1,0)[12]",
"SARIMA(0,1,1)(1,1,0)[12]",
"SARIMA(1,1,0)(0,1,1)[12]"
),
AIC = c(
AIC(m1), AIC(m2), AIC(m3), AIC(m4), AIC(m5)
),
BIC = c(
BIC(m1), BIC(m2), BIC(m3), BIC(m4), BIC(m5)
),
AICc = c(
AICc = c(m1$aicc, m2$aicc, m3$aicc, m4$aicc, m5$aicc)
)
)
# tampilkan hasil
perbandingan_sarima## Model AIC BIC AICc
## AICc1 SARIMA(1,1,0)(1,1,0)[12] 4577.344 4589.052 4577.410
## AICc2 SARIMA(0,1,1)(0,1,1)[12] 4475.285 4486.992 4475.351
## AICc3 SARIMA(1,1,1)(1,1,0)[12] 4579.343 4594.953 4579.453
## AICc4 SARIMA(0,1,1)(1,1,0)[12] 4577.341 4589.049 4577.408
## AICc5 SARIMA(1,1,0)(0,1,1)[12] 4475.290 4486.998 4475.356Model terbaik yang diperoleh adalah:
SARIMA(0,1,1)(0,1,1)[12]
sehingga:
\[p = 0, d = 1, q = 1\] \[P = 0, D = 1, Q = 1\] Karena tidak terdapat komponen AR maupun Seasonal AR, maka:
\[ϕ(B)=1\]
dan
\[Φ(B^{12})=1\] Komponen MA menjadi:
\[θ(B)=1+θ_1B\] dan Seasonal MA:
\[Θ(B^{12})=1+Θ_1B^{12}\] Substitusi ke model umum menghasilkan:
\[(1−B)(1−B^{12})Yt=(1+θ_1B)(1+Θ_1B^{12})ε_t\] \[Y_t-Y_{t-1}-Y_{t-12}+Y_{t-13}= ε_t+θ_1 ε_{t-1}+Θ_1 ε_{t-12}+θ_1 Θ_1 ε_{t-13}\] \[Y_t=Y_{t-1}+Y_{t-12}-Y_{t-13}+ε_t+θ_1 ε_{t-1}+Θ_1 ε_{t-12}+θ_1 Θ_1 ε_{t-13}\] Dari hasil estimasi diperoleh:
\[θ_1=-0,0161\]
\[Θ_1=-1,0000\]
\[θ_1 Θ_1=0,0161\]
Sehingga:
\[Y_t=Y_{t-1}+Y_{t-12}-Y_{t-13}+ε_t-0,0161ε_{t-1}-1,0000ε_{t-12}+0,0161ε_{t-13}\] Model SARIMA menunjukkan bahwa data dipengaruhi oleh pola non-musiman dan musiman secara bersamaan. Komponen musiman pada lag 12 sangat dominan, sehingga model ini mampu menangkap pola berulang setiap periode 12.
DIAGNOSTIC MODEL
Diagnostic model dilakukan untuk memeriksa apakah residual telah memenuhi asumsi white noise.
# Model terbaik ARIMA
best_arima <- Arima(train, order = c(0,1,1))
# Diagnostik residual
checkresiduals(best_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)
## Q* = 8.2263, df = 9, p-value = 0.5115
##
## Model df: 1. Total lags used: 10# Model terbaik SARIMA
best_sarima <- Arima(
train,
order = c(0,1,1),
seasonal = list(order = c(0,1,1), period = 12)
)
# Diagnostik residual
checkresiduals(best_sarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 8.2673, df = 8, p-value = 0.4078
##
## Model df: 2. Total lags used: 10Berdasarkan hasil diagnostic model, plot ACF residual menunjukkan bahwa sebagian besar nilai autokorelasi berada dalam batas kepercayaan 95%, sehingga tidak terdapat autokorelasi yang signifikan.
FORECASTING
Kode ini digunakan untuk menghasilkan prediksi nilai JST pada periode mendatang sebanyak h periode.
Forecast diperoleh berdasarkan parameter model terbaik yang telah dibentuk sebelumnya.
h <- length(test)
forecast_arima <- forecast(
best_arima,
h = h
)
forecast_sarima <- forecast(
best_sarima,
h = h
)EVALUASI MODEL
Semakin kecil nilai MAE, RMSE, dan MAPE, maka semakin baik performa model forecasting yang dihasilkan.
# ==============================
# Akurasi model
# ==============================
akurasi_arima <- accuracy(forecast_arima, test)
akurasi_sarima <- accuracy(forecast_sarima, test)
# ==============================
# Model evaluasi gabungan
# ==============================
hasil_evaluasi <- data.frame(
Model = c(
"ARIMA(0,1,1)",
"SARIMA(0,1,1)(0,1,1)[12]"
),
RMSE = c(
akurasi_arima["Test set", "RMSE"],
akurasi_sarima["Test set", "RMSE"]
),
MAE = c(
akurasi_arima["Test set", "MAE"],
akurasi_sarima["Test set", "MAE"]
),
MAPE = c(
akurasi_arima["Test set", "MAPE"],
akurasi_sarima["Test set", "MAPE"]
),
AIC = c(
AIC(best_arima),
AIC(best_sarima)
),
BIC = c(
BIC(best_arima),
BIC(best_sarima)
),
AICc = c(
best_arima$aicc,
best_sarima$aicc
)
)
# tampilkan hasil
hasil_evaluasi## Model RMSE MAE MAPE AIC BIC
## 1 ARIMA(0,1,1) 303.2581 249.5739 3.360734 4573.809 4581.679
## 2 SARIMA(0,1,1)(0,1,1)[12] 314.9210 264.6136 3.565866 4475.285 4486.992
## AICc
## 1 4573.841
## 2 4475.351Hasil ini menunjukkan bahwa model ARIMA(0,1,1) memiliki tingkat kesalahan prediksi yang lebih rendah sehingga memberikan performa peramalan yang lebih baik dan dipilih sebagai model terbaik untuk forecasting data Jakarta Stock Exchange Composite (JST).