1. Packages

library(rio)
library(ggplot2)
library(tsibble)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(MASS)
library(readr)
library(TSA)
## 
## Attaching package: 'TSA'
## The following object is masked from 'package:readr':
## 
##     spec
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(forecast)
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA

2. Import Data

df <- read_csv("https://raw.githubusercontent.com/rafiadabhi/MPDW/main/Pertemuan%201/Data/crude_oil_clean_period_price.csv")
## Rows: 510 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (2): periode, price
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
data <- df$price[255:382]
print(data)
##   [1]  39.88  37.05  43.80  42.12  49.64  51.76  49.13  43.45  48.20  51.75
##  [11]  55.40  49.72  51.97  56.50  60.57  68.94  66.24  59.76  57.32  61.04
##  [21]  67.92  61.41  66.63  71.88  71.29  73.93  74.40  70.26  62.91  58.73
##  [31]  63.13  61.05  58.14  61.79  65.87  65.71  64.01  70.68  78.21  74.04
##  [41]  81.66  94.53  88.71  95.98  91.75 101.84 101.58 113.46 127.35 140.00
##  [51] 124.08 115.46 100.64  67.81  54.43  44.60  41.68  44.76  49.66  51.12
##  [61]  66.31  69.89  69.45  69.96  70.61  77.00  77.28  79.36  72.89  79.66
##  [71]  83.76  86.15  73.97  75.63  78.95  71.92  79.97  81.43  84.11  91.38
##  [81]  92.19  96.97 106.72 113.93 102.70  95.42  95.70  88.81  79.20  93.19
##  [91] 100.36  98.83  98.48 107.07 103.02 104.87  86.53  84.96  88.06  96.47
## [101]  92.19  86.24  88.91  91.82  97.49  92.05  97.23  93.46  91.97  96.56
## [111] 105.03 107.65 102.33  96.38  92.72  98.42  97.49 102.59 101.58  99.74
## [121] 102.71 105.37  98.17  95.96  91.16  80.54  66.15  53.27

2.1 Karakteristik Data

length(data) 
## [1] 128
head(data)     
## [1] 39.88 37.05 43.80 42.12 49.64 51.76
tail(data)     
## [1] 98.17 95.96 91.16 80.54 66.15 53.27
is.numeric(data)  
## [1] TRUE

2.2 Data Spliting

n <- length(data)          # total = 128
n_train <- floor(0.8 * n)  # 80% = 102

train_data <- data[1:n_train]
test_data  <- data[(n_train+1):n]

length(train_data)  # 102
## [1] 102
length(test_data)   # 26
## [1] 26

3. Eksplorasi Data

3.1 Plot Data Penuh

plot.ts(data, 
        lty = 1, 
        xlab = "Waktu", 
        ylab = "Harga", 
        main = "Plot Data Harga Crude Oil")

3.2 Plot Data Latih

plot.ts(train_data, lty=1, xlab="waktu", ylab="harga", main="Plot Harga Crude Oil Train")

3.3 Plot Data Uji

plot.ts(test_data, lty=1, xlab="waktu", ylab="Harga", main="Plot Harga Crude Oil Test")

4. Uji Stasioneritas Data

4.1 Plot ACF

acf(train_data)

Berdasarkan plot ACF, terlihat bahwa plot ACF data menurun secara perlahan (tails of slowly). Hal ini juga menjadi indikasi bahwa data tidak stasioner dalam rataan

4.2 Uji ADF

tseries::adf.test(train_data)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_data
## Dickey-Fuller = -3.7827, Lag order = 4, p-value = 0.02244
## alternative hypothesis: stationary

\(H_0\) : Data tidak stasioner dalam rataan \(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.0224 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa data stasioner dalam rataan.Namun jika menggunakan taraf nyata 1% hipotesis nya tak tolak \(H_0\), sehingga ada indikasi rataan tidak stasioner

4.3 Plot Box-Cox

# Uji Box-Cox
index <- seq(1:102)
bc <- boxcox(
  train_data - min(train_data) + 1 ~ index,
  lambda = seq(0, 5, by = 0.001)
)

# Nilai lambda optimum (dibulatkan)
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] 0.592
# Selang kepercayaan 95% untuk lambda
ci_lambda <- bc$x[bc$y > max(bc$y) - 0.5 * qchisq(0.95, 1)]

