Pertemuan 7 MPDW - Pendugaan Parameter, Diagnostik Model, dan Peramalan

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

library(tsibble)

## Warning: package 'tsibble' was built under R version 4.3.3

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

## Warning: package 'tseries' was built under R version 4.3.2

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

library(forecast)

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

library(MASS)

## Warning: package 'MASS' was built under R version 4.3.3

library(TSA)

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

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

library(TTR)

## Warning: package 'TTR' was built under R version 4.3.2

library(aTSA)

## Warning: package 'aTSA' was built under R version 4.3.2

## 
## Attaching package: 'aTSA'

## The following object is masked from 'package:forecast':
## 
##     forecast

## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test

## The following object is masked from 'package:graphics':
## 
##     identify

library(graphics)

I. Data Keseluruhan

Data Time Series

data <- read.csv("C:\\Users\\user\\Downloads\\dt_pt2.csv")
data

##        Tanggal Penumpang.Berangkat
## 1   2023-04-26                7538
## 2   2023-04-27                6954
## 3   2023-04-28                7420
## 4   2023-04-29                8118
## 5   2023-04-30                8970
## 6   2023-05-01                9472
## 7   2023-05-02                7850
## 8   2023-05-03                7657
## 9   2023-05-04                7515
## 10  2023-05-05                7646
## 11  2023-05-06                7336
## 12  2023-05-07                7525
## 13  2023-05-08                6769
## 14  2023-05-09                6509
## 15  2023-05-10                7562
## 16  2023-05-11                6471
## 17  2023-05-12                6810
## 18  2023-05-13                6686
## 19  2023-05-14                6932
## 20  2023-05-15                6952
## 21  2023-05-16                6071
## 22  2023-05-17                8257
## 23  2023-05-18                7842
## 24  2023-05-19                5854
## 25  2023-05-20                5879
## 26  2023-05-21                7127
## 27  2023-05-22                6446
## 28  2023-05-23                6292
## 29  2023-05-24                6786
## 30  2023-05-25                6979
## 31  2023-05-26                7588
## 32  2023-05-27                6390
## 33  2023-05-28                7530
## 34  2023-05-29                6593
## 35  2023-05-30                7177
## 36  2023-05-31                8949
## 37  2023-06-01                8128
## 38  2023-06-02                5582
## 39  2023-06-03                5132
## 40  2023-06-04                6595
## 41  2023-06-05                6304
## 42  2023-06-06                5984
## 43  2023-06-07                6616
## 44  2023-06-08                6899
## 45  2023-06-09                6989
## 46  2023-06-10                6830
## 47  2023-06-11                7631
## 48  2023-06-12                7022
## 49  2023-06-13                7091
## 50  2023-06-14                8598
## 51  2023-06-15                9241
## 52  2023-06-16                8803
## 53  2023-06-17                9322
## 54  2023-06-18               10038
## 55  2023-06-19                9265
## 56  2023-06-20                9051
## 57  2023-06-21                9895
## 58  2023-06-22                9556
## 59  2023-06-23                9250
## 60  2023-06-24                9446
## 61  2023-06-25                9704
## 62  2023-06-26                8872
## 63  2023-06-27                3673
## 64  2023-06-28                8940
## 65  2023-06-29                5570
## 66  2023-06-30                6883
## 67  2023-07-01                6791
## 68  2023-07-02                7998
## 69  2023-07-03                7496
## 70  2023-07-04                7273
## 71  2023-07-05                8557
## 72  2023-07-06                8148
## 73  2023-07-07                7973
## 74  2023-07-08                7397
## 75  2023-07-09                7762
## 76  2023-07-10                7674
## 77  2023-07-11                7265
## 78  2023-07-12                8403
## 79  2023-07-13                7763
## 80  2023-07-14                8054
## 81  2023-07-15                7800
## 82  2023-07-16                7758
## 83  2023-07-17                6046
## 84  2023-07-18                6562
## 85  2023-07-19                7811
## 86  2023-07-20                7045
## 87  2023-07-21                7294
## 88  2023-07-22                7001
## 89  2023-07-23                7560
## 90  2023-07-24                7610
## 91  2023-07-25                6927
## 92  2023-07-26                7188
## 93  2023-07-27                6824
## 94  2023-07-28                7521
## 95  2023-07-29                6204
## 96  2023-07-30                5918
## 97  2023-07-31                6013
## 98  2023-08-01                6121
## 99  2023-08-02                7177
## 100 2023-08-03                6217
## 101 2023-08-04                6922
## 102 2023-08-05                6672
## 103 2023-08-06                6539
## 104 2023-08-07                6457
## 105 2023-08-08                5438
## 106 2023-08-09                7027
## 107 2023-08-10                7257
## 108 2023-08-11                7154
## 109 2023-08-12                7083
## 110 2023-08-13                7122
## 111 2023-08-14                6481
## 112 2023-08-15                5898
## 113 2023-08-16                6762
## 114 2023-08-17                6139
## 115 2023-08-18                5795

