Background

Time Series adalah sekumpulan data yang dikumpulkan berdasarkan urutan waktu. Data Time Series dapat kita gunakan untuk memahami pola yang terjadi di masa lalu dan melakukan suatu forecasting atau peramalan suatu nilai di masa depan.

Salah satu bentuk dari time series adalah Multivariate Time Series yaitu data time series yang memiliki lebih dari satu variabel dan setiap variabel tersebut tidak hanya bergantung pada nilai di masa lalunya saja namun juga bergantung dengan nilai variabel lainnya. Salah satu contohnya adalah temperatur dan tekanan udara. Nilai temperatur tidak hanya dipengaruhi oleh nilainya di masa lampau namun dapat dipengaruhi juga oleh nilai tekanan udara. Hal yang sama dapat terjadi pula untuk sebaliknya. Pada artikel ini akan dibahas bagaimana membangun model Multivariate Time Series menggunakan Vector Autoregressive (VAR).

Vector Autoregressive (VAR)

Model VAR adalah generalisasi dari model univariate autoregressive (model AR dengan satu variabel). Dalam model VAR kita akan membangun persamaan untuk setiap variabel yang akan membentuk suatu sistem persamaan. Kita notasikan \(y_{1,t}\) adalah nilai variabel \(y_{1}\) pada observasi ke-t dan \(y_{2,t}\) adalah nilai variabel \(y_{2}\) pada observasi ke-t, dan seterusnya. Sebelum membangun model VAR, kita perlu menentukan terlebih dahulu banyaknya variabel yang akan digunakan (dinotasikan sebagai \(K\)) dan banyaknya lag (dinotasikan sebagai \(p\)) yang akan terlibat dalam sistem VAR. Misalkan kita memiliki VAR yang terdiri dari dua buah variabel dengan satu lag. Maka kita peroleh \(VAR(1)\) dengan persamaan:

\[ y_{1,t} = c_{1} + \phi_{11,1}y_{1,t-1} + \phi_{12,1}y_{2,t-1} + e_{1,t}\\ y_{2,t} = c_{2} + \phi_{21,1}y_{1,t-1} + \phi_{22,1}y_{2,t-1} + e_{2,t} \]

dengan

\(c_{i}\) : konstanta untuk variabel \(i\)

\(\phi_{ij,l}\) : Pengaruh variabel \(j\) pada lag ke-\(l\) terhadap variabel \(i\)

\(e_{i,t}\) : Residual/error variabel \(i\)

Sistem persamaan di atas dapat kita tuliskan juga dalam notasi matriks berikut

\[ \begin{bmatrix} y_{1,t} \\ y_{2,t} \end{bmatrix} = \begin{bmatrix} c_{1} \\ c_{2} \end{bmatrix} + \begin{bmatrix} \phi_{11} & \phi_{12} \\ \phi_{21} & \phi_{22} \end{bmatrix} \begin{bmatrix} y_{1,t-1} \\ y_{2,t-1} \end{bmatrix} + \begin{bmatrix} e_{1,t} \\ e_{2,t} \end{bmatrix} \]

Import Library

Untuk mengimplementasikan VAR di R, kita akan menggunakan package vars. Selain itu, kita akan menggunakan beberapa package lainnya untuk melakukan preprocessing data, visualisasi data, serta evaluasi model.

library(dplyr)
library(padr)
library(lubridate)
library(forecast)
library(ggplot2)
library(MLmetrics)
library(lmtest)
library(vars)

Read Data

Kita akan menggunakan data dongsi.csv berikut yang merupakan data terkait kualitas udara di subdistrict Dongsi, Beijing. Deskripsi lebih lanjut dari data dapat dilihat pada link berikut

dongsi <- read.csv("dataset/PRSA_Data_20130301-20170228/PRSA_Data_Dongsi_20130301-20170228.csv")
head(dongsi)

Keterangan:

year : Tahun pengamatan
month : Bulan pengamatan
day : Tanggal pengamatan
hour : Jam pengamatan
PM2.5 : Konsentrasi PM2.5 (\(\mu g/m^3\))
PM10 : Konsentrasi PM10 (\(\mu g/m^3\))
SO2 : Konsentrasi \(SO_{2}\) (\(\mu g/m^3\))
NO2 : Konsentrasi \(NO_{2}\) (\(\mu g/m^3\))
CO : Konsentrasi \(CO\) (\(\mu g/m^3\))
O3 : Konsentrasi \(O_3\) (\(\mu g/m^3\))
TEMP : Temperatur (Celsius)
PRES : Tekanan (hPa)
DEWP : dew point temperature (Celsius)
RAIN : Curah Hujan (mm)
wd : Arah angin
WSPM : Kecepatan angin (m/s)
station : Stasiun pengamatan

Data Preprocessing

Tahapan preprocessing yang dilakukan adalah:

Pada dataset yang kita miliki, observasi diambil setiap jam (hourly). Karena kita akan melakukan pemodelan time series dalam rentang waktu harian (daily) maka kita hanya akan mengambil data pada jam 12:00 saja
Membuat kolom Date yang berisi tanggal pengambilan observasi
Cek NA values dan lakukan imputasi

