“Praktikum 04 - Pembangkitan Arma”

Angga Fathan Rofiqy

21 September, 2023

1 Rpubs

1.1 White Noise

Pembangkitan data berpola AR, MA, ARMA, dan banyak proses deret waktu lainnya diawali pembangkitan white noise. White noise merupakan sederet nilai dari peubah bebas stokastik identik. Oleh karena itu, white noise memiliki dua karakteristik penting:

White noise tidak memiliki autokorelasi (karena saling bebas)
Nilai harapan dan ragam white noise sama (karena berasal dari peubah acak bebas stokastik identik)

White noise dibangkitkan dari suatu peubah acak, umumnya peubah acak normal.

wn <- rnorm(300)
ts.plot(wn)

Dapat terlihat bahwa white noise tidak memiliki autokorelasi dari ACF. Perhatikan bahwa lag ke-0 adalah korelasi observasi ke-t dengan dirinya sendiri. Nilai korelasi tersebut pasti 1. Sebagai alternatif, lag pertama di plot ACF dapat ditetapkan sebagai 1 (alih-alih 0) dengan menambahkan argumen xlim(1, lag akhir). Plot tersebut dapat disandingkan bersamaan dengan membuat matriks \(1 \times 2\) dengan par(mfrow = c(1,2)).

par(mfrow = c(1, 2)) 
acf(wn)
acf(wn, xlim = c(1, 20))

2 Proses MA

Proses MA dapat dituliskan sebagai berikut:

\[ y_{t} = c + e_t + \theta_{1}e_{t-1} + \theta_{2}e_{t-2} + \dots + \theta_{q}e_{t-q} = c+{e_t+\sum_{i=1}^p \theta_ie_{t-i}} \] Terlihat bahwa \(e_t\), atau white noise, berperan penting dalam pembangkitan proses MA.

2.1 Pembangkitan Proses MA(2)

Akan dicoba membangkitkan proses MA paling sederhana, yaitu MA(1) dengan \(\theta = 0.5\) sebanyak 200 observasi dan \(c=0\). Karena diperlukan satu nilai awal untuk \(e_{t-1}\), masukkan nilai pertama white noise sebagai nilai awal tersebut.

set.seed(006)
ma <- wn[c(1,2)]

Nilai-nilai selanjutnya dapat dicari melalui loop. Bentuk loop dapat dilihat dari rumus MA(1) yang hendak dibangkitkan:

\[ y_t = e_t + 0.4e_{t-1} + 0.6e_{t-2} \]

for(i in 3:300){
   ma[i] <- wn[i] + 0.4 * wn[i - 1] + 0.6 * wn[i - 2]  
}

install_load('DT')
datatable(as.matrix(ma), filter = 'top', 
          options = list(pageLength = 5))

Selain menggunakan cara di atas, pembangkitan proses MA(2) dapat dilakukan dengan fungsi arima.sim() sebagai berikut.

ma2 <- arima.sim(list(order=c(0,0,2), ma=c(0.4, 0.6)), n=300)

datatable(as.matrix(ma2), filter = 'top', 
          options = list(pageLength = 5))

2.2 Karakteristik MA(1)

2.2.1 Plot Time Series

ts.plot(ma)

Berdasarkan plot time series, terlihat bahwa data MA(2) yang dibangkitkan stasioner dalam rataan.

2.2.2 Plot ACF

acf(ma,lag.max = 20)

Berdasarkan plot AFC tersebut, terlihat bahwa plot ACF cuts off di lag kedua

2.2.3 Plot PACF

pacf(ma)

Berdasarkan plot PACF tersebut, terlihat bahwa plot PACF cenderung tails off dan membentuk gelombang sinus

2.2.4 Plot EACF

TSA::eacf(ma)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o  o  x  o 
## 1 x x x o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  x  o 
## 3 o x o x o o o o o o o  o  o  o 
## 4 o x x x o o o o o o o  o  o  o 
## 5 x x x o o o o o o o o  o  o  o 
## 6 x x x x o o o o o o o  o  o  o 
## 7 x x x x o o o o o o o  o  o  o

Berdasarkan pola segitiga nol pada plot EACF, terlihat bahwa segitiga nol berada pada ordo AR(0) dan ordo MA(2)

Dari keempat plot tersebut, terlihat bahwa data stasioner.

2.3 Scatterplot Antar Lag

2.3.1 Korelasi antara \(Y_t\) dengan \(Y_{t-1}\)

#Yt
yt_ma <- ma[-1]
datatable(as.matrix(yt_ma), filter = 'top', 
          options = list(pageLength = 5))

