Membaca Data Curah Hujan

data_hujan <- read.csv("C:/Users/ADVAN/OneDrive/文件/Tugas Apie/Topik Dalam Statistika 1/data curah hujan.csv")
head(data_hujan)

##    Tanggal Data_Curah_Hujan
## 1 1/4/2024              3.2
## 2 2/4/2024              0.0
## 3 3/4/2024              7.8
## 4 4/4/2024              2.4
## 5 5/4/2024             14.6
## 6 6/4/2024              1.2

Mengubah data tidak valid menjadi NA

data_hujan$Data_Curah_Hujan[data_hujan$Data_Curah_Hujan == 8888 | data_hujan$Data_Curah_Hujan == ""] <- NA

Mengubah data tidak kosong menjadi NA

data_hujan$Data_Curah_Hujan[data_hujan$Data_Curah_Hujan == "" ] <- NA

Mengecek Jumlah Missing Value

sum(is.na(data_hujan$Data_Curah_Hujan))

## [1] 68

Melakukan Interpolasi Missing Value

library(zoo)

## Warning: package 'zoo' was built under R version 4.5.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

data_hujan$Data_Curah_Hujan <- na.approx(data_hujan$Data_Curah_Hujan)

sum(is.na(data_hujan$Data_Curah_Hujan))

## [1] 0

names(data_hujan)

## [1] "Tanggal"          "Data_Curah_Hujan"

Mengubah Format Tanggal Menjadi Format Date

# Ubah menjadi character
data_hujan$Tanggal <- as.character(data_hujan$Tanggal)

# Membuat kolom tanggal baru
tanggal_baru <- rep(NA_character_, nrow(data_hujan))

# Format dd/mm/yyyy
idx1 <- grepl("/", data_hujan$Tanggal)

tanggal_baru[idx1] <- as.character(
  as.Date(data_hujan$Tanggal[idx1],
          format = "%d/%m/%Y")
)

# Format dd-mm-yyyy
idx2 <- grepl("-", data_hujan$Tanggal)

tanggal_baru[idx2] <- as.character(
  as.Date(data_hujan$Tanggal[idx2],
          format = "%d-%m-%Y")
)

# Konversi menjadi Date
data_hujan$Tanggal <- as.Date(tanggal_baru)

Membentuk Variabel Tahun dan Bulan

data_hujan$Tahun <- format(data_hujan$Tanggal, "%Y")
data_hujan$Bulan <- format(data_hujan$Tanggal, "%m")

#data_hujan$Data_Curah_Hujan <- as.numeric(data_hujan$Data_Curah_Hujan)

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.5.3

## 
## 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(lubridate)

## Warning: package 'lubridate' was built under R version 4.5.3

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# Membuat variabel minggu
data_hujan$Minggu <- floor_date(data_hujan$Tanggal,
                                unit = "week",
                                week_start = 1)

# Mengubah data harian menjadi mingguan
data_mingguan <- data_hujan %>%
  group_by(Minggu) %>%
  summarise(
    Curah_Hujan_Mingguan = sum(Data_Curah_Hujan,
                               na.rm = TRUE)
  )

# Melihat hasil data mingguan
head(data_mingguan)

## # A tibble: 6 × 2
##   Minggu     Curah_Hujan_Mingguan
##   <date>                    <dbl>
## 1 2024-04-01                 40.8
## 2 2024-04-08                 42  
## 3 2024-04-15                 68.2
## 4 2024-04-22                 52.5
## 5 2024-04-29                 71.4
## 6 2024-05-06                109.

data_mingguan <- data_hujan %>%
  group_by(Minggu) %>%
  summarise(
    Curah_Hujan_Mingguan = round(
      mean(Data_Curah_Hujan, na.rm = TRUE),
      1
    )
  )

# Menampilkan hasil
print(data_mingguan)