# Batas bawah dan batas atas
ci_lower <- min(ci_lambda)
ci_upper <- max(ci_lambda)

c(Batas_Bawah = ci_lower, Batas_Atas = ci_upper)
## Batas_Bawah  Batas_Atas 
##       0.390       0.813

Hasil uji Box–Cox menunjukkan bahwa nilai λ optimum adalah 0,592 dengan interval kepercayaan 95% sebesar (0,390–0,813). Karena selang tersebut tidak memuat nilai λ = 1, maka dapat disimpulkan bahwa data tidak stasioner dalam ragam sehingga diperlukan transformasi.

5. Penanganan Stasioneritas Data

5.1 Transformasi Box-Cox

# Data digeser agar positif
y_shift <- train_data - min(train_data) + 1

# Transformasi Box-Cox dengan lambda ~ 0.5 (akar kuadrat)
y_bc <- y_shift^0.5
# Uji Box-Cox
bc <- boxcox(
  y_bc ~ index,
  lambda = seq(0, 5, by = 0.001)
)

# Nilai lambda optimum (dibulatkan)
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] 1.183
# Selang kepercayaan 95% untuk lambda
ci_lambda <- bc$x[bc$y > max(bc$y) - 0.5 * qchisq(0.95, 1)]

# Batas bawah dan batas atas
ci_lower <- min(ci_lambda)
ci_upper <- max(ci_lambda)

c(Batas_Bawah = ci_lower, Batas_Atas = ci_upper)
## Batas_Bawah  Batas_Atas 
##       0.780       1.626

6. Identifikasi Model

6.1 Plot ACF untuk MA(q)

acf(y_bc)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung tails off pada lag ke 8

6.2 Plot PACF untuk AR(p=)

pacf(y_bc)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 2

eacf(y_bc)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x o o o  o  o  o 
## 1 x x o o o o o o o o o  o  o  o 
## 2 x o o o o o o o o o o  o  o  o 
## 3 x x 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 x x x o o o o o o o  o  o  o 
## 6 o x x o o o o o o o o  o  o  o 
## 7 o x x o o o o o o o o  o  o  o

Berdasarkan plot tersebut, terlihat bahwa plot EACF menunjukan kandidat model ARMA(1,1) atau ARMA(2,1)

7. Pendugaan Parameter Model ARMA

ARIMA(1,0,2)

model102.da <- Arima(y_bc, order = c(1, 0, 2))  
summary(model102.da) #AIC 200.63
## Series: y_bc 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1     ma2    mean
##       0.9072  0.2137  0.1558  5.8088
## s.e.  0.0480  0.1106  0.0946  0.8191
## 
## sigma^2 = 0.3835:  log likelihood = -95
## AIC=200   AICc=200.63   BIC=213.12
## 
## Training set error measures:
##                      ME      RMSE      MAE       MPE    MAPE      MASE
## Training set 0.04238235 0.6070134 0.491925 -1.242403 10.6233 0.9870068
##                     ACF1
## Training set 0.004864779
lmtest::coeftest(model102.da) #terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       0.907229   0.048018 18.8934 < 2.2e-16 ***
## ma1       0.213687   0.110639  1.9314   0.05344 .  
## ma2       0.155765   0.094619  1.6462   0.09972 .  
## intercept 5.808795   0.819130  7.0914 1.327e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (2,0,1)

model201.da <- Arima(y_bc, order = c(2, 0, 1))  
summary(model201.da) #AIC 198.71
## Series: y_bc 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1    mean
##       1.5894  -0.6399  -0.4361  5.9091
## s.e.  0.1877   0.1778   0.2128  0.6425
## 
## sigma^2 = 0.3787:  log likelihood = -94.36
## AIC=198.71   AICc=199.34   BIC=211.84
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.03075141 0.6032026 0.4882564 -1.470673 10.59297 0.9796459
##                     ACF1
## Training set -0.03636957
lmtest::coeftest(model201.da) # semua parameter signifikan
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1        1.58944    0.18771  8.4676 < 2.2e-16 ***
## ar2       -0.63988    0.17783 -3.5982 0.0003205 ***
## ma1       -0.43605    0.21282 -2.0490 0.0404648 *  
## intercept  5.90914    0.64248  9.1974 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (0,0,8)