#Yt-1
yt_1_ma <- ma[-300]

datatable(as.matrix(yt_1_ma), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_ma,x=yt_1_ma)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-1}\). Hal ini sesuai dengan teori yang ada

cat("Nilai Korelasi : ", cor(yt_ma,yt_1_ma))

## Nilai Korelasi :  0.3213018

Korelasi antara \(Y_t\) dengan \(Y_{t-1}\) dari hasil simulasi mendekati perhitungan teoritis.

2.3.2 Korelasi antara \(Y_t\) dengan \(Y_{t-2}\)

#Yt
yt_ma2 <- ma[-c(1,2)]
datatable(as.matrix(yt_ma2), filter = 'top', 
          options = list(pageLength = 5))

#Yt-2
yt_2_ma <- ma[-c(299,300)]

datatable(as.matrix(yt_2_ma), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_ma2, x=yt_2_ma)

Berdasarkan scatterplot tersebut, terlihat bahwa cenderung tidak terdapat hubungan antara \(Y_t\) dengan \(Y_{t-2}\).

cat("Nilai Korelasi : ", cor(yt_ma2,yt_2_ma))

## Nilai Korelasi :  0.364761

Korelasi antara \(Y_t\) dengan \(Y_{t-2}\) hasil simulasi yaitu \(0.4017548\).

2.3.3 Korelasi antara \(Y_t\) dengan \(Y_{t-3}\)

#Yt
yt_ma3 <- ma[-c(1,2,3)]
datatable(as.matrix(yt_ma3), filter = 'top', 
          options = list(pageLength = 5))

#Yt-2
yt_3_ma <- ma[-c(298,299,300)]

datatable(as.matrix(yt_3_ma), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_ma3, x=yt_3_ma)

Berdasarkan scatterplot tersebut, terlihat bahwa cenderung tidak terdapat hubungan antara \(Y_t\) dengan \(Y_{t-3}\).

cat("Nilai Korelasi : ", cor(yt_ma3, yt_3_ma))

## Nilai Korelasi :  -0.02723429

Korelasi antara \(Y_t\) dengan \(Y_{t-3}\) hasil simulasi mendekati teori yang ada yaitu \(0\).

3 Proses AR

Proses AR dapat dituliskan sebagai berikut:

\[ y_{t} = c + e_t + \phi_{1}Y_{t-1} + \phi_{2}Y_{t-2} + \dots + \phi_{q}Y_{t-q} = c+{e_t+\sum_{i=1}^p \phi_iY_{t-i}} \] Terlihat bahwa \(Y_t\) berperan penting dalam pembangkitan proses AR.

3.1 Pembangkitan Proses AR

Akan dicoba membangkitkan proses AR paling sederhana, yaitu AR(1) dengan \(\phi = 0.7\) sebanyak 200 observasi dan \(c=0\).

set.seed(006)

Nilai-nilai selanjutnya dapat dicari melalui loop. Bentuk loop dapat dilihat dari rumus AR(1) yang hendak dibangkitkan:

\[ Y_t = e_t+ 0.5Y_{t-1} + 0.2Y_{t-2} \]

n<-length(wn)
cat("n : ", n)

## n :  300

ar <- c(1:n) 
for (i in 3:n) {ar[i]<-wn[i]+0.5*ar[i-1]+0.2*ar[i-2]}
datatable(as.matrix(ar), filter = 'top', 
          options = list(pageLength = 5))

Selain menggunakan cara di atas, pembangkitan proses AR dapat dilakukan dengan fungsi arima.sim() sebagai berikut.

ar2 <- arima.sim(list(order=c(2,0,0), ar=c(0.5,0.2)), n=300)
datatable(as.matrix(ar2), filter = 'top', 
          options = list(pageLength = 5))

3.2 Karakteristik AR(2)

3.2.1 Plot Time Series

ts.plot(ar)

Berdasarkan plot time series tersebut terlihat bahwa data cenderung stasioner pada rataan

3.2.2 Plot ACF

acf(ar)

Berdasarkan plot ACF tersebut terlihat bahwa plot ACF cenderung tails off dan cenderung membentuk pola grafik sinus

3.2.3 Plot PACF

pacf(ar)

Berdasarkan plot PACF tersebut, terlihat bahwa plot PACF cuts off pada lag kedua, sejalan dengan teori yang ada

3.2.4 Plot EACF

TSA::eacf(ar)

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

Berdasarkan pola segitiga nol pada plot EACF, terlihat bahwa segitiga nol berada pada ordo AR(2) dan ordo MA(0)

3.3 Scatterplot Antar Lag