data$Tanggal <- as.Date(data$Tanggal, format ="%Y-%m-%d")
summary(data)

##     Tanggal           Penumpang.Berangkat
##  Min.   :2023-04-26   Min.   : 3673      
##  1st Qu.:2023-05-24   1st Qu.: 6594      
##  Median :2023-06-22   Median : 7154      
##  Mean   :2023-06-22   Mean   : 7298      
##  3rd Qu.:2023-07-20   3rd Qu.: 7806      
##  Max.   :2023-08-18   Max.   :10038

Plot Time Series

ggplot(data) + geom_line(aes(x=Tanggal, y=Penumpang.Berangkat), size=1)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

data_ts <- ts(data$Penumpang.Berangkat)
plot(data_ts)

Plot deret waktu di atas menunjukkan bahwa data tidak stasioner dalam rataan, ditandai dengan adanya beberapa perubahan tren, terutama setelah titik sekitar lag 60, di mana terdapat lonjakan yang cukup besar diikuti oleh penurunan. Ini menandakan bahwa rataan data tidak konstan sepanjang waktu. Selain itu, lebar pita memiliki lebar yang sekilas hampir sama sehingga kemungkinan data stasioner dalam ragam.

Plot ACF

acf(data_ts)

Berdasarkan plot AFC tersebut, terlihat bahwa plot ACF tails off dan menurun secara perlahan seiring bertambahnya lag.

Uji ADF

