# Memuat paket yang diperlukan
library(readxl)

# Membaca data dari file Excel
data <- read_excel("databaru.xlsx")

# Memeriksa struktur data
str(data)
## tibble [60 Ă— 3] (S3: tbl_df/tbl/data.frame)
##  $ Tahun      : num [1:60] 2019 2019 2019 2019 2019 ...
##  $ Bulan      : num [1:60] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Curah Hujan: num [1:60] 637 236 445 346 60 61 2 179 179 179 ...
# Mengambil jumlah baris dan kolom data
row_data <- nrow(data)
col_data <- ncol(data)

# Transformasi logaritma pada kolom 'Curah Hujan'
data$`Curah Hujan` <- log(data$`Curah Hujan`)

# Menentukan proporsi data untuk pelatihan (misalnya 80%)
proporsi_pelatihan <- 0.6
n_pelatihan <- floor(proporsi_pelatihan * row_data)

# Membagi data menjadi set pelatihan dan pengujian
data_pelatihan <- data[1:n_pelatihan, ]
data_pengujian <- data[(n_pelatihan + 1):row_data, ]

# Menghitung nilai awal L0, T0, dan S0 menggunakan data pelatihan

# Menghitung L0 sebagai rata-rata dari 12 bulan pertama data pelatihan
data_l0 <- data_pelatihan[1:12, ]
l0 <- mean(data_l0$`Curah Hujan`)
print(paste("L0:", l0))
## [1] "L0: 4.86005841986022"
# Menghitung T0 sebagai rata-rata perbedaan antara tahun kedua dan pertama
data_y1 <- data_pelatihan[1:12, ]
data_y2 <- data_pelatihan[13:24, ]
y2_y1 <- data_y2$`Curah Hujan` - data_y1$`Curah Hujan`
T0 <- mean(y2_y1) / 12
print(paste("T0:", T0))
## [1] "T0: -0.0107226704295978"
# Menghitung S0 sebagai perbedaan antara nilai aktual dan L0 untuk 12 bulan pertama
S0 <- data_l0$`Curah Hujan` - l0
print("S0:")
## [1] "S0:"
print(S0)
##  [1]  1.5967112  0.6037734  1.2380159  0.9863804 -0.7657139 -0.7491846
##  [7] -4.1669112  0.3273274  0.3273274  0.3273274 -0.5033496  0.7782962
# Inisialisasi vektor Lt, Tt, St, dan Prediksi
Lt <- c(l0)
Tt <- c(T0)
St <- as.numeric(S0)
Prediksi <- numeric()

# Menentukan parameter smoothing
alpha <- 0.07494387
beta <- 0
gamma <- 0.4682306

# Menerapkan metode Holt-Winters pada data pelatihan
for (i in 1:nrow(data_pelatihan)) {
  if (i <= 12) {
    # Untuk 12 bulan pertama, gunakan nilai S0 yang telah dihitung
    st <- St[i]
  } else {
    # Untuk bulan berikutnya, hitung st berdasarkan data sebelumnya
    st <- gamma * (data_pelatihan$`Curah Hujan`[i] - Lt[i]) + (1 - gamma) * St[i - 12]
    St <- c(St, st)
  }
  
  lt <- alpha * (data_pelatihan$`Curah Hujan`[i] - st) + (1 - alpha) * (Lt[i] + Tt[i])
  Lt <- c(Lt, lt)
  
  tt <- beta * (Lt[i + 1] - Lt[i]) + (1 - beta) * Tt[i]
  Tt <- c(Tt, tt)
  
  prediksi <- (Lt[i + 1] + Tt[i + 1]) + st
  Prediksi <- c(Prediksi, prediksi)
}

# Menghapus nilai pertama dari Lt dan Tt karena merupakan nilai awal
Lt <- Lt[-1]
Tt <- Tt[-1]

