Model SARIMA

Membaca Data Curah Hujan

data_hujan <- read.csv("C:/Users/User/OneDrive/Dokumen/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)

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

# sum(is.na(data_hujan$Tanggal))

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)

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

## 
## 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
print(data_mingguan, n = Inf)

## # A tibble: 105 × 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. 
##   7 2024-05-13                 64  
##   8 2024-05-20                101. 
##   9 2024-05-27                151. 
##  10 2024-06-03                 74.2
##  11 2024-06-10                 22.4
##  12 2024-06-17                  9.3
##  13 2024-06-24                 14.1
##  14 2024-07-01                 48.2
##  15 2024-07-08                  3  
##  16 2024-07-15                 27.4
##  17 2024-07-22                  0  
##  18 2024-07-29                  6  
##  19 2024-08-05                  1  
##  20 2024-08-12                 14  
##  21 2024-08-19                 87.2
##  22 2024-08-26                 50.4
##  23 2024-09-02                 39.4
##  24 2024-09-09                 24.6
##  25 2024-09-16                  0  
##  26 2024-09-23                 75  
##  27 2024-09-30                 54.2
##  28 2024-10-07                 48.8
##  29 2024-10-14                135. 
##  30 2024-10-21                  3.6
##  31 2024-10-28                 37  
##  32 2024-11-04                 48  
##  33 2024-11-11                145. 
##  34 2024-11-18                111. 
##  35 2024-11-25                 79.2
##  36 2024-12-02                126. 
##  37 2024-12-09                 72.6
##  38 2024-12-16                106. 
##  39 2024-12-23                 66.4
##  40 2024-12-30                 97.2
##  41 2025-01-06                 65.2
##  42 2025-01-13                111. 
##  43 2025-01-20                216. 
##  44 2025-01-27                159. 
##  45 2025-02-03                 16.6
##  46 2025-02-10                 39.2
##  47 2025-02-17                119. 
##  48 2025-02-24                 12.2
##  49 2025-03-03                103. 
##  50 2025-03-10                180. 
##  51 2025-03-17                 50.5
##  52 2025-03-24                 38.1
##  53 2025-03-31                 95.2
##  54 2025-04-07                189. 
##  55 2025-04-14                 66.2
##  56 2025-04-21                 48.4
##  57 2025-04-28                 42  
##  58 2025-05-05                 15.2
##  59 2025-05-12                125. 
##  60 2025-05-19                 52  
##  61 2025-05-26                228. 
##  62 2025-06-02                  6.3
##  63 2025-06-09                 35.9
##  64 2025-06-16                 73.2
##  65 2025-06-23                 20.8
##  66 2025-06-30                  8  
##  67 2025-07-07                 17.8
##  68 2025-07-14                 13.4
##  69 2025-07-21                  5.8
##  70 2025-07-28                  1.8
##  71 2025-08-04                 93.4
##  72 2025-08-11                 70.2
##  73 2025-08-18                 80.6
##  74 2025-08-25                 11  
##  75 2025-09-01                 23.6
##  76 2025-09-08                 73.1
##  77 2025-09-15                 56.8
##  78 2025-09-22                 10.2
##  79 2025-09-29                 43.2
##  80 2025-10-06                 30  
##  81 2025-10-13                 27.6
##  82 2025-10-20                224  
##  83 2025-10-27                178. 
##  84 2025-11-03                 84.9
##  85 2025-11-10                 61.8
##  86 2025-11-17                 89.3
##  87 2025-11-24                 72  
##  88 2025-12-01                 74.8
##  89 2025-12-08                  7.2
##  90 2025-12-15                196. 
##  91 2025-12-22                 78.6
##  92 2025-12-29                 65.2
##  93 2026-01-05                 63.5
##  94 2026-01-12                  0.8
##  95 2026-01-19                  0  
##  96 2026-01-26                138. 
##  97 2026-02-02                 68.3
##  98 2026-02-09                 26.6
##  99 2026-02-16                 47.1
## 100 2026-02-23                 31.6
## 101 2026-03-02                 35.6
## 102 2026-03-09                  6.2
## 103 2026-03-16                 45.6
## 104 2026-03-23                 26  
## 105 2026-03-30                 51.2

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