tseries::adf.test(data_ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_ts
## Dickey-Fuller = -2.8291, Lag order = 4, p-value = 0.2324
## 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.2324 yang lebih besar dari taraf nyata 5% sehingga terima \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF.

Plot Box-Cox

index <- seq(1:115)
bc = boxcox(data_ts~index, lambda = seq(-2,3,by=0.01))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda

## [1] 0.65

#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]

##   [1] -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01
##  [13]  0.00  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09  0.10  0.11
##  [25]  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.20  0.21  0.22  0.23
##  [37]  0.24  0.25  0.26  0.27  0.28  0.29  0.30  0.31  0.32  0.33  0.34  0.35
##  [49]  0.36  0.37  0.38  0.39  0.40  0.41  0.42  0.43  0.44  0.45  0.46  0.47
##  [61]  0.48  0.49  0.50  0.51  0.52  0.53  0.54  0.55  0.56  0.57  0.58  0.59
##  [73]  0.60  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68  0.69  0.70  0.71
##  [85]  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.80  0.81  0.82  0.83
##  [97]  0.84  0.85  0.86  0.87  0.88  0.89  0.90  0.91  0.92  0.93  0.94  0.95
## [109]  0.96  0.97  0.98  0.99  1.00  1.01  1.02  1.03  1.04  1.05  1.06  1.07
## [121]  1.08  1.09  1.10  1.11  1.12  1.13  1.14  1.15  1.16  1.17  1.18  1.19
## [133]  1.20  1.21  1.22  1.23  1.24  1.25  1.26  1.27  1.28  1.29  1.30  1.31
## [145]  1.32  1.33  1.34  1.35  1.36  1.37  1.38  1.39  1.40  1.41  1.42  1.43
## [157]  1.44  1.45  1.46  1.47  1.48  1.49

Gambar di atas menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 0,65 dan pada selang kepercayaan 95% nilai memiliki batas bawah -0.12 dan batas atas 1.49. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data bangkitan stasioner dalam ragam.

II. Partisi Data - Data Awal Sebesar 50%

Data Time Series

dt_stas1 <- data_ts[1:58] |> ts()
head(dt_stas1)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 7538 6954 7420 8118 8970 9472

mean(dt_stas1)

## [1] 7431.448

var(dt_stas1)

## [1] 1307058

Plot Time Series

dt_stas1 |> as_tsibble() |> 
  ggplot(aes(x = index, y = value)) +
  geom_line() + theme_bw() +
  xlab("Obs") + ylab("Nilai")

Plot deret waktu di atas menunjukkan bahwa data tidak stasioner dalam rataan, ditandai dengan adanya beberapa perubahan tren, terutama setelah titik sekitar lag 30 dan 40, di mana terdapat penurunan yang cukup besar diikuti oleh lonjakan. Ini menandakan bahwa rataan data tidak konstan sepanjang waktu. Selain itu, lebar pita memiliki lebar yang sekilas hampir sama sehingga kemungkinan data stasioner dalam ragam.

Plot ACF

acf(dt_stas1)

Berdasarkan plot ACF, terlihat bahwa plot ACF pada data tersebut cenderung tails off dan menurun secara perlahan seiring bertambahnya lag.

Uji ADF

tseries::adf.test(dt_stas1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dt_stas1
## Dickey-Fuller = -1.3173, Lag order = 3, p-value = 0.8503
## 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.8503 yang lebih besar dari taraf nyata 5% sehingga terima \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF.

Plot Box-Cox

index <- seq(1:58)
bc = boxcox(dt_stas1~index, lambda = seq(-4,6,by=1))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda

## [1] -0.2626263

#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]

##  [1] -1.77777778 -1.67676768 -1.57575758 -1.47474747 -1.37373737 -1.27272727
##  [7] -1.17171717 -1.07070707 -0.96969697 -0.86868687 -0.76767677 -0.66666667
## [13] -0.56565657 -0.46464646 -0.36363636 -0.26262626 -0.16161616 -0.06060606
## [19]  0.04040404  0.14141414  0.24242424  0.34343434  0.44444444  0.54545455
## [25]  0.64646465  0.74747475  0.84848485  0.94949495  1.05050505  1.15151515
## [31]  1.25252525

Gambar di atas menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar -0.2626263 dan pada selang kepercayaan 95% nilai memiliki batas bawah -1.77777778 dan batas atas 1.25252525. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data bangkitan stasioner dalam ragam.

III. Partisi Data - Data Akhir Sebesar 50%

Data Time Series

dt_stas2 <- data_ts[59:115] |> ts()
head(dt_stas2)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 9250 9446 9704 8872 3673 8940

mean(dt_stas2)

## [1] 7162.018

var(dt_stas2)

## [1] 1132565

Plot Time Series

dt_stas2 |> as_tsibble() |> 
  ggplot(aes(x = index, y = value)) +
  geom_line() + theme_bw() +
  xlab("Obs") + ylab("Nilai")

Plot deret waktu di atas menunjukkan bahwa data kemungkinan stasioner dalam rataan ditandai dengan data yang menyebar di sekitar nilai tengahnya (7162.018). Namun, plot juga mengindikasikan adanya beberapa perubahan tren. Ini menandakan bahwa rataan data mungkin saja tidak konstan sepanjang waktu. Selain itu, lebar pita memiliki lebar yang sekilas hampir sama sehingga kemungkinan data stasioner dalam ragam. Sehingga, perlu dilakukan analisis lanjut.

Plot ACF

acf(dt_stas2)

Berdasarkan plot ACF tersebut, terlihat bahwa data stasioner dalam rataan ditandai dengan plot ACF yang cuts off pada lag ke-2.

Uji ADF

tseries::adf.test(dt_stas2)

## Warning in tseries::adf.test(dt_stas2): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dt_stas2
## Dickey-Fuller = -4.8825, Lag order = 3, p-value = 0.01
## 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.01 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF.

Plot Box-Cox

index <- seq(59:115)
bc = boxcox(dt_stas2~index, lambda = seq(0,4,by=1))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda

## [1] 2.020202

#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]

##  [1] 0.969697 1.010101 1.050505 1.090909 1.131313 1.171717 1.212121 1.252525
##  [9] 1.292929 1.333333 1.373737 1.414141 1.454545 1.494949 1.535354 1.575758
## [17] 1.616162 1.656566 1.696970 1.737374 1.777778 1.818182 1.858586 1.898990
## [25] 1.939394 1.979798 2.020202 2.060606 2.101010 2.141414 2.181818 2.222222
## [33] 2.262626 2.303030 2.343434 2.383838 2.424242 2.464646 2.505051 2.545455
## [41] 2.585859 2.626263 2.666667 2.707071 2.747475 2.787879 2.828283 2.868687
## [49] 2.909091 2.949495 2.989899 3.030303 3.070707 3.111111 3.151515

Gambar di atas menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 2.020202 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0.969697 dan batas atas 3.151515. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data bangkitan stasioner dalam ragam.

IV. Penanganan Ketidakstasioneran Data

train.diff<-diff(dt_stas1,differences = 1) 
plot.ts(train.diff, lty=1, main="Plot Difference")

Berdasarkan plot data deret waktu, terlihat bahwa data sudah stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data).

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cuts off pada lag ke 2. Hal ini menandakan data sudah stasioner dalam rataan dan ketidakstasioneran data telah berhasil tertangani.

tseries::adf.test(train.diff)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff
## Dickey-Fuller = -5.0654, Lag order = 3, p-value = 0.01
## 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.01 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) atau data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan

V. Identifikasi Model

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 2, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,2).

