Analisis Deret Waktu Praktikum 6

Khusnia Nurul Khikmah (G1501211049)

3/7/2022

Soal

Data yang digunakan pada tugas praktikum 6 ini adalah data time series yang sudah Anda pilih/tentukan pada tugas praktikum 1.

Lakukan identifikasi model tentatif ARIMA dengan menggunakan ACF, PACF dan EACF. Tuliskan beberapa model tentatifnya, berikan penjelasan Anda!

Silahkan upload tugas praktikum 6 di sini. Aturan penamaan folder yang sudah di-compressed adalah NIM_Nama_Prakt6 (7Z atau ZIP) Folder tersebut terdiri dari file data, script R beserta hasil.

Jawaban

Jawaban ini tersedia di : rpubs

Data

Data yang digunakan adalah data Daily Delhi climate data yang diambil dari Kaggle. Data ini memiliki beberapa variabel seperti mean temperature, humidity, wind speed, and mean pressure. Namun, pada tugas ini variabel yang digunakan adalah wind speed.

Membaca data

DailyDelhiClimateTest <- read.csv("D:/S2/IPB/STA542 Analisis Deret Waktu/1/Res/DailyDelhiClimateTest.csv", header=TRUE)
head(DailyDelhiClimateTest)
##         date meantemp humidity wind_speed meanpressure
## 1 2017-01-01 15.91304 85.86957   2.743478       59.000
## 2 2017-01-02 18.50000 77.22222   2.894444     1018.278
## 3 2017-01-03 17.11111 81.88889   4.016667     1018.333
## 4 2017-01-04 18.70000 70.05000   4.545000     1015.700
## 5 2017-01-05 18.38889 74.94444   3.300000     1014.333
## 6 2017-01-06 19.31818 79.31818   8.681818     1011.773

Visualisasi data

Sebelum dilakukan visualisasi data, variabel yang digunakan yaitu wind speed akan diubah menjadi time-series object dengan fungsits. Data yang digunakan terdiri dari 114 observasi atau yang diukur dalam 114 hari. Hasil dari visualisasi data tersaji pada plot berikut:

DailyDelhiClimateTest.ts <- ts(DailyDelhiClimateTest$wind_speed, start=1) #form a time-series object
head(DailyDelhiClimateTest.ts)
## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 2.743478 2.894444 4.016667 4.545000 3.300000 8.681818
ts.plot(DailyDelhiClimateTest.ts, col="blue", ylab = "Wind speed", xlab = "Day")
title(main = "Time Series Plot of Wind Speed", sub = "Source : Kaggle (https://www.kaggle.com/)", 
      cex.sub = 0.8)
points(DailyDelhiClimateTest.ts, pch = 20, col = "blue")

Split data

Data yang akan digunakan untuk menentukan model terbaik berdasarkan kandidat model akan dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian bagian pertama yaitu sebagai data training sebesar 80% dan bagian kedua sebagai data testing sebesar 20%. Splitting data ini menggunakan fungsi subset. Rasio splitting data yang digunakan 80:20 karena data wind speed dari Daily Delhi Climate memiliki range yang masuk akal yaitu antara 0 hingga 20 dan jika dilihat dari nilai rata-ratanya sebesar 8.1439241. Selain itu rasio splitting data 80:20 ratio adalah the best splitting ratio yang sebelumnya telah diteliti oleh [Rácz, A. et al.].

options(digits = 8)

# Training set
train.DailyDelhiClimateTest.ts = subset(DailyDelhiClimateTest.ts, start=1,end=94)

# Testing set
test.DailyDelhiClimateTest.ts = subset(DailyDelhiClimateTest.ts, start=95,end=114)

Visualisasi data training yang akan digunakan untuk mencari model terbaik

ts.plot(train.DailyDelhiClimateTest.ts, col="blue", ylab = "Wind speed", xlab = "Day")
title(main = "Time Series Plot of Wind Speed", sub = "Source : Kaggle (https://www.kaggle.com/)", cex.sub = 0.8)
points(train.DailyDelhiClimateTest.ts, pch = 20, col = "blue")

Model ARIMA