## # A tibble: 105 × 2
##    Minggu     Curah_Hujan_Mingguan
##    <date>                    <dbl>
##  1 2024-04-01                  5.8
##  2 2024-04-08                  6  
##  3 2024-04-15                  9.7
##  4 2024-04-22                  7.5
##  5 2024-04-29                 10.2
##  6 2024-05-06                 15.6
##  7 2024-05-13                  9.1
##  8 2024-05-20                 14.4
##  9 2024-05-27                 21.6
## 10 2024-06-03                 10.6
## # ℹ 95 more rows

Analisis Statistika Deskriptif

# Analisis Statistika Deskriptif
rata_rata <- mean(data_mingguan$Curah_Hujan_Mingguan, na.rm = TRUE)

kuartil <- quantile(
  data_mingguan$Curah_Hujan_Mingguan,
  probs = c(0.25, 0.50, 0.75),
  na.rm = TRUE
)

q1 <- kuartil[1]
median <- kuartil[2]
q3 <- kuartil[3]

varians <- var(data_mingguan$Curah_Hujan_Mingguan, na.rm = TRUE)
simpangan_baku <- sd(data_mingguan$Curah_Hujan_Mingguan, na.rm = TRUE)

maksimum <- max(data_mingguan$Curah_Hujan_Mingguan, na.rm = TRUE)
minimum <- min(data_mingguan$Curah_Hujan_Mingguan, na.rm = TRUE)

# Menampilkan hasil
statistika_deskriptif <- data.frame(
  Statistik = c(
    "Minimum",
    "Kuartil 1 (Q1)",
    "Median (Q2)",
    "Kuartil 3 (Q3)",
    "Maksimum",
    "Rata-rata",
    "Varians",
    "Simpangan Baku"
  ),
  Nilai = c(
    minimum,
    q1,
    median,
    q3,
    maksimum,
    rata_rata,
    varians,
    simpangan_baku
  )
)

# Membulatkan 2 angka di belakang koma
statistika_deskriptif$Nilai <- round(
  statistika_deskriptif$Nilai,
  2
)

print(statistika_deskriptif)

##        Statistik Nilai
## 1        Minimum  0.00
## 2 Kuartil 1 (Q1)  3.40
## 3    Median (Q2)  7.40
## 4 Kuartil 3 (Q3) 12.50
## 5       Maksimum 32.60
## 6      Rata-rata  9.25
## 7        Varians 60.88
## 8 Simpangan Baku  7.80

Mengubah data menjadi Time Series

data_mingguan.ts <- ts(data_mingguan$Curah_Hujan_Mingguan, start=1)
head(data_mingguan.ts)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1]  5.8  6.0  9.7  7.5 10.2 15.6

Visualisasi Data Time Series

ts.plot(data_mingguan.ts, col="blue", ylab = "Curah Hujan", xlab = "Minggu")
title(main = "Time Series Plot Data Curah Hujan", sub =  "Sumber : Data Online BMKG Stasiun Meteorologi Tjilik Riwut ", 
      cex.sub = 0.8)
points(data_mingguan.ts, pch = 20, col = "blue")

spliting data

# Training set
train.Curah_Hujan_Mingguan.ts <- data_mingguan.ts[1:84]

# Testing set
test.Curah_Hujan_Mingguan.ts <- data_mingguan.ts[85:105]

Visualisasi data train yang akan digunakan untuk mencari model terbaik

ts.plot(train.Curah_Hujan_Mingguan.ts, col="blue", ylab = "Curah_Hujan_Mingguan", xlab = "Minggu")
title(main = "Train Time Series Plot Curah Hujan", cex.sub = 0.8)
points(train.Curah_Hujan_Mingguan.ts, pch = 20, col = "blue")

Visualisasi data testing yang akan digunakan untuk mencari model terbaik

ts.plot(test.Curah_Hujan_Mingguan.ts, col="blue", ylab = "Curah_Hujan_Mingguan", xlab = "Minggu")
title(main = "Test Time Series Plot Curah Hujan", cex.sub = 0.8)
points(test.Curah_Hujan_Mingguan.ts, pch = 20, col = "blue")

