MPDW Project 6

Packages

library(ggplot2)

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

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(MASS)
library(forecast)

## Warning: package 'forecast' 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.3

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)

Data

Import Data

data <- read.csv ("C:/Users/Muhammad Rizqa Salas/Downloads/Semester 5/MPDW/Project/Data Air MPDW.csv")
data

##     Periode Air
## 1         1 185
## 2         2 194
## 3         3 192
## 4         4 198
## 5         5 210
## 6         6 208
## 7         7 214
## 8         8 220
## 9         9 225
## 10       10 219
## 11       11 213
## 12       12 197
## 13       13 193
## 14       14 167
## 15       15 155
## 16       16 144
## 17       17 140
## 18       18 143
## 19       19 140
## 20       20 157
## 21       21 170
## 22       22 178
## 23       23 184
## 24       24 186
## 25       25 190
## 26       26 195
## 27       27 201
## 28       28 208
## 29       29 215
## 30       30 222
## 31       31 225
## 32       32 225
## 33       33 230
## 34       34 225
## 35       35 223
## 36       36 203
## 37       37 193
## 38       38 167
## 39       39 155
## 40       40 144
## 41       41 140
## 42       42 143
## 43       43 140
## 44       44 157
## 45       45 170
## 46       46 178
## 47       47 184
## 48       48 186
## 49       49 190
## 50       50 195
## 51       51 201
## 52       52 208
## 53       53 215
## 54       54 222
## 55       55 225
## 56       56 225
## 57       57 230
## 58       58 225
## 59       59 223
## 60       60 203
## 61       61 199
## 62       62 186
## 63       63 170
## 64       64 157
## 65       65 148
## 66       66 136
## 67       67 135
## 68       68 133
## 69       69 147
## 70       70 155
## 71       71 163
## 72       72 169
## 73       73 184
## 74       74 185
## 75       75 190
## 76       76 194
## 77       77 202
## 78       78 207
## 79       79 215
## 80       80 225
## 81       81 227
## 82       82 227
## 83       83 223
## 84       84 214
## 85       85 200
## 86       86 188
## 87       87 170
## 88       88 155
## 89       89 142
## 90       90 135
## 91       91 137
## 92       92 135
## 93       93 180
## 94       94 158
## 95       95 165
## 96       96 175
## 97       97 178
## 98       98 185
## 99       99 188
## 100     100 190

dim(data)

## [1] 100   2

Eksplorasi Data

Plot Data Penuh

data.ts <- ts(data$Air) 
ts.plot(data.ts, xlab="Time Period", ylab="Kurs", main="Time Series Plot")

Split Data

train <- data[1:80,]
test <- data[81:100,]
train.ts <- ts(train$Air)
test.ts <- ts(test$Air)

Data Latih

Plot

plot(train.ts, col="blue",main="Plot data training")

Berdasarkan plot tersebut, terdapat pola naik turun yang berulang, mengindikasikan adanya tren yang kuat dalam data, yang mana rataan data tidak konstan sepanjang waktu. Maka dari itu, data training tersebut cenderung tidak stasioner dalam rataan.

Data Uji

Plot

plot(test.ts, col="blue",main="Plot data testing")

Identifikasi Stasioneritas Data

Plot ACF

acf(train.ts)

Berdasarkan plot ACF, terlihat bahwa plot ACF data menurun tails off.

Uji ADF

