Pembangkitan Data

# Import Library
library(tseries)
## Warning: package 'tseries' was built under R version 4.4.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
## Warning: package 'forecast' was built under R version 4.4.3
library(TSA)
## Warning: package 'TSA' was built under R version 4.4.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
# Parameter

# Panjang Data
n <- 400

# Koefisien Ar dan MA
Ar <- 0.8
Ma <- 0.9
set.seed(2)

# Generate Data Arima (1,1,1)
data <- arima.sim(model = list(order = c(1,1,1), ar = Ar, ma = Ma), n = n)

Analisis ARIMA Dengan Data Simulasi

1. Plot data

# Visualisasi Data Simulasi
ts.plot(data, ylab = "Data",
              main = "Data Simulasi ARIMA(1,1,1)")

Dari visualisasi di atas, pola data simulasi cenderung naik seiring bertambahanya waktu. Hal ini menandakan bahwa plot memiliki trend naik dan tidak stasioner.

# Plot ACF dan PACF data

# Acf
Acf(data, type = 'correlation', plot = T, main = "Data Simulasi ARIMA(1,1,1)")

# PACF
Acf(data, type = 'partial', plot = T, main = "Data Simulasi ARIMA(1,1,1)")

  • Plot ACF

Pada plot ACF dari data simulasi terlihat ACF beberapa lag mencapai batas signifkansi (garis biru putus-putus), hal ini menendakan bahwa ada autokorelasi yang kuat. Pola plot ACF juga turun secara perlahan (tailing off) yang menandakan bahwa data simulasi tidak stasioner.

  • Plot PACF

Pada plot PACF terlihat spike signifikan pada lag 1, sedangkan lag lainnya relatif tidak signifikan. Hal ini menunjukkan adanya pengaruh kuat dari lag pertama

2. Uji Kestasioneran Data Simulasi

Uji kestasioneran data dapat diuji dengan uji Augmented Duckey Fuller ataupun uji Kpss.

# Uji Stationary

# Uji Augmented Duckey Fuller
adf.test(data, alternative = 'stationary')
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data
## Dickey-Fuller = -1.9589, Lag order = 7, p-value = 0.5946
## alternative hypothesis: stationary

\(\alpha: 0.05\)

\(H_0:\) Data simulasi tidak stasioner

\(H_1:\) Data simulasi stasioner

P-value yang didapatkan dalam uji ADF adalah 0.5946. \(p-value > \alpha\) maka gagal tolak \(H_0\) sehingga dapat disimpulkan bahwa data simulasi tidak stasioner pada taraf signifikansi 5%.

# Uji Kpss
kpss.test(data, null = 'Trend')
## Warning in kpss.test(data, null = "Trend"): p-value smaller than printed
## p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  data
## KPSS Trend = 1.3964, Truncation lag parameter = 5, p-value = 0.01

\(\alpha: 0.05\)

\(H_0:\) Data simulasi stasioner

\(H_1:\) Data simulasi tidak stasioner

P-value yang didapatkan dalam uji Kpss adalah 0.01. \(p-value < \alpha\) maka tolak \(H_0\) sehingga dapat disimpulkan bahwa data simulasi tidak stasioner pada taraf signifikansi 5%.

Dari dua uji di atas didapatkan kesimpulan bahwa data simulasi tidak stasioner. ARIMA bekerja pada data yang stasioner sehingga data perlu di-differencing.

3. Differencing

# Differencing 1
dif1 <- diff(data, differences = 1)

# Visualisasi
ts.plot(dif1)

Setelah data simulasi di-differencing 1 terlihat bahwa pola data cenderung menyebar di sekitar rata-rata dan tidak terlihat adanya tren. Selain itu, variansi data relatif stabil dari waktu ke waktu. Hal ini menandakan bahwa data differencing 1 telah stasioner dan siap digunakan untuk identifikasi parameter ARIMA p dan q. Agar bisa terbukti stasioner secara statistik, perlu dilakukan uji formal.

# Uji Stationary

# Uji Augmented Duckey Fuller
adf.test(dif1, alternative = 'stationary')
## Warning in adf.test(dif1, alternative = "stationary"): p-value smaller than
## printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dif1
## Dickey-Fuller = -5.71, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

\(\alpha: 0.05\)

\(H_0:\) Data differencing tidak stasioner

\(H_1:\) Data differencing stasioner

P-value yang didapatkan dalam uji ADF adalah 0.01. \(p-value < \alpha\) maka tolak \(H_0\) sehingga dapat disimpulkan bahwa data differencing telah stasioner pada taraf signifikansi 5%.