# Model ARIMA

1. Cek Cek kestasioneran data

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

adf.test(train.Curah_Hujan_Mingguan.ts, alternative = "stationary", k=trunc((length(train.Curah_Hujan_Mingguan.ts)-1)^(1/3)))

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.Curah_Hujan_Mingguan.ts
## Dickey-Fuller = -2.8635, Lag order = 4, p-value = 0.2217
## alternative hypothesis: stationary

Karena nilai \(p-value = 0.2217 > \alpha = 0.05\) sehingag dilakukan differencing hingga data train stasioner

# Differencing Ordo 1
Curah_Hujan_Mingguan.diff1 <- diff(train.Curah_Hujan_Mingguan.ts, difference=1)

library(tseries)
adf.test(Curah_Hujan_Mingguan.diff1, alternative = "stationary", k=trunc((length(Curah_Hujan_Mingguan.diff1)-1)^(1/3)))

## Warning in adf.test(Curah_Hujan_Mingguan.diff1, alternative = "stationary", :
## p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Curah_Hujan_Mingguan.diff1
## Dickey-Fuller = -6.3133, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

2. Identifikasi kandidat model

acf(Curah_Hujan_Mingguan.diff1, lag.max = 30, col = "blue")

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

pacf(Curah_Hujan_Mingguan.diff1, lag.max = 30, col = "blue")

Berdasarkan plot PACF di atas kandidat model yang diperoleh adalah \(ARIMA(2,1,0)\)

library(TSA)

## Warning: package 'TSA' was built under R version 4.5.3

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

eacf(Curah_Hujan_Mingguan.diff1)

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

Berdasarkan hasil EACF di atas kandidat model yang diperoleh adalah \(ARIMA(1,1,0)\),\(ARIMA(2,1,1)\),\(ARIMA(1,1,3)\)

Karena data yang digunakan adalah data yang hasil differencing orde pertama maka \(d=1\) dan dapat ditentukan beberapa kandidat model yaitu: \(ARIMA(0,1,1),ARIMA(2,1,0),ARIMA(1,1,0),ARIMA(2,1,1),ARIMA(1,1,3)\)

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

## 
## Call:
## arima(x = Curah_Hujan_Mingguan.diff1, order = c(0, 0, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ma1  intercept
##       -0.8095     0.0622
## s.e.   0.0836     0.1749
## 
## sigma^2 estimated as 61.47:  log likelihood = -289.23,  aic = 582.45

# ARIMA(2,1,0)
arima210 <- arima(Curah_Hujan_Mingguan.diff1, order = c(2,0,0), include.mean = TRUE, method = "ML")
arima210

## 
## Call:
## arima(x = Curah_Hujan_Mingguan.diff1, order = c(2, 0, 0), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.5087  -0.3546     0.1585
## s.e.   0.1043   0.1035     0.4924
## 
## sigma^2 estimated as 68.73:  log likelihood = -293.54,  aic = 593.07

# ARIMA(1,1,0)
arima110 <- arima(Curah_Hujan_Mingguan.diff1, order = c(1,0,0), include.mean = TRUE, method = "ML")
arima110

## 
## Call:
## arima(x = Curah_Hujan_Mingguan.diff1, order = c(1, 0, 0), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1  intercept
##       -0.3699     0.1197
## s.e.   0.1025     0.7129
## 
## sigma^2 estimated as 78.63:  log likelihood = -298.98,  aic = 601.97

# ARIMA(2,1,1)
arima211 <- arima(Curah_Hujan_Mingguan.diff1, order = c(2,0,1), include.mean = TRUE, method = "ML")
arima211

## 
## Call:
## arima(x = Curah_Hujan_Mingguan.diff1, order = c(2, 0, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##          ar1      ar2      ma1  intercept
##       0.1682  -0.0899  -0.8550     0.0454
## s.e.  0.1756   0.1524   0.1517     0.1481
## 
## sigma^2 estimated as 59.65:  log likelihood = -288.04,  aic = 584.08