# Menampilkan Lt, Tt, St, dan Prediksi pada data pelatihan
print("Lt:")
## [1] "Lt:"
print(Lt)
##  [1] 4.850139 4.840964 4.832476 4.824624 4.817360 4.810641 4.804426 4.798676
##  [9] 4.793357 4.788437 4.783885 4.779675 4.768670 4.795307 4.767160 4.711475
## [17] 4.752567 4.754089 4.808335 4.704241 4.618758 4.570695 4.618567 4.665711
## [25] 4.686571 4.690728 4.713108 4.739597 4.714392 4.714204 4.801386 4.796913
## [33] 4.807258 4.797899 4.821001 4.835089
print("Tt:")
## [1] "Tt:"
print(Tt)
##  [1] -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267
##  [7] -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267
## [13] -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267
## [19] -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267
## [25] -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267
## [31] -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267 -0.01072267
print("St:")
## [1] "St:"
print(St)
##  [1]  1.5967112  0.6037734  1.2380159  0.9863804 -0.7657139 -0.7491846
##  [7] -4.1669112  0.3273274  0.3273274  0.3273274 -0.5033496  0.7782962
## [13]  1.5839533  1.0332683  1.0238595  0.4486782 -0.1663856 -0.6147664
## [19] -3.4130381 -0.7791360 -0.5604724 -0.1208165  0.1756350  1.4487316
## [25]  1.9455733  1.1986497  1.4033332  0.8764351 -0.3459783 -0.5004338
## [31] -2.2721924 -0.7151604 -0.3223809 -0.1142432  0.5636037  1.7307916
print("Prediksi pada data pelatihan:")
## [1] "Prediksi pada data pelatihan:"
print(Prediksi)
##  [1] 6.4361279 5.4340144 6.0597688 5.8002814 4.0409237 4.0507339 0.6267917
##  [8] 5.1152806 5.1099617 5.1050415 4.2698131 5.5472486 6.3419006 5.8178525
## [15] 5.7802970 5.1494308 4.5754590 4.1285999 1.3845742 3.9143820 4.0475625
## [22] 4.4391561 4.7834794 6.1037203 6.6214217 5.8786553 6.1057181 5.6053089
## [29] 4.3576906 4.2030474 2.5184713 4.0710295 4.4741548 4.6729328 5.3738821
## [36] 6.5551581
# Menghitung error pada data pelatihan
error_pelatihan <- data_pelatihan$`Curah Hujan` - Prediksi
print("Error pada data pelatihan:")
## [1] "Error pada data pelatihan:"
print(error_pelatihan)
##  [1]  0.020641742  0.029817441  0.038305477  0.046157387  0.053420844
##  [6]  0.060139949  0.066355499  0.072105232  0.077424057  0.082344269
## [11]  0.086895741  0.091106109  0.007238381  0.471863091 -0.204347914
## [16] -0.544260618  0.650287680  0.161859546  0.812650407 -1.141793233
## [21] -0.912068307 -0.450172024  0.733973466  0.724991803  0.400554721
## [26]  0.194389213  0.419311513  0.470037122 -0.168035890  0.140758039
## [31]  1.219198288  0.087853603  0.270777354  0.027547560  0.428236255
## [36]  0.316969988
absolute_error_pelatihan <- abs(error_pelatihan)
print("Absolute Error pada data pelatihan:")
## [1] "Absolute Error pada data pelatihan:"
print(absolute_error_pelatihan)
##  [1] 0.020641742 0.029817441 0.038305477 0.046157387 0.053420844 0.060139949
##  [7] 0.066355499 0.072105232 0.077424057 0.082344269 0.086895741 0.091106109
## [13] 0.007238381 0.471863091 0.204347914 0.544260618 0.650287680 0.161859546
## [19] 0.812650407 1.141793233 0.912068307 0.450172024 0.733973466 0.724991803
## [25] 0.400554721 0.194389213 0.419311513 0.470037122 0.168035890 0.140758039
## [31] 1.219198288 0.087853603 0.270777354 0.027547560 0.428236255 0.316969988
MAE_pelatihan <- mean(absolute_error_pelatihan)
print(paste("MAE Pelatihan:", MAE_pelatihan))
## [1] "MAE Pelatihan: 0.324552493402081"
MAPE_pelatihan <- 100 * MAE_pelatihan / mean(data_pelatihan$`Curah Hujan`)
print(paste("MAPE Pelatihan:", MAPE_pelatihan, "%"))
## [1] "MAPE Pelatihan: 6.51485057100715 %"
# Menggunakan model yang telah dilatih untuk memprediksi data pengujian
# Mengambil nilai Lt, Tt, dan St terakhir dari data pelatihan
Lt_pengujian <- tail(Lt, 1)
Tt_pengujian <- tail(Tt, 1)
St_pengujian <- tail(St, 12)  # Mengambil 12 nilai musiman terakhir