3.3.1 Korelasi antara \(Y_t\) dengan \(Y_{t-1}\)

#Yt
yt_ar <- ar[-1]

datatable(as.matrix(yt_ar), filter = 'top', 
          options = list(pageLength = 5))

#Yt-1
yt_1_ar <- ar[-300]
datatable(as.matrix(yt_1_ar), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_ar,x=yt_1_ar)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-1}\). Hal ini sesuai dengan teori yang ada

cat("Nilai Korelasi : ", cor(yt_ar,yt_1_ar))

## Nilai Korelasi :  0.5571268

Korelasi antara \(Y_t\) dengan \(Y_{t-1}\) dari hasil simulasi mendekati perhitungan teoritis yaitu \(\rho_1=\phi^1=0.64\)

3.3.2 Korelasi antara \(Y_t\) dengan \(Y_{t-2}\)

#Yt
yt_ar2 <- ar[-c(1,2)]

datatable(as.matrix(yt_ar2), filter = 'top', 
          options = list(pageLength = 5))

#Yt-2
yt_2_ar <- ar[-c(299,300)]

datatable(as.matrix(yt_2_ar), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_ar2,x=yt_2_ar)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-2}\). Hal ini sesuai dengan teori yang ada

cat("Nilai Korelasi : ", cor(yt_ar2,yt_2_ar))

## Nilai Korelasi :  0.440635

Korelasi antara \(Y_t\) dengan \(Y_{t-2}\) dari hasil simulasi mendekati perhitungan teoritis yaitu \(\rho_2=\phi^2=0.48\).

3.3.3 Korelasi antara \(Y_t\) dengan \(Y_{t-3}\)

#Yt
yt_ar3 <- ar[-c(1,2,3)]
datatable(as.matrix(yt_ar3), filter = 'top', 
          options = list(pageLength = 5))

#Yt-3
yt_3_ar <- ar[-c(298,299,300)]
datatable(as.matrix(yt_3_ar), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_ar3,x=yt_3_ar)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-3}\). Hal ini sesuai dengan teori yang ada

cat("Nilai Korelasi : ", cor(yt_ar3,yt_3_ar))

## Nilai Korelasi :  0.3513907

Korelasi antara \(Y_t\) dengan \(Y_{t-3}\) dari hasil simulasi mendekati perhitungan teoritis yaitu \(\rho_2=\phi^2=0.35\).

4 Fungsi pembangkitan ARMA

Setelah mengetahui cara membangkitkan data berpola AR, MA, dan ARMA sederhana, bagaimana cara melakukan pembangkitan data berpola tersebut yang lebih kompleks? Apakah dapat dibuat suatu fungsi yang fleksibel yang memungkinan pembangkitan dengan berapapun jumlah koefisien?

Pertama, lihat kembali bentuk umum data berpola ARMA.

\[ y_{t} = c + \sum_{i=1}^p \phi_{i}y_{t-i} + \sum_{j=1}^q e_{t-j}+ e_{t} \]

Dari prinsip ini, dapat dibuat fungsi umum untuk membangkitkan data ARMA. Input dari fungsi adalah jumlah data yang hendak dibangkitkan, koefisien MA, dan koefisien AR.

arma.sim <- function(n, macoef, arcoef){
  manum <- length(macoef)
  arnum <- length(arcoef)
  stopifnot(manum < n & arnum < n)
  
  wn <- rnorm(n, sd = 0.5)
  init <- max(manum, arnum)

  arma <- wn[1:init]
  for(i in {init+1}:n){
   mastart <- i - manum
   maend <- i-1
   arstart <- i - arnum
   arend <- i-1
   arma[i] <- sum(arcoef * arma[arstart:arend]) + sum(macoef * wn[mastart:maend])  + wn[i]
   }
  return(arma)
}

4.0.0.1 Pembangkitan Proses ARMA(2,2) manual

set.seed(006)
n = length(wn)
phi1 = 0.5
phi2 = 0.2

theta1 = 0.4
theta2 = 0.6

y.arma=c(1:n)
for (i in 5:n){y.arma[i] = phi1*y.arma[i-1] +phi2*y.arma[i-2]+ theta1*wn[i-1] + theta2*wn[i-2] +wn[i]}
datatable(as.matrix(y.arma), filter = 'top', 
          options = list(pageLength = 5))

4.0.0.2 Pembangkitan Proses ARMA(2,2) manual dengan fungsi arima.sim

Pembangkitan ARMA(p,q) juga dapat dilakukan dengan fungsi arima.sim sebagai berikut.

