—- LANGKAH 0. Paket & pengaturan ———————————–

install.packages(c(“rugarch”, “forecast”, “tseries”, “xts”, “readxl”)) # sekali saja

library(rugarch) # model ARCH/GARCH & variasinya library(forecast) # auto.arima untuk persamaan rata-rata library(tseries) # uji ADF & KPSS library(xts) # format data runtun waktu library(readxl) # membaca file Excel

FOLDER <- file.path(Sys.getenv(ifelse(.Platform$OS.type == “windows”, “USERPROFILE”, “HOME”)), “Downloads”) # <– GANTI jika file di folder lain FILE_DATA <- file.path(FOLDER, “data uas no 1 arw.xlsx”) # <– GANTI jika nama/ekstensi beda if (!file.exists(FILE_DATA)) { # file tidak ketemu -> pilih manual message(“File tidak ditemukan di:”, FILE_DATA, “pilih file secara manual.”) FILE_DATA <- file.choose() } SHEET <- 1 # sheet ke-1 (“Sheet1”) KODE_SAHAM <- “PTBA” # emiten: PT Bukit Asam Tbk KOL_TANGGAL <- “Date” KOL_HARGA <- “Close”

DIST <- “std” # distribusi galat: “std” = Student-t (cocok untuk fat tail), # “norm” = normal. Coba keduanya, bandingkan AIC/BIC. HORIZON <- 30 # jumlah hari ke depan untuk peramalan volatilitas HARI_THN <- 252 # hari bursa per tahun (untuk anualisasi volatilitas)

Kosongkan (NULL) agar model dipilih otomatis, atau isi nama model, mis. “GARCH(1,1)”

MODEL_PILIHAN <- NULL # <– GANTI jika ingin memaksa model tertentu

— Fungsi bantu —————————————————–