# Uji Kpss
kpss.test(dif1, null = 'Trend')
## Warning in kpss.test(dif1, null = "Trend"): p-value greater than printed
## p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  dif1
## KPSS Trend = 0.083421, Truncation lag parameter = 5, p-value = 0.1

\(\alpha: 0.05\)

\(H_0:\) Data differencing stasioner

\(H_1:\) Data differencing tidak stasioner

P-value yang didapatkan dalam uji Kpss adalah 0.1. \(p-value > \alpha\) maka gagal tolak \(H_0\) sehingga dapat disimpulkan bahwa data differencing telah stasioner pada taraf signifikansi 5%.

4. Plot ACF dan PACF diff1

Plot ACF dan PACF differencing yang telah stsioner dapat kita jadikan salah satu cara menetapkan parameter p dan q.

# Plot ACF dan PACF diff1

# Acf
Acf(dif1, type = 'correlation', plot = T)

# PACF
Acf(dif1, type = 'partial', plot = T)

  • Plot ACF

Berdasarkan plot ACF, terlihat bahwa nilai autokorelasi menurun secara bertahap tanpa adanya cut-off yang jelas. Hal ini menunjukkan bahwa model yang sesuai bukan murni Moving Average, melainkan model campuran dengan orde yang relatif kecil. Oleh karena itu, nilai q yang dipertimbangkan adalah dalam rentang kecil, seperti q = 0 sampai q = 2. Namun demikian, untuk memastikan model terbaik, dilakukan eksplorasi terhadap beberapa nilai q yang lebih besar

  • Plot PACF

Meskipun terdapat beberapa spike signifikan pada beberapa lag di plot PACF, spike tersebut tidak menunjukkan pola cut-off yang jelas dan cenderung menurun secara bertahap. Hal ini mengindikasikan bahwa spike pada lag yang lebih tinggi kemungkinan merupakan efek noise. Oleh karena itu, komponen autoregressive yang lebih sesuai adalah dengan orde kecil, seperti p = 1 atau p = 2. Namun, untuk memastikan model terbaik, dilakukan eksplorasi terhadap beberapa kombinasi parameter yang lebih luas.

5. Pemilihan Model

# Membuat data frame
Tabel <- data.frame(Mod = character(),
                    AIC = numeric())
# Pembuatan tabel Model dan AIC-nya
for (i in 0:4){
  for (a in 0:4){
    Model <- arima(data, order = c(i,1,a), method = 'ML')
    Mod <- paste("ARIMA(", i, ",1,", a, ")", sep = "")
    AIC <- Model$aic
    
    # Masukan hasil ke tabel
    Tabel <- rbind(Tabel, data.frame(Mod = Mod,
                            AIC =AIC))
  }
}
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg = xreg,
## : possible convergence problem: optim gave code = 1
# Penguurtan Model Berdasarkan AIC
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
## 
## 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
Tabel_Terurut <- Tabel %>% arrange(AIC)
Tabel_Terurut
##             Mod      AIC
## 1  ARIMA(4,1,3) 1164.196
## 2  ARIMA(1,1,1) 1167.709
## 3  ARIMA(2,1,1) 1167.710
## 4  ARIMA(1,1,2) 1167.860
## 5  ARIMA(1,1,3) 1169.366
## 6  ARIMA(3,1,1) 1169.503
## 7  ARIMA(2,1,2) 1169.602
## 8  ARIMA(4,1,1) 1171.215
## 9  ARIMA(4,1,4) 1171.266
## 10 ARIMA(1,1,4) 1171.306
## 11 ARIMA(3,1,2) 1171.320
## 12 ARIMA(2,1,3) 1171.376
## 13 ARIMA(4,1,2) 1173.318
## 14 ARIMA(3,1,3) 1173.320
## 15 ARIMA(2,1,4) 1173.458
## 16 ARIMA(3,1,4) 1175.301
## 17 ARIMA(0,1,4) 1207.263
## 18 ARIMA(4,1,0) 1209.325
## 19 ARIMA(3,1,0) 1223.569
## 20 ARIMA(0,1,3) 1253.856
## 21 ARIMA(2,1,0) 1267.671
## 22 ARIMA(0,1,2) 1353.747
## 23 ARIMA(1,1,0) 1391.953
## 24 ARIMA(0,1,1) 1598.519
## 25 ARIMA(0,1,0) 2089.854