# ARIMA(1,1,3)
arima113 <- arima(Curah_Hujan_Mingguan.diff1, order = c(1,0,3), include.mean = TRUE, method = "ML")
arima113

## 
## Call:
## arima(x = Curah_Hujan_Mingguan.diff1, order = c(1, 0, 3), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1     ma1      ma2      ma3  intercept
##       -0.8136  0.1356  -0.7865  -0.1091     0.0403
## s.e.   0.2867  0.3065   0.2095   0.1817     0.1264
## 
## sigma^2 estimated as 58.94:  log likelihood = -287.63,  aic = 585.25

Model Terbaik

AICKandidatModel <- c(arima011$aic, arima210$aic, arima110$aic, arima211$aic, arima113$aic)
KandidatModelARIMA <- c("ARIMA(0,1,1)", "ARIMA(2,1,0)","ARIMA(1,1,0)","ARIMA(2,1,1)", "ARIMA(1,1,3)")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA

##   Kandidat Model        Nilai AIC
## 1   ARIMA(0,1,1) 582.450422970504
## 2   ARIMA(2,1,0) 593.074312771881
## 3   ARIMA(1,1,0) 601.965328533654
## 4   ARIMA(2,1,1) 584.080230574667
## 5   ARIMA(1,1,3) 585.250321746231

Model terbaik diperoleh berdasarkan nilai AIC terkecil. Oleh karena itu, model terbaik yang diperoleh adalah \(ARIMA(0,1,1)\)

Model terbaik ini selanjutnya dilakukan uji hipotesis untuk mengetahui signifikansi parameternya, dimana:

intercept_mu <- 0.0622/0.1749
intercept_mu

## [1] 0.3556318

ma_theta1 <- -0.8095/0.0836 
ma_theta1

## [1] -9.683014

prmtrdugaanarima011 <- c("mu", "theta(1)")
thitungprmtrarima011 <- c(abs(intercept_mu), abs(ma_theta1))
ttabel <- c("1.661961084", "1.661961084")
keputusan <- c("Tidak Signifikan","Signifikan")
tksignimodelarima011 <- cbind(prmtrdugaanarima011, thitungprmtrarima011, ttabel, keputusan)

colnames(tksignimodelarima011) <- c("Parameter dugaan", "T-Hitung", "T-Tabel", "Keputusan")
tksignimodelarima011 <- as.data.frame(tksignimodelarima011)
tksignimodelarima011

##   Parameter dugaan          T-Hitung     T-Tabel        Keputusan
## 1               mu 0.355631789594054 1.661961084 Tidak Signifikan
## 2         theta(1)  9.68301435406699 1.661961084       Signifikan

lmtest::coeftest(arima011)

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)    
## ma1       -0.809516   0.083585 -9.6849   <2e-16 ***
## intercept  0.062194   0.174945  0.3555   0.7222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dibandingkan Arima

library(forecast)

## Warning: package 'forecast' was built under R version 4.5.3

## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA

model1 <- Arima(Curah_Hujan_Mingguan.diff1,order=c(0,0,1))
summary(model1)

## Series: Curah_Hujan_Mingguan.diff1 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##           ma1    mean
##       -0.8096  0.0622
## s.e.   0.0836  0.1749
## 
## sigma^2 = 62.99:  log likelihood = -289.23
## AIC=584.45   AICc=584.75   BIC=591.71
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 0.06691802 7.840362 6.033166 76.36537 162.4157 0.5120261 0.1074457

model2 <- Arima(Curah_Hujan_Mingguan.diff1,order=c(2,0,0))
summary(model2)