model008.da <- Arima(y_bc, order = c(0, 0, 8))  
summary(model008.da) #AIC 196.02
## Series: y_bc 
## ARIMA(0,0,8) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7    ma8    mean
##       1.0659  1.2090  1.3131  1.5416  1.3346  1.1772  0.6991  0.401  5.9391
## s.e.  0.1099  0.1357  0.1339  0.1785  0.1999  0.1930  0.1557  0.127  0.5128
## 
## sigma^2 = 0.3317:  log likelihood = -88.01
## AIC=196.02   AICc=198.43   BIC=222.27
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.02272521 0.5499788 0.4466986 -1.348351 9.336127 0.8962638
##                     ACF1
## Training set 0.005214646
lmtest::coeftest(model008.da) # semua parameter signifikan
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1        1.06591    0.10986  9.7025 < 2.2e-16 ***
## ma2        1.20904    0.13569  8.9101 < 2.2e-16 ***
## ma3        1.31311    0.13392  9.8053 < 2.2e-16 ***
## ma4        1.54157    0.17846  8.6383 < 2.2e-16 ***
## ma5        1.33460    0.19994  6.6749 2.475e-11 ***
## ma6        1.17716    0.19300  6.0992 1.066e-09 ***
## ma7        0.69914    0.15573  4.4895 7.137e-06 ***
## ma8        0.40097    0.12702  3.1568  0.001595 ** 
## intercept  5.93913    0.51277 11.5823 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA (3,0,2)

model302.da <- Arima(y_bc, order = c(3, 0, 2))  
summary(model302.da) #AIC 201.34
## Series: y_bc 
## ARIMA(3,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ar3     ma1      ma2    mean
##       0.6220  0.8410  -0.5683  0.5147  -0.3041  5.8882
## s.e.  0.2537  0.2662   0.2179  0.2845   0.2514  0.6547
## 
## sigma^2 = 0.3808:  log likelihood = -93.67
## AIC=201.34   AICc=202.53   BIC=219.71
## 
## Training set error measures:
##                      ME      RMSE       MAE      MPE    MAPE      MASE
## Training set 0.03161691 0.5986676 0.4871075 -1.38691 10.4275 0.9773408
##                      ACF1
## Training set -0.007871802
lmtest::coeftest(model302.da) #terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.62199    0.25367  2.4519  0.014209 *  
## ar2        0.84096    0.26623  3.1588  0.001584 ** 
## ar3       -0.56832    0.21795 -2.6076  0.009118 ** 
## ma1        0.51466    0.28453  1.8088  0.070478 .  
## ma2       -0.30410    0.25143 -1.2095  0.226487    
## intercept  5.88819    0.65471  8.9936 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8. Ringkasan AIC

model_tentatif <- data.frame(
  Model = c("ARIMA(1,0,2)","ARIMA(2,0,1)","ARIMA(0,0,8)", "ARIMA(3,0,2)"),
  AIC   = c(AIC(model102.da),AIC(model201.da), AIC(model008.da), AIC(model302.da)),
  BIC   = c(BIC(model102.da),BIC(model201.da), BIC(model008.da), BIC(model302.da))
)

model_tentatif
##          Model      AIC      BIC
## 1 ARIMA(1,0,2) 200.0000 213.1249
## 2 ARIMA(2,0,1) 198.7148 211.8397
## 3 ARIMA(0,0,8) 196.0157 222.2654
## 4 ARIMA(3,0,2) 201.3390 219.7138

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(0,0,8) dan parameter model ARIMA(0,0,8) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(0,0,8).

9. Analisis Sisaan

Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan sisaan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.

9.1 Eksplorasi Sisaan

sisaan.da <- model008.da$residuals

par(mfrow = c(2,2),    
    mar = c(4,4,2,1),  
    oma = c(0,0,2,0))  

# Q-Q Plot
qqnorm(sisaan.da, main = "Normal Q-Q Plot")
qqline(sisaan.da, col = "blue", lwd = 2)

# Plot Sisaan vs Waktu
plot(sisaan.da, type = "p",
     main = "Plot Sisaan vs Waktu",
     xlab = "Waktu", ylab = "Sisaan")