Berdasarkan nilai AIC, model yang memiliki nilai AIC terkecil adalah Model ARIMA(4,1,3). Model ini berbeda dengan model untuk membuat data simulasi, ARIMA(1,1,1). Model terbaik berdasarkan AIC adalah model ARIMA(4,1,3). Perbedaan model terbaik yang didapatkan dari data simulasi kemungkinana disebabkan:

  • Pengaruh Komponen Acak

Data yang dibangkitkan dari model ARIMA mengandung unsur acak, sehingga pola yang terbentuk pada data sampel tidak selalu mencerminkan struktur model teoritis secara sempurna.

  • Keterbatasan Pada Data Sampel

Meskipun jumlah data cukup besar, tetap terdapat kemungkinan bahwa sampel yang diperoleh belum sepenuhnya merepresentasikan karakteristik populasi, sehingga estimasi parameter menjadi kurang akurat.

  • Kecenderungan AIC Memilih Model Yang Kompleks

Kriteria AIC cenderung memilih model dengan jumlah parameter lebih banyak karena model tersebut memiliki fleksibilitas yang lebih tinggi dalam menyesuaikan data, meskipun peningkatan kualitas model tidak selalu signifikan.

  • Kemungkinan Overfitting

Model ARIMA(4,1,3) memiliki jumlah parameter yang lebih banyak dibandingkan model pembangkitan. Hal ini memungkinkan model menangkap noise dalam data sebagai bagian dari pola, sehingga menghasilkan nilai AIC yang lebih kecil namun kurang representatif terhadap struktur sebenarnya.

Namun, selisih AIC dengan model ARIMA(1,1,1) sekitar 3.512346. Selisih dengan model dengan parameter yang lebih sederhana, ARIMA(1,1,1), tidak terlalu tinggi sehingga model ini masih layak dipilih tanpa kehilangan banyak informasi. Berdasarkan prinsip parsimony, model ARIMA(1,1,1) tetap dapat dianggap sebagai model yang lebih representatif karena lebih sederhana dan sesuai dengan proses pembangkitan data.

Tabel_Terurut$AIC[2]-Tabel_Terurut$AIC[1]
## [1] 3.512346

6. Diagnostik Residual

# Diagnostik Residual
Model_ARIMA <- arima(data, order = c(4,1,3), method = 'ML')
checkresiduals(Model_ARIMA)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,1,3)
## Q* = 3.483, df = 3, p-value = 0.323
## 
## Model df: 7.   Total lags used: 10
  • Uji Ljung-Box

\(\alpha: 0.05\)

\(H_0:\) Residuals ARIMA(4,1,3) tidak memiliki autokorelasi

\(H_1:\) Residuals ARIMA(4,1,3) memiliki autokorelasi

P-value yang didapatkan adalah 0.323. \(p-value > \alpha\) maka gagal tolak \(H_0\) sehingga dapat disimpulkan bahwa Residuals ARIMA(4,1,3) tidak memiliki autokorelasi.

  • Plot Residual

Dari plot di atas terlihat residual tersebar di sekitar rata-rata dan tidak membentuk pola (acak).Hal ini mengindikasikan bahwa model telah mampu menangkap pola data dengan baik.

  • Plot ACF

Dari plot ACF residual terlihat bahwa tidak ada spike yang signifikan. Hal ini menandakan bahwa tidak ada autokorelasi yang tersisa pada residual ARIMA(4,1,3).

  • Histogram + kurva normal

Dari visualisasi di atas terlihat bahwa distribusi residual cenderung simetris dan kurva yang dihasilkan berbentuk menyerupai lonceng. Hal ini menandakan bahwa secara visual, Residual ARIMA(4,1,3) mendekati bentuk berdistribusi normal. Namun, kita bisa uji normalitas sebagai berikut.

library(nortest)
lillie.test(Model_ARIMA$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Model_ARIMA$residuals
## D = 0.030571, p-value = 0.4815

\(\alpha: 0.05\)

\(H_0:\) Residuals ARIMA(4,1,3) berdistribusi normal

\(H_1:\) Residuals ARIMA(4,1,3) tidak berdistribusi normal

p-value yang didapatkan adalah 0.4815. \(P-value > \alpha\) maka gagal tolak \(H_0\) sehingga dapat disimpulkan bahwa residuals ARIMA(4,1,3) dianggap berdistribusi normal pada taraf signifikansi 5%.

Secara keseluruhan, residual model ARIMA(4,1,3) telah memenuhi asumsi white noise, yaitu tidak terdapat autokorelasi dan residual dapat dianggap berdistribusi normal.