UAS TOPIK STATISTIKA I

SAYYID RAFIF

Data

# Membaca data
data <- read_excel("~/topik statistika/UAS TOPIK STAT/data harian BMKG.xlsx")
data$TANGGAL <- as.Date(data$TANGGAL, format = "%d-%m-%Y")
summary(data)

##     TANGGAL           jumlah_curah_hujan
##  Min.   :2025-05-01   Min.   :  0.000   
##  1st Qu.:2025-07-31   1st Qu.:  0.000   
##  Median :2025-10-30   Median :  0.200   
##  Mean   :2025-10-30   Mean   :  9.128   
##  3rd Qu.:2026-01-29   3rd Qu.:  7.800   
##  Max.   :2026-04-30   Max.   :213.400

Ubah data ke time series

# Mengubah kolom 'jumlah_curah_hujan' menjadi objek time series
data.ts <- ts(data$jumlah_curah_hujan, start = 1)
# Menampilkan hasil time series
print(head(data.ts, 35))

## Time Series:
## Start = 1 
## End = 35 
## Frequency = 1 
##  [1]  0.0  0.0  3.0 42.2  0.0  0.0  3.0  3.7  0.0  0.0  0.8 27.0  0.0 55.2  0.0
## [16] 37.2  2.1  5.7  1.0  0.0  0.0 40.4  4.9  7.5  0.5 33.0 21.3  0.0  0.0 20.0
## [31]  0.0  1.6  0.0  0.0  0.0

Uji Pola Musiman

plot(data$TANGGAL, data.ts, 
     type = "l",           
     xaxt = "n",          
     xlab = "Bulan", 
     ylab = "Jumlah Curah Hujan (Yt)", 
     col = "blue",         
     lwd = 1.5,
     main = "Time Series Plot Curah Hujan Harian (Mei 2025 - April 2026)")
titik_bulan <- seq(from = min(data$TANGGAL), to = max(data$TANGGAL), by = "1 month")
axis.Date(side = 1, at = titik_bulan, format = "%b %Y", las = 2, cex.axis = 0.8)

data$Bulan_Tahun <- format(data$TANGGAL, "%Y-%m")
data_bulanan <- aggregate(jumlah_curah_hujan ~ Bulan_Tahun, data = data, FUN = sum)
label_bulan <- format(as.Date(paste0(data_bulanan$Bulan_Tahun, "-01")), "%b %Y")
barplot(data_bulanan$jumlah_curah_hujan, 
        names.arg = label_bulan, 
        las = 2,              
        col = "skyblue", 
        main = "Distribusi Musiman Jumlah Curah Hujan (Mei 2025 - April 2026)",
        ylab = "Jumlah Curah Hujan (Yt)", 
        cex.names = 0.8)

berdasarkan plot time series dan diagram distribusi musiman jumlah curah hujan tidak terlihat pola musiman pada data sehingga tidak representatif jika menggunakan model SARIMA maka analisis selajutnya akan dilakukan menggunakan model AR,MA,ARMA, dan ARIMA

Splitting data

80% Data Train = 292 Hari, 20% Data Test = 73 Hari

train.ts <- window(data.ts, end = 292)
test.ts  <- window(data.ts, start = 293)

Uji kestasioneran varians data

# Stasioneritas Varians (Box-Cox)
lambda_opt <- BoxCox.lambda(train.ts)

## Warning in guerrero(x, lower, upper): Guerrero's method for selecting a Box-Cox
## parameter (lambda) is given for strictly positive data.

cat("Nilai Lambda Box-Cox:", round(lambda_opt, 4), "\n")

## Nilai Lambda Box-Cox: 0.0283

Uji kestasioneran rataan data

Mengidentifikasi kestasioneran data dilakukan dengan uji ADF. Uji Augmented Dickey-Fuller (ADF) adalah salah satu alat identifikasi untuk menguji kestasioneran data dan memiliki hipotesis sebagai berikut:

\[H_0: \text{Data tidak stasioner (mengandung unit root)}\] \[H_1: \text{Data stasioner (tidak mengandung unit root)}\]

cat("\n--- Uji ADF Data Train ---\n")

## 
## --- Uji ADF Data Train ---

print(adf.test(train.ts))

## Warning in adf.test(train.ts): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -4.6802, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5\%\) di atas dapat diketahui bersama bahwa \(p-value = 0.01 < \alpha = 0.05\) sehingga keputusan yang diambil adalah Tolak \(H_0\) artinya data yang digunakan merupakan data yang stasioner.

Identifikasi kandidat model ARIMA

Identifikasi kandidat model diperoleh berdasarkan nilai p dan q dimana nilai d=0

acf(train.ts, main="ACF Data Train")

Berdasarkan plot ACF terlihat bahwa plot cuts off setelah lag ke-10 sehingga kandidat model yang diperoleh adalah ARMA(0,0,10)

pacf(train.ts, main="PACF Data Train")

Berdasarkan plot PACF di atas kandidat model yang diperoleh adalah ARMA(8,0,0)

cat("\n--- Extended ACF (EACF) ---\n")

## 
## --- Extended ACF (EACF) ---

eacf(train.ts)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o x o o x o x o  o  o  o 
## 1 x o o o x o o x o x o  o  o  o 
## 2 x x o o x o x o o x o  o  o  o 
## 3 x o o o x o x o o x o  o  o  o 
## 4 x x x x x o x o o x o  o  o  o 
## 5 o o x o x o o x o o o  o  o  o 
## 6 x o x o x o o x o o o  o  o  o 
## 7 x o x x x o o x x o o  o  o  o

Berdasarkan hasil EACF di atas kandidat model yang diperoleh adalah ARMA(1,0,1),ARMA(0,0,1), dan ARMA(2,0,2)

Dibandingkan dengan auto.arima

auto.arima(train.ts, trace = TRUE)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2) with drift         : 2605.718
##  ARIMA(0,1,0) with drift         : 2779.127
##  ARIMA(1,1,0) with drift         : 2701.139
##  ARIMA(0,1,1) with drift         : 2599.393
##  ARIMA(0,1,0)                    : 2777.099
##  ARIMA(1,1,1) with drift         : 2602.462
##  ARIMA(0,1,2) with drift         : 2601.448
##  ARIMA(1,1,2) with drift         : Inf
##  ARIMA(0,1,1)                    : 2597.572
##  ARIMA(1,1,1)                    : 2600.601
##  ARIMA(0,1,2)                    : 2599.598
##  ARIMA(1,1,0)                    : 2699.097
##  ARIMA(1,1,2)                    : Inf
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,1,1)                    : 2605.125
## 
##  Best model: ARIMA(0,1,1)

## Series: train.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.9280
## s.e.   0.0334
## 
## sigma^2 = 444.6:  log likelihood = -1300.54
## AIC=2605.08   AICc=2605.12   BIC=2612.43

Berdasarkan hasil dari auto arima diperoleh best model yaitu ARIMA(0,1,1)

Pemilihan kandidat model terbaik

# ARMA(0,0,10)
arma0010 <- arima(train.ts, order = c(0,0,10), include.mean = TRUE, method = "ML")
arma0010

## 
## Call:
## arima(x = train.ts, order = c(0, 0, 10), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ma1     ma2      ma3      ma4     ma5     ma6      ma7     ma8     ma9
##       0.0660  0.0485  -0.0130  -0.0110  0.3067  0.0310  -0.0317  0.1366  0.0492
## s.e.  0.0584  0.0583   0.0588   0.0582  0.0578  0.0559   0.0571  0.0532  0.0654
##         ma10  intercept
##       0.1457     9.2496
## s.e.  0.0561     1.9861
## 
## sigma^2 estimated as 392.6:  log likelihood = -1286.78,  aic = 2595.56