# ACF
acf(sisaan.da, main = "ACF Sisaan")

# PACF
pacf(sisaan.da, main = "PACF Sisaan")

Berdasarkan plot kuantil-kuantil normal, secara eksploratif terlihat bahwa sisaan cenderung menyebar normal karena. Selain itu, lebar pita sebaran sisaan tampak seragam, yang mengindikasikan adanya homogenitas ragam. Sementara itu, plot ACF dan PACF sisaan dari model ARIMA(0,0,8) menunjukkan tidak adanya spike signifikan pada 20 lag pertama, sehingga secara visual sisaan dapat dianggap saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

10. Uji Formal

#1) Sisaan Menyebar Normal 
ks.test(sisaan.da,"pnorm")  #gagal tolak H0 > sisaan menyebar normal
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.da
## D = 0.16543, p-value = 0.007523
## alternative hypothesis: two-sided

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar \(0.007523\) yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.

#2) Sisaan saling bebas/tidak ada autokorelasi 
Box.test(sisaan.da, type = "Ljung")  #gagal tolak H0 > sisaan saling bebas
## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 0.002856, df = 1, p-value = 0.9574

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.9574 yang lebih besar dari taraf nyata 5% sehingga gagal tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini berbeda dengan eksplorasi.

#3) Sisaan homogen 
Box.test((sisaan.da)^2, type = "Ljung")  #gagal tolak H0 > sisaan homogen
## 
##  Box-Ljung test
## 
## data:  (sisaan.da)^2
## X-squared = 0.18846, df = 1, p-value = 0.6642

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_0\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.6642 yang lebih dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen.

#4) Nilai tengah sisaan sama dengan nol 
t.test(sisaan.da, mu = 0, conf.level = 0.95)  #gagal tolak h0 > nilai tengah sisaan sama dengan 0
## 
##  One Sample t-test
## 
## data:  sisaan.da
## t = 0.41562, df = 101, p-value = 0.6786
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.08574163  0.13119205
## sample estimates:
##  mean of x 
## 0.02272521

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_0\) : nilai tengah sisaan tidak sama dengan 0

Berdasarkan uji-ttersebut, didapat p-value sebesar 0.6786 yang lebih besar dari taraf nyata 5% sehingga gagal tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.

11. Overfitting

Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(0,0,8) Kandidat model overfitting adalah ARIMA(0,0,9) dan ARIMA(1,0,8). Karena model hasil overfitting sudah termasuk di dalam model tentatif dan hasilnya tidak lebih baik, maka model terbaik yang dipilih adalah ARIMA(0,1,1) dilanjutkan dengan peramalan.

