Analisis ARDL Sederhana: Tinggi Muka Air (TMA) dengan Curah Hujan (PC) dan Kelembapan Relatif (RH)
Disederhanakan dari sintaks Fajryanti Kusuma Wardani (Analisis ARDL SKRIPSI)
Rio Faudzan Pritajaya
30-09-2026
Peubah yang digunakan
- TMA = Tinggi Muka Air Pintu Air Manggarai (peubah respon)
- PC = Curah Hujan Kota Bogor (peubah penjelas 1)
- RH = Kelembapan Relatif Kota Bogor (peubah penjelas 2)
Tujuan: mengidentifikasi hubungan jangka pendek dan jangka panjang antara PC dan RH terhadap TMA, serta meramal TMA harian periode Januari-Maret 2025.
Catatan: peubah TM dan WS dari analisis asli tidak dipakai agar model lebih sederhana. PC dan RH dipilih karena keduanya memiliki korelasi paling kuat dengan TMA pada matriks korelasi analisis asli (PC = 0,39 dan RH = 0,53).
1 Library
library(readxl)
library(dplyr)
library(ggplot2)
library(moments)
library(corrplot)
library(forecast)
library(tseries)
library(dynlm)
library(dLagM)
library(lmtest)
library(car)
library(ARDL)
library(tidyr)
library(patchwork)
library(dynamac)2 Import Data
Letakkan file Data_Harian_2022_2024.xlsx dan
Data_Jan25_Mar25.xlsx di folder yang sama dengan file Rmd
ini (atau ubah path di bawah).
dataset <- read_excel("Data_Harian_2022_2024.xlsx")
dataset <- dataset[, c("Date", "TMA", "PC", "RH")] # hanya peubah yang dipakai
str(dataset)## tibble [1,096 × 4] (S3: tbl_df/tbl/data.frame)
## $ Date: POSIXct[1:1096], format: "2022-01-01" "2022-01-02" ...
## $ TMA : num [1:1096] 595 588 580 576 574 ...
## $ PC : num [1:1096] 0.85 1.79 0.4 0.22 5.7 ...
## $ RH : num [1:1096] 86.4 88.1 83.2 80.3 83.2 ...Data harian periode 1 Jan 2022 - 31 Des 2024.
3 Eksplorasi Data
3.1 Statistika Deskriptif
summary(dataset)## Date TMA PC
## Min. :2022-01-01 00:00:00 Min. :555.4 Min. : 0.0000
## 1st Qu.:2022-10-01 18:00:00 1st Qu.:597.7 1st Qu.: 0.7375
## Median :2023-07-02 12:00:00 Median :620.0 Median : 4.0000
## Mean :2023-07-02 12:00:00 Mean :625.3 Mean : 6.5700
## 3rd Qu.:2024-04-01 06:00:00 3rd Qu.:650.1 3rd Qu.: 9.1325
## Max. :2024-12-31 00:00:00 Max. :781.2 Max. :145.0300
## RH
## Min. :67.03
## 1st Qu.:82.92
## Median :86.29
## Mean :85.29
## 3rd Qu.:88.53
## Max. :94.31deskriptif <- data.frame(
Peubah = c("TMA", "PC", "RH"),
SD = c(sd(dataset$TMA), sd(dataset$PC), sd(dataset$RH)),
Skewness = c(skewness(dataset$TMA), skewness(dataset$PC), skewness(dataset$RH)),
Kurtosis = c(kurtosis(dataset$TMA), kurtosis(dataset$PC), kurtosis(dataset$RH))
)
deskriptif## Peubah SD Skewness Kurtosis
## 1 TMA 39.449620 0.5304754 2.984404
## 2 PC 8.838259 5.0708708 62.196788
## 3 RH 4.710888 -0.9881591 4.0250323.2 Plot Time Series
data_long <- dataset %>%
pivot_longer(cols = c(TMA, PC, RH),
names_to = "Variabel",
values_to = "Nilai")
ggplot(data_long, aes(x = Date, y = Nilai)) +
geom_line() +
facet_wrap(~ Variabel, scales = "free_y", ncol = 1) +
labs(title = "Time Series Plot of TMA, PC, dan RH",
x = "Date", y = "Value") +
theme_minimal()
3.3 Matriks Korelasi
cor_matrix <- cor(dataset[, c("TMA", "PC", "RH")], use = "complete.obs")
cor_matrix## TMA PC RH
## TMA 1.0000000 0.3945932 0.5250838
## PC 0.3945932 1.0000000 0.4196228
## RH 0.5250838 0.4196228 1.0000000corrplot(cor_matrix, method = "color", addCoef.col = "black",
tl.col = "black", tl.srt = 50, number.cex = 0.8,
title = "Matriks Korelasi",
mar = c(0, 0, 1, 0))
4 Data Splitting
Data latih 936 hari pertama dan data uji 160 hari terakhir (25 Jul 2024 - 31 Des 2024), sama seperti analisis asli.
total_rows <- nrow(dataset)
test_size <- 160
train_index <- 1:(total_rows - test_size)
test_index <- (total_rows - test_size + 1):total_rows
train_2225 <- dataset[train_index, ]
test_2225 <- dataset[test_index, ]
str(train_2225)## tibble [936 × 4] (S3: tbl_df/tbl/data.frame)
## $ Date: POSIXct[1:936], format: "2022-01-01" "2022-01-02" ...
## $ TMA : num [1:936] 595 588 580 576 574 ...
## $ PC : num [1:936] 0.85 1.79 0.4 0.22 5.7 ...
## $ RH : num [1:936] 86.4 88.1 83.2 80.3 83.2 ...str(test_2225)## tibble [160 × 4] (S3: tbl_df/tbl/data.frame)
## $ Date: POSIXct[1:160], format: "2024-07-25" "2024-07-26" ...
## $ TMA : num [1:160] 600 600 600 600 600 ...
## $ PC : num [1:160] 0.01 0.01 0 0 0 ...
## $ RH : num [1:160] 85.3 84.4 82.3 82.3 81.6 ...train_plot <- train_2225 %>% mutate(DataType = "Train")
test_plot <- test_2225 %>% mutate(DataType = "Test")
combined_data <- bind_rows(train_plot, test_plot)
data_long <- combined_data %>%
pivot_longer(cols = c(TMA, PC, RH),
names_to = "Variabel",
values_to = "Nilai")
ggplot(data_long, aes(x = Date, y = Nilai, color = DataType)) +
geom_line() +
facet_wrap(~ Variabel, scales = "free_y", ncol = 1) +
scale_color_manual(values = c("Train" = "#00BFC4", "Test" = "#F8766D")) +
labs(title = "Time Series (Train vs. Test)",
x = "Date", y = "Value", color = "Data Type") +
theme_minimal()
5 Kestasioneran Data
5.1 Plot ACF dan PACF
acf_TMA <- ggAcf(train_2225$TMA) + labs(title = "ACF TMA")
acf_PC <- ggAcf(train_2225$PC) + labs(title = "ACF PC")
acf_RH <- ggAcf(train_2225$RH) + labs(title = "ACF RH")
pacf_TMA <- ggPacf(train_2225$TMA) + labs(title = "PACF TMA")
pacf_PC <- ggPacf(train_2225$PC) + labs(title = "PACF PC")
pacf_RH <- ggPacf(train_2225$RH) + labs(title = "PACF RH")
(acf_TMA + acf_PC + acf_RH) / (pacf_TMA + pacf_PC + pacf_RH)
5.2 Uji ADF
H0: data tidak stasioner; H1: data stasioner.
adf.test(train_2225$TMA)##
## Augmented Dickey-Fuller Test
##
## data: train_2225$TMA
## Dickey-Fuller = -4.0015, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationaryadf.test(train_2225$PC)##
## Augmented Dickey-Fuller Test
##
## data: train_2225$PC
## Dickey-Fuller = -7.5589, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationaryadf.test(train_2225$RH)##
## Augmented Dickey-Fuller Test
##
## data: train_2225$RH
## Dickey-Fuller = -3.3561, Lag order = 9, p-value = 0.06101
## alternative hypothesis: stationary5.3 Plot CCF (data level)
ccf_result1 <- ggCcf(train_2225$TMA, train_2225$PC, lag.max = 9) +
ggtitle("CCF antara TMA dan PC")
ccf_result2 <- ggCcf(train_2225$TMA, train_2225$RH, lag.max = 9) +
ggtitle("CCF antara TMA dan RH")
ccf_result1 + ccf_result2
6 First Difference
6.1 Differencing
DTMA <- diff(train_2225$TMA, differences = 1)
DPC <- diff(train_2225$PC, differences = 1)
DRH <- diff(train_2225$RH, differences = 1)
df <- data.frame(DTMA, DPC, DRH)
df.ts <- ts(df)
head(df)## DTMA DPC DRH
## 1 -7.708333 0.94 1.75
## 2 -7.500000 -1.39 -4.93
## 3 -4.375000 -0.18 -2.91
## 4 -1.458333 5.48 2.87
## 5 9.583333 3.75 2.39
## 6 -10.208333 6.72 4.066.2 Plot ACF dan PACF (data diferensi)
PACF DTMA digunakan untuk menentukan lag optimum peubah respon.
acf_DTMA <- ggAcf(DTMA) + labs(title = "ACF DTMA")
acf_DPC <- ggAcf(DPC) + labs(title = "ACF DPC")
acf_DRH <- ggAcf(DRH) + labs(title = "ACF DRH")
pacf_DTMA <- ggPacf(DTMA) + labs(title = "PACF DTMA")
pacf_DPC <- ggPacf(DPC) + labs(title = "PACF DPC")
pacf_DRH <- ggPacf(DRH) + labs(title = "PACF DRH")
(acf_DTMA + acf_DPC + acf_DRH) / (pacf_DTMA + pacf_DPC + pacf_DRH)
6.3 Uji ADF (data diferensi)
adf.test(DTMA)##
## Augmented Dickey-Fuller Test
##
## data: DTMA
## Dickey-Fuller = -12.518, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationaryadf.test(DPC)##
## Augmented Dickey-Fuller Test
##
## data: DPC
## Dickey-Fuller = -14.984, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationaryadf.test(DRH)##
## Augmented Dickey-Fuller Test
##
## data: DRH
## Dickey-Fuller = -13.91, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary6.4 Plot CCF (data diferensi)
Digunakan untuk menentukan lag optimum peubah penjelas.
ccf_d1 <- ggCcf(DTMA, DPC, lag.max = 9) + ggtitle("CCF antara DTMA dan DPC")
ccf_d2 <- ggCcf(DTMA, DRH, lag.max = 9) + ggtitle("CCF antara DTMA dan DRH")
ccf_d1 + ccf_d2
ccf_r1 <- ccf(DTMA, DPC, lag.max = 9, plot = FALSE)
ccf_r2 <- ccf(DTMA, DRH, lag.max = 9, plot = FALSE)
ccf_all_df <- rbind(
data.frame(Lag = ccf_r1$lag, Correlation = ccf_r1$acf, Pair = "DTMA vs DPC"),
data.frame(Lag = ccf_r2$lag, Correlation = ccf_r2$acf, Pair = "DTMA vs DRH")
)
ccf_all_df## Lag Correlation Pair
## 1 -9 0.0529586274 DTMA vs DPC
## 2 -8 -0.0550552403 DTMA vs DPC
## 3 -7 -0.0075649185 DTMA vs DPC
## 4 -6 0.0272775490 DTMA vs DPC
## 5 -5 -0.0260687707 DTMA vs DPC
## 6 -4 0.0095200923 DTMA vs DPC
## 7 -3 -0.0985237382 DTMA vs DPC
## 8 -2 0.1056947588 DTMA vs DPC
## 9 -1 -0.2068571316 DTMA vs DPC
## 10 0 0.1086869199 DTMA vs DPC
## 11 1 0.2580985831 DTMA vs DPC
## 12 2 -0.1309036731 DTMA vs DPC
## 13 3 -0.0120556745 DTMA vs DPC
## 14 4 0.0379023799 DTMA vs DPC
## 15 5 0.0102910918 DTMA vs DPC
## 16 6 -0.0779871914 DTMA vs DPC
## 17 7 0.0489950877 DTMA vs DPC
## 18 8 -0.0192893364 DTMA vs DPC
## 19 9 0.0083648908 DTMA vs DPC
## 20 -9 0.0273578093 DTMA vs DRH
## 21 -8 -0.0421464810 DTMA vs DRH
## 22 -7 0.0127058736 DTMA vs DRH
## 23 -6 0.0392240962 DTMA vs DRH
## 24 -5 -0.0190854520 DTMA vs DRH
## 25 -4 -0.0046866215 DTMA vs DRH
## 26 -3 -0.0759843923 DTMA vs DRH
## 27 -2 0.0320354704 DTMA vs DRH
## 28 -1 -0.0338906315 DTMA vs DRH
## 29 0 0.0936044467 DTMA vs DRH
## 30 1 0.0960394721 DTMA vs DRH
## 31 2 -0.0636944028 DTMA vs DRH
## 32 3 0.0459606309 DTMA vs DRH
## 33 4 0.0052088145 DTMA vs DRH
## 34 5 -0.0575980410 DTMA vs DRH
## 35 6 0.0126319156 DTMA vs DRH
## 36 7 -0.0261257009 DTMA vs DRH
## 37 8 -0.0003282647 DTMA vs DRH
## 38 9 0.0054043507 DTMA vs DRH7 Dynardl (dynamac)
7.1 Uji Kointegrasi (Bounds Test)
Lag maksimum mengikuti analisis asli (max.p = 1,
max.q = 9).
modelbound <- ardlBound(data = train_2225, formula = TMA ~ PC + RH,
case = 2, max.p = 1, max.q = 9)##
## Orders being calculated with max.p = 1 and max.q = 9 ...
##
## Autoregressive order: 10 and p-orders: 2 1
## ------------------------------------------------------
##
## Breusch-Godfrey Test for the autocorrelation in residuals:
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelFull$model
## LM test = 1.5157, df1 = 1, df2 = 909, p-value = 0.2186
##
## ------------------------------------------------------
##
## Ljung-Box Test for the autocorrelation in residuals:
##
## Box-Ljung test
##
## data: res
## X-squared = 0.023902, df = 1, p-value = 0.8771
##
## ------------------------------------------------------
##
## Breusch-Pagan Test for the homoskedasticity of residuals:
##
## studentized Breusch-Pagan test
##
## data: modelFull$model
## BP = 199.46, df = 15, p-value < 2.2e-16
##
## The p-value of Breusch-Pagan test for the homoskedasticity of residuals: 2.734809e-34 < 0.05!
## ------------------------------------------------------
##
## Shapiro-Wilk test of normality of residuals:
##
## Shapiro-Wilk normality test
##
## data: modelFull$model$residual
## W = 0.88703, p-value < 2.2e-16
##
## The p-value of Shapiro-Wilk test normality of residuals: 1.908888e-25 < 0.05!
## ------------------------------------------------------
##
## PESARAN, SHIN AND SMITH (2001) COINTEGRATION TEST
##
## Observations: 935
## Number of Regressors (k): 2
## Case: 2
##
## ------------------------------------------------------
## - F-test -
## ------------------------------------------------------
## <------- I(0) ------------ I(1) ----->
## 10% critical value 2.63 3.35
## 5% critical value 3.1 3.87
## 1% critical value 4.13 5
##
##
## F-statistic = 31.6583766474097
##
## ------------------------------------------------------
## F-statistic note: Asymptotic critical values used.
##
## ------------------------------------------------------
##
## Ramsey's RESET Test for model specification:
##
## RESET test
##
## data: modelECM$model
## RESET = 14.337, df1 = 1, df2 = 912, p-value = 0.0001628
##
## the p-value of RESET test: 0.0001628165 < 0.05!
## ------------------------------------------------------## ------------------------------------------------------
## Error Correction Model Output:
##
## Time series regression with "ts" data:
## Start = 10, End = 935
##
## Call:
## dynlm(formula = as.formula(model.text), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -130.086 -11.717 -2.635 6.745 136.100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ec.1 -0.17775 0.01577 -11.272 < 2e-16 ***
## dPC.t 0.67150 0.09149 7.340 4.73e-13 ***
## dPC.1 0.35361 0.09905 3.570 0.000376 ***
## dRH.t 0.71637 0.27437 2.611 0.009176 **
## dTMA.1 -0.42205 0.03010 -14.023 < 2e-16 ***
## dTMA.2 -0.27043 0.03320 -8.146 1.23e-15 ***
## dTMA.3 -0.19856 0.03398 -5.844 7.10e-09 ***
## dTMA.4 -0.12458 0.03440 -3.621 0.000309 ***
## dTMA.5 -0.14675 0.03396 -4.322 1.72e-05 ***
## dTMA.6 -0.12230 0.03394 -3.603 0.000332 ***
## dTMA.7 -0.09617 0.03356 -2.865 0.004262 **
## dTMA.8 -0.04125 0.03250 -1.269 0.204753
## dTMA.9 0.02274 0.02935 0.775 0.438574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.23 on 913 degrees of freedom
## Multiple R-squared: 0.3707, Adjusted R-squared: 0.3618
## F-statistic: 41.37 on 13 and 913 DF, p-value: < 2.2e-16
##
## ------------------------------------------------------
## Long-run coefficients:
## TMA.1 PC.1 RH.1
## -0.1777542 1.2022135 0.3960311
## 7.2 Model
7.2.1 Jangka Pendek
modelardl <- dynardl(
DTMA ~ DPC + DRH,
data = df,
lags = list(
"DTMA" = 1:9,
"DPC" = 1,
"DRH" = 1
),
ec = TRUE
)## [1] "Error correction (EC) specified; dependent variable to be run in differences."summary(modelardl)##
## Call:
## lm(formula = as.formula(paste(paste(dvnamelist), "~", paste(colnames(IVs),
## collapse = "+"), collapse = " ")))
##
## Residuals:
## Min 1Q Median 3Q Max
## -82.905 -11.899 -3.051 7.310 140.121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07954 0.75371 0.106 0.915980
## l.1.DTMA -1.48549 0.03170 -46.861 < 2e-16 ***
## l.2.DTMA -0.26216 0.03553 -7.378 3.60e-13 ***
## l.3.DTMA -0.20926 0.03615 -5.788 9.76e-09 ***
## l.4.DTMA -0.12204 0.03682 -3.315 0.000954 ***
## l.5.DTMA -0.16120 0.03649 -4.418 1.12e-05 ***
## l.6.DTMA -0.12655 0.03656 -3.461 0.000563 ***
## l.7.DTMA -0.10256 0.03618 -2.835 0.004685 **
## l.8.DTMA -0.05487 0.03493 -1.571 0.116531
## l.9.DTMA 0.02288 0.03169 0.722 0.470456
## l.1.DPC 0.72043 0.08740 8.243 5.79e-16 ***
## l.1.DRH 0.62171 0.29720 2.092 0.036725 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.94 on 914 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.726, Adjusted R-squared: 0.7227
## F-statistic: 220.1 on 11 and 914 DF, p-value: < 2.2e-16AIC(modelardl$model)## [1] 8443.5117.2.2 Jangka Panjang
modelardll <- dynardl(
formula = TMA ~ PC + RH,
data = train_2225,
levels = c("PC", "RH"),
ec = FALSE
)## [1] "Dependent variable to be run in levels."summary(modelardll)##
## Call:
## lm(formula = as.formula(paste(paste(dvnamelist), "~", paste(colnames(IVs),
## collapse = "+"), collapse = " ")))
##
## Residuals:
## Min 1Q Median 3Q Max
## -81.035 -12.390 -3.133 8.906 138.920
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 99.84097 16.21708 6.157 1.10e-09 ***
## l.1.TMA 0.69406 0.02133 32.540 < 2e-16 ***
## PC 0.83239 0.09610 8.662 < 2e-16 ***
## RH 1.00643 0.19987 5.035 5.73e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.34 on 931 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6785, Adjusted R-squared: 0.6775
## F-statistic: 655 on 3 and 931 DF, p-value: < 2.2e-16AIC(modelardll$model)## [1] 8550.5057.3 Uji Asumsi
7.3.1 Autokorelasi residual (Box-Ljung Test)
H0: tidak ada autokorelasi pada residual.
res1 <- residuals(modelardl$model)
Box.test(res1, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: res1
## X-squared = 19.038, df = 20, p-value = 0.51937.3.2 Normalitas (Jarque-Bera Test)
H0: residual berdistribusi normal.
jarque.bera.test(res1)##
## Jarque Bera Test
##
## data: res1
## X-squared = 1281.6, df = 2, p-value < 2.2e-167.3.3 Heteroskedastisitas (Breusch-Pagan Test)
H0: ragam residual konstan (homoskedastis).
fitted_val1 <- fitted(modelardl$model)
bptest(res1 ~ fitted_val1)##
## studentized Breusch-Pagan test
##
## data: res1 ~ fitted_val1
## BP = 1.3632, df = 1, p-value = 0.2437.3.4 Multikolinearitas (VIF)
vif(modelardl$model)## l.1.DTMA l.2.DTMA l.3.DTMA l.4.DTMA l.5.DTMA l.6.DTMA l.7.DTMA l.8.DTMA
## 1.249782 1.570812 1.626521 1.687545 1.657613 1.664197 1.629196 1.518988
## l.9.DTMA l.1.DPC l.1.DRH
## 1.250541 1.153269 1.0971027.4 Fungsi bantu evaluasi
evaluasi <- function(aktual, prediksi) {
MAE <- mean(abs(prediksi - aktual), na.rm = TRUE)
RMSE <- sqrt(mean((prediksi - aktual)^2, na.rm = TRUE))
MAPE <- mean(abs((aktual - prediksi) / aktual), na.rm = TRUE) * 100
data.frame(MAE = MAE, RMSE = RMSE, MAPE_persen = MAPE)
}
# Koefisien model jangka pendek (dipakai untuk peramalan, seperti analisis asli)
coefs <- summary(modelardl)$coefficients
intercept <- coefs["(Intercept)", "Estimate"]
b1 <- coefs["l.1.DPC", "Estimate"]
b2 <- coefs["l.1.DRH", "Estimate"]
coefs## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07953787 0.75371316 0.1055280 9.159800e-01
## l.1.DTMA -1.48549061 0.03170006 -46.8608170 2.895372e-245
## l.2.DTMA -0.26216448 0.03553161 -7.3783460 3.597338e-13
## l.3.DTMA -0.20926129 0.03615189 -5.7883922 9.762272e-09
## l.4.DTMA -0.12204373 0.03682087 -3.3145262 9.540847e-04
## l.5.DTMA -0.16120030 0.03649032 -4.4176182 1.117423e-05
## l.6.DTMA -0.12654844 0.03656265 -3.4611401 5.627992e-04
## l.7.DTMA -0.10255553 0.03617559 -2.8349377 4.684534e-03
## l.8.DTMA -0.05487372 0.03492939 -1.5709897 1.165311e-01
## l.9.DTMA 0.02288324 0.03169259 0.7220373 4.704562e-01
## l.1.DPC 0.72043361 0.08739820 8.2431169 5.793166e-16
## l.1.DRH 0.62171161 0.29720344 2.0918722 3.672540e-028 Forecasting
8.1 Data Training
train_baru <- data.frame(DPC, DRH)
train_baru$DTMA_pred <- intercept + b1 * train_baru$DPC + b2 * train_baru$DRH
# TMA ke-n = TMA ke-(n-1) + delta TMA ke-n
train_baru$TMA_pred <- NA
train_baru$TMA_pred[1] <- train_2225$TMA[1] + train_baru$DTMA_pred[1]
for (i in 2:nrow(train_baru)) {
train_baru$TMA_pred[i] <- train_baru$TMA_pred[i - 1] + train_baru$DTMA_pred[i]
}
TMA_actual <- train_2225$TMA[-1] # karena ada differencing
TMA_pred <- train_baru$TMA_pred
plot(TMA_actual, type = "l", col = "black", main = "Train: Actual vs Predicted TMA")
lines(TMA_pred, col = "blue")
legend("topleft", legend = c("Actual", "Predicted"), col = c("black", "blue"), lty = 1)
evaluasi(TMA_actual, TMA_pred)## MAE RMSE MAPE_persen
## 1 31.66414 41.07022 5.1303718.2 Data Testing
TDPC <- diff(test_2225$PC, differences = 1)
TDRH <- diff(test_2225$RH, differences = 1)
test_baru <- data.frame(DPC = TDPC, DRH = TDRH)
test_baru$DTMA_pred <- intercept + b1 * test_baru$DPC + b2 * test_baru$DRH
# Nilai awal: TMA pertama pada data uji
test_baru$TMA_pred <- NA
test_baru$TMA_pred[1] <- test_2225$TMA[1] + test_baru$DTMA_pred[1]
for (i in 2:nrow(test_baru)) {
test_baru$TMA_pred[i] <- test_baru$TMA_pred[i - 1] + test_baru$DTMA_pred[i]
}
TMA_actual <- test_2225$TMA[-1]
TMA_pred <- test_baru$TMA_pred
date_seq <- seq.Date(from = as.Date(test_2225$Date[2]), by = "day",
length.out = length(TMA_pred))
plot(date_seq, TMA_actual, type = "l", col = "black",
main = "Actual vs Predicted TMA (Data Uji)", xlab = "Tanggal", ylab = "TMA",
ylim = range(c(TMA_actual, TMA_pred)))
lines(date_seq, TMA_pred, col = "blue")
legend("topleft", legend = c("Actual", "Predicted"), col = c("black", "blue"), lty = 1)
evaluasi(TMA_actual, TMA_pred)## MAE RMSE MAPE_persen
## 1 18.02884 24.49108 2.8194288.3 Prediksi Januari - Maret 2025
data_prediksi <- read_excel("Data_Jan25_Mar25.xlsx")
data_prediksi$Date <- as.Date(data_prediksi$Date)
str(data_prediksi)## tibble [90 × 6] (S3: tbl_df/tbl/data.frame)
## $ Date: Date[1:90], format: "2025-01-01" "2025-01-02" ...
## $ TMA : num [1:90] 641 640 640 640 640 ...
## $ PC : num [1:90] 5.27 2.45 6.22 6.03 2.63 ...
## $ TM : num [1:90] 24.6 25 25.3 25.1 24.9 ...
## $ RH : num [1:90] 88 85.2 85.5 86.4 84.6 ...
## $ WS : num [1:90] 2.56 1.85 1.75 2.56 2.38 1.82 1.87 2.22 1.64 1.64 ...DPC_25 <- diff(data_prediksi$PC, differences = 1)
DRH_25 <- diff(data_prediksi$RH, differences = 1)
test_25 <- data.frame(DPC = DPC_25, DRH = DRH_25)
test_25$DTMA_pred <- intercept + b1 * test_25$DPC + b2 * test_25$DRH
test_25$TMA_pred <- NA
test_25$TMA_pred[1] <- data_prediksi$TMA[1] + test_25$DTMA_pred[1]
for (i in 2:nrow(test_25)) {
test_25$TMA_pred[i] <- test_25$TMA_pred[i - 1] + test_25$DTMA_pred[i]
}
hasil_prediksi <- data.frame(
date = data_prediksi$Date[-1], # karena differencing
TMA_actual = data_prediksi$TMA[-1],
TMA_pred = test_25$TMA_pred
)
head(hasil_prediksi, 10)## date TMA_actual TMA_pred
## 1 2025-01-02 640.0000 637.0512
## 2 2025-01-03 640.0000 640.0706
## 3 2025-01-04 640.0000 640.5603
## 4 2025-01-05 640.0000 637.0962
## 5 2025-01-06 640.0000 653.5259
## 6 2025-01-07 640.0000 646.5268
## 7 2025-01-08 641.4583 655.8792
## 8 2025-01-09 666.2500 654.3403
## 9 2025-01-10 642.5000 643.4396
## 10 2025-01-11 637.5000 648.2091plot(hasil_prediksi$date, hasil_prediksi$TMA_actual, type = "l", col = "black",
ylab = "TMA", xlab = "Tanggal", main = "TMA Aktual vs Prediksi (Jan-Mar 2025)",
ylim = range(c(hasil_prediksi$TMA_actual, hasil_prediksi$TMA_pred)))
lines(hasil_prediksi$date, hasil_prediksi$TMA_pred, col = "blue")
legend("topleft", legend = c("Aktual", "Prediksi"), col = c("black", "blue"), lty = 1)
evaluasi(hasil_prediksi$TMA_actual, hasil_prediksi$TMA_pred)## MAE RMSE MAPE_persen
## 1 16.18396 25.94375 2.386003