Prediksi_pengujian <- numeric()
for (i in 1:nrow(data_pengujian)) {
  st <- gamma * (data_pengujian$`Curah Hujan`[i] - Lt_pengujian) + (1 - gamma) * St_pengujian[(i - 1) %% 12 + 1]
  St_pengujian <- c(St_pengujian, st)
  
  lt <- alpha * (data_pengujian$`Curah Hujan`[i] - st) + (1 - alpha) * (Lt_pengujian + Tt_pengujian)
  Lt_pengujian <- lt
  
  tt <- beta * (Lt_pengujian - lt) + (1 - beta) * Tt_pengujian
  Tt_pengujian <- tt
  
  prediksi <- (Lt_pengujian + Tt_pengujian) + st
  Prediksi_pengujian <- c(Prediksi_pengujian, prediksi)
}

# Menampilkan Prediksi pada data pengujian
print("Prediksi pada data pengujian:")
## [1] "Prediksi pada data pengujian:"
print(Prediksi_pengujian)
##  [1] 6.610389 6.235401 5.822393 4.876237 5.048733 4.447990 2.194189 3.854886
##  [9] 4.094875 5.256911 5.547672 6.524506 6.585790 6.452941 5.595522 5.496508
## [17] 4.272674 3.905698 3.281874 3.112449 4.884516 3.274918 4.863088 5.640387
# Menghitung error pada data pengujian
error_pengujian <- data_pengujian$`Curah Hujan` - Prediksi_pengujian
print("Error pada data pengujian:")
## [1] "Error pada data pengujian:"
print(error_pengujian)
##  [1] -0.12422853  0.25683875 -0.35012191 -0.71735373  0.70066036  0.22483851
##  [7] -0.24827900 -0.11721607 -0.24472764  0.70923608  0.29005877  0.09423308
## [13] -0.05054898  0.55917380 -0.50177160 -0.03692222 -0.02417875 -0.21681824
## [19]  0.92281898 -0.80986397  0.60028132 -1.19547661 -0.20912718 -0.55898313
absolute_error_pengujian <- abs(error_pengujian)
print("Absolute Error pada data pengujian:")
## [1] "Absolute Error pada data pengujian:"
print(absolute_error_pengujian)
##  [1] 0.12422853 0.25683875 0.35012191 0.71735373 0.70066036 0.22483851
##  [7] 0.24827900 0.11721607 0.24472764 0.70923608 0.29005877 0.09423308
## [13] 0.05054898 0.55917380 0.50177160 0.03692222 0.02417875 0.21681824
## [19] 0.92281898 0.80986397 0.60028132 1.19547661 0.20912718 0.55898313
MAE_pengujian <- mean(absolute_error_pengujian)
print(paste("MAE Pengujian:", MAE_pengujian))
## [1] "MAE Pengujian: 0.406823216974582"
# Menghitung MAPE pada data pengujian
MAPE_pengujian <- mean(abs(error_pengujian) / abs(data_pengujian$`Curah Hujan`)) * 100
print(paste("MAPE Pengujian:", MAPE_pengujian, "%"))
## [1] "MAPE Pengujian: 10.5757396820473 %"
# Mengambil nilai terakhir dari Lt, Tt, dan St untuk prediksi
Lt_prediksi <- tail(Lt, 1)
Tt_prediksi <- tail(Tt, 1)
St_prediksi <- tail(St, 12)  # Mengambil 12 nilai musiman terakhir

# Inisialisasi vektor untuk menyimpan prediksi 12 bulan ke depan
Prediksi_tahunan <- numeric()

# Melakukan prediksi untuk 12 bulan ke depan
for (i in 1:12) {
  # Menggunakan indeks musiman yang berulang untuk setiap bulan
  st <- St_prediksi[(i - 1) %% 12 + 1]
  
  # Menghitung prediksi berdasarkan Lt, Tt, dan St terakhir
  prediksi <- (Lt_prediksi + (i * Tt_prediksi)) + st
  Prediksi_tahunan <- c(Prediksi_tahunan, prediksi)
}

# Menampilkan prediksi 12 bulan ke depan
print("Prediksi untuk 12 bulan ke depan:")
## [1] "Prediksi untuk 12 bulan ke depan:"
print(Prediksi_tahunan)
##  [1] 6.769940 6.012294 6.206254 5.668634 4.435498 4.270319 2.487838 4.034147
##  [9] 4.416204 4.613619 5.280744 6.437209
# Memuat paket yang diperlukan untuk plot
library(ggplot2)