Model ARIMA terbaik dari data wind speed dari Daily Delhi Climate sebagai berikut:

1. Cek kestasioneran data

Mengidentifikasi kestasioneran data dilakukan dengan uji ADF.

Uji Augmented Dickey-Fuller (ADF) adalah salah satu alat identifikasi untuk menguji kestasioneran data dan memiliki hipotesis sebagai berikut:

\(H_{0}\) : Data tidak stasioner (unit roots mendekati 1)

\(H_{1}\) : Data stasioner (unit roots tidak mendekati 1)

adf.test(train.DailyDelhiClimateTest.ts, alternative = "stationary", k=trunc((length(train.DailyDelhiClimateTest.ts)-1)^(1/3)))
## Warning in adf.test(train.DailyDelhiClimateTest.ts, alternative =
## "stationary", : p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.DailyDelhiClimateTest.ts
## Dickey-Fuller = -4.96428, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5% = 0.05\) di atas dapat diketahui bersama bahwa, \(p-value = 0.01 < \alpha = 0.05\) sehingga keputusan yang diambil adalah Tolak \(H_{0}\) artinya data yang digunakan merupakan data yang stasioner

2. Identifikasi kandidat model

Identifikasi kandidat model diperoleh berdasarkan nilai \(p\) dan \(q\) dimana nilai \(d=0\)

acf(train.DailyDelhiClimateTest.ts, lag.max = 20, col = "blue")

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

pacf(train.DailyDelhiClimateTest.ts, lag.max = 20, col = "blue")

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

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

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

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

# ARIMA(0,0,1)
arima001 <- arima(train.DailyDelhiClimateTest.ts, order = c(0,0,1), include.mean = TRUE, method = "ML")
arima001
## 
## Call:
## arima(x = train.DailyDelhiClimateTest.ts, order = c(0, 0, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ma1  intercept
##       0.30050    7.97552
## s.e.  0.08785    0.41696
## 
## sigma^2 estimated as 9.7094:  log likelihood = -240.26,  aic = 484.53
# ARIMA(1,0,0)
arima100 <- arima(train.DailyDelhiClimateTest.ts, order = c(1,0,0), include.mean = TRUE, method = "ML")
arima100
## 
## Call:
## arima(x = train.DailyDelhiClimateTest.ts, order = c(1, 0, 0), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1  intercept
##       0.37489    7.97396
## s.e.  0.09854    0.50271
## 
## sigma^2 estimated as 9.4007:  log likelihood = -238.77,  aic = 481.55
# ARIMA(1,0,1)
arima101 <- arima(train.DailyDelhiClimateTest.ts, order = c(1,0,1), include.mean = TRUE, method = "ML")
arima101
## 
## Call:
## arima(x = train.DailyDelhiClimateTest.ts, order = c(1, 0, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1       ma1  intercept
##       0.51209  -0.15761    7.96742
## s.e.  0.21032   0.23360    0.53995
## 
## sigma^2 estimated as 9.3572:  log likelihood = -238.56,  aic = 483.12
# ARIMA(2,0,1)
arima201 <- arima(train.DailyDelhiClimateTest.ts, order = c(2,0,1), include.mean = TRUE, method = "ML")
arima201
## 
## Call:
## arima(x = train.DailyDelhiClimateTest.ts, order = c(2, 0, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1      ar2      ma1  intercept
##       0.09171  0.17727  0.25468    7.96839
## s.e.  0.79485  0.29812  0.79974    0.53636
## 
## sigma^2 estimated as 9.3355:  log likelihood = -238.45,  aic = 484.91
# ARIMA(3,0,1)
arima301 <- arima(train.DailyDelhiClimateTest.ts, order = c(3,0,1), include.mean = TRUE, method = "ML")
arima301
## 
## Call:
## arima(x = train.DailyDelhiClimateTest.ts, order = c(3, 0, 1), include.mean = TRUE, 
##     method = "ML")
## 
## Coefficients:
##           ar1       ar2       ar3       ma1  intercept
##       0.70722  -0.02716  -0.08113  -0.35959    7.97629
## s.e.  0.73026   0.28721   0.11474   0.72723    0.50064
## 
## sigma^2 estimated as 9.3156:  log likelihood = -238.36,  aic = 486.71
AICKandidatModel <- c(484.53, 481.55, 483.12, 484.91, 486.71)
KndidatModelARIMA <- c("ARIMA(0,0,1)", "ARIMA(1,0,0)", "ARIMA(1,0,1)", "ARIMA(2,0,1)", "ARIMA(3,0,1)")
compmodelARIMA <- cbind(KndidatModelARIMA, AICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA
##   Kandidat Model Nilai AIC
## 1   ARIMA(0,0,1)    484.53
## 2   ARIMA(1,0,0)    481.55
## 3   ARIMA(1,0,1)    483.12
## 4   ARIMA(2,0,1)    484.91
## 5   ARIMA(3,0,1)    486.71

Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu \(ARIMA(1,0,0)\)

Model \(ARIMA(1,0,0)\), maka \(Y_{t}\) diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(1,0,0)\): \[\phi_{p} (1-B)^{d} Y_{t}=\mu+\theta_{q}(B)e_{t}\] \[(1-\phi_{1}(B))(1-B)^{0} Y_{t}=\mu+\theta_{0}(B)e_{t}\] \[(1-\phi_{1}(B))(1) Y_{t}=\mu+(1)e_{t}\] \[(Y_{t}-\phi_{1}Y_{t-1})=\mu+e_{t}\] \[Y_{t}=\mu+\phi_{1}Y_{t-1}+e_{t}\]