statdes <- function(x) { # statistik deskriptif ala EViews x <- as.numeric(x); n <- length(x); m <- mean(x); s0 <- sqrt(mean((x - m)^2)) S <- mean((x - m)^3) / s0^3; K <- mean((x - m)^4) / s0^4 JB <- n / 6 * (S^2 + (K - 3)^2 / 4) c(Mean = m, Median = median(x), Maximum = max(x), Minimum = min(x), Std.Dev = sd(x), Skewness = S, Kurtosis = K, Jarque.Bera = JB, Probability = pchisq(JB, 2, lower.tail = FALSE), Observations = n) } uji_arch_lm <- function(e, lags = c(1, 5, 10)) { # ARCH-LM: N x R^2 ~ Chi2(lag) e2 <- as.numeric(e)^2 do.call(rbind, lapply(lags, function(lag) { X <- embed(e2, lag + 1) NR2 <- nrow(X) * summary(lm(X[, 1] ~ X[, -1, drop = FALSE]))$r.squared data.frame(Lag = lag, ARCH_LM = round(NR2, 3), Prob = round(pchisq(NR2, lag, lower.tail = FALSE), 4)) })) } uji_lb <- function(x, lags = c(10, 20)) { # Ljung-Box pada beberapa lag do.call(rbind, lapply(lags, function(l) { t <- Box.test(as.numeric(x), lag = l, type = "Ljung-Box") data.frame(Lag = l, Q = round(unname(t$statistic), 3), Prob = round(t$p.value, 4)) })) } bintang <- function(p) ifelse(p < 0.001, “”, ifelse(p < 0.01, ””, ifelse(p < 0.05, ””, ifelse(p < 0.1, “.”, ““))))

—- LANGKAH 1. Baca data harga penutupan —————————-

Baca mentah (tanpa nama kolom), cari sel header “Date” lalu ambil data di bawahnya

mentah <- suppressMessages(read_excel(FILE_DATA, sheet = SHEET, col_names = FALSE, col_types = “text”)) pos <- which(as.matrix(mentah) == KOL_TANGGAL, arr.ind = TRUE) if (nrow(pos) == 0) stop(“Header ‘“, KOL_TANGGAL,”’ tidak ditemukan. Cek nama kolom / sheet.”) br <- pos[1, “row”]; kt <- pos[1, “col”] kh <- which(as.character(unlist(mentah[br, ])) == KOL_HARGA)[1] if (is.na(kh)) stop(“Header ‘“, KOL_HARGA,”’ tidak ditemukan di baris yang sama.”) dat <- data.frame(Tanggal = as.character(unlist(mentah[-(1:br), kt])), Harga = as.character(unlist(mentah[-(1:br), kh])), stringsAsFactors = FALSE) head(dat) # CEK: kolom tanggal & harga harus terlihat benar (mis. 2021-10-01 / 2720)

Ubah tanggal ke format Date. Mendukung:

- angka serial Excel (mis. 44470.6667)

- teks “dd/mm/yyyy hh:mm:ss” (format Indonesia, seperti di Google Sheets Anda)

- teks “yyyy-mm-dd hh:mm:ss”

t_raw <- trimws(dat$Tanggal) tgl <- rep(as.Date(NA), length(t_raw)) angka <- suppressWarnings(as.numeric(t_raw)) i_ang <- !is.na(angka) & !grepl("[/-]", t_raw) tgl[i_ang] <- as.Date(floor(angka[i_ang]), origin = "1899-12-30") i_dmy <- grepl("^[0-9]{1,2}/[0-9]{1,2}/[0-9]{4}", t_raw) tgl[i_dmy] <- as.Date(substr(t_raw[i_dmy], 1, 10), format = "%d/%m/%Y") i_iso <- grepl("^[0-9]{4}-[0-9]{2}-[0-9]{2}", t_raw) tgl[i_iso] <- as.Date(substr(t_raw[i_iso], 1, 10), format = "%Y-%m-%d") hrg <- suppressWarnings(as.numeric(gsub(",", ".", dat$Harga)))

ok <- !is.na(tgl) & !is.na(hrg) & hrg > 0 harga <- xts(hrg[ok], order.by = tgl[ok]) harga <- harga[!duplicated(index(harga))] # buang tanggal ganda colnames(harga) <- KODE_SAHAM

cat(“Saham:”, KODE_SAHAM, “| Data harian:”, nrow(harga), “hari |”, format(start(harga)), “s.d.”, format(end(harga)), “”) summary(as.numeric(harga))

—- LANGKAH 2. Return log harian & grafik —————————

R_t = ln(P_t / P_t-1) x 100 (return log dalam persen)

r <- na.omit(diff(log(harga))) * 100 colnames(r) <- “Return”

par(mfrow = c(2, 1), mar = c(3, 4, 2, 1)) plot(index(harga), as.numeric(harga), type = “l”, main = paste(“Harga penutupan”, KODE_SAHAM), xlab = ““, ylab =”Harga (Rp)“) plot(index(r), as.numeric(r), type =”l”, main = paste(“Return harian”, KODE_SAHAM, “(%)”), xlab = ““, ylab =”Return (%)“) abline(h = 0, col =”grey50”) par(mfrow = c(1, 1))

—- LANGKAH 3. Statistik deskriptif & uji stasioneritas ————-

print(round(as.matrix(statdes(r)), 4))

Uji ADF (H0: ada akar unit / tidak stasioner) dan KPSS (H0: stasioner)

adf_p <- suppressWarnings(adf.test(as.numeric(log(harga)))) adf_r <- suppressWarnings(adf.test(as.numeric(r))) kps_r <- suppressWarnings(kpss.test(as.numeric(r), null = “Level”)) cat(“log-harga : p =”, round(adf_p$p.value, 4), ifelse(adf_p$p.value < 0.05, “-> stasioner”, “-> TIDAK stasioner”), “”) cat(“ADF return : p =”, round(adf_r$p.value, 4), ifelse(adf_r$p.value < 0.05, “-> stasioner”, “-> TIDAK stasioner”), “”) cat(“KPSS return : p =”, round(kps_r$p.value, 4), ifelse(kps_r$p.value > 0.05, “-> stasioner”, “-> TIDAK stasioner”), “”) # BACA: yang dimodelkan adalah RETURN, jadi return harus stasioner # (ADF p < 0,05 dan KPSS p > 0,05). Harga/log-harga umumnya tidak stasioner.

—- LANGKAH 4. Persamaan rata-rata & UJI EFEK ARCH ——————

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1)) Acf(as.numeric(r), main = “ACF return”, lag.max = 30) Acf(as.numeric(r)^2, main = “ACF return kuadrat”, lag.max = 30) par(mfrow = c(1, 1)) # BACA: ACF return biasanya kecil (hampir tidak ada autokorelasi), tetapi ACF # return KUADRAT nyata -> varians bergantung pada masa lalu = efek ARCH.