## Series: Curah_Hujan_Mingguan.diff1 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2    mean
##       -0.5087  -0.3546  0.1585
## s.e.   0.1043   0.1035  0.4924
## 
## sigma^2 = 71.31:  log likelihood = -293.54
## AIC=595.07   AICc=595.59   BIC=604.75
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE    MAPE      MASE
## Training set 0.01264219 8.290571 6.089607 70.99019 120.366 0.5168161
##                     ACF1
## Training set -0.07471012

model3 <- Arima(Curah_Hujan_Mingguan.diff1,order=c(1,0,0))
summary(model3)

## Series: Curah_Hujan_Mingguan.diff1 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1    mean
##       -0.3698  0.1212
## s.e.   0.1024  0.7130
## 
## sigma^2 = 80.57:  log likelihood = -298.98
## AIC=603.97   AICc=604.27   BIC=611.22
## 
## Training set error measures:
##                        ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.001841727 8.867355 6.479491 87.36389 104.6858 0.5499051
##                    ACF1
## Training set -0.1325953

model4 <- Arima(Curah_Hujan_Mingguan.diff1,order=c(2,0,1))
summary(model4)

## Series: Curah_Hujan_Mingguan.diff1 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1    mean
##       0.1682  -0.0898  -0.8550  0.0454
## s.e.  0.1757   0.1524   0.1518  0.1481
## 
## sigma^2 = 62.67:  log likelihood = -288.04
## AIC=586.08   AICc=586.86   BIC=598.17
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.04578903 7.723374 5.913358 85.77726 157.9522 0.5018582
##                      ACF1
## Training set -0.006465298

model5 <- Arima(Curah_Hujan_Mingguan.diff1,order=c(1,0,3))
summary(model5)

## Series: Curah_Hujan_Mingguan.diff1 
## ARIMA(1,0,3) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1      ma2      ma3    mean
##       -0.8136  0.1356  -0.7865  -0.1091  0.0403
## s.e.   0.2867  0.3065   0.2095   0.1817  0.1264
## 
## sigma^2 = 62.71:  log likelihood = -287.63
## AIC=587.25   AICc=588.36   BIC=601.76
## 
## Training set error measures:
##                      ME    RMSE      MAE      MPE     MAPE      MASE
## Training set 0.02340331 7.67704 5.883528 95.68344 161.8555 0.4993265
##                      ACF1
## Training set -0.006080369

AICKandidatModel <- c(model1$aic, model2$aic, model3$aic, model4$aic, model5$aic)
AICcKandidatModel <- c(model1$aicc, model2$aicc, model3$aicc, model4$aicc, model5$aicc)
BICKandidatModel <- c(model1$bic, model2$bic, model3$bic, model4$bic, model5$bic)
KandidatModelARIMA <- c("ARIMA(0,1,1)", "ARIMA(2,1,0)","ARIMA(1,1,0)","ARIMA(2,1,1)", "ARIMA(1,1,3)")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, BICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC", "Nilai AICc", "Nilai BIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA

##   Kandidat Model        Nilai AIC       Nilai AICc        Nilai BIC
## 1   ARIMA(0,1,1) 584.450422709736 584.754220178091 591.706944533126
## 2   ARIMA(2,1,0) 595.074312771582 595.587133284403 604.749675202769
## 3   ARIMA(1,1,0) 603.965335012339 604.269132480693 611.221856835728
## 4   ARIMA(2,1,1) 586.080230568305 586.859451347526 598.174433607288
## 5   ARIMA(1,1,3) 587.250321746231 588.355584904125  601.76336539301

Model terbaik yang didapatkan adalah \(ARIMA(0,1,1)\) dengan fungsi Arima di mana nilai AIC sebesar \(584.45\) merupakan nilai AIC terkecil.

Oleh karena itu, model yang disarankan untuk dianalisis lebih lanjut yaitu dengan menggunakan model \(ARIMA(0,1,1)\) yang merupakan model terbaik berdasarkan nilai AIC terkecil dari kandidat model berdasarkan ACF, PACF, dan EACF.