# Membuat data frame untuk semua hasil (data aktual, prediksi pelatihan, pengujian, dan ramalan tahunan)
# Gabungkan data aktual dan prediksi dalam satu frame untuk visualisasi
data_aktual <- c(data_pelatihan$`Curah Hujan`, data_pengujian$`Curah Hujan`)
prediksi_akhir <- c(Prediksi, Prediksi_pengujian, Prediksi_tahunan)

# Menentukan rentang waktu untuk seluruh data (pelatihan, pengujian, dan prediksi)
waktu_aktual <- 1:length(data_aktual)
waktu_prediksi <- 1:length(prediksi_akhir)

# Membuat data frame untuk ggplot
plot_data <- data.frame(
  Waktu = c(waktu_aktual, (length(waktu_aktual) + 1):(length(waktu_aktual) + 12)),
  Curah_Hujan = c(data_aktual, rep(NA, 12)),  # Tambahkan NA untuk data aktual pada prediksi tahunan
  Prediksi = prediksi_akhir
)

# Plot data aktual dan prediksi
ggplot(plot_data, aes(x = Waktu)) +
  geom_line(aes(y = Curah_Hujan, color = "Data Aktual"), size = 1) +
  geom_line(aes(y = Prediksi, color = "Prediksi"), size = 1, linetype = "dashed") +
  labs(title = "Hasil Ramalan Curah Hujan",
       x = "Waktu (bulan)",
       y = "Curah Hujan",
       color = "Legenda") +
  scale_color_manual(values = c("Data Aktual" = "blue", "Prediksi" = "red")) +
  theme_minimal() +
  theme(legend.position = "bottom")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Invers transformasi logaritma pada data aktual, prediksi pengujian, dan prediksi tahunan
data_aktual_invers <- exp(c(data_pelatihan$`Curah Hujan`, data_pengujian$`Curah Hujan`))
Prediksi_invers <- exp(Prediksi)
Prediksi_pengujian_invers <- exp(Prediksi_pengujian)

# Prediksi tahunan
Lt_prediksi <- tail(Lt, 1)
Tt_prediksi <- tail(Tt, 1)
St_prediksi <- tail(St, 12)
Prediksi_tahunan <- numeric()

for (i in 1:12) {
  st <- St_prediksi[(i - 1) %% 12 + 1]
  prediksi <- (Lt_prediksi + (i * Tt_prediksi)) + st
  Prediksi_tahunan <- c(Prediksi_tahunan, prediksi)
}

# Invers transformasi logaritma pada prediksi tahunan
Prediksi_tahunan_invers <- exp(Prediksi_tahunan)

print("Prediksi untuk 12 bulan ke depan (Invers Transformasi):")
## [1] "Prediksi untuk 12 bulan ke depan (Invers Transformasi):"
print(Prediksi_tahunan_invers)
##  [1] 871.25948 408.41901 495.84056 289.63850  84.39410  71.54448  12.03523
##  [8]  56.49474  82.78148 100.84849 196.51593 624.66077
# Gabungkan data aktual dan prediksi dalam satu frame untuk visualisasi
prediksi_akhir <- c(Prediksi_invers, Prediksi_pengujian_invers, Prediksi_tahunan_invers)
waktu_aktual <- 1:length(data_aktual_invers)

# Membuat data frame untuk ggplot
plot_data <- data.frame(
  Waktu = c(waktu_aktual, (length(waktu_aktual) + 1):(length(waktu_aktual) + 12)),
  Curah_Hujan = c(data_aktual_invers, rep(NA, 12)),
  Prediksi = prediksi_akhir
)

# Plot data aktual dan prediksi
ggplot(plot_data, aes(x = Waktu)) +
  geom_line(aes(y = Curah_Hujan, color = "Data Aktual"), size = 1) +
  geom_line(aes(y = Prediksi, color = "Prediksi"), size = 1, linetype = "dashed") +
  labs(title = "Hasil Ramalan Curah Hujan",
       x = "Waktu (bulan)",
       y = "Curah Hujan",
       color = "Legenda") +
  scale_color_manual(values = c("Data Aktual" = "blue", "Prediksi" = "red")) +
  theme_minimal() +
  theme(legend.position = "bottom")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).