Plot PACF

pacf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 2, sehingga jika plot ACF dianggap tails of, maka model tentatifnya adalah ARIMA(2,1,0).

Jika baik plot ACF maupun plot PACF keduanya dianggap tails off, maka model yang terbentuk adalah ARIMA(2,1,2)

Plot EACF

eacf(train.diff)

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

Identifikasi model menggunakan plot EACF dilakukan dengan melihat ujung segitiga pada pola segitiga nol. Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,2), ARIMA(1,1,2), ARIMA(2,1,0), ARIMA(2,1,1), dan ARIMA(2,1,2).

VI. Pendugaan Parameter Model Tentatif

#ARIMA(0,1,2)
model1.da=Arima(train.diff, order=c(0,0,2),method="ML")
summary(model1.da) #AIC=929.41

## Series: train.diff 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##           ma1      ma2     mean
##       -0.1695  -0.4076  30.6923
## s.e.   0.1352   0.1383  45.7179
## 
## sigma^2 = 642727:  log likelihood = -460.71
## AIC=929.41   AICc=930.18   BIC=937.59
## 
## Training set error measures:
##                    ME     RMSE    MAE      MPE     MAPE      MASE        ACF1
## Training set 3.765875 780.3197 632.82 58.06482 158.2342 0.5729933 -0.06310222

lmtest::coeftest(model1.da) #terdapat parameter yang tidak signifikan

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value Pr(>|z|)   
## ma1       -0.16951    0.13517 -1.2541 0.209812   
## ma2       -0.40756    0.13831 -2.9467 0.003212 **
## intercept 30.69229   45.71785  0.6713 0.502003   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(1,1,2)
model2.da=Arima(train.diff, order=c(1,0,2),method="ML")
summary(model2.da) #AIC=929.54

## Series: train.diff 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1      ma2     mean
##       -0.4244  0.2099  -0.5325  34.0596
## s.e.   0.2204  0.1931   0.1107  49.7125
## 
## sigma^2 = 631075:  log likelihood = -459.77
## AIC=929.54   AICc=930.72   BIC=939.75
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 1.765858 766.0215 617.7286 33.44867 193.6986 0.5593287
##                      ACF1
## Training set -0.007168734

lmtest::coeftest(model2.da) #terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       -0.42439    0.22035 -1.9259   0.05411 .  
## ma1        0.20994    0.19312  1.0871   0.27700    
## ma2       -0.53254    0.11070 -4.8106 1.505e-06 ***
## intercept 34.05959   49.71254  0.6851   0.49326    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(2,1,0)
model3.da=Arima(train.diff, order=c(2,0,0),method="ML")
summary(model3.da) #AIC=927.92

## Series: train.diff 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2     mean
##       -0.2124  -0.4482  36.4984
## s.e.   0.1179   0.1165  62.1634
## 
## sigma^2 = 625894:  log likelihood = -459.96
## AIC=927.92   AICc=928.69   BIC=936.1
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -2.470455 770.0341 621.9286 58.38714 167.1116 0.5631317
##                     ACF1
## Training set -0.02837839

