# ==============================
# PRAKTIKUM 13 - TIME SERIES
# Forecasting Curah Hujan
# ==============================
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.3
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(writexl)
library(tidyr)
# 1. Import data
data <- read_excel("D:/S2 IPB/Semester 4/Praktikum Stat for Sosial/PRAKTIKUM SESI UAS/P14/Data.xlsx")
# Cek data
head(data)
## # A tibble: 6 × 2
## Tanggal `Curah Hujan`
## <dttm> <dbl>
## 1 2017-01-01 00:00:00 282
## 2 2017-02-01 00:00:00 183
## 3 2017-03-01 00:00:00 166.
## 4 2017-04-01 00:00:00 211.
## 5 2017-05-01 00:00:00 185.
## 6 2017-06-01 00:00:00 135.
str(data)
## tibble [83 × 2] (S3: tbl_df/tbl/data.frame)
## $ Tanggal : POSIXct[1:83], format: "2017-01-01" "2017-02-01" ...
## $ Curah Hujan: num [1:83] 282 183 166 211 185 ...
# 2. Ubah menjadi data time series bulanan
hujan_ts <- ts(data$`Curah Hujan`,
start = c(2017, 1),
frequency = 12)
# 3. Plot data time series
autoplot(hujan_ts) +
labs(title = "Plot Time Series Curah Hujan",
x = "Tahun",
y = "Curah Hujan") +
theme_minimal()
# 4. Membuat model smoothing
# =========================================================
# MEMBUAT SEMUA MODEL
# =========================================================
# 1. Simple Exponential Smoothing
model_simple <- ets(hujan_ts, model = "ANN")
# 2. Holt's Linear Trend
model_holt <- ets(hujan_ts, model = "AAN")
# 3. Brown's Linear Trend
model_brown <- holt(hujan_ts,
h = 12,
exponential = FALSE,
damped = FALSE)
# 4. Damped Trend
model_damped <- ets(hujan_ts,
model = "AAN",
damped = TRUE)
# 5. Simple Seasonal
model_seasonal <- ets(hujan_ts, model = "ANA")
# 6. Winters Additive
model_win_add <- ets(hujan_ts, model = "AAA")
# 7. Winters Multiplicative
# Karena model multiplicative tidak boleh ada nilai 0
hujan_baru <- hujan_ts + 1
model_win_mult <- ets(hujan_baru, model = "MAM")
# 5. Forecast 12 bulan ke depan
fc_simple <- forecast(model_simple, h = 12)
fc_holt <- forecast(model_holt, h = 12)
fc_brown <- model_brown
fc_damped <- forecast(model_damped, h = 12)
fc_seasonal <- forecast(model_seasonal, h = 12)
fc_win_add <- forecast(model_win_add, h = 12)
fc_win_mult <- forecast(model_win_mult, h = 12)
# =========================================================
# 5. FUNGSI AKURASI
# =========================================================
hitung_akurasi <- function(model, aktual){
fitted_value <- as.numeric(fitted(model))
MAD <- mean(abs(aktual - fitted_value), na.rm = TRUE)
MSD <- mean((aktual - fitted_value)^2, na.rm = TRUE)
MAPE <- mean(abs((aktual - fitted_value) / aktual), na.rm = TRUE) * 100
return(c(MAD, MSD, MAPE))
}
hitung_akurasi_brown <- function(model, aktual){
fitted_value <- as.numeric(model$fitted)
MAD <- mean(abs(aktual - fitted_value), na.rm = TRUE)
MSD <- mean((aktual - fitted_value)^2, na.rm = TRUE)
MAPE <- mean(abs((aktual - fitted_value) / aktual), na.rm = TRUE) * 100
return(c(MAD, MSD, MAPE))
}
# =========================================================
# 6. TABEL AKURASI MODEL
# =========================================================
aktual <- as.numeric(hujan_ts)
akurasi <- data.frame(
Model = c("Simple",
"Holt's Linear Trends",
"Brown's linear trend",
"Damped trend",
"Simple Seasonal",
"Winters' additive",
"Winters' multiplicative"),
rbind(
hitung_akurasi(model_simple, aktual),
hitung_akurasi(model_holt, aktual),
hitung_akurasi_brown(model_brown, aktual),
hitung_akurasi(model_damped, aktual),
hitung_akurasi(model_seasonal, aktual),
hitung_akurasi(model_win_add, aktual),
hitung_akurasi(model_win_mult, aktual + 1)
)
)
colnames(akurasi) <- c("Model", "MAD", "MSD", "MAPE")
akurasi <- akurasi %>%
arrange(MAPE)
akurasi
## Model MAD MSD MAPE
## 1 Winters' multiplicative 193.2617 232453.66 5542.932
## 2 Holt's Linear Trends 137.2623 33634.61 995088.680
## 3 Brown's linear trend 137.2663 33634.61 995488.072
## 4 Damped trend 137.3083 33599.97 1009667.820
## 5 Simple 137.4444 33605.81 1011463.389
## 6 Simple Seasonal 116.7898 25022.18 1016290.203
## 7 Winters' additive 118.1566 25422.86 1058433.042
cat("Model terbaik berdasarkan MAPE adalah:",
akurasi$Model[1])
## Model terbaik berdasarkan MAPE adalah: Winters' multiplicative
# =========================================================
# 7. TABEL HASIL FORECAST 12 BULAN
# =========================================================
hasil_forecast <- data.frame(
Bulan = seq(as.Date("2023-12-01"),
by = "month",
length.out = 12),
Simple = as.numeric(fc_simple$mean),
Holt = as.numeric(fc_holt$mean),
Brown = as.numeric(fc_brown$mean),
Damped = as.numeric(fc_damped$mean),
Simple_Seasonal = as.numeric(fc_seasonal$mean),
Winters_Additive = as.numeric(fc_win_add$mean),
Winters_Multiplicative = as.numeric(fc_win_mult$mean - 1)
)
hasil_forecast
## Bulan Simple Holt Brown Damped Simple_Seasonal
## 1 2023-12-01 526.2805 523.8660 523.7625 526.3737 398.4404
## 2 2024-01-01 526.2805 521.3435 521.2407 526.4028 311.9300
## 3 2024-02-01 526.2805 518.8209 518.7190 526.4287 358.4291
## 4 2024-03-01 526.2805 516.2984 516.1973 526.4517 290.5097
## 5 2024-04-01 526.2805 513.7758 513.6755 526.4722 358.7827
## 6 2024-05-01 526.2805 511.2533 511.1538 526.4904 450.3155
## 7 2024-06-01 526.2805 508.7307 508.6321 526.5065 369.2043
## 8 2024-07-01 526.2805 506.2082 506.1103 526.5209 252.0884
## 9 2024-08-01 526.2805 503.6856 503.5886 526.5337 209.3446
## 10 2024-09-01 526.2805 501.1631 501.0669 526.5450 228.3205
## 11 2024-10-01 526.2805 498.6405 498.5451 526.5551 444.8381
## 12 2024-11-01 526.2805 496.1180 496.0234 526.5641 478.1367
## Winters_Additive Winters_Multiplicative
## 1 408.8182 263.71819
## 2 328.0666 182.33103
## 3 390.1592 167.98569
## 4 335.7617 117.56665
## 5 398.0604 124.07468
## 6 502.7473 153.45375
## 7 431.7009 116.84293
## 8 305.3906 49.91876
## 9 285.1676 36.85410
## 10 311.9151 43.25642
## 11 533.9425 261.59579
## 12 562.9853 381.00743
# =========================================================
# 8. DATA AKTUAL DAN FITTED VALUE
# =========================================================
tanggal_aktual <- seq(as.Date("2017-01-01"),
by = "month",
length.out = length(hujan_ts))
data_fitted <- data.frame(
Bulan = tanggal_aktual,
Aktual = as.numeric(hujan_ts),
Simple = as.numeric(fitted(model_simple)),
Holt = as.numeric(fitted(model_holt)),
Brown = as.numeric(fc_brown$fitted),
Damped = as.numeric(fitted(model_damped)),
Simple_Seasonal = as.numeric(fitted(model_seasonal)),
Winters_Additive = as.numeric(fitted(model_win_add)),
Winters_Multiplicative = as.numeric(fitted(model_win_mult) - 1)
)
head(data_fitted)
## Bulan Aktual Simple Holt Brown Damped Simple_Seasonal
## 1 2017-01-01 282.00 274.7301 281.3341 281.3341 277.9523 195.4398
## 2 2017-02-01 183.00 281.5339 279.3866 279.3873 278.6778 309.8532
## 3 2017-03-01 165.50 189.3177 186.5331 186.5583 186.3892 142.4285
## 4 2017-04-01 210.75 167.0271 164.2513 164.2590 164.4068 228.8284
## 5 2017-05-01 184.75 207.9466 205.2224 205.2119 205.6303 306.1845
## 6 2017-06-01 135.25 186.2373 183.4683 183.4736 184.1682 129.8534
## Winters_Additive Winters_Multiplicative
## 1 96.11674 152.7525
## 2 303.03930 192.3220
## 3 155.22001 133.6812
## 4 225.58127 156.6410
## 5 318.78563 223.2270
## 6 143.41943 159.5508
# =========================================================
# 9. GABUNG DATA AKTUAL, FITTED, DAN FORECAST
# =========================================================
data_fitted_long <- data_fitted %>%
pivot_longer(cols = -Bulan,
names_to = "Model",
values_to = "Nilai") %>%
mutate(Tipe = ifelse(Model == "Aktual", "Aktual", "Prediksi/Fitted"))
data_forecast_long <- hasil_forecast %>%
pivot_longer(cols = -Bulan,
names_to = "Model",
values_to = "Nilai") %>%
mutate(Tipe = "Forecast")
data_plot_all <- bind_rows(data_fitted_long, data_forecast_long)
# =========================================================
# 10. VISUALISASI AKTUAL, FITTED, DAN FORECAST SEMUA MODEL
# =========================================================
ggplot(data_plot_all,
aes(x = Bulan,
y = Nilai,
color = Tipe,
linetype = Tipe)) +
geom_line(linewidth = 0.8, na.rm = TRUE) +
facet_wrap(~ Model, scales = "free_y", ncol = 2) +
labs(title = "Visualisasi Nilai Aktual, Prediksi, dan Forecast",
x = "Bulan",
y = "Curah Hujan",
color = "Keterangan",
linetype = "Keterangan") +
theme_minimal() +
theme(legend.position = "bottom")
# =========================================================
# VISUALISASI AKTUAL, PREDIKSI/FITTED, DAN FORECAST
# UNTUK SEMUA MODEL
# Merah = Aktual, Biru = Prediksi/Fitted, Hijau = Forecast
# =========================================================
library(ggplot2)
library(dplyr)
library(tidyr)
tanggal_aktual <- seq(as.Date("2017-01-01"),
by = "month",
length.out = length(hujan_ts))
tanggal_forecast <- seq(as.Date("2023-12-01"),
by = "month",
length.out = 12)
# Data aktual dan fitted semua model
data_plot <- data.frame(
Bulan = rep(tanggal_aktual, 8),
Nilai = c(
as.numeric(hujan_ts),
as.numeric(fitted(model_simple)),
as.numeric(fitted(model_holt)),
as.numeric(fc_brown$fitted),
as.numeric(fitted(model_damped)),
as.numeric(fitted(model_seasonal)),
as.numeric(fitted(model_win_add)),
as.numeric(fitted(model_win_mult) - 1)
),
Model = c(
rep("Aktual", length(hujan_ts)),
rep("Simple", length(hujan_ts)),
rep("Holt", length(hujan_ts)),
rep("Brown", length(hujan_ts)),
rep("Damped", length(hujan_ts)),
rep("Simple Seasonal", length(hujan_ts)),
rep("Winters Additive", length(hujan_ts)),
rep("Winters Multiplicative", length(hujan_ts))
),
Tipe = c(
rep("Aktual", length(hujan_ts)),
rep("Prediksi", length(hujan_ts) * 7)
)
)
# Data forecast semua model
data_forecast <- data.frame(
Bulan = rep(tanggal_forecast, 7),
Nilai = c(
as.numeric(fc_simple$mean),
as.numeric(fc_holt$mean),
as.numeric(fc_brown$mean),
as.numeric(fc_damped$mean),
as.numeric(fc_seasonal$mean),
as.numeric(fc_win_add$mean),
as.numeric(fc_win_mult$mean - 1)
),
Model = c(
rep("Simple", 12),
rep("Holt", 12),
rep("Brown", 12),
rep("Damped", 12),
rep("Simple Seasonal", 12),
rep("Winters Additive", 12),
rep("Winters Multiplicative", 12)
),
Tipe = rep("Forecast", 12 * 7)
)
# Gabungkan aktual ke setiap model agar setiap panel punya garis merah
data_aktual_per_model <- data.frame(
Bulan = rep(tanggal_aktual, 7),
Nilai = rep(as.numeric(hujan_ts), 7),
Model = rep(c("Simple",
"Holt",
"Brown",
"Damped",
"Simple Seasonal",
"Winters Additive",
"Winters Multiplicative"),
each = length(hujan_ts)),
Tipe = "Aktual"
)
# Data fitted tanpa panel aktual sendiri
data_fitted_model <- data_plot %>%
filter(Model != "Aktual")
# Gabungkan semua
data_visualisasi <- bind_rows(
data_aktual_per_model,
data_fitted_model,
data_forecast
)
# Visualisasi semua model
ggplot(data_visualisasi,
aes(x = Bulan,
y = Nilai,
color = Tipe,
linetype = Tipe)) +
geom_line(linewidth = 0.8, na.rm = TRUE) +
facet_wrap(~ Model, scales = "free_y", ncol = 2) +
scale_color_manual(values = c(
"Aktual" = "red",
"Prediksi" = "blue",
"Forecast" = "green"
)) +
scale_linetype_manual(values = c(
"Aktual" = "solid",
"Prediksi" = "dashed",
"Forecast" = "dotted"
)) +
labs(title = "Visualisasi Nilai Aktual, Prediksi, dan Forecast Semua Model",
x = "Bulan",
y = "Curah Hujan",
color = "Keterangan",
linetype = "Keterangan") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "bottom"
)
# =========================================================
# VISUALISASI AKTUAL (MERAH)
# FITTED/PREDIKSI (BIRU)
# FORECAST (HIJAU)
# =========================================================
library(ggplot2)
# =========================================================
# CONTOH UNTUK SIMPLE SEASONAL
# =========================================================
# Data aktual
tanggal_aktual <- seq(as.Date("2017-01-01"),
by = "month",
length.out = length(hujan_ts))
# Data forecast
tanggal_forecast <- seq(as.Date("2024-01-01"),
by = "month",
length.out = 12)
# Membuat dataframe
plot_data <- data.frame(
Bulan = c(tanggal_aktual,
tanggal_aktual,
tanggal_forecast),
Nilai = c(as.numeric(hujan_ts),
as.numeric(fitted(model_seasonal)),
as.numeric(fc_seasonal$mean)),
Tipe = c(rep("Aktual", length(hujan_ts)),
rep("Prediksi", length(hujan_ts)),
rep("Forecast", 12))
)
# =========================================================
# VISUALISASI
# =========================================================
ggplot(plot_data,
aes(x = Bulan,
y = Nilai,
color = Tipe,
linetype = Tipe)) +
geom_line(linewidth = 1) +
scale_color_manual(values = c(
"Aktual" = "red",
"Prediksi" = "blue",
"Forecast" = "green"
)) +
labs(title = "Aktual vs Prediksi vs Forecast",
x = "Bulan",
y = "Curah Hujan") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5,
face = "bold"),
legend.position = "bottom"
)
# =========================================================
# 11. VISUALISASI ACTUAL VS FORECAST MODEL SIMPLE SEASONAL
# =========================================================
plot_seasonal <- data.frame(
Bulan = c(tanggal_aktual, hasil_forecast$Bulan),
Nilai = c(as.numeric(hujan_ts),
hasil_forecast$Simple_Seasonal),
Tipe = c(rep("Aktual", length(hujan_ts)),
rep("Forecast", 12))
)
ggplot(plot_seasonal,
aes(x = Bulan,
y = Nilai,
color = Tipe)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
labs(title = "Actual vs Forecast - Simple Seasonal",
x = "Bulan",
y = "Curah Hujan",
color = "Keterangan") +
theme_minimal()
# =========================================================
# 12. VISUALISASI FORECAST SEMUA MODEL SAJA
# =========================================================
forecast_long <- hasil_forecast %>%
pivot_longer(cols = -Bulan,
names_to = "Model",
values_to = "Forecast")
ggplot(forecast_long,
aes(x = Bulan,
y = Forecast,
color = Model)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
labs(title = "Perbandingan Forecast Semua Model",
x = "Bulan",
y = "Curah Hujan") +
theme_minimal() +
theme(legend.position = "bottom")
# =========================================================
# 13. SIMPAN HASIL KE EXCEL
# =========================================================
write_xlsx(
list(
Akurasi_Model = akurasi,
Forecast_12_Bulan = hasil_forecast,
Aktual_dan_Fitted = data_fitted
),
"Hasil_Praktikum_13_Time_Series.xlsx"
)