# 1. INSTALL PACKAGE
packages <- c(
"readxl",
"forecast",
"tseries",
"openxlsx"
)
installed <- rownames(installed.packages())
for (p in packages) {
if (!(p %in% installed)) {
install.packages(p)
}
}
# 2. LOAD PACKAGE
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.5.3
# 3. MEMBACA DATA EXCEL
data_raw <- read_excel(
"C:/Users/ASUS/Downloads/Analisis Deret Waktu/Data UTS ADW.xlsx",
col_names = FALSE
)
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
# 4. MENGAMBIL DATA YANG DIPERLUKAN
data <- data_raw[4:nrow(data_raw), c(2, 3)]
colnames(data) <- c(
"Tanggal",
"Y"
)
# 5. MEMBERSIHKAN DATA
data$Tanggal <- as.Date(
data$Tanggal,
format = "%d-%b-%y"
)
data$Y <- as.numeric(data$Y)
data <- data[
!is.na(data$Tanggal) &
!is.na(data$Y),
]
rownames(data) <- NULL
cat(
"Jumlah observasi =",
nrow(data),
"\n"
)
## Jumlah observasi = 220
# 6. MEMBAGI DATA TRAINING DAN TESTING
TRAIN_SIZE <- 176
train <- data[
1:TRAIN_SIZE,
]
test <- data[
(TRAIN_SIZE + 1):nrow(data),
]
cat(
"Jumlah data training =",
nrow(train),
"\n"
)
## Jumlah data training = 176
cat(
"Jumlah data testing =",
nrow(test),
"\n"
)
## Jumlah data testing = 44
# 7. GRAFIK DATA AKTUAL
plot(
data$Tanggal,
data$Y,
type = "l",
xlab = "Tanggal",
ylab = "Y",
main = "Data Aktual",
lwd = 2
)
abline(
v = train$Tanggal[nrow(train)],
lty = 2
)
# 8. TRANSFORMASI PANGKAT
LAMBDA <- 0.256958404
train$Transformasi <- (
train$Y ^ LAMBDA
)
print(
head(
train[
,
c(
"Tanggal",
"Y",
"Transformasi"
)
],
10
)
)
## # A tibble: 10 × 3
## Tanggal Y Transformasi
## <date> <dbl> <dbl>
## 1 2021-01-02 6 1.58
## 2 2021-01-03 5 1.51
## 3 2021-01-04 2 1.19
## 4 2021-01-05 2 1.19
## 5 2021-01-06 1 1
## 6 2021-01-08 6 1.58
## 7 2021-01-09 7 1.65
## 8 2021-01-10 3 1.33
## 9 2021-01-13 9 1.76
## 10 2021-01-16 8 1.71
# 9. DIFFERENCING
train$Differencing <- c(
NA,
diff(train$Transformasi)
)
diff_train <- na.omit(
train$Differencing
)
print(
head(
train[
,
c(
"Tanggal",
"Y",
"Transformasi",
"Differencing"
)
],
15
)
)
## # A tibble: 15 × 4
## Tanggal Y Transformasi Differencing
## <date> <dbl> <dbl> <dbl>
## 1 2021-01-02 6 1.58 NA
## 2 2021-01-03 5 1.51 -0.0725
## 3 2021-01-04 2 1.19 -0.317
## 4 2021-01-05 2 1.19 0
## 5 2021-01-06 1 1 -0.195
## 6 2021-01-08 6 1.58 0.585
## 7 2021-01-09 7 1.65 0.0640
## 8 2021-01-10 3 1.33 -0.323
## 9 2021-01-13 9 1.76 0.433
## 10 2021-01-16 8 1.71 -0.0524
## 11 2021-01-17 1 1 -0.706
## 12 2021-01-18 6 1.58 0.585
## 13 2021-01-20 9 1.76 0.174
## 14 2021-01-21 5 1.51 -0.247
## 15 2021-01-23 8 1.71 0.194
# 10. GRAFIK DIFFERENCING
plot(
train$Tanggal[-1],
diff_train,
type = "l",
xlab = "Tanggal",
ylab = "Differencing",
main = "Data Setelah Differencing",
lwd = 2
)
abline(
h = 0,
lty = 2
)
# 11. UJI AUGMENTED DICKEY-FULLER (ADF)
adf_test <- adf.test(
diff_train
)
## Warning in adf.test(diff_train): p-value smaller than printed p-value
print(adf_test)
##
## Augmented Dickey-Fuller Test
##
## data: diff_train
## Dickey-Fuller = -8.4233, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
if (
adf_test$p.value < 0.05
) {
cat(
"Kesimpulan ADF: Data stasioner\n"
)
} else {
cat(
"Kesimpulan ADF: Data belum stasioner\n"
)
}
## Kesimpulan ADF: Data stasioner
# 12. UJI KPSS
kpss_test <- kpss.test(
diff_train,
null = "Level"
)
## Warning in kpss.test(diff_train, null = "Level"): p-value greater than printed
## p-value
print(kpss_test)
##
## KPSS Test for Level Stationarity
##
## data: diff_train
## KPSS Level = 0.061574, Truncation lag parameter = 4, p-value = 0.1
if (
kpss_test$p.value > 0.05
) {
cat(
"Kesimpulan KPSS: Data stasioner\n"
)
} else {
cat(
"Kesimpulan KPSS: Data belum stasioner\n"
)
}
## Kesimpulan KPSS: Data stasioner
# 13. ACF
acf(
diff_train,
lag.max = 30,
main = "ACF Data Setelah Differencing"
)
# 14. PACF
pacf(
diff_train,
lag.max = 30,
main = "PACF Data Setelah Differencing"
)
# 15. MEMBENTUK MODEL ARIMA(0,1,1)
model_arima <- Arima(
train$Transformasi,
order = c(
0,
1,
1
),
include.drift = FALSE
)
print(
summary(model_arima)
)
## Series: train$Transformasi
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.7682
## s.e. 0.0425
##
## sigma^2 = 0.04982: log likelihood = 14.18
## AIC=-24.36 AICc=-24.29 BIC=-18.03
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02882704 0.2219338 0.1748295 -0.08402798 10.08546 0.7884469
## ACF1
## Training set -0.0280261
# 16. KOEFISIEN MODEL
cat(
"Koefisien model:\n"
)
## Koefisien model:
print(
coef(model_arima)
)
## ma1
## -0.768157
# 17. AIC DAN BIC
cat(
"AIC =",
AIC(model_arima),
"\n"
)
## AIC = -24.36424
cat(
"BIC =",
BIC(model_arima),
"\n"
)
## BIC = -18.03467
# 18. RESIDUAL MODEL
residual <- residuals(
model_arima
)
print(
head(
residual,
10
)
)
## Time Series:
## Start = 1
## End = 10
## Frequency = 1
## [1] 0.001584719 -0.057518195 -0.319066274 -0.211085059 -0.339665572
## [6] 0.325777569 0.307652346 -0.087899336 0.364280789 0.226136887
# 19. GRAFIK RESIDUAL
plot(
residual,
type = "l",
xlab = "Observasi",
ylab = "Residual",
main = "Residual ARIMA(0,1,1)"
)
abline(
h = 0,
lty = 2
)
# 20. UJI LJUNG-BOX
lb6 <- Box.test(
residual,
lag = 6,
type = "Ljung-Box"
)
lb12 <- Box.test(
residual,
lag = 12,
type = "Ljung-Box"
)
lb18 <- Box.test(
residual,
lag = 18,
type = "Ljung-Box"
)
lb24 <- Box.test(
residual,
lag = 24,
type = "Ljung-Box"
)
cat(
"Ljung-Box lag 6:\n"
)
## Ljung-Box lag 6:
print(lb6)
##
## Box-Ljung test
##
## data: residual
## X-squared = 1.8715, df = 6, p-value = 0.9311
cat(
"Ljung-Box lag 12:\n"
)
## Ljung-Box lag 12:
print(lb12)
##
## Box-Ljung test
##
## data: residual
## X-squared = 3.67, df = 12, p-value = 0.9887
cat(
"Ljung-Box lag 18:\n"
)
## Ljung-Box lag 18:
print(lb18)
##
## Box-Ljung test
##
## data: residual
## X-squared = 6.3297, df = 18, p-value = 0.9947
cat(
"Ljung-Box lag 24:\n"
)
## Ljung-Box lag 24:
print(lb24)
##
## Box-Ljung test
##
## data: residual
## X-squared = 17.599, df = 24, p-value = 0.822
# 21. FORECAST DATA TESTING
jumlah_testing <- nrow(test)
forecast_model <- forecast(
model_arima,
h = jumlah_testing
)
forecast_transformasi <- as.numeric(
forecast_model$mean
)
# 22. MENGEMBALIKAN FORECAST KE SKALA ASLI
forecast_y <- (
forecast_transformasi ^
(1 / LAMBDA)
)
# 23. MEMBUAT TABEL AKTUAL VS FORECAST
hasil_testing <- data.frame(
Tanggal = test$Tanggal,
Aktual = test$Y,
Forecast = forecast_y
)
hasil_testing$Error <- (
hasil_testing$Aktual -
hasil_testing$Forecast
)
hasil_testing$Absolute_Error <- abs(
hasil_testing$Error
)
hasil_testing$APE <- (
abs(
hasil_testing$Error /
hasil_testing$Aktual
) * 100
)
print(
hasil_testing
)
## Tanggal Aktual Forecast Error Absolute_Error APE
## 1 2021-08-11 55 45.02159 9.97841 9.97841 18.142564
## 2 2021-08-12 53 45.02159 7.97841 7.97841 15.053605
## 3 2021-08-13 37 45.02159 -8.02159 8.02159 21.679972
## 4 2021-08-14 52 45.02159 6.97841 6.97841 13.420020
## 5 2021-08-15 56 45.02159 10.97841 10.97841 19.604304
## 6 2021-08-16 55 45.02159 9.97841 9.97841 18.142564
## 7 2021-08-17 35 45.02159 -10.02159 10.02159 28.633113
## 8 2021-08-18 40 45.02159 -5.02159 5.02159 12.553974
## 9 2021-08-19 31 45.02159 -14.02159 14.02159 45.230934
## 10 2021-08-20 52 45.02159 6.97841 6.97841 13.420020
## 11 2021-08-21 42 45.02159 -3.02159 3.02159 7.194261
## 12 2021-08-22 33 45.02159 -12.02159 12.02159 36.429059
## 13 2021-08-23 56 45.02159 10.97841 10.97841 19.604304
## 14 2021-08-24 42 45.02159 -3.02159 3.02159 7.194261
## 15 2021-08-25 37 45.02159 -8.02159 8.02159 21.679972
## 16 2021-08-26 59 45.02159 13.97841 13.97841 23.692221
## 17 2021-08-27 56 45.02159 10.97841 10.97841 19.604304
## 18 2021-08-28 40 45.02159 -5.02159 5.02159 12.553974
## 19 2021-08-29 43 45.02159 -2.02159 2.02159 4.701371
## 20 2021-08-30 42 45.02159 -3.02159 3.02159 7.194261
## 21 2021-08-31 58 45.02159 12.97841 12.97841 22.376570
## 22 2021-09-01 44 45.02159 -1.02159 1.02159 2.321794
## 23 2021-09-02 40 45.02159 -5.02159 5.02159 12.553974
## 24 2021-09-03 56 45.02159 10.97841 10.97841 19.604304
## 25 2021-09-04 55 45.02159 9.97841 9.97841 18.142564
## 26 2021-09-05 55 45.02159 9.97841 9.97841 18.142564
## 27 2021-09-06 38 45.02159 -7.02159 7.02159 18.477867
## 28 2021-09-07 52 45.02159 6.97841 6.97841 13.420020
## 29 2021-09-08 33 45.02159 -12.02159 12.02159 36.429059
## 30 2021-09-09 30 45.02159 -15.02159 15.02159 50.071965
## 31 2021-09-10 43 45.02159 -2.02159 2.02159 4.701371
## 32 2021-09-11 51 45.02159 5.97841 5.97841 11.722373
## 33 2021-09-12 30 45.02159 -15.02159 15.02159 50.071965
## 34 2021-09-13 53 45.02159 7.97841 7.97841 15.053605
## 35 2021-09-14 25 45.02159 -20.02159 20.02159 80.086358
## 36 2021-09-15 30 45.02159 -15.02159 15.02159 50.071965
## 37 2021-09-16 28 45.02159 -17.02159 17.02159 60.791391
## 38 2021-09-17 30 45.02159 -15.02159 15.02159 50.071965
## 39 2021-09-18 47 45.02159 1.97841 1.97841 4.209384
## 40 2021-09-19 32 45.02159 -13.02159 13.02159 40.692467
## 41 2021-09-20 38 45.02159 -7.02159 7.02159 18.477867
## 42 2021-09-21 28 45.02159 -17.02159 17.02159 60.791391
## 43 2021-09-22 42 45.02159 -3.02159 3.02159 7.194261
## 44 2021-09-23 38 45.02159 -7.02159 7.02159 18.477867
# 24. MENGHITUNG MAE
MAE <- mean(
abs(
hasil_testing$Aktual -
hasil_testing$Forecast
)
)
# 25. MENGHITUNG MSE
MSE <- mean(
(
hasil_testing$Aktual -
hasil_testing$Forecast
)^2
)
# 26. MENGHITUNG RMSE
RMSE <- sqrt(
MSE
)
# 27. MENGHITUNG MAPE
MAPE <- mean(
abs(
(
hasil_testing$Aktual -
hasil_testing$Forecast
) /
hasil_testing$Aktual
)
) * 100
# 28. HASIL EVALUASI MODEL
cat(
"MAE =",
round(MAE, 4),
"\n"
)
## MAE = 9.0958
cat(
"MSE =",
round(MSE, 4),
"\n"
)
## MSE = 104.405
cat(
"RMSE =",
round(RMSE, 4),
"\n"
)
## RMSE = 10.2179
cat(
"MAPE =",
round(MAPE, 4),
"%\n"
)
## MAPE = 23.8565 %
# 29. INTERPRETASI MAPE
if (
MAPE < 10
) {
cat(
"Interpretasi MAPE: Sangat baik\n"
)
} else if (
MAPE < 20
) {
cat(
"Interpretasi MAPE: Baik\n"
)
} else if (
MAPE < 50
) {
cat(
"Interpretasi MAPE: Layak / reasonable\n"
)
} else {
cat(
"Interpretasi MAPE: Kurang baik\n"
)
}
## Interpretasi MAPE: Layak / reasonable
# 30. GRAFIK AKTUAL VS FORECAST
plot(
train$Tanggal,
train$Y,
type = "l",
xlim = range(data$Tanggal),
ylim = range(
c(
data$Y,
forecast_y
)
),
xlab = "Tanggal",
ylab = "Y",
main = "Aktual vs Forecast ARIMA(0,1,1)",
lwd = 2
)
lines(
test$Tanggal,
test$Y,
lwd = 2
)
lines(
test$Tanggal,
forecast_y,
lty = 2,
lwd = 2
)
abline(
v = train$Tanggal[nrow(train)],
lty = 2
)
legend(
"topright",
legend = c(
"Training",
"Aktual Testing",
"Forecast"
),
lty = c(
1,
1,
2
),
lwd = c(
2,
2,
2
)
)
# 31. PERBANDINGAN BEBERAPA MODEL ARIMA
model_011 <- Arima(
train$Transformasi,
order = c(
0,
1,
1
),
include.drift = FALSE
)
model_110 <- Arima(
train$Transformasi,
order = c(
1,
1,
0
),
include.drift = FALSE
)
model_111 <- Arima(
train$Transformasi,
order = c(
1,
1,
1
),
include.drift = FALSE
)
model_210 <- Arima(
train$Transformasi,
order = c(
2,
1,
0
),
include.drift = FALSE
)
model_012 <- Arima(
train$Transformasi,
order = c(
0,
1,
2
),
include.drift = FALSE
)
perbandingan_model <- data.frame(
Model = c(
"ARIMA(0,1,1)",
"ARIMA(1,1,0)",
"ARIMA(1,1,1)",
"ARIMA(2,1,0)",
"ARIMA(0,1,2)"
),
AIC = c(
AIC(model_011),
AIC(model_110),
AIC(model_111),
AIC(model_210),
AIC(model_012)
),
BIC = c(
BIC(model_011),
BIC(model_110),
BIC(model_111),
BIC(model_210),
BIC(model_012)
)
)
perbandingan_model <- perbandingan_model[
order(
perbandingan_model$AIC
),
]
print(
perbandingan_model
)
## Model AIC BIC
## 1 ARIMA(0,1,1) -24.364237 -18.034665
## 5 ARIMA(0,1,2) -22.400545 -12.906188
## 3 ARIMA(1,1,1) -22.396172 -12.901815
## 4 ARIMA(2,1,0) -2.804853 6.689505
## 2 ARIMA(1,1,0) 13.804190 20.133762
# 32. FORECAST 10 PERIODE KE DEPAN
forecast_10 <- forecast(
model_arima,
h = 10
)
forecast_10_transformasi <- as.numeric(
forecast_10$mean
)
forecast_10_y <- (
forecast_10_transformasi ^
(1 / LAMBDA)
)
# 33. TABEL FORECAST 10 PERIODE
forecast_masa_depan <- data.frame(
Periode = 1:10,
Forecast = forecast_10_y
)
print(
forecast_masa_depan
)
## Periode Forecast
## 1 1 45.02159
## 2 2 45.02159
## 3 3 45.02159
## 4 4 45.02159
## 5 5 45.02159
## 6 6 45.02159
## 7 7 45.02159
## 8 8 45.02159
## 9 9 45.02159
## 10 10 45.02159
# 34. GRAFIK FORECAST MASA DEPAN
plot(
forecast_10,
main = "Forecast 10 Periode ke Depan"
)
# 35. MENYIAPKAN HASIL EVALUASI
hasil_ringkas <- data.frame(
Ukuran = c(
"MAE",
"MSE",
"RMSE",
"MAPE"
),
Nilai = c(
MAE,
MSE,
RMSE,
MAPE
)
)
# 36. MENYIAPKAN HASIL LJUNG-BOX
ljung_box <- data.frame(
Lag = c(
6,
12,
18,
24
),
P_Value = c(
lb6$p.value,
lb12$p.value,
lb18$p.value,
lb24$p.value
)
)
# 37. MENYIMPAN HASIL KE EXCEL
wb <- createWorkbook()
addWorksheet(
wb,
"Data"
)
writeData(
wb,
"Data",
data
)
addWorksheet(
wb,
"Training"
)
writeData(
wb,
"Training",
train
)
addWorksheet(
wb,
"Testing"
)
writeData(
wb,
"Testing",
test
)
addWorksheet(
wb,
"Forecast Testing"
)
writeData(
wb,
"Forecast Testing",
hasil_testing
)
addWorksheet(
wb,
"Evaluasi"
)
writeData(
wb,
"Evaluasi",
hasil_ringkas
)
addWorksheet(
wb,
"Perbandingan Model"
)
writeData(
wb,
"Perbandingan Model",
perbandingan_model
)
addWorksheet(
wb,
"Forecast 10 Periode"
)
writeData(
wb,
"Forecast 10 Periode",
forecast_masa_depan
)
addWorksheet(
wb,
"Ljung Box"
)
writeData(
wb,
"Ljung Box",
ljung_box
)
saveWorkbook(
wb,
"Hasil_Pengujian_ARIMA_011.xlsx",
overwrite = TRUE
)
cat(
"Hasil berhasil disimpan sebagai Hasil_Pengujian_ARIMA_011.xlsx\n"
)
## Hasil berhasil disimpan sebagai Hasil_Pengujian_ARIMA_011.xlsx