# Step 1
dongsi <- dongsi[dongsi$hour == 12,]

# Step 2
dongsi$Date <- paste(dongsi$year, "-", dongsi$month, "-", dongsi$day) %>% ymd()

head(dongsi, 10)

# Step 3
dongsi %>% is.na() %>% colSums()

##      No    year   month     day    hour   PM2.5    PM10     SO2     NO2      CO 
##       0       0       0       0       0      53      36      32      69     143 
##      O3    TEMP    PRES    DEWP    RAIN      wd    WSPM station    Date 
##      31       0       0       0       0       0       0       0       0

Terdapat beberapa kolom yang memiliki nilai missing value. Maka kita lakukan imputasi dengan metode forward filling (mengisi NA values berdasarkan nilai terakhir sebelum NA)

dongsi <- dongsi %>%
  pad(start_val = min(dongsi$Date), end_val = max(dongsi$Date)) %>% 
  mutate_if(is.numeric, na.locf)

## pad applied on the interval: day

colSums(is.na(dongsi))

##      No    year   month     day    hour   PM2.5    PM10     SO2     NO2      CO 
##       0       0       0       0       0       0       0       0       0       0 
##      O3    TEMP    PRES    DEWP    RAIN      wd    WSPM station    Date 
##       0       0       0       0       0       0       0       0       0

Exploratory Data Analysis

Pertama, kita lihat terlebih dahulu visualisasi pada masing-masing kolom.

ts(dongsi[,c(6:15)], start = c(2013,60), freq = 365) %>% 
  autoplot(facets = TRUE)

Perhatikan bahwa pada kolom TEMP dan PRES terjadi hubungan yang cukup kuat. Ketika nilai TEMP meningkat maka nilai PRES cenderung menurun dan begitu pula untuk sebaliknya. Untuk mengetahui hubungan apakah nilai dari suatu variabel di masa lampau dapat memengaruhi variabel lainnya di masa depan, kita dapat menggunakan uji Granger Causality Test. Berikut adalah hipotesis untuk uji Granger Causality:

\(H_0\) : X tidak dapat digunakan untuk memprediksi Y

\(H_1\) : X dapat digunakan untuk memprediksi Y

Pada kasus multivariate time series, kita perlu melakukan pengujian secara dua arah yaitu apakah X dapat digunakan untuk memprediksi Y dan apakah Y dapat digunakan untuk memprediksi X. Untuk melakukan uji tersebut, kita dapat menggunakan fungsi grangertest dari package lmtest

grangertest(x = dongsi$TEMP, y = dongsi$PRES)

Diperoleh bahwa nilai p-values < 0.05 maka dapat disimpulkan bahwa variabel TEMP di masa lampau dapat digunakan untuk memprediksi nilai PRES di masa depan

Kita juga perlu menguji apakah hal sebaliknya dapat berlaku, yaitu apakah variabel PRES dapat digunakan untuk memprediksi nilai TEMP di masa depan

grangertest(x = dongsi$PRES, y = dongsi$TEMP)

Diperoleh bahwa nilai p-values < 0.05 maka dapat disimpulkan variabel PRES dapat digunakan untuk memprediksi nilai TEMP di masa depan

Berdasarkan hasil di atas, maka kita dapat membangun model multivariate time series yang terdiri dari dua buah variabel yaitu TEMP dan PRES.

Modeling

Pertama, kita buat kolom TEMP dan PRES menjadi object time series.

data_ts <- ts(dongsi[,c("TEMP", "PRES")], start = c(2013,60), freq=365)
data_ts %>% autoplot(facets = TRUE)

Lakukan spliting data train dan data test

n_train <- as.integer(0.8*nrow(data_ts))
train_data <- head(data_ts, n_train)
test_data <- tail(data_ts, -n_train)

Pada package vars, terdapat fungsi VARselect yang dapat digunakan untuk mencari lag yang menghasilkan model dengan performa yang paling baik. Parameter lag.max menunjukkan nilai maksimum lag yang akan dicari dan season menunjukkan seasonality pada data. Fungsi tersebut akan mencoba membuat model VAR dari lag 1,2, sampai lag.max. Pada masing-masing lag tersebut akan dihitung seberapa baik model tersebut berdasarkan kriteria AIC, HQ, SC, dan FPE. Nilai-nilai tersebut dapat diinterpretasikan sebagai banyaknya informasi yang hilang dari model yang telah dibuat. Semakin kecil nilai tersebut maka model semakin baik. Fungsi VARselect akan memberikan output lag mana yang memberikan nilai-nilai paling kecil untuk masing-masing kriteria tersebut.

VARselect(train_data, lag.max = 10, season = 365)

## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      5      3      2      5 
## 
## $criteria
##                 1          2          3          4          5          6
## AIC(n)   5.799175   5.773589   5.759686   5.756928   5.751587   5.756739
## HQ(n)    7.008131   6.989133   6.981819   6.985649   6.986896   6.998636
## SC(n)    9.002945   8.994819   8.998375   9.013076   9.025194   9.047805
## FPE(n) 345.286776 336.825767 332.436318 331.784221 330.282564 332.258893
##                 7          8          9         10
## AIC(n)   5.758723   5.763874   5.765059   5.755199
## HQ(n)    7.007209   7.018948   7.026721   7.023449
## SC(n)    9.067249   9.089859   9.108503   9.116102
## FPE(n) 333.193405 335.193426 335.873903 332.862435