Penduga parameter \(\widehat{\mu}=7.97396\), \(\widehat{\phi}_{1}=0.37489\), dengan sigma^2 estimated as 9.4007 adalah nilai dugaan bagi \(\sigma_{e}^2\). Model terbaik yang diperoleh yaitu model \(ARIMA(1,0,0)\).

\[Y_{t}=\mu+\phi_{1}Y_{t-1}+e_{t}=7.97396+0.37489Y_{t-1}+e_{t}\]

Oleh karena itu, model terbaik tersebut selanjutnya bisa digunakan untuk peramalan.

Model terbaik ini selanjutnya dilakukan uji hipotesis untuk mengetahui signifikansi sebuah parameter.

\(H_{0}\) : parameter \(= 0\)

\(H_{1}\) : parameter \(\neq\) \(0\)

Pengambilan keputusan yaitu Tolak \(H_{0}\) jika \(|t_{hitung}| > t_{\frac{\alpha}{2},(n-1)}\)

T-hitung parameter dugaan model \(ARIMA(1,0,0)\) \[t_{hitung}=\frac{parameter estimasi}{SE parameter estimasi}\]

intercept_mu <- 7.97396/0.50271
intercept_mu
## [1] 15.861948
ar_phi1 <- 0.37489/0.09854
ar_phi1
## [1] 3.8044449
prmtrdugaanarima100 <- c("mu", "phi(1)")
thitungprmtrarima100 <- c(abs(intercept_mu), abs(ar_phi1))
ttabel <- c("1.981180359", "1.981180359")
keputusan <- c("Signifikan", "Signifikan")
tksignimodelarima100 <- cbind(prmtrdugaanarima100, thitungprmtrarima100, ttabel, keputusan)
colnames(tksignimodelarima100) <- c("Parameter dugaan", "T-Hitung", "T-Tabel", "Keputusan")
tksignimodelarima100 <- as.data.frame(tksignimodelarima100)
tksignimodelarima100
##   Parameter dugaan         T-Hitung     T-Tabel  Keputusan
## 1               mu 15.8619482405363 1.981180359 Signifikan
## 2           phi(1) 3.80444489547392 1.981180359 Signifikan
lmtest::coeftest(arima100)
## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value   Pr(>|z|)    
## ar1       0.3748900  0.0985368  3.80457 0.00014205 ***
## intercept 7.9739594  0.5027069 15.86204 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan tabel di atas maka pengambilan keputusan dari semua paramater dugaan dari model \(ARIMA(1,0,0)\) semua parameternya signifikan.

Dibandingkan dengan autoarima