length(data_mingguan$Curah_Hujan_Mingguan)

## [1] 105

Ubah data menjadi time series

Dilakukan perubahan data Curah_Hujan_Mingguan menjadi data time series dengan fungsi ts

data_mingguan.ts <- ts(data_mingguan, start=c(2024,1),frequency=52)

Visualisasi data time series

ts.plot(data_mingguan.ts[,"Curah_Hujan_Mingguan"], type="l", ylab="Curah_Hujan_Mingguan", col="blue")
title(main = "Time Series Plot Curah Hujan Mingguan", cex.sub = 0.8)
points(data_mingguan.ts, pch = 20, col = "blue")

Plot di atas memperlihatkan pola musiman dengan $s = 52$

Visualisasi per Musim

library(forecast)
seasonplot(data_mingguan.ts[,"Curah_Hujan_Mingguan"],52,main="Seasonal Plot Curah Hujan", ylab="Minggu",year.labels = TRUE, 
col=rainbow(18))

# Visualisasi deskriptif

monthplot(data_mingguan.ts[,"Curah_Hujan_Mingguan"],ylab="Curah_Hujan_Mingguan", col="blue")

Spliting data

train.ts <- subset(data_mingguan.ts[,"Curah_Hujan_Mingguan"],start=1,end=84)
test.ts <- subset(data_mingguan.ts[,"Curah_Hujan_Mingguan"],start=85,end=105)

Visualisasi data train yang akan digunakan untuk mencari model terbaik

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

Visualisasi data testing yang akan digunakan untuk mencari model terbaik

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

Uji stasioneritas data

library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

adf.test(train.ts)

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

par (mfrow=c(1,2))
Acf(train.ts, lag.max = 104)
Pacf(train.ts, lag.max = 104)

Differensi Musiman Orde 1

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

ts.plot(Curah_Hujan_Mingguan.diff1, col="red", main="Curah Hujan Mingguan diff1")

adf.test(Curah_Hujan_Mingguan.diff1)

## Warning in adf.test(Curah_Hujan_Mingguan.diff1): 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

Differensi Non-Musiman Orde 1

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

ts.plot(Curah_Hujan_Mingguan.ddiff1, col="red", main="Curah Hujan Mingguan ddiff1")

### Uji stationeritas data differencing musiman dan non-musiman ###
adf.test(Curah_Hujan_Mingguan.ddiff1)

## Warning in adf.test(Curah_Hujan_Mingguan.ddiff1): p-value smaller than printed
## p-value

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

Identifikasi Model SARIMA

acf3 <- acf(Curah_Hujan_Mingguan.ddiff1,lag.max=104,xaxt="n", main="ACF Curah_Hujan_Mingguan.ddiff1", col="blue")
axis(1, at=0:104/52, labels=0:104)

pacf3 <- pacf(Curah_Hujan_Mingguan.ddiff1,lag.max=104,xaxt="n", main="PACF Curah_Hujan_Mingguan.ddiff1", col="blue")
axis(1, at=0:104/52, labels=0:104)

library(TSA)

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

## 
## 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.ddiff1)

## 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 o o o o o x o x o o  o  o  o 
## 2 x o o o o o x x o o o  o  o  o 
## 3 x o o o o o o x o o o  x  o  o 
## 4 x o o o o o o o o o o  o  o  o 
## 5 x x x o x o x o o o o  o  x  o 
## 6 x o o x x o o o o o o  o  o  o 
## 7 o o o o x o o o o o o  o  o  o

Karena, kedua komponen telah stasioner. Identifikasi komponen non-seasonal adalah $ARIMA(0,1,1), ARIMA(2,1,1), ARIMA(1,1,0), ARIMA(1,1,1),ARIMA(1,1,3)$, dan $ARIMA(2,1,0)$. identifikasi komponen seasonal adalah $ARIMA(0,1,1)_{52}$, Sehingga model yang diperoleh adalah: $ARIMA(0,1,1) \times ARIMA(0,1,1)_{52}$ $ARIMA(2,1,1) \times ARIMA(0,1,1)_{52}$ $ARIMA(1,1,0) \times ARIMA(0,1,1)_{52}$ $ARIMA(1,1,1) \times ARIMA(0,1,1)_{52}$ $ARIMA(1,1,3) \times ARIMA(0,1,1)_{52}$ $ARIMA(2,1,0) \times ARIMA(0,1,1)_{52}$