Model \(ARIMA(0,1,1)\), maka \(Y_t\) diperoleh dari penjabaran operator backshift untuk model \(ARIMA(0,1,1)\) sebagai berikut:

\[ \phi_p(B)(1-B)^dY_t = \mu + \theta_q(B)e_t \] \[ \phi_0(B)(1-B)^dY_t = \mu + \theta_q(B)e_t \] \[ (1-B)Y_t = \mu + \theta_1(B)e_t \] \[ Y_t-Y_{t-1} = \mu + e_t-\theta_1e_{t-1} \] \[ Y_t = \mu + Y_{t-1} + e_t - \theta_1e_{t-1} \]

Pendugaan parameter \(\hat{\mu}\) = \(0.0622\), \(\theta_1 = -0.8096\). Model terbaik yang diperoleh yaitu model ARIMA(0,1,1).

] \[ Y_t = 0.0622 + Y_{t-1} + e_t + 0.8096e_{t-1} \]

Diagnostic model

ARIMA011diag <- stats::arima(Curah_Hujan_Mingguan.diff1 , order = c(0,0,1), method = "ML")
checkresiduals(ARIMA011diag)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,1) with non-zero mean
## Q* = 12.003, df = 9, p-value = 0.2131
## 
## Model df: 1.   Total lags used: 10

Berdasarkan plot di atas terlihat bahwa sisaan tidak mengikuti sebaran normal. Selanjutnya, dari plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal:

sisaan <- arima011$residuals

Uji normalitas data \(H_0\) : sisaan menyebar normal \(H_1\) : sisaan tidak mengikuti sebaran normal

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

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.44682, p-value = 8.882e-16
## alternative hypothesis: two-sided

Berdasarkan hasil \(p-value = 8.882e-16 < \alpha = 0.05\) yang berarti tolak \(H_0\), artinya sisaan tidak menyebar normal.

Uji nilai tengah sisaan \(H_0\) : \(\mu = 0\) \(H_1\) : \(\mu \ne 0\)

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

## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.077401, df = 82, p-value = 0.9385
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.655323  1.789351
## sample estimates:
##  mean of x 
## 0.06701357

Berdasarkan hasil \(p-value = 0.9385 > \alpha = 0.05\) yang berarti terima \(H_0\), artinya nilai tengah sisaan sama dengan \(0\)

Uji autokorelasi \(H_0\) : tidak ada autokorelasi \(H_1\) : terdapat autokorelasi

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

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

Hasil \(p-value = 0.08067 > \alpha = 0.05\) yang berarti terima \(H_0\) atau tidak terdapat gejala autokorelasi

Kesimpulan: Asumsi terpenuhi, kecuali sisaan tidak menyebar normal.

Peramalan

Peramalan hasil dari pemodelan dengan menggunakan data training sebagai berikut:

Fitted model \(ARIMA(0,1,1)\) dengan data training

predictARIMA011 <- fitted(arima011)
fittedARIMA011 <- cbind(Curah_Hujan_Mingguan.diff1, predictARIMA011)

Visualisasi hasil fitted model \(ARIMA(0,1,1)\) dengan data training

ts.plot(Curah_Hujan_Mingguan.diff1, col="blue", ylab = "Curah Hujan", xlab = "Minggu")
title(main = "Fitted Time Series Plot Curah Hujan Mingguan", cex.sub = 0.8)
points(Curah_Hujan_Mingguan.diff1, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA011)

forecasting <- predict(arima011, n.ahead = 5)
hasilforcasting <- as.data.frame(forecasting)
hasilforcasting

##         pred        se
## 1 1.27838137  7.840371
## 2 0.06219376 10.087351
## 3 0.06219376 10.087351
## 4 0.06219376 10.087351
## 5 0.06219376 10.087351

# Membentuk model ARIMA(0,1,1) pada data training
library(forecast)

model_arima011 <- Arima(
  train.Curah_Hujan_Mingguan.ts,
  order = c(0,1,1)
)