am <- auto.arima(as.numeric(r), max.p = 3, max.q = 3, max.d = 0, seasonal = FALSE, stationary = TRUE, ic = “bic”, stepwise = FALSE, approximation = FALSE) ord <- arimaorder(am); ORDE_ARMA <- unname(c(ord[“p”], ord[“q”])) cat(“rata-rata terpilih: ARMA(”, ORDE_ARMA[1], “,”, ORDE_ARMA[2], “)”) # ORDE_ARMA <- c(0, 0) # <– (opsional) aktifkan untuk memaksa hanya konstanta

e_mean <- residuals(am) cat(“ARCH-LM pada residual persamaan rata-rata:”); print(uji_arch_lm(e_mean)) cat(“-Box pada residual KUADRAT:”); print(uji_lb(e_mean^2))

—- LANGKAH 5. Estimasi banyak model & pemilihan model terbaik ——

daftar_model <- list( “ARCH(1)” = list(vm = “sGARCH”, ord = c(1, 0)), “ARCH(2)” = list(vm = “sGARCH”, ord = c(2, 0)), “ARCH(5)” = list(vm = “sGARCH”, ord = c(5, 0)), “GARCH(1,1)” = list(vm = “sGARCH”, ord = c(1, 1)), “GARCH(1,2)” = list(vm = “sGARCH”, ord = c(1, 2)), “GARCH(2,1)” = list(vm = “sGARCH”, ord = c(2, 1)), “EGARCH(1,1)” = list(vm = “eGARCH”, ord = c(1, 1)), “GJR-GARCH(1,1)” = list(vm = “gjrGARCH”, ord = c(1, 1)), “IGARCH(1,1)” = list(vm = “iGARCH”, ord = c(1, 1)), “GARCH-M(1,1)” = list(vm = “sGARCH”, ord = c(1, 1), archm = TRUE))

buat_spec <- function(m) ugarchspec( mean.model = list(armaOrder = ORDE_ARMA, include.mean = TRUE, archm = isTRUE(m$archm), archpow = 2), variance.model = list(model = m$vm, garchOrder = m$ord), distribution.model = DIST)