Estimasi Parameter Model SARIMA

library(forecast)
model1 <- Arima(train.ts,order=c(0,1,1),seasonal=c(0,1,1))
summary(model1)

## Series: train.ts 
## ARIMA(0,1,1)(0,1,1)[52] 
## 
## Coefficients:
##           ma1      sma1
##       -1.0000    0.0000
## s.e.   0.1188  112.9645
## 
## sigma^2 = 108.5:  log likelihood = -117.33
## AIC=240.65   AICc=241.54   BIC=244.96
## 
## Training set error measures:
##                     ME    RMSE      MAE  MPE MAPE      MASE       ACF1
## Training set -0.189563 6.11958 2.919111 -Inf  Inf 0.3926504 0.04112259

model2 <- Arima(train.ts,order=c(2,1,1),seasonal=c(0,1,1))
summary(model2)

## Series: train.ts 
## ARIMA(2,1,1)(0,1,1)[52] 
## 
## Coefficients:

## Warning in sqrt(diag(x$var.coef)): NaNs produced

##          ar1      ar2      ma1  sma1
##       0.0509  -0.0853  -1.0000     0
## s.e.  0.1803   0.2012   0.1188   NaN
## 
## sigma^2 = 115.2:  log likelihood = -117.19
## AIC=244.38   AICc=246.78   BIC=251.55
## 
## Training set error measures:
##                     ME    RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set -0.157822 6.08508 2.909548 -Inf  Inf 0.3913642 -0.02724735

model3 <- Arima(train.ts,order=c(1,1,0),seasonal=c(0,1,1))
summary(model3)

## Series: train.ts 
## ARIMA(1,1,0)(0,1,1)[52] 
## 
## Coefficients:
##           ar1     sma1
##       -0.4737   0.0000
## s.e.   0.1595  75.0509
## 
## sigma^2 = 163.3:  log likelihood = -122.06
## AIC=250.12   AICc=251.01   BIC=254.42
## 
## Training set error measures:
##                      ME    RMSE      MAE  MPE MAPE      MASE       ACF1
## Training set 0.02096236 7.50786 3.468423 -Inf  Inf 0.4665386 -0.1568929

model4 <- Arima(train.ts,order=c(1,1,1),seasonal=c(0,1,1))
summary(model4)

## Series: train.ts 
## ARIMA(1,1,1)(0,1,1)[52] 
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.0560  -1.0000   0.0000
## s.e.  0.1811   0.1113  55.0615
## 
## sigma^2 = 112.4:  log likelihood = -117.28
## AIC=242.56   AICc=244.1   BIC=248.29
## 
## Training set error measures:
##                      ME     RMSE    MAE  MPE MAPE      MASE        ACF1
## Training set -0.1816832 6.120846 2.9115 -Inf  Inf 0.3916267 -0.01182591

model5 <- Arima(train.ts,order=c(1,1,3),seasonal=c(0,1,1))
summary(model5)

## Series: train.ts 
## ARIMA(1,1,3)(0,1,1)[52] 
## 
## Coefficients:
##           ar1     ma1      ma2      ma3    sma1
##       -0.9516  0.0404  -0.9990  -0.0405  0.7796
## s.e.   0.0945  0.2389   0.1371   0.1962  3.6246
## 
## sigma^2 = 73.25:  log likelihood = -117.05
## AIC=246.11   AICc=249.61   BIC=254.71
## 
## Training set error measures:
##                      ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set -0.1474401 4.761443 2.277768 -Inf  Inf 0.3063833 -0.02122717

model6 <- Arima(train.ts,order=c(2,1,0),seasonal=c(0,1,1))
summary(model6)