# ARMA(8,0,0)
arma800 <- arima(train.ts, order = c(8,0,0), include.mean = TRUE, method = "ML")
arma800

## 
## Call:
## arima(x = train.ts, order = c(8, 0, 0), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1     ar2      ar3     ar4     ar5     ar6      ar7     ar8
##       0.0665  0.0579  -0.0221  0.0094  0.3154  0.0181  -0.0414  0.1295
## s.e.  0.0579  0.0579   0.0579  0.0552  0.0552  0.0583   0.0581  0.0579
##       intercept
##          9.2212
## s.e.     2.4461
## 
## sigma^2 estimated as 394.9:  log likelihood = -1287.56,  aic = 2593.13

# ARMA(1,0,1)
arma101 <- arima(train.ts, order = c(1,0,1), include.mean = TRUE, method = "ML")
arma101

## 
## Call:
## arima(x = train.ts, order = c(1, 0, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.9407  -0.8654     9.2524
## s.e.  0.0381   0.0539     2.6845
## 
## sigma^2 estimated as 434.1:  log likelihood = -1301.1,  aic = 2608.2

# ARMA(0,0,1)
arma001 <- arima(train.ts, order = c(0,0,1), include.mean = TRUE, method = "ML")
arma001

## 
## Call:
## arima(x = train.ts, order = c(0, 0, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ma1  intercept
##       0.0820     9.2491
## s.e.  0.0547     1.3475
## 
## sigma^2 estimated as 453.2:  log likelihood = -1307.3,  aic = 2618.61

# ARMA(2,0,2)
arma202 <- arima(train.ts, order = c(2,0,2), include.mean = TRUE, method = "ML")
arma202

## 
## Call:
## arima(x = train.ts, order = c(2, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1     ar2      ma1      ma2  intercept
##       0.1706  0.7271  -0.0777  -0.6888     9.2467
## s.e.  0.2680  0.2500   0.2615   0.2229     2.6955
## 
## sigma^2 estimated as 433.6:  log likelihood = -1300.96,  aic = 2611.92

# ARIMA(0,1,1)
arima011 <- arima(train.ts, order = c(0,1,1), include.mean = TRUE, method = "ML")
arima011

## 
## Call:
## arima(x = train.ts, order = c(0, 1, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##           ma1
##       -0.9280
## s.e.   0.0334
## 
## sigma^2 estimated as 443.1:  log likelihood = -1300.54,  aic = 2603.08

KndidatModelARIMA <- c("ARMA(0,0,10)", "ARMA(8,0,0)", "ARMA(1,0,1)", 
                       "ARMA(0,0,1)", "ARMA(2,0,2)", "ARIMA(0,1,1)")

AICKandidatModel <- c(arma0010$aic, arma800$aic, arma101$aic, 
                      arma001$aic, arma202$aic, arima011$aic)

compmodelARIMA <- data.frame(
  "Kandidat Model" = KndidatModelARIMA,
  "Nilai AIC" = round(AICKandidatModel, 2),
  check.names = FALSE
)
compmodelARIMA

##   Kandidat Model Nilai AIC
## 1   ARMA(0,0,10)   2595.56
## 2    ARMA(8,0,0)   2593.13
## 3    ARMA(1,0,1)   2608.20
## 4    ARMA(0,0,1)   2618.61
## 5    ARMA(2,0,2)   2611.92
## 6   ARIMA(0,1,1)   2603.08

Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu ARMA(8,0,0)