lmtest::coeftest(model3.da) #terdapat parameter yang tidak signifikan

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value Pr(>|z|)    
## ar1       -0.21240    0.11790 -1.8015 0.071630 .  
## ar2       -0.44821    0.11647 -3.8482 0.000119 ***
## intercept 36.49845   62.16337  0.5871 0.557111    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#model terpilih

#ARIMA(2,1,1)
model4.da=Arima(train.diff, order=c(2,0,1),method="ML")
summary(model4.da) #AIC=929.73

## Series: train.diff 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2      ma1     mean
##       -0.1018  -0.4312  -0.1397  36.1618
## s.e.   0.2833   0.1294   0.3218  57.9669
## 
## sigma^2 = 635431:  log likelihood = -459.86
## AIC=929.73   AICc=930.91   BIC=939.95
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE    MAPE      MASE
## Training set -1.791302 768.6606 621.6785 61.63982 161.007 0.5629052
##                      ACF1
## Training set -0.002583302

lmtest::coeftest(model4.da) #terdapat parameter tidak signifikan

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       -0.10180    0.28328 -0.3593 0.7193352    
## ar2       -0.43117    0.12940 -3.3319 0.0008625 ***
## ma1       -0.13971    0.32183 -0.4341 0.6641995    
## intercept 36.16176   57.96694  0.6238 0.5327364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(2,1,2)
model5.da=Arima(train.diff, order=c(2,0,2),method="ML")
summary(model5.da) #AIC=931.65

## Series: train.diff 
## ARIMA(2,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2      ma1      ma2     mean
##       -0.1306  -0.3601  -0.1119  -0.0972  35.4339
## s.e.   0.3035   0.2986   0.3231   0.3600  54.9929
## 
## sigma^2 = 646606:  log likelihood = -459.82
## AIC=931.65   AICc=933.33   BIC=943.91
## 
## Training set error measures:
##                      ME     RMSE     MAE      MPE     MAPE     MASE
## Training set -0.1565411 768.0402 621.664 58.94574 166.1434 0.562892
##                      ACF1
## Training set -0.003933697

lmtest::coeftest(model5.da) #tidak ada parameter yang signifikan

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)
## ar1       -0.130621   0.303520 -0.4304   0.6669
## ar2       -0.360125   0.298649 -1.2058   0.2279
## ma1       -0.111903   0.323072 -0.3464   0.7291
## ma2       -0.097203   0.359968 -0.2700   0.7871
## intercept 35.433870  54.992949  0.6443   0.5194

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(0,1,2) namun ada parameter yang tidak signifikan. Namun, model ARIMA(2,1,0) memiliki nilai AIC ketiga terkecil dan parameter model ARIMA(2,1,0) juga seluruhnya lebih banyak signifikan daripada model lain jika intercept dihilangkan. Sehingga, model yang dipilih adalah model ARIMA(2,1,0).

VII. Overfitting

Overfitting dilakukan dengan menambahkan nilai MA yang semula dari ARIMA(2,1,0) menjadi ARIMA(2,1,1). Dan juga mencoba menambahkan nilai AR yang semula dari ARIMA(2,1,0) menjadi ARIMA(3,1,0)

model3a.da<- forecast::Arima(train.diff, order=c(2,0,1),method="ML")
summary(model3a.da) #AIC= 929.73

## Series: train.diff 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2      ma1     mean
##       -0.1018  -0.4312  -0.1397  36.1618
## s.e.   0.2833   0.1294   0.3218  57.9669
## 
## sigma^2 = 635431:  log likelihood = -459.86
## AIC=929.73   AICc=930.91   BIC=939.95
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE    MAPE      MASE
## Training set -1.791302 768.6606 621.6785 61.63982 161.007 0.5629052
##                      ACF1
## Training set -0.002583302

lmtest::coeftest(model3a.da) #terdapat model tidak signifikan

## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       -0.10180    0.28328 -0.3593 0.7193352    
## ar2       -0.43117    0.12940 -3.3319 0.0008625 ***
## ma1       -0.13971    0.32183 -0.4341 0.6641995    
## intercept 36.16176   57.96694  0.6238 0.5327364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model3b.da<- forecast::Arima(train.diff, order=c(3,0,0),method="ML")
summary(model3b.da) #AIC= 929.76