model009.da <- Arima(y_bc, order = c(0, 0, 9))  
summary(model009.da) #AIC 197.94
## Series: y_bc 
## ARIMA(0,0,9) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8     ma9
##       1.0774  1.2271  1.3482  1.5791  1.3782  1.2174  0.7409  0.4327  0.0319
## s.e.  0.1164  0.1531  0.1879  0.2232  0.2529  0.2415  0.2199  0.1721  0.1172
##         mean
##       5.9387
## s.e.  0.5277
## 
## sigma^2 = 0.3351:  log likelihood = -87.97
## AIC=197.94   AICc=200.87   BIC=226.82
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.02218081 0.5498078 0.4475429 -1.307644 9.359995 0.8979576
##                      ACF1
## Training set -0.004110879
lmtest::coeftest(model009.da) #terdapat ada parameter signifikan
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1       1.077391   0.116441  9.2527 < 2.2e-16 ***
## ma2       1.227137   0.153076  8.0165 1.088e-15 ***
## ma3       1.348234   0.187879  7.1761 7.174e-13 ***
## ma4       1.579060   0.223153  7.0761 1.482e-12 ***
## ma5       1.378223   0.252926  5.4491 5.062e-08 ***
## ma6       1.217422   0.241502  5.0410 4.630e-07 ***
## ma7       0.740925   0.219936  3.3688 0.0007549 ***
## ma8       0.432654   0.172119  2.5137 0.0119477 *  
## ma9       0.031864   0.117232  0.2718 0.7857754    
## intercept 5.938748   0.527674 11.2546 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model108.da <- Arima(y_bc, order = c(1, 0, 8))  
summary(model108.da) #AIC 197.96
## Series: y_bc 
## ARIMA(1,0,8) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
##       0.0574  1.0167  1.1611  1.2692  1.4935  1.2780  1.1302  0.6620  0.3839
## s.e.  0.2462  0.2386  0.2499  0.2359  0.2774  0.3147  0.2771  0.2276  0.1516
##         mean
##       5.9386
## s.e.  0.5243
## 
## sigma^2 = 0.3352:  log likelihood = -87.98
## AIC=197.96   AICc=200.9   BIC=226.84
## 
## Training set error measures:
##                     ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.0223316 0.5498171 0.4473381 -1.315821 9.352584 0.8975469
##                      ACF1
## Training set -0.001502436
lmtest::coeftest(model108.da) #parameter signifikan
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       0.057413   0.246167  0.2332  0.815585    
## ma1       1.016695   0.238551  4.2620 2.026e-05 ***
## ma2       1.161078   0.249859  4.6469 3.369e-06 ***
## ma3       1.269219   0.235930  5.3796 7.464e-08 ***
## ma4       1.493511   0.277411  5.3837 7.295e-08 ***
## ma5       1.278042   0.314743  4.0606 4.895e-05 ***
## ma6       1.130159   0.277147  4.0778 4.546e-05 ***
## ma7       0.661950   0.227577  2.9087  0.003629 ** 
## ma8       0.383901   0.151601  2.5323  0.011332 *  
## intercept 5.938605   0.524268 11.3274 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Karena model hasil overfitting sudah termasuk hasilnya tidak lebih baik, maka model terbaik yang dipilih adalah ARIMA(0,0,8) dilanjutkan dengan peramalan.

12. Model Terbaik : ARIMA(0,0,8)

bestmodel.da=Arima(train_data, order=c(0,0,8),method="ML")
summary(bestmodel.da)
## Series: train_data 
## ARIMA(0,0,8) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8     mean
##       1.1636  1.2218  1.2532  1.2988  0.9763  0.7201  0.3533  0.1468  75.0702
## s.e.  0.0971  0.1473  0.1776  0.1941  0.1884  0.1647  0.1457  0.1201   5.3257
## 
## sigma^2 = 50.82:  log likelihood = -341.94
## AIC=703.88   AICc=706.3   BIC=730.13
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.213402 6.807347 5.475297 -0.7255524 7.375013 0.9439523
##                      ACF1
## Training set 0.0002475575
lmtest::coeftest(bestmodel.da)
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ma1        1.163644   0.097147 11.9782 < 2.2e-16 ***
## ma2        1.221777   0.147262  8.2966 < 2.2e-16 ***
## ma3        1.253240   0.177581  7.0573 1.698e-12 ***
## ma4        1.298809   0.194129  6.6904 2.225e-11 ***
## ma5        0.976320   0.188355  5.1834 2.179e-07 ***
## ma6        0.720088   0.164685  4.3725 1.228e-05 ***
## ma7        0.353317   0.145658  2.4257   0.01528 *  
## ma8        0.146775   0.120105  1.2221   0.22169    
## intercept 75.070179   5.325694 14.0958 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13. Peramalan

Peramalan dilakukan menggunakan fungsi forecast(). Contoh peramalan berikut ini dilakukan untuk 30 hari ke depan.

# ---FORECAST---
ramalan.da <- forecast::forecast(bestmodel.da, h = length(test_data))
data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

# 3. Kembalikan hasil forecast ke skala asli
ramalan.da$mean <- InvBoxCox(ramalan.da$mean, lambda)
ramalan.da$lower <- InvBoxCox(ramalan.da$lower, lambda)
ramalan.da$upper <- InvBoxCox(ramalan.da$upper, lambda)
autoplot(ramalan.da) +
  ggtitle("Hasil Forecast Model ARIMA (Skala Asli)") +
  xlab("Waktu") +
  ylab("Nilai (Skala Asli)") +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    panel.grid.minor = element_blank()
  )

hasil_akurasi <- accuracy(ramalan.da)
hasil_akurasi
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.213402 6.807347 5.475297 -0.7255524 7.375013 0.9439523
##                      ACF1
## Training set 0.0002475575