Model ARMA(8,0,0), maka \(Y_t\) diperoleh dari penjabaran operator backshift sebagai berikut: \[ \begin{aligned} \phi_p(B)(1 - B)^d Y_t &= \mu + \theta_q(B) e_t \\ (1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3 - \phi_4 B^4 - \phi_5 B^5 - \phi_6 B^6 - \phi_7 B^7 - \phi_8 B^8)(1 - B)^0 Y_t &= \mu + (1) e_t \\ (1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3 - \phi_4 B^4 - \phi_5 B^5 - \phi_6 B^6 - \phi_7 B^7 - \phi_8 B^8) Y_t &= \mu + e_t \\ Y_t - \phi_1 B Y_t - \phi_2 B^2 Y_t - \phi_3 B^3 Y_t - \phi_4 B^4 Y_t - \phi_5 B^5 Y_t - \phi_6 B^6 Y_t - \phi_7 B^7 Y_t - \phi_8 B^8 Y_t &= \mu + e_t \\ Y_t - \phi_1 Y_{t-1} - \phi_2 Y_{t-2} - \phi_3 Y_{t-3} - \phi_4 Y_{t-4} - \phi_5 Y_{t-5} - \phi_6 Y_{t-6} - \phi_7 Y_{t-7} - \phi_8 Y_{t-8} &= \mu + e_t \\ Y_t &= \mu + \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \phi_3 Y_{t-3} + \phi_4 Y_{t-4} + \phi_5 Y_{t-5} + \phi_6 Y_{t-6} + \phi_7 Y_{t-7} + \phi_8 Y_{t-8} + e_t \end{aligned} \]

Diagnostic model

ARMA800diag <- stats::arima(train.ts, order = c(8,0,0), method = "ML")
checkresiduals(ARMA800diag)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(8,0,0) with non-zero mean
## Q* = 1.3013, df = 3, p-value = 0.7288
## 
## Model df: 8.   Total lags used: 11

sisaan <- arma800$residuals

Berdasarkan plot di atas terlihat bahwa sisaan tidak mengikuti sebaran normal. Selanjutnya, ditinjau dari plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan.

Uji normalitas data

\[H_0: \text{sisaan menyebar normal}\]

\[H_1: \text{sisaan tidak mengikuti sebaran normal}\]

# Uji normalitas data
ks.test(sisaan,"pnorm")

## Warning in ks.test.default(sisaan, "pnorm"): ties should not be present for the
## one-sample Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.71594, p-value < 2.2e-16
## alternative hypothesis: two-sided

Berdasarkan hasil Kolmogorov-Smirnov Test di atas, didapatkan nilai p-value < \(\alpha = 0.05\), maka keputusan yang diambil adalah tolak \(H_0\). Secara matematis, hal ini membuktikan bahwa sisaan dari model tidak mengikuti sebaran normal.

Uji nilai tengah sisaan

\[H_0: \mu = 0\] \[H_1: \mu \neq 0\]

# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.011325, df = 291, p-value = 0.991
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.279553  2.305938
## sample estimates:
##  mean of x 
## 0.01319249

Berdasarkan hasil uji-t nilai p-value > \(\alpha = 0.05\), maka terima \(H_0\). Ini menunjukkan bahwa secara statistik, nilai harapan (expected value) dari galat terbukti sama dengan nol.

Uji autokorelasi

\[H_0: \text{ tidak ada autokorelasi}\] \[H_1: \text{terdapat autokorelasi}\]

# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 6.7262, df = 23, p-value = 0.9996

Hasil pengujian memberikan nilai p-value > \(\alpha = 0.05\), maka diambil keputusan terima \(H_0\). Hal ini menyimpulkan bahwa sisaan bersifat saling bebas / tidak terdapat autokorelasi (white noise)

Peramalan

h_test <- 73

forecasting_arma0010 <- predict(arma0010, n.ahead = h_test)
forecasting_arma800  <- predict(arma800, n.ahead = h_test)
forecasting_arma101  <- predict(arma101, n.ahead = h_test)
forecasting_arma001  <- predict(arma001, n.ahead = h_test)
forecasting_arma202  <- predict(arma202, n.ahead = h_test)
forecasting_arima011 <- predict(arima011, n.ahead = h_test)

# Menampilkan hasil forecast model terbaik (ARMA 8,0,0) sebagai data frame
hasilforecasting <- as.data.frame(forecasting_arma800)
head(hasilforecasting)