## Series: train.diff 
## ARIMA(3,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2      ar3     mean
##       -0.2372  -0.4595  -0.0536  36.3056
## s.e.   0.1324   0.1195   0.1312  58.9774
## 
## sigma^2 = 635745:  log likelihood = -459.88
## AIC=929.76   AICc=930.93   BIC=939.97
## 
## Training set error measures:
##                     ME     RMSE     MAE      MPE     MAPE      MASE
## Training set -2.127473 768.8504 621.448 61.60605 160.8978 0.5626965
##                      ACF1
## Training set -0.005931505

lmtest::coeftest(model3b.da) ##terdapat model tidak signifikan

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1       -0.237183   0.132423 -1.7911 0.0732773 .  
## ar2       -0.459524   0.119506 -3.8452 0.0001205 ***
## ar3       -0.053646   0.131223 -0.4088 0.6826739    
## intercept 36.305619  58.977398  0.6156 0.5381683    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh nilai AIC overfitting lebih besar sehingga tetap dipakai model sebelum overfitting.

VIII. 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.

Eksplorasi Sisaan

#Eksplorasi 
sisaan.da <- model3.da$residuals 
par(mfrow=c(2,2)) 
qqnorm(sisaan.da) 
qqline(sisaan.da, col = "blue", lwd = 2) 
plot(c(1:length(sisaan.da)),sisaan.da) 
acf(sisaan.da) 
pacf(sisaan.da)

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan cenderung menyebar normal ditandai dengan titik titik yang cenderung mengikuti garis \(45^\circ\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Pada ARIMA(2,1,0), plot ACF signifikan pada lag 14 sedangkan plot PACF sisaan signifikan pada lag 12 dan 14 pada 17 lag awal. Kondisi ini akan diuji lebih lanjut dengan uji formal.

Uji Formal

#1) Sisaan Menyebar Normal 
ks.test(sisaan.da,"pnorm")  #tolak H0 >> sisaan tidak menyebar normal

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.da
## D = 0.52632, p-value = 3.775e-15
## 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 Jarque Bera. Hipotesis pada uji Jarque Bera adalah sebagai berikut.

\(H_0\): Sisaan menyebar normal

\(H_1\): Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 3.775e-15 yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal.

#2) Sisaan saling bebas/tidak ada autokorelasi 
Box.test(sisaan.da, type = "Ljung")  #tak tolak H0 >> sisaan saling bebas

## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 0.048363, df = 1, p-value = 0.8259

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.8259 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini sesuai dengan hasil eksplorasi.

#3) Sisaan homogen 
Box.test((sisaan.da)^2, type = "Ljung")  #tak tolak H0 >> sisaan homogen

## 
##  Box-Ljung test
## 
## data:  (sisaan.da)^2
## X-squared = 0.97205, df = 1, p-value = 0.3242

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\): Ragam sisaan homogen

\(H_1\): Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.3436 yang lebih dari taraf nyata 5% sehingga tak 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)  #tak tolak h0 >> nilai tengah sisaan sama dengan 0

## 
##  One Sample t-test
## 
## data:  sisaan.da
## t = -0.024008, df = 56, p-value = 0.9809
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -208.6031  203.6622
## sample estimates:
## mean of x 
## -2.470455

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_1\): nilai tengah sisaan tidak sama dengan 0

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

IX. Peramalan

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

#---FORECAST---#
ramalan.da <- forecast::forecast(model3.da, h = 30) 
ramalan.da