## Series: train.ts 
## ARIMA(2,1,0)(0,1,1)[52] 
## 
## Coefficients:
##           ar1      ar2     sma1
##       -0.6495  -0.3617   0.0000
## s.e.   0.1735   0.1774  63.3452
## 
## sigma^2 = 148.2:  log likelihood = -120.15
## AIC=248.31   AICc=249.85   BIC=254.04
## 
## Training set error measures:
##                     ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 0.1239895 7.027961 3.351554 -Inf  Inf 0.4508185 -0.03830072

Model Terbaik

AICKandidatModel <- c(model1$aic, model2$aic, model3$aic, model4$aic, model5$aic,model6$aic)
AICcKandidatModel <- c(model1$aicc, model2$aicc, model3$aicc, model4$aicc, model5$aicc,model6$aicc)
BICKandidatModel <- c(model1$bic, model2$bic, model3$bic, model4$bic, model5$bic)
KandidatModelARIMA <- c("ARIMA(0,1,1)(0,1,1)52", "ARIMA(2,1,1)(0,1,1)52", "ARIMA(1,1,0)(0,1,1)52", "ARIMA(1,1,1)(0,1,1)52", "ARIMA(1,1,3)(0,1,1)52","ARIMA(2,1,0)(0,1,1)52")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, BICKandidatModel)

## Warning in cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, :
## number of rows of result is not a multiple of vector length (arg 4)

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)(0,1,1)52 240.654151114995 241.543040003884  244.95611272845
## 2 ARIMA(2,1,1)(0,1,1)52 244.380582569051 246.780582569051 251.550518591477
## 3 ARIMA(1,1,0)(0,1,1)52 250.118965094095 251.007853982984 254.420926707551
## 4 ARIMA(1,1,1)(0,1,1)52 242.558292467139 244.096754005601  248.29424128508
## 5 ARIMA(1,1,3)(0,1,1)52 246.107939001978 249.607939001978 254.711862228889
## 6 ARIMA(2,1,0)(0,1,1)52 248.307958132539    249.846419671  244.95611272845

Model terbaik diperoleh berdasarkan nilai AIC, AICc, dan BIC dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu $ARIMA(0,1,1)(0,1,1)_{52}$

Pengujian parameter model

printstatarima <- function (x, digits = 4,se=TRUE){
       if (length(x$coef) > 0) {
         cat("\nCoefficients:\n")
         coef <- round(x$coef, digits = digits)
         if (se && nrow(x$var.coef)) {
           ses <- rep(0, length(coef))
           ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
           coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
           coef <- rbind(coef, s.e. = ses)
           statt <- coef[1,]/ses
           pval  <- 2*pt(abs(statt), df=length(x$residuals)-1, lower.tail = FALSE)
           coef <- rbind(coef, t=round(statt,digits=digits),sign.=round(pval,digits=digits))
           coef <- t(coef)
         }
         print.default(coef, print.gap = 2)
       }
     }

printstatarima(model1)

## 
## Coefficients:
##               s.e.        t  sign.
## ma1   -1    0.1188  -8.4175      0
## sma1   0  112.9645   0.0000      1

Model terbaik adalah model $ARIMA(0,1,1)(0,1,1)_{52}$ karena semua dugaan parameter berpengaruh nya.

Pendugaan parameter $\theta_1 = -1, \Theta_1 = 0$ dan sigma^2 estimated $108.5$ adalah nilai dugaan.

Model $ARIMA(0,1,1)(0,1,1)_{52}$, maka $X_t$ diperoleh dari penjabaran operator backshift sehingga untuk model$ARIMA(0,1,1)(0,1,1)_{52}$: $ p = 0, d = 1, q = 1, P = 0, D = 1, Q = 1$, dan $s = 52$ \[ \phi_p(B)\Phi_P(B)^s(1-B)^d(1-B^s)^DX_t = \mu + \theta_q(B)\Theta_Q(B^s)e_t \] \[ \phi_0(B)\Phi_0(B)^{52}(1-B)^1(1-B^{52})^1X_t = \mu + \theta_1(B)\Theta_1(B^{52})e_t \] \[ (1-B)(1-B^{52})X_t = \mu + (1-\theta_1B)(1-\Theta_1B^{52})e_t \] \[ (1-B^{52}-B+B^{53})X_t = \mu + (1-\Theta_1B^{52}-\theta_1B+\theta_1\Theta_1B^{53})e_t \] \[ X_t-X_{t-52}-X_{t-1}+X_{t-53} = \mu + e_t-\Theta_1e_{t-52}-\theta_1e_{t-1}+\theta_1\Theta_1e_{t-53} \] \[ X_t = \mu + X_{t-52} + X_{t-1} - X_{t-53}+ e_t-\Theta_1e_{t-52}-\theta_1e_{t-1}+\theta_1\Theta_1e_{t-53} \] \[ X_t = \mu + X_{t-52} + X_{t-1} - X_{t-53}+ e_t+e_{t-1} \]