auto.arima(train.DailyDelhiClimateTest.ts, trace = TRUE)
## 
##  ARIMA(2,1,2) with drift         : Inf
##  ARIMA(0,1,0) with drift         : 509.6631
##  ARIMA(1,1,0) with drift         : 498.68014
##  ARIMA(0,1,1) with drift         : Inf
##  ARIMA(0,1,0)                    : 507.68091
##  ARIMA(2,1,0) with drift         : 497.92096
##  ARIMA(3,1,0) with drift         : 497.96692
##  ARIMA(2,1,1) with drift         : Inf
##  ARIMA(1,1,1) with drift         : Inf
##  ARIMA(3,1,1) with drift         : Inf
##  ARIMA(2,1,0)                    : 495.93179
##  ARIMA(1,1,0)                    : 496.70215
##  ARIMA(3,1,0)                    : 495.96707
##  ARIMA(2,1,1)                    : 482.9352
##  ARIMA(1,1,1)                    : 480.95136
##  ARIMA(0,1,1)                    : 486.38975
##  ARIMA(1,1,2)                    : 483.00556
##  ARIMA(0,1,2)                    : 482.44406
##  ARIMA(2,1,2)                    : Inf
## 
##  Best model: ARIMA(1,1,1)
## Series: train.DailyDelhiClimateTest.ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1       ma1
##       0.32134  -0.94195
## s.e.  0.10897   0.04151
## 
## sigma^2 estimated as 9.6898:  log likelihood=-237.34
## AIC=480.68   AICc=480.95   BIC=488.28

Diagnostic model

ARIMA100diag <- stats::arima(train.DailyDelhiClimateTest.ts, order = c(1,0,0), method = "ML")
checkresiduals(ARIMA100diag)

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

Berdasarkan pl0t di atas terlihat bahwa sisaan tidak mengikuti sebaran normal. Selanjutnya, ditinjau 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 <- arima100$residuals
# Uji formal normalitas data
ks.test(sisaan,"pnorm")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.303924, p-value = 3.3471e-08
## alternative hypothesis: two-sided

Hasil pengujian, Normalitas Data:

\(H_{0}\) : sisaan menyebar normal

\(H_{1}\) : sisaan tidak mengikuti sebaran normal

Hasil : \(p-value=0.000000033471 < \alpha=0.05\) yang berarti TOLAK \(H_{0}\). Sisaan tidak menyebar normal

# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.0783562, df = 93, p-value = 0.93771
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.60642519  0.65624802
## sample estimates:
##   mean of x 
## 0.024911413

Nilai Tengah Sisaan

\(H_{0}:\mu = 0\)

\(H_{1}:\mu \neq 0\)