arma22 <- arima.sim(list(order=c(2,0,2), ar = c(0.5,0.2), ma =c(0.4,0.6)), n=300)
datatable(as.matrix(arma22), filter = 'top', 
          options = list(pageLength = 5))

4.1 Karakteristik ARMA(2,2)

4.1.1 Plot Time Series

par(mfrow = c(1, 2))
ts.plot(y.arma)
ts.plot(arma22)

par(mfrow = c(1, 1))

Berdasarkan plot time series tersebut, terlihat bahwa model ARMA(2,2) cenderung stasioner dalam rataan

4.1.2 Plot ACF

par(mfrow = c(1, 2))
acf(y.arma)
acf(arma22)

par(mfrow = c(1, 1))

Berdasarkan plot ACF tersebut, terlihat bahwa model ARMA(2,2) hasil simulasi memiliki plot ACF yang tails off, sesuai dengan teori yang ada.

4.1.3 Plot PACF

par(mfrow = c(1, 2))
pacf(y.arma)
pacf(arma22)

par(mfrow = c(1, 1))

Berdasarkan plot PACF tersebut, terlihat bahwa model ARMA(2,2) hasil simulasi memiliki plot PACF yang tails off, sesuai dengan teori.

4.1.4 Plot EACF

TSA::eacf(y.arma)

## 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 x o x  x  x  x 
## 1 x x x o o o o o o o o  o  o  o 
## 2 x x o x o o o o o o o  o  o  o 
## 3 x x x x o o o o o o o  o  x  o 
## 4 x x o x o o o o o o o  o  o  o 
## 5 x x o x o o o o o o o  o  o  o 
## 6 x x x x o o o o o o o  o  o  o 
## 7 x x x x o o o o o o o  o  o  o

TSA::eacf(arma22)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x o o o o o o  o  o  o 
## 1 o x x o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x x o o o o x o o o o  o  o  o 
## 4 x x o o o o o o o o o  o  o  o 
## 5 x x x x x x o o o o o  o  o  o 
## 6 x x o x o x o o o o o  o  o  o 
## 7 x x x o o o x x o o o  o  o  o

Berdasarkan pola segitiga nol pada plot EACF, terlihat bahwa segitiga nol berada pada ordo AR(2) dan ordo MA(2).

4.2 Scatterplot Antar Lag

4.2.1 Korelasi antara \(Y_t\) dengan \(Y_{t-1}\)

#Yt
yt_arma <- arma22[-1]
datatable(as.matrix(yt_arma), filter = 'top', 
          options = list(pageLength = 5))

#Yt-1
yt_1_arma <- arma22[-300]
datatable(as.matrix(yt_1_arma), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_arma,x=yt_1_arma)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-1}\). Hal ini sesuai dengan teori yang ada.

cat("Nilai Korelasi : ", cor(yt_arma,yt_1_arma))

## Nilai Korelasi :  0.8097838

4.2.2 Korelasi antara \(Y_t\) dengan \(Y_{t-2}\)

#Yt
yt_arma <- arma22[-c(1,2)]
datatable(as.matrix(yt_arma), filter = 'top', 
          options = list(pageLength = 5))

#Yt-1
yt_2_arma <- arma22[-c(299,300)]
datatable(as.matrix(yt_2_arma), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_arma,x=yt_2_arma)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-2}\). Hal ini sesuai dengan teori yang ada.

cat("Nilai Korelasi : ", cor(yt_arma,yt_2_arma))

## Nilai Korelasi :  0.678883

4.2.3 Korelasi antara \(Y_t\) dengan \(Y_{t-3}\)

#Yt
yt_arma <- arma22[-c(1,2,3)]
datatable(as.matrix(yt_arma), filter = 'top', 
          options = list(pageLength = 5))

#Yt-1
yt_3_arma <- arma22[-c(298,299,300)]
datatable(as.matrix(yt_3_arma), filter = 'top', 
          options = list(pageLength = 5))

plot(y=yt_arma,x=yt_3_arma)

Berdasarkan scatterplot tersebut, terlihat bahwa terdapat hubungan positif antara \(Y_t\) dengan \(Y_{t-3}\). Hal ini sesuai dengan teori yang ada.

cat("Nilai Korelasi : ", cor(yt_arma,yt_3_arma))

## Nilai Korelasi :  0.4147366

5 Referensi

Rpubs : https://rpubs.com/Prass_Code/911845
Rmd : Read This Nerd!
DT : DT
Theme : Rmd Themes
Hide Code : Hide Code In RMD
Kable : Kable1, Kable2
Other Rmd settings : Set1, Set2, Set3
Notes : Hanya Orang Gabut yang baca ini
https://packer.john-coene.com/#/guide/style