Dengan mempertimbangkan nilai AIC pada output di atas, diperoleh bahwa model VAR(5) merupakan model yang terbaik dibandingkan dengan lag lainnya. Maka langkah selanjutnya adalah membangun model tersebut dengan fungsi VAR

model <- VAR(train_data, p=5, season = 365)

Evaluation

Setelah memperoleh model, maka kita akan menggunakan model tersebut untuk melakukan evaluasi menggunakan test_data. Pertama akan kita lakukan forecasting sebanyak baris pada data tersebut.

prediction <- forecast(model, h = nrow(test_data))

Kemudian, mari kita visualisasikan hasil forecasting dan membandingkannya dengan nilai pada test_data

# plot TEMP
train_data[,"TEMP"] %>% 
  autoplot(main = "TEMPERATURE") +
  autolayer(test_data[,"TEMP"], series = "Test Data") +
  autolayer(prediction$forecast$TEMP$mean, series = "Forecast")

# Plot PRES
train_data[,"PRES"] %>% autoplot(main = "PRESSURE") +
  autolayer(test_data[,"PRES"], series = "Test Data") +
  autolayer(prediction$forecast$PRES$mean, series = "Forecast")

Dari visualisasi di atas, dapat kita lihat bahwa hasil forecast sudah mengikuti pola dari nilai aktualnya. Untuk memeriksa lebih lanjut, mari kita lihat nilai MAE dari kedua variabel tersebut untuk mengetahui seberapa jauh rata-rata penyimpangan prediksi kita.

mae_temp <- MAE(prediction$forecast$TEMP$mean, test_data[,"TEMP"])
mae_pres <- MAE(prediction$forecast$PRES$mean, test_data[,"PRES"])
paste("MAE untuk temperature: ", mae_temp)

## [1] "MAE untuk temperature:  2.7086742102155"

paste("MAE untuk pressure: ", mae_pres)

## [1] "MAE untuk pressure:  4.74033499797873"

Dari output di atas, diperoleh bahwa nilai MAE untuk temperatur dan pressure adalah 2.70 dan 4.74. Hal ini menunjukkan bahwa secara rata-rata hasil prediksi kita akan menyimpang terhadap nilai aktual sebesar 2.70 derajat celcius untuk temperatur dan 4.74 hPa untuk pressure. Dengan nilai tersebut maka dapat disimpulkan bahwa model yang telah dibuat sudah cukup baik untuk melakukan forecasting data.

Forecasting

Langkah selanjutnya adalah kita akan coba lakukan forecasting pada waktu satu tahun ke depan dari data yang kita punya. Untuk melakukan hal tersebut, kita akan coba lakukan pemodelan ulang dengan lag dan parameter yang sama seperti sebelumnya namun data yang kita gunakan adalah keseluruhan data (data sebelum dilakukan spliting), yaitu data_ts.

# pemodelan
final_model <- VAR(data_ts, p = 5, season = 365)

# forecasting
forecasting <- forecast(final_model, h = 365)

Kemudian, kita coba lakukan visualisasi hasil forecasting

# Plot temperature
data_ts[,"TEMP"] %>% 
  autoplot(main = "TEMPERATURE") +
  autolayer(forecasting$forecast$TEMP$mean, series = "Forecast")

# Plot Pressure
data_ts[,"PRES"] %>% autoplot(main = "PRESSURE") +
  autolayer(forecasting$forecast$PRES$mean, series = "Forecast")

Kesimpulan

Dari artikel ini, kita telah mempelajari terkait multivariate time series yaitu time series yang terdiri dari beberapa variabel/prediktor. Salah satu model yang dapat digunakan untuk memodelkan tipe time series ini adalah Vector Autoregressive (VAR). Model ini dapat kita manfaatkan untuk menggambarkan pengaruh dari satu variabel ke variabel lainnya. Tentunya, model ini dapat digunakan di berbagai bidang seperti makroekonomi, meteorologi, dan kesehatan.

Pada model VAR, kita perlu menentukan variabel apa yang dapat dimasukkan ke dalam model dan seberapa jauh lag yang akan diambil.

Untuk menentukan variabel yang dapat digunakan, kita dapat melakukan uji Granger Causality Test. Uji tersebut digunakan untuk melihat apakah suatu variabel dapat digunakan untuk memprediksi variabel lainnya. Uji ini perlu dilakukan secara dua arah.
Untuk menentukan lag yang akan diambil, kita dapat memanfaatkan fungsi VARselect() yang akan membuat beberapa model dengan lag yang berbeda dan akan memberikan output lag yang memberikan performa terbaik berdasarkan kriteria AIC, HQ, SC, atau FPE.

Introduction to Vector Autoregressive (VAR)

Tria Rahmat Mauludin

5/25/2022