--- title: "Revisi MPDW" author: "Dalilah Husna - G1401221003" date: "2025-04-13" output: html_document --- ## LOAD LIBRARY ```{r} library(ggplot2) library(tsibble) library(tseries) library(MASS) library(forecast) library(TSA) library(TTR) library(aTSA) library(graphics) ``` ## READ DATA ```{r} library(readxl) data <- read_xlsx("C:\\Users\\USER\\OneDrive - apps.ipb.ac.id\\SEM 5 Husna\\Metode Peramalan Deret Waktu\\Tugas Kelompok\\mpdwlengkap.xlsx") data.ts <- ts(data$jumlah) data ``` ## PLOT ```{r} plot(data.ts, xlab = "Tanggal", ylab = "Jumlah Kunjungan") ``` ## DATA LATIH ```{r} datatrain<-data$jumlah[1:400] train.ts<-ts(datatrain) plot.ts(train.ts, lty=1, xlab="Waktu", ylab="Jumlah Kunjungan", main="Plot Data Train") ``` ## DATA TEST ```{r} datatest<-data$jumlah[401:500] test.ts<-ts(datatest) plot.ts(test.ts, lty=1, xlab="Waktu", ylab="Jumlah Kunjungan", main="Plot Data Test") ``` ## UJI STASIONERITAS DATA ### PLOT ACF ```{r} acf(train.ts) ``` ### PLOT BOXCOX ```{r} index <- seq(1:400) bc = boxcox(train.ts~index, lambda = seq(-2,2,by=0.01)) lambda <- bc$x[which.max(bc$y)] lambda bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] ``` ## TRANSFORMASI BOXCOX ```{r} transformed <- BoxCox(data.ts, 0.14) transformed.ts <- ts(transformed) ``` ```{r} trt <- transformed.ts[1:400] trt.ts <- ts(trt) tet <- transformed.ts[401:500] tet.ts <- ts(tet) ``` ```{r} acf(trt.ts) ``` ## DIFFERENCING ```{r} trtdiff <- diff(trt.ts, differences = 1) plot.ts(trtdiff) ``` ```{r} acf(trtdiff) ``` ```{r} pacf(trtdiff) ``` ```{r} eacf(trtdiff) ``` MA(1) ARIMA(0,1,1) ARIMA(0,1,2) ARIMA(1,1,1) ARIMA(1,1,2) ## PENDUGAAN PARAMETER MODEL TENTATIF ### MA(1) ```{r} model1 <- Arima(trtdiff, order = c(0,0,1), method = "ML") summary(model1) lmtest::coeftest(model1) ``` ### ARIMA(0,1,1) ```{r} model2 <- Arima(trtdiff, order = c(0,0,1), method = "ML") summary(model2) lmtest::coeftest(model2) ``` ### ARIMA(0,1,2) ```{r} model3 <- Arima(trtdiff, order = c(0,0,2), method = "ML") summary(model3) lmtest::coeftest(model3) ``` ### ARIMA(1,1,2) ```{r} model4 <- Arima(trtdiff, order = c(1,0,2), method = "ML") summary(model4) lmtest::coeftest(model4) ``` ### ARIMA(1,1,1) ```{r} model5 <- Arima(trtdiff, order = c(1,0,1), method = "ML") summary(model5) lmtest::coeftest(model5) ``` ## UJI DIAGNOSTIK ### EKSPLORASI SISAAN ```{r} sisaanmodel <- model4$residuals par(mfrow=c(2,2)) qqnorm(sisaanmodel) qqline(sisaanmodel, col = "blue", lwd = 2) plot(c(1:length(sisaanmodel)),sisaanmodel) acf(sisaanmodel) pacf(sisaanmodel) ``` ## UJI FORMAL ```{r} #1) Sisaan Menyebar Normal ks.test(sisaanmodel,"pnorm") ``` ```{r} #2) Sisaan saling bebas/tidak ada autokorelasi Box.test(sisaanmodel, type = "Ljung") #tak tolak H0 > sisaan saling bebas ``` ```{r} #3) Sisaan homogen Box.test((sisaanmodel)^2, type = "Ljung") #tak tolak H0 > sisaan homogen ``` ```{r} #4) Nilai tengah sisaan sama dengan nol t.test(sisaanmodel, mu = 0, conf.level = 0.95) #tak tolak h0 > nilai tengah sisaan sama dengan 0 ``` ## OVERFITTING ARIMA(1,1,2) jadi ARIMA(1,1,3) ARIMA(2,1,2) ```{r} model1a=Arima(trtdiff, order=c(1,0,3),method="ML") summary(model1a) lmtest::coeftest(model1a) ``` ```{r} model1b=Arima(trtdiff, order=c(2,0,2),method="ML") summary(model1b) lmtest::coeftest(model1b) ``` ARIMA(2,1,2) ## SISAAN MODEL TERBAIK ```{r} sisaanmodell <- model1b$residuals par(mfrow=c(2,2)) qqnorm(sisaanmodell) qqline(sisaanmodell, col = "black", lwd = 2) plot(c(1:length(sisaanmodell)), sisaanmodell, type = "p", # Plot tipe point main = "Residual vs Order", xlab = "Order", ylab = "Residual", col = "black") # Warna residual acf(sisaanmodell) pacf(sisaanmodell) ``` ```{r} #1) Sisaan Menyebar Normal ks.test(sisaanmodell,"pnorm") ``` ```{r} #2) Sisaan saling bebas/tidak ada autokorelasi Box.test(sisaanmodell, type = "Ljung") #tak tolak H0 > sisaan saling bebas ``` ```{r} #3) Sisaan homogen Box.test((sisaanmodell)^2, type = "Ljung") #tak tolak H0 > sisaan homogen ``` ```{r} #4) Nilai tengah sisaan sama dengan nol t.test(sisaanmodell, mu = 0, conf.level = 0.95) #tak tolak h0 > nilai tengah sisaan sama dengan 0 ``` ## PERAMALAN ```{r} ramalan <- forecast::forecast(model1b, h = 100) ramalan dataramalan <- ramalan$mean plot(ramalan) ``` ```{r} tra <- InvBoxCox(ramalan$mean, 0.14) ``` ```{r} traa <- diffinv(tra, differences = 1) ``` ```{r} perbandingan.da<-matrix(data=c(head(test.ts, n=100), tra[-1]), nrow = 100, ncol = 2) colnames(perbandingan.da)<-c("Aktual","Hasil Forecast") perbandingan.da ``` ```{r} accuracy(ts(traa[-1]), head(test.ts, n=100)) ``` ## UJI EFEK ARCH GARCH ```{r} library(FinTS) ``` ```{r} cobaefek <- residuals(model1b) arch_test <- ArchTest(cobaefek, lags = 15) print(arch_test) ``` GARCH ## TAMBAH EFEK ARCH GARCH ```{r} library(rugarch) ``` ```{r} spec1 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)), mean.model = list(armaOrder = c(2,2), include.mean = TRUE), distribution.model = "norm") garch_fit1 <- ugarchfit(spec = spec1, data = trtdiff) garch_fit1 ``` ARCH(1) : tidk signifikan ARCH(2) : 2.5658, signifikan semua GARCH(1,1) : 2.4897, signifikan semua GARCH(1,2) : tidk signifikan GARCH(2,1) : tidk signifikan GARCH(2,2) : tidk signifikan GARCH(3,1) -(3,5) : tidk signifikan GARCH(3,2) : tidk signifikan GARCH(3,3) : ```{r} # Ambil residual standar dari model GARCH garch_residuals1 <- residuals(garch_fit1, standardize = TRUE) # Melakukan uji ARCH pada residual standar ArchTest(garch_residuals1) ``` ```{r} forc <- ugarchforecast(garch_fit1, data = trtdiff, n.ahead = 100) ramalan <- forc@forecast$seriesFor[,1] pt_1 <- train.ts #nilai akhir data latih hasil.forc <- diffinv(ramalan, differences = 1) + pt_1 hasil.forcc <- InvBoxCox(hasil.forc, 0.14) perbandingan <- data.frame("Aktual"= train.ts, "Ramalan" = hasil.forcc) perbandingan ``` ```{r} forc <- ugarchforecast(garch_fit1, data = trtdiff, n.ahead = 100, n.roll = 0) forc ``` ```{r} plot(forc, which= 1) ``` ```{r} forc <- ugarchforecast(garch_fit1, data = trtdiff, n.ahead = 100) ramalan <- forc@forecast$seriesFor[,1] pt_1 <- train.ts #nilai akhir data latih hasil.forc <- diffinv(ramalan, differences = 1) + pt_1 hasil.forcc <- InvBoxCox(hasil.forc, 0.14) perbandingan <- data.frame("Aktual"= train.ts, "Ramalan" = hasil.forcc) perbandingan ``` ```{r} forecast_and_mape <- function(horizon, model, data_test, last_value, lambda) { forc <- ugarchforecast(model, data = trtdiff, n.ahead = horizon) ramalan <- forc@forecast$seriesFor[,1] # Transformasi balik differencing (dengan nilai awal dari trt.ts) hasil_diffinv <- diffinv(ramalan, differences = 1, xi = tail(last_value, 1))[-1] # Transformasi balik Box-Cox hasil_invboxcox <- InvBoxCox(hasil_diffinv, lambda) # Ambil aktual dari test data sebanyak 'horizon' actual <- head(data_test, horizon) # Hitung MAPE mape_val <- mape(actual, hasil_invboxcox[1:horizon]) # Tampilkan hasil list( horizon = horizon, forecast = hasil_invboxcox[1:horizon], actual = actual, mape = mape_val ) } ``` ```{r} # Nilai terakhir dari trt.ts (bukan train.ts) last_value <- trt.ts lambda <- 0.14 hasil_7 <- forecast_and_mape(7, garch_fit1, test.ts, last_value, lambda) hasil_30 <- forecast_and_mape(30, garch_fit1, test.ts, last_value, lambda) hasil_60 <- forecast_and_mape(60, garch_fit1, test.ts, last_value, lambda) hasil_100 <- forecast_and_mape(100, garch_fit1, test.ts, last_value, lambda) # Tampilkan MAPEs cat("MAPE 7 periode:", hasil_7$mape, "\n") cat("MAPE 30 periode:", hasil_30$mape, "\n") cat("MAPE 60 periode:", hasil_60$mape, "\n") cat("MAPE 100 periode:", hasil_100$mape, "\n") ``` ```{r} hasil_7 ``` ```{r} pred1 = c(14.96, 15.47, 15.68, 15.83, 15.97, 16.1, 16.23) pred2 = c(15.79,16.84,17.56,18.04,18.39,18.67,18.92) aktual7 = c(12, 14, 18, 17, 14, 19, 18) errorgarch = abs(aktual7 - pred1) errorgarch errorlstm = abs(aktual7- pred2) errorlstm ``` ```{r} pred2 = c(15.79,16.84,17.56,18.04,18.39,18.67,18.92) ``` ```{r} t.test(errorlstm, errorgarch, paired = TRUE) ```