Hasil : \(p-value=0.93771 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Nilai tengah sisaan sama dengan \(0\)

# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 14.7715, df = 23, p-value = 0.90264

Autokorelasi

\(H_{0}\) : tidak ada autokorelasi

\(H_{1}\) : terdapat autokorelasi

Hasil : \(p-value=0.90264 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Tidak terdapat gejala autokorelasi

Kesimpulan : Asumsi terpenuhi, kecuali sisaan tidak menyebar normal.

Peramalan

Selanjutnya yaitu melakukan peramalan hasil dari pemodelan dengan menggunakan data training, yaitu sebagai berikut:

ARIMAdata <- Arima(train.DailyDelhiClimateTest.ts, order=c(1,0,0), method="ML")
forecastARIMA100 <- forecast(ARIMAdata, h=20)
hasil_forecastARIMA100 <- as.data.frame(forecastARIMA100)
hasil_forecastARIMA100
##     Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
##  95     10.3530266 6.3812303 14.324823 4.2786883 16.427365
##  96      8.8658479 4.6241209 13.107575 2.3786863 15.353010
##  97      8.3083195 4.0300205 12.586618 1.7652259 14.851413
##  98      8.0993077 3.8158938 12.382722 1.5483915 14.650224
##  99      8.0209512 3.7368190 12.305083 1.4689364 14.572966
## 100      7.9915762 3.7073430 12.275809 1.4394069 14.543745
## 101      7.9805638 3.6963164 12.264811 1.4283728 14.532755
## 102      7.9764353 3.6921859 12.260685 1.4242413 14.528629
## 103      7.9748876 3.6906379 12.259137 1.4226932 14.527082
## 104      7.9743074 3.6900577 12.258557 1.4221129 14.526502
## 105      7.9740899 3.6898402 12.258340 1.4218954 14.526284
## 106      7.9740083 3.6897586 12.258258 1.4218138 14.526203
## 107      7.9739777 3.6897280 12.258227 1.4217832 14.526172
## 108      7.9739663 3.6897166 12.258216 1.4217718 14.526161
## 109      7.9739620 3.6897123 12.258212 1.4217675 14.526156
## 110      7.9739604 3.6897107 12.258210 1.4217659 14.526155
## 111      7.9739598 3.6897101 12.258209 1.4217653 14.526154
## 112      7.9739595 3.6897098 12.258209 1.4217650 14.526154
## 113      7.9739595 3.6897098 12.258209 1.4217650 14.526154
## 114      7.9739594 3.6897097 12.258209 1.4217649 14.526154

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

plot(forecastARIMA100)

data_forecastARIMA100 <- hasil_forecastARIMA100$`Point Forecast`
data_forecastARIMA100 <- as.data.frame(data_forecastARIMA100)
head(data_forecastARIMA100)
##   data_forecastARIMA100
## 1            10.3530266
## 2             8.8658479
## 3             8.3083195
## 4             8.0993077
## 5             8.0209512
## 6             7.9915762

Akurasi hasil peramalan

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

dataaktual <- as.data.frame(test.DailyDelhiClimateTest.ts)

compmodelforecast <- cbind(dataaktual, data_forecastARIMA100)
colnames(compmodelforecast) <- c("Data Aktual", "Data hasil peramalan ARIMA(1,0,0)")
compmodelforecast <- as.data.frame(compmodelforecast)
compmodelforecast
##    Data Aktual Data hasil peramalan ARIMA(1,0,0)
## 1   14.3846154                        10.3530266
## 2   13.5777778                         8.8658479
## 3    4.6500000                         8.3083195
## 4    8.3375000                         8.0993077
## 5   14.1250000                         8.0209512
## 6   19.3142857                         7.9915762
## 7   15.5125000                         7.9805638
## 8    9.4875000                         7.9764353
## 9    4.9444444                         7.9748876
## 10   1.3875000                         7.9743074
## 11   5.9666667                         7.9740899
## 12   2.1000000                         7.9740083
## 13   5.3666667                         7.9739777
## 14   7.8111111                         7.9739663
## 15   9.0250000                         7.9739620
## 16   5.5625000                         7.9739604
## 17   6.9625000                         7.9739598
## 18   8.8900000                         7.9739595
## 19   9.9625000                         7.9739595
## 20  12.1571429                         7.9739594
# Akurasi
akurasi.arima100 <- accuracy(data_forecastARIMA100$data_forecastARIMA100,dataaktual$x)
akurasi.arima100
##                  ME      RMSE      MAE        MPE      MAPE       ACF1
## Test set 0.81200923 4.4974205 3.547018 -36.418448 67.241012 0.57755553
##           Theil's U
## Test set 0.70311985
modelakursi <- c("ME", "RMSE", "MAE", "MPE", "MAPE")
ModelARIMA100 <- c(0.81200923, 4.497420, 3.547018, -36.418448, 67.241012)
Akurasiall <- cbind(modelakursi, ModelARIMA100)
colnames(Akurasiall) <- c("Akurasi", "Model ARIMA(1,0,0)")
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall
##   Akurasi Model ARIMA(1,0,0)
## 1      ME         0.81200923
## 2    RMSE            4.49742
## 3     MAE           3.547018
## 4     MPE         -36.418448
## 5    MAPE          67.241012

Berdasarkan hasil peramalan dengan menggunakan metode ARIMA yang diperoleh model terbaik yaitu \(ARIMA(1,0,0)\). Pemilihan model terbaik ini didapatkan berdasarkan nilai akurasi AIC yang terkecil dari model tentatif dan diperoleh akurasi dari peramalan berdasarkan nilai MAPE sebesar \(67.241012\).