# Peramalan 5 minggu ke depan
forecast_arima011 <- forecast(
  model_arima011,
  h = 5
)

# Menampilkan hasil peramalan
forecast_arima011

##    Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 85       13.04968 2.932917 23.16645 -2.422575 28.52194
## 86       13.04968 2.751175 23.34819 -2.700525 28.79989
## 87       13.04968 2.572586 23.52678 -2.973655 29.07302
## 88       13.04968 2.396989 23.70238 -3.242206 29.34157
## 89       13.04968 2.224241 23.87513 -3.506402 29.60577

plot(
  forecast_arima011,
  main = "Peramalan Curah Hujan Mingguan dengan ARIMA(0,1,1)",
  xlab = "Minggu",
  ylab = "Curah Hujan (mm)",
  col = "blue"
)

lines(train.Curah_Hujan_Mingguan.ts,
      col = "black",
      lwd = 2)

legend(
  "topleft",
  legend = c("Data Aktual", "Ramalan"),
  col = c("black", "blue"),
  lty = 1,
  lwd = 2
)

Visualisasi hasil peramalan model \(ARIMA(0,1,1) 5\) waktu ke depan

Fitting \(ARIMA(0,1,1)\) dengan data testing

arima011test <- arima(test.Curah_Hujan_Mingguan.ts, order = c(0,1,1), include.mean = TRUE, method = "ML")
arima011test

## 
## Call:
## arima(x = test.Curah_Hujan_Mingguan.ts, order = c(0, 1, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ma1
##       -1.0000
## s.e.   0.2043
## 
## sigma^2 estimated as 58.38:  log likelihood = -70.57,  aic = 143.14

model7 <- Arima(test.Curah_Hujan_Mingguan.ts,order=c(0,0,1))
summary(model7)

## Series: test.Curah_Hujan_Mingguan.ts 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##           ma1    mean
##       -0.1475  8.8170
## s.e.   0.2662  1.4036
## 
## sigma^2 = 60.51:  log likelihood = -71.84
## AIC=149.67   AICc=151.09   BIC=152.81
## 
## Training set error measures:
##                       ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 0.004951582 7.399039 5.782727 -Inf  Inf 0.7846305 0.003978717

predictARIMA011test <- fitted(arima011test)
fittedARIMA011test <- cbind(test.Curah_Hujan_Mingguan.ts, predictARIMA011test)

Visualisasi hasil fitted model \(ARIMA(0,1,1)\) dengan testing

ts.plot(test.Curah_Hujan_Mingguan.ts, col="blue", ylab = "Curah_Hujan_Mingguan", xlab = "Minggu")
title(main = "Fitted Time Series Plot Curah Hujan Mingguan", cex.sub = 0.8)
points(test.Curah_Hujan_Mingguan.ts, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA011test)

forecastingtest <- predict(arima011test, n.ahead = 5)
hasilforcastingtest <- as.data.frame(forecastingtest)

Akurasi hasil peramalan

Akurasi hasil peramalan model \(ARIMA(0,1,1)\). Akurasi hasil peramalan model \(ARIMA(0,1,1)\) dibandingkan dengan data aktual adalah Data testing.

library(forecast)

# Model terbaik dari data training
model_arima <- Arima(train.Curah_Hujan_Mingguan.ts,
                     order = c(0,1,1))

# Meramalkan sebanyak panjang data testing
forecast_test <- forecast(model_arima,
                          h = length(test.Curah_Hujan_Mingguan.ts))

# Akurasi terhadap data testing
akurasi_testing <- accuracy(
  forecast_test,
  test.Curah_Hujan_Mingguan.ts
)

akurasi_testing

##                      ME     RMSE      MAE  MPE MAPE      MASE      ACF1
## Training set  0.3669447 7.799610 5.929500 -Inf  Inf 0.8408482 0.1069515
## Test set     -4.1163499 8.517147 7.368821 -Inf  Inf 1.0449550        NA