tseries::adf.test(train.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -4.0619, Lag order = 4, p-value = 0.01112
## 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.01112 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, sehingga perlu dilakukannya penanganan terhadap ketidakstasioneran dalam rataan.

Plot Boxcox

index <- seq(1:80)
bc = boxcox(train.ts~index, lambda = seq(0,10,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.5050505 0.6060606 0.7070707 0.8080808 0.9090909 1.0101010 1.1111111
##  [8] 1.2121212 1.3131313 1.4141414 1.5151515 1.6161616 1.7171717 1.8181818
## [15] 1.9191919 2.0202020 2.1212121 2.2222222 2.3232323 2.4242424 2.5252525
## [22] 2.6262626 2.7272727 2.8282828 2.9292929 3.0303030 3.1313131 3.2323232
## [29] 3.3333333 3.4343434 3.5353535

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 2,02 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0,5 dan batas atas 3,53. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data bangkitan stasioner dalam ragam.

Penanganan Ketidakstasioneran Data

Plot

train.diff<-diff(train.ts,differences = 1) 
plot.ts(train.diff, lty=1, xlab="waktu", ylab="Data Difference 1 Air", main="Plot Difference Air")

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

Plot ACF

acf(train.diff)

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

Uji ADF

tseries::adf.test(train.diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff
## Dickey-Fuller = -3.7856, Lag order = 4, p-value = 0.02387
## alternative hypothesis: stationary

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.02387 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

Identifikasi Model

Plot ACF

acf(train.diff)

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

Plot PACF

pacf(train.diff)

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

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

Plot EACF

eacf(train.diff)

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

Pendugaan Parameter Model Tentatif

ARIMA(0,1,3)

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

## Series: train.diff 
## ARIMA(0,1,3) 
## 
## Coefficients:
##           ma1     ma2     ma3
##       -0.5253  0.4216  0.0590
## s.e.   0.1084  0.1093  0.1039
## 
## sigma^2 = 49.1:  log likelihood = -261.32
## AIC=530.64   AICc=531.19   BIC=540.07
## 
## Training set error measures:
##                      ME     RMSE     MAE  MPE MAPE      MASE       ACF1
## Training set 0.01411674 6.827184 5.57133 -Inf  Inf 0.9072311 0.00988099

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

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.525338   0.108415 -4.8456 1.262e-06 ***
## ma2  0.421569   0.109311  3.8566  0.000115 ***
## ma3  0.058997   0.103902  0.5678  0.570159    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,2)

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

## Series: train.diff 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.1550  -0.6732  0.5009
## s.e.  0.2379   0.2094  0.1258
## 
## sigma^2 = 49.05:  log likelihood = -261.28
## AIC=530.56   AICc=531.11   BIC=539.99
## 
## Training set error measures:
##                      ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 0.01421396 6.823736 5.575637 -Inf  Inf 0.9079325 0.004119044

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

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.15495    0.23791  0.6513  0.514851    
## ma1 -0.67324    0.20942 -3.2148  0.001305 ** 
## ma2  0.50092    0.12584  3.9806 6.875e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,2)

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

## Series: train.diff 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1     ar2      ma1     ma2
##       0.1536  0.0270  -0.6730  0.4825
## s.e.  0.2383  0.2215   0.2129  0.2060
## 
## sigma^2 = 49.7:  log likelihood = -261.27
## AIC=532.55   AICc=533.38   BIC=544.33
## 
## Training set error measures:
##                      ME     RMSE     MAE  MPE MAPE      MASE        ACF1
## Training set 0.01445081 6.823303 5.58623 -Inf  Inf 0.9096575 0.008623022

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

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ar1  0.153608   0.238252  0.6447 0.519102   
## ar2  0.027028   0.221546  0.1220 0.902902   
## ma1 -0.673049   0.212886 -3.1615 0.001569 **
## ma2  0.482534   0.206041  2.3419 0.019184 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(3,1,3)

model4.da=Arima(train.diff, order=c(3,1,3),method="ML")
summary(model4.da) #AIC=531.48

## Series: train.diff 
## ARIMA(3,1,3) 
## 
## Coefficients:
##           ar1      ar2      ar3      ma1    ma2     ma3
##       -0.1882  -0.4145  -0.4826  -0.3519  0.778  0.2790
## s.e.   0.2462   0.1316   0.1473   0.2610  0.204  0.2226
## 
## sigma^2 = 46.79:  log likelihood = -258.74
## AIC=531.48   AICc=533.08   BIC=547.98
## 
## Training set error measures:
##                      ME     RMSE      MAE  MPE MAPE     MASE       ACF1
## Training set 0.02017721 6.530118 5.300154 -Inf  Inf 0.863073 0.02044634

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

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)    
## ar1 -0.18819    0.24615 -0.7645 0.444566    
## ar2 -0.41445    0.13162 -3.1488 0.001640 ** 
## ar3 -0.48260    0.14728 -3.2767 0.001050 ** 
## ma1 -0.35189    0.26101 -1.3482 0.177591    
## ma2  0.77803    0.20402  3.8135 0.000137 ***
## ma3  0.27899    0.22265  1.2531 0.210179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(3,1,1)

model5.da=Arima(train.diff, order=c(1,1,3),method="ML")
summary(model5.da) #AIC=530.64

## Series: train.diff 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1     ma2      ma3
##       0.1809  -0.6992  0.5157  -0.0128
## s.e.  0.5920   0.5857  0.3362   0.2730
## 
## sigma^2 = 49.71:  log likelihood = -261.28
## AIC=532.56   AICc=533.4   BIC=544.35
## 
## Training set error measures:
##                      ME     RMSE      MAE  MPE MAPE     MASE        ACF1
## Training set 0.01424017 6.823663 5.576751 -Inf  Inf 0.908114 0.004707977

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

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.180852   0.592012  0.3055   0.7600
## ma1 -0.699245   0.585747 -1.1938   0.2326
## ma2  0.515677   0.336213  1.5338   0.1251
## ma3 -0.012814   0.273042 -0.0469   0.9626