##    Point Forecast      Lo 80     Hi 80     Lo 95    Hi 95
## 59     -245.67984 -1259.5599  768.2002 -1796.276 1304.916
## 60      264.73625  -771.7611 1301.2336 -1320.450 1849.922
## 61      114.49740  -999.6656 1228.6605 -1589.468 1818.463
## 62      -82.36776 -1211.5127 1046.7772 -1809.246 1644.510
## 63       26.78519 -1111.5359 1165.1063 -1714.127 1767.697
## 64       91.83902 -1052.0580 1235.7360 -1657.600 1841.279
## 65       29.09783 -1115.5227 1173.7184 -1721.448 1779.644
## 66       13.26591 -1132.8848 1159.4166 -1739.620 1766.152
## 67       44.75006 -1101.4146 1190.9147 -1708.158 1797.658
## 68       45.15897 -1101.3414 1191.6593 -1708.262 1798.580
## 69       30.96048 -1115.5447 1177.4657 -1722.468 1784.389
## 70       33.79293 -1112.7722 1180.3581 -1719.727 1787.313
## 71       39.55529 -1107.0168 1186.1274 -1713.975 1793.086
## 72       37.06183 -1109.5187 1183.6424 -1716.482 1790.605
## 73       35.00867 -1111.5751 1181.5925 -1718.540 1788.557
## 74       36.56236 -1110.0223 1183.1470 -1716.988 1790.112
## 75       37.15261 -1109.4330 1183.7383 -1716.399 1790.704
## 76       36.33086 -1110.2548 1182.9165 -1717.221 1789.882
## 77       36.24084 -1110.3451 1182.8268 -1717.311 1789.793
## 78       36.62828 -1109.9576 1183.2142 -1716.924 1790.180
## 79       36.58633 -1109.9996 1183.1723 -1716.966 1790.138
## 80       36.42159 -1110.1644 1183.0076 -1717.130 1789.974
## 81       36.47538 -1110.1106 1183.0614 -1717.077 1790.027
## 82       36.53780 -1110.0482 1183.1238 -1717.014 1790.090
## 83       36.50043 -1110.0856 1183.0864 -1717.052 1790.052
## 84       36.48039 -1110.1056 1183.0664 -1717.072 1790.032
## 85       36.50139 -1110.0846 1183.0874 -1717.051 1790.053
## 86       36.50591 -1110.0801 1183.0919 -1717.046 1790.058
## 87       36.49554 -1110.0904 1183.0815 -1717.056 1790.047
## 88       36.49572 -1110.0903 1183.0817 -1717.056 1790.048

data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(2,1,0) cenderung stabil hingga akhir periode. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.

pt_1 <- dt_stas1[58] #nilai akhir data latih
hasil.forc.Diff <- data.ramalan.da
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1

#has.1 sama hasilnta dengan: cumsum(c(pt_1,hasil.forc.Diff))
ts.plot(dt_stas1,hasil)

perbandingan.da<-matrix(data=c(head(dt_stas2, n=30), hasil[-1]),
                     nrow = 30, ncol = 2)
colnames(perbandingan.da)<-c("Aktual","Hasil Forecast")
perbandingan.da

##       Aktual Hasil Forecast
##  [1,]   9250       9310.320
##  [2,]   9446       9575.056
##  [3,]   9704       9689.554
##  [4,]   8872       9607.186
##  [5,]   3673       9633.971
##  [6,]   8940       9725.810
##  [7,]   5570       9754.908
##  [8,]   6883       9768.174
##  [9,]   6791       9812.924
## [10,]   7998       9858.083
## [11,]   7496       9889.043
## [12,]   7273       9922.836
## [13,]   8557       9962.392
## [14,]   8148       9999.454
## [15,]   7973      10034.462
## [16,]   7397      10071.025
## [17,]   7762      10108.177
## [18,]   7674      10144.508
## [19,]   7265      10180.749
## [20,]   8403      10217.377
## [21,]   7763      10253.963
## [22,]   8054      10290.385
## [23,]   7800      10326.860
## [24,]   7758      10363.398
## [25,]   6046      10399.899
## [26,]   6562      10436.379
## [27,]   7811      10472.880
## [28,]   7045      10509.386
## [29,]   7294      10545.882
## [30,]   7001      10582.378

accuracy(ts(hasil[-1]), head(dt_stas2, n=30))

##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -2441.281 2762.423 2442.244 -36.51939 36.52931 0.1375262  1.199523

Diperoleh bahwa nilai MAPE mencapai 36.529%, yang mengindikasikan bahwa model prediksi rata-rata mengalami deviasi sekitar 37% dari nilai sebenarnya. Saat dilakukan pengecekan stasioneritas, ditemukan bahwa data tidak menunjukkan stasioneritas dalam rataan. Dalam analisis ini, masalah ketidakstasioneran dalam rataan telah diatasi. Namun, pada plot sisaan, acf dan pacf signifikan di lag tertentu. Ketika dilakukan uji formal, sisaan tidak menyebar normal. Hal ini dapat berpotensi menjadi salah satu penyebab mengapa nilai MAPE dari model masih relatif tinggi.