r_num <- as.numeric(r) fit <- list() for (nm in names(daftar_model)) { f <- tryCatch(ugarchfit(buat_spec(daftar_model[[nm]]), data = r, solver = “hybrid”), error = function(e) NULL) if (is.null(f) || f@fit$convergence != 0) { f <- tryCatch(ugarchfit(buat_spec(daftar_model[[nm]]), data = r, solver = "gosolnp"), error = function(e) NULL) # coba solver cadangan } if (!is.null(f) && f@fit$convergence == 0) fit[[nm]] <- f else cat(“Estimasi gagal konvergen (dilewati):”, nm, “”) }

Tabel perbandingan model

tab_model <- do.call(rbind, lapply(names(fit), function(nm) { f <- fit[[nm]]; ic <- infocriteria(f) z <- as.numeric(residuals(f, standardize = TRUE)) data.frame(Model = nm, LogLik = round(likelihood(f), 2), AIC = round(ic[“Akaike”, 1], 4), BIC = round(ic[“Bayes”, 1], 4), Persistensi = round(persistence(f), 4), p_LB_z2 = round(Box.test(z^2, lag = 10, type = “Ljung-Box”)$p.value, 4), p_ARCHLM_z = round(uji_arch_lm(z, lags = 10)$Prob, 4)) })) tab_model <- tab_model[order(tab_model$BIC), ] print(tab_model, row.names = FALSE)

pers <- persistence(fit_b) hl <- tryCatch(halflife(fit_b), error = function(e) NA) vol_tak_bersyarat <- tryCatch(sqrt(uncvariance(fit_b)), error = function(e) NA) cat(“:”, round(pers, 4), “”) cat(“Half-life shock (hari) :”, round(hl, 2), ” -> waktu agar dampak shock tinggal separuh“) cat(”Volatilitas jangka panjang :“, round(vol_tak_bersyarat, 4),”% per hari (“, round(vol_tak_bersyarat * sqrt(HARI_THN), 2),”% per tahun )“)

—- LANGKAH 7. Diagnostik residual standar ————————–

z <- as.numeric(residuals(fit_b, standardize = TRUE)) cat(“-Box residual standar (z):”); print(uji_lb(z)) cat(“-Box residual standar KUADRAT (z^2):”); print(uji_lb(z^2)) cat(“-LM residual standar:”); print(uji_arch_lm(z)) cat(“-Bera residual standar: p =”, round(statdes(z)[“Probability”], 4), “| Kurtosis =”, round(statdes(z)[“Kurtosis”], 3), “”) cat(“Sign Bias (asimetri sisa):”); print(signbias(fit_b)) # BACA: model memadai bila # - Ljung-Box z dan z^2 : Prob > 0,05 (tidak tersisa autokorelasi / efek ARCH) # - ARCH-LM z : Prob > 0,05 (efek ARCH sudah tertangkap) # - Sign Bias Joint Effect: Prob > 0,05 (tidak ada asimetri yang terlewat). # Jika < 0,05, coba EGARCH atau GJR-GARCH. # - Jarque-Bera z boleh tetap signifikan; itu sebabnya dipakai Student-t.

—- LANGKAH 8. Volatilitas bersyarat & news impact curve ————

sig <- sigma(fit_b) # simpangan baku bersyarat (% per hari) par(mfrow = c(2, 1), mar = c(3, 4, 2, 1)) plot(index(r), abs(as.numeric(r)), type = “h”, col = “grey70”, main = paste(“Volatilitas bersyarat”, KODE_SAHAM, “-”, MODEL_PILIHAN), xlab = ““, ylab =”% per hari”) lines(index(sig), as.numeric(sig), col = “firebrick”, lwd = 1.5) legend(“topright”, c(“|Return|”, “Sigma bersyarat”), col = c(“grey70”, “firebrick”), lty = 1, bty = “n”, cex = 0.8) ni <- newsimpact(fit_b) plot(ni$zx, ni$zy, type = “l”, lwd = 2, xlab = “Shock (z)”, ylab = “Varians bersyarat”, main = “News impact curve”) abline(v = 0, col = “grey60”, lty = 2) par(mfrow = c(1, 1))

Lima hari paling bergejolak

top <- head(sig[order(-as.numeric(sig))], 5) cat(“hari dengan volatilitas bersyarat tertinggi:”) print(data.frame(Tanggal = format(index(top)), Sigma = round(as.numeric(top), 4)), row.names = FALSE) # TIPS: cocokkan tanggal-tanggal ini dengan peristiwa (COVID-19 2020, lonjakan/ # anjloknya harga batu bara/minyak, kebijakan DMO, dsb.) untuk pembahasan.

—- LANGKAH 9. Peramalan volatilitas ke depan ———————–

fc <- ugarchforecast(fit_b, n.ahead = HORIZON) sig_fc <- as.numeric(sigma(fc)) tab_fc <- data.frame(Hari_ke = 1:HORIZON, Sigma_harian = round(sig_fc, 4), Sigma_tahunan = round(sig_fc * sqrt(HARI_THN), 2)) print(head(tab_fc, 10), row.names = FALSE) plot(1:HORIZON, sig_fc, type = “b”, pch = 16, xlab = “Hari ke depan”, ylab = “Sigma (% per hari)”, main = paste(“Ramalan volatilitas”, KODE_SAHAM, “-”, MODEL_PILIHAN)) abline(h = vol_tak_bersyarat, col = “grey50”, lty = 2) # BACA: ramalan volatilitas bergerak menuju volatilitas jangka panjang (garis # putus-putus). Makin tinggi persistensi, makin lambat konvergen. # (IGARCH tidak konvergen ke rata-rata karena shock permanen.)

—- LANGKAH 10. Ringkasan interpretasi otomatis ———————

p_ <- function(nama) if (nama %in% rownames(mc)) mc[nama, “Prob”] else NA k_ <- function(nama) if (nama %in% rownames(mc)) mc[nama, “Koefisien”] else NA cat(“================ RINGKASAN INTERPRETASI ================”) cat(“Saham :”, KODE_SAHAM, “| Sampel:”, format(start(r)), “s.d.”, format(end(r)), “(”, length(r), “return harian )”) cat(“Efek ARCH :”, ifelse(uji_arch_lm(e_mean, 5)\(Prob < 0.05, "ADA (varians berubah menurut waktu) -> ARCH/GARCH layak dipakai", "TIDAK terdeteksi"), "\n", sep = "") cat("Model terpilih :", MODEL_PILIHAN, "\n") if (!is.na(p_("alpha1"))) cat("alpha1 (ARCH) :", round(k_("alpha1"), 4), bintang(p_("alpha1")), "-> respon varians terhadap shock kemarin", ifelse(p_("alpha1") < 0.05, "(signifikan)", "(tidak signifikan)"), "\n") if (!is.na(p_("beta1"))) cat("beta1 (GARCH) :", round(k_("beta1"), 4), bintang(p_("beta1")), "-> pewarisan volatilitas kemarin", ifelse(p_("beta1") < 0.05, "(signifikan)", "(tidak signifikan)"), "\n") cat("Persistensi :", round(pers, 4), ifelse(pers >= 0.95, "-> sangat persisten (gejolak lama mereda)", ifelse(pers >= 0.8, "-> persisten sedang", "-> cepat mereda")), "\n") if (!is.na(hl) && is.finite(hl)) cat("Half-life :", round(hl, 1), "hari\n") if (!is.na(p_("gamma1"))) { lev <- if (grepl("EGARCH", MODEL_PILIHAN)) k_("alpha1") < 0 && p_("alpha1") < 0.05 else k_("gamma1") > 0 && p_("gamma1") < 0.05 cat("Leverage effect :", ifelse(lev, "ADA (kabar buruk > kabar baik)", "tidak terbukti"), "\n") } if (!is.na(p_("archm"))) cat("Premi risiko : archm =", round(k_("archm"), 4), bintang(p_("archm")), ifelse(p_("archm") < 0.05, "(signifikan)", "(tidak signifikan)"), "\n") cat("Diagnostik : p LB z^2 =", uji_lb(z^2, 10)\)Prob, “| p ARCH-LM z =”, uji_arch_lm(z, 10)$Prob, ifelse(uji_lb(z^2, 10)$Prob > 0.05 && uji_arch_lm(z, 10)$Prob > 0.05, “-> model memadai”, “-> masih ada sisa efek ARCH, coba model lain”), “”) cat(“========================================================”)

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

Script R UAS Rifqy Aunur No. 1

2026-09-26