Diagnostik Model

par(mar = c(3, 3, 2, 1))
tsdisplay(residuals(model1), lag.max=52, main='ARIMA(0,1,1)(0,1,1)52 Model Residuals', col="blue")

library(portes)
ljbtest2 <- LjungBox(residuals(model1),lags=seq(5,30,5),fitdf = 2)
ljbtest2

##  lags statistic df   p-value
##     5  3.636727  3 0.3034572
##    10  7.720647  8 0.4612229
##    15 13.767737 13 0.3904022
##    20 19.001207 18 0.3917491
##    25 26.782089 23 0.2653534
##    30 35.063144 28 0.1680367

Berdasarkan hasil uji LjungBox di atas terdapat autokorelasi pada sisaan, karena nilai $p-value$ ada lag yang tidak signifikan atau $p-value>\alpha=0.05$

Selanjutnya dilakukan uji asumsi formal terhadap kenormalan sisan dengan menggunakan uji Jarque Bera.

library(tseries)
jarque.bera.test(residuals(model1))

## 
##  Jarque Bera Test
## 
## data:  residuals(model1)
## X-squared = 222.41, df = 2, p-value < 2.2e-16

Berdasarkan hasil uji kenormalan dengan uji arque Bera sisaan tidak menyebar normal, karena nilai $p-value = 2.2e-16 < \alpha = 0.05$

Forecasting

Validasi model

ramalan_arima2 = forecast(model1, 5)
ramalan_arima2

##          Point Forecast      Lo 80    Hi 80      Lo 95   Hi 95
## 2025.615       23.06564  9.5112367 36.62004  2.3359709 43.7953
## 2025.635       18.26564  4.7112367 31.82004 -2.4640291 38.9953
## 2025.654       13.66564  0.1112367 27.22004 -7.0640291 34.3953
## 2025.673       20.36564  6.8112367 33.92004 -0.3640291 41.0953
## 2025.692       12.76564 -0.7887633 26.32004 -7.9640291 33.4953

accuracy(ramalan_arima2,test.ts)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.189563 6.119580 2.919111      -Inf      Inf 0.3926504
## Test set     -8.905637 9.762826 8.905637 -300.8765 300.8765 1.1978998
##                     ACF1 Theil's U
## Training set  0.04112259        NA
## Test set     -0.01771734  1.543505

plot(ramalan_arima2, col="blue")

forecast_arima2 <- cbind(ramalan_arima2$mean,ramalan_arima2$lower,ramalan_arima2$upper)
forecast_arima2

## Time Series:
## Start = c(2025, 33) 
## End = c(2025, 37) 
## Frequency = 52 
##          ramalan_arima2$mean ramalan_arima2$lower.80% ramalan_arima2$lower.95%
## 2025.615            23.06564                9.5112367                2.3359709
## 2025.635            18.26564                4.7112367               -2.4640291
## 2025.654            13.66564                0.1112367               -7.0640291
## 2025.673            20.36564                6.8112367               -0.3640291
## 2025.692            12.76564               -0.7887633               -7.9640291
##          ramalan_arima2$upper.80% ramalan_arima2$upper.95%
## 2025.615                 36.62004                  43.7953
## 2025.635                 31.82004                  38.9953
## 2025.654                 27.22004                  34.3953
## 2025.673                 33.92004                  41.0953
## 2025.692                 26.32004                  33.4953