##        pred       se
## 1  7.118789 19.87213
## 2 17.684006 19.91595
## 3  6.694200 19.95447
## 4  5.068215 19.95643
## 5 10.280688 19.95753
## 6  7.442357 20.91339

# 2. Visualisasi hasil peramalan model terbaik (ARMA 8,0,0)
ts.plot(train.ts, forecasting_arma800$pred, 
        ylab = "Jumlah Curah Hujan", xlab = "Waktu (Indeks Hari)", 
        col = c("blue", "red"),
        main = "Visualisasi Hasil Peramalan ARMA(8,0,0)")

# Menambahkan batas kepercayaan (Upper 95% dan Lower 95%) menggunakan Standard Error 
lines(forecasting_arma800$pred + 1.96 * forecasting_arma800$se, col = "green", lty = 2)
lines(forecasting_arma800$pred - 1.96 * forecasting_arma800$se, col = "green", lty = 2)

legend("topleft", 
       legend = c("Data Aktual (Train)", "Forecast", "Upper/Lower 95%"), 
       col = c("blue", "red", "green"), 
       lty = c(1, 1, 2), 
       cex = 0.7, inset = 0.02)

Evaluasi Akurasi hasil peramalan

akurasi.arma0010 <- accuracy(forecasting_arma0010$pred, test.ts)
akurasi.arma800  <- accuracy(forecasting_arma800$pred, test.ts)
akurasi.arma101  <- accuracy(forecasting_arma101$pred, test.ts)
akurasi.arma001  <- accuracy(forecasting_arma001$pred, test.ts)
akurasi.arma202  <- accuracy(forecasting_arma202$pred, test.ts)
akurasi.arima011 <- accuracy(forecasting_arima011$pred, test.ts)

KandidatModel <- c("ARMA(0,0,10)", "ARMA(8,0,0)", "ARMA(1,0,1)", 
                   "ARMA(0,0,1)", "ARMA(2,0,2)", "ARIMA(0,1,1)")

RMSE_Testing <- c(akurasi.arma0010[,"RMSE"], akurasi.arma800[,"RMSE"], 
                  akurasi.arma101[,"RMSE"], akurasi.arma001[,"RMSE"], 
                  akurasi.arma202[,"RMSE"], akurasi.arima011[,"RMSE"])


# Menggabungkan semuanya ke dalam Data Frame
EvaluasiAkurasi <- data.frame(
  "Kandidat Model" = KandidatModel,
  "RMSE" = round(RMSE_Testing, 2),
  check.names = FALSE
)

# Menampilkan Tabel Komparasi Akhir
cat("\n--- TABEL EVALUASI AKURASI MODEL ---\n")

## 
## --- TABEL EVALUASI AKURASI MODEL ---

print(EvaluasiAkurasi)

##   Kandidat Model  RMSE
## 1   ARMA(0,0,10) 21.20
## 2    ARMA(8,0,0) 21.23
## 3    ARMA(1,0,1) 20.85
## 4    ARMA(0,0,1) 20.94
## 5    ARMA(2,0,2) 20.83
## 6   ARIMA(0,1,1) 22.21

Hasil evaluasi menunjukkan bahwa peramalan pada set data pengujian memposisikan model ARMA(2,0,2) dengan nilai simpangan galat (RMSE) paling minimum sebesar 20,83. Namun, filosofi pemodelan dalam pendekatan Box-Jenkins menekankan pada pemilihan arsitektur model berdasarkan keseimbangan terbaik antara parsimony (efisiensi parameter) pada data training dan ketepatan prediksi pada data testing. Mengacu pada signifikansi nilai efisiensi informasi dari kriteria AIC, model ARMA(8,0,0) terpilih sebagai instrumen model dasar yang paling optimal secara struktural, meskipun memiliki angka galat (RMSE sebesar 21.23) yang sedikit berada di bawah model ARMA(2,0,2) pada fase pengujian. Oleh karenanya, ARMA(8,0,0) ditetapkan sebagai model paling representatif dalam kerangka analitik ini.