Membuat Model Time Series ARIMA Menggunakan Python

ARIMA atau Auto Regressive Integrated Moving Average adalah model yang banyak dipakai dalam peramalan data time series univariat. Sesuai namanya, model ARIMA terdiri 3 komponen yaitu Auto Regressive (AR), Integrated (I) dan Moving Average (MA) dan dinotasikan sebagai ARIMA(p, d, q).

Komponen p menunjukkan jumlah orde AR pada model. Model AR dapat dikatakan sebagai model regresi dengan peubah penjelasnya adalah data periode sebelumnya. Contoh, nilai p=2 menunjukkan untuk memprediksi data pada waktu ke-t, model menggunakan peubah penjelas yaitu data pada waktu ke-(t-1) dan (t-2).

Komponen d menunjukkan jumlah proses diferensiasi (d). Proses ini perlu dilakukan untuk menangani kondisi non-stasioner (misalnya trend) hingga mencapai stasioner. Nilai d=0 menunjukkan kondisi daata yang sudah stasioner tanpa perlu melakukan proses diferensiasi atau disebut juga model ARMA(p, q).

Komponen terakhir q menunjukkan orde atau lag dari error pada model MA. Nilai q menunjukkan jumlah periode yang digunakan untuk memprediksi error saat ini berdasarkan error periode sebelumnya.

Pemilihan kombinasi yang tepat untuk p, d dan q merupakan bagian penting untuk memperoleh model ARIMA yang optimal. Nilai p dan q umumnya dapat ditentukan dengan melihat plot Autocorrelation function (ACF) dan Partial Autocorrelation Function (PACF). Sementara itu untuk d terdapat uji yang dapat digunakan yaitu uji Augmented Dickey–Fuller (ADF).

Pendekatan lain yang banyak digunakan saat ini khususnya dalam konteks machine learning adalah grid search. Grid search atau dalam penerapan pada beberapa software menggunakan terminologi auto ARIMA adalah proses pencarian model secara iteratif dari berbagai kombinasi nilai p, d dan q. Nilai p, d dan q terbaik ditentukan berdasarkan metrik tertentu seperti AIC atau BIC.

Pada tutorial ini kita akan membuat model ARIMA menggunakan bahasa pemrograman python dengan library statsmodel dan memfokuskan pada penggunaan grid search dalam menentukan orde p, d dan q yang optimal.

Contoh yang akan digunakan adalah data ekspor Indonesia periode bulanan selama tahun 2016-2019 (48 bulan). download

Data terdiri dari dua kolom yaitu periode dan nilai_jt_usd. Kita perlu merubah kolom periode sebagai indeks bagi dataframe. Langkah ini dilakukan dengan mengisi nilai parameter parse_dates. Selain itu, kita juga perlu menentukan frekuensi data yang digunakan yaitu bulanan.

Untuk eksplorasi data, kita dapat memvisualisasikan secara langsung menggunakan metode plot pada objek dataframe ekspor seperti kode berikut:

Python

import pandas as pd

# Memuat dataset ekspor-idn.csv
ekspor = pd.read_csv("ekspor-idn.csv", index_col="periode", parse_dates=True)

# mengatur frekuensi data bulanan
ekspor.index.freq = 'MS'

# print data
print(ekspor)

# membuat plot data
ekspor.plot(figsize=(12, 6))

Output

# Output

            nilai_jt_usd
periode                 
2016-01-01       39479.3
2016-02-01       38596.9
2016-03-01       42908.5
...
2019-10-01       61413.0
2019-11-01       56022.2
2019-12-01       55720.1

Pembuatan model ARIMA dilakukan dengan fungsi ARIMA dari library statsmodel. Pada bagian ini, kita akan membuat model ARIMA dengan nilai p,d,q yang sudah ditentukan yaitu p=1, d=1 dan q=1. Selanjutnya, dari model tersebut kita dapat melakukan forecasting atau prediksi nilai ekspor Indonesia untuk periode-periode berikutnya. Misalkan seperti contoh berikut ini, kita melakukan forecasting untuk 3 bulan yaitu Januari-Maret 2020.

# Memuat fungsi ARIMA
from statsmodels.tsa.arima.model import ARIMA

# Membagun model ARIMA dengan p=1, d=1 dan q=1
model = ARIMA(data, order=(1, 1, 1))
model_fit = model.fit()

# Melakukan forecasting
# 3 bulan (Jan-Mar 2020)
print(model_fit.forecast(3))

# alternatif prediksi
# model_fit.predict(start=len(data), end=len(data)+2, typ="levels")

# menampilkan summary model
print(model_fit.summary())

# forecast

2020-01-01    55472.142703
2020-02-01    55586.231825
2020-03-01    55533.737594
Freq: MS, Name: predicted_mean, dtype: float64


# summary

SARIMAX Results                                
==============================================================================
Dep. Variable:           nilai_jt_usd   No. Observations:                   48
Model:                 ARIMA(1, 1, 1)   Log Likelihood                -461.005
Date:                Sat, 21 Jan 2023   AIC                            928.009
Time:                        14:38:14   BIC                            933.560
Sample:                    01-01-2016   HQIC                           930.098
                         - 12-01-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4601      0.245     -1.878      0.060      -0.940       0.020
ma.L1          0.1835      0.272      0.675      0.500      -0.350       0.717
sigma2      1.218e+07   7.16e-09    1.7e+15      0.000    1.22e+07    1.22e+07
==============================================================================
Ljung-Box (L1) (Q):              9.50   Jarque-Bera (JB):                 0.38
Prob(Q):                         0.00   Prob(JB):                         0.83
Heteroskedasticity (H):          1.27   Skew:                             0.14
Prob(H) (two-sided):             0.63   Kurtosis:                         2.65
==============================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 3.29e+31. 
    Standard errors may be unstable.

Seperti yang telah disampaikan, model ARIMA memerlukan penentuan nilai p,d dan q. Salah satu teknik yang dapat digunakan adalah grid search. Langkah pertama yang perlu dilakukan adalah menyiapkan nilai-nilai p, d dan q yang akan dicobakan. Untuk memudahkan pekerjaan ini, kita dapat menggunakan fungsi product pada library itertools. Contoh yang akan kita gunakan yaitu nilai p (0-4), d (0-2) dan q (0-4). Secara total, pengaturan ini akan membandingkan 5x3x5 atau 75 model dengan kombinasi nilai p, d dan q yang berbeda.

from itertools import product

# Membuat list untuk p, d, dan q
# ingat : indeks pada python dimulai dari 0
p = range(0, 5)  # 0-4
d = range(0, 3)  # 0-2
q = range(0, 5)  # 0-4

# Menggunakan fungsi product dari itertools 
# untuk membuat kombinasi p, d, dan q
pdq = list(product(p, d, q))

print(pdq)

Output:

[(0, 0, 0),
 (0, 0, 1),
 (0, 0, 2),
 ...
 (4, 2, 2),
 (4, 2, 3),
 (4, 2, 4)]

Berikutnya, kita lakukan iterasi untuk setiap kombinasi nilai pada variabel pdq. Pada setiap iterasi kita membuat model ARIMA dengan nilai p, d dan q tersebut. Dari setiap model yang dibuat, kita juga menyimpan metrik aic yang nantinya digunakan dalam penentuan model terbaik

from statsmodels.tsa.arima.model import ARIMA

# Membuat list untuk menyimpan Skor AIC
aic_scores = []

# Melakukan grid search manual untuk menemukan p, d, q yang optimal
for param in pdq:
    # Melakukan fitting ARIMA model
    model = ARIMA(data, order=param)
    model_fit = model.fit()
    # Menambahkan aic score ke list
    aic_scores.append({'par': param, 'aic': model_fit.aic})
    
# Memilih aic score terkecil
best_aic = min(aic_scores, key=lambda x: x['aic'])

print(best_aic)

Output:

{'par': (2, 1, 2), 'aic': 911.7672265060809}

Dari proses grid search diperoleh nilai p, d dan q yang optimal adalah p=2, d=1 dan q=2 dengan skor AIC terkecil yaitu 911,767.

Langkah berikutnya, buat model ARIMA menggunakan nilai p, d dan q tersebut. Model kemudian dapat digunakan untuk forecasting sama seperti sebelumnya

# Membuat model dengan p, d, dan q hasil grid search
model = ARIMA(data, order=best_aic['par']) # order=(2, 1, 2)
model_fit = model.fit()

# forecasting 3 periode waktu ke depan
model_fit.forecast(3)

# alternatif prediksi
# model_fit.predict(start=len(data), end=len(data)+2, typ="levels")

2020-01-01    51357.694104
2020-02-01    54836.830177
2020-03-01    52991.158003
Freq: MS, Name: predicted_mean, dtype: float64

Pada bagian sebelumnya, kita membangun model ARIMA menggunakan SELURUH DATA yang ada sehingga model tidak dapat dievaluasi dengan data baru. Cara tersebut terkadang dapat menghasilkan model yang overfitting. Praktik yang disarankan adalah membagi data menjadi setidaknya 2 bagian yaitu data training dan testing. Model akan dilatih dan dibentuk hanya berdasarkan data training saja, kemudian dilakukan evaluasi performa model menggunakan data testing (out-of-sample forecasting).

Tidak seperti validasi pada data cross section, pada data time series pembagian data tidak dapat dilakukan secara acak. Data harus dibagi 2, tepat pada suatu titik waktu tertentu. Setiap data yang berada pada periode sebelum titik tersebut digunakan sebagai data training, sementara data-data setelahnya menjadi menjadi data testing.

Dalam contoh ini, kita akan menggunakan data 5 bulan terakhir (Agustus-Desember 2019) sebagai data testing sehingga data sebelum periode tersebut menjadi data training.

Untuk menentukan nilai p,d dan q kita akan menggunakan cara yang sama seperti sebelumnya, namun data untuk fitting model adalah data training saja, bukan seluruh data.

from statsmodels.tsa.arima.model import ARIMA

# Membagi data training dan testing
data_train = data[:-5]
data_test = data[-5:]

# Membuat list untuk menyimpan Skor AIC
aic_scores = []

# Melakukan grid search manual untuk menemukan p, d, q yang optimal
for param in pdq:
    # Melakukan fitting model ARIMA
    model = ARIMA(data_train, order=param)
    model_fit = model.fit()
    # Menambahkan aic score ke list
    aic_scores.append({'par': param, 'aic': model_fit.aic})
    
# Mencari skor aic terkecil
best_aic = min(aic_scores, key=lambda x: x['aic'])

print(best_aic)

# Membuat model ARIMA dengan p, d, dan q hasil grid search
model = ARIMA(data_train, order=best_aic['par'])
model_fit = model.fit()

# Melakukan prediksi 5 periode ke depan (data testing)
# (Agustus - Desember 2019)
preds = model_fit.forecast(5)

print(preds)

# alternatif prediksi periode data testing
# model_fit.predict(start=data_test.index[0], end=data_test.index[-1], typ="levels")

Output:

{'par': (3, 1, 3), 'aic': 819.6467518173176}

2019-08-01    53694.007736
2019-09-01    54841.950539
2019-10-01    56074.166643
2019-11-01    53314.075703
2019-12-01    57545.297465
Freq: MS, Name: predicted_mean, dtype: float64

Hasil grid search menggunakan data training memberikan nilai yang sedikit berbeda dari sebelumnya. Nilai AIC terkecil diperoleh pada p=3, d=1 dan q=3.

Selanjutnya, untuk mengevaluasi model, kita melakukan prediksi pada periode data testing. Hasil prediksi model kemudian dibandingkan dengan nilai sebenarnya. Di sini kita akan menggunakan metrik RMSE. Untuk menghitung nilai RMSE dapat menggunakan fungsi rmse pada library statsmodels.

Dari hasil penghitungan, diperoleh nilai RMSE data testing dari model ARIMA(3, 1, 3) adalah sebesar 2861,661. Jika melihat kisaran nilai pada data testing yang berada pada rentang 50.000 – 60.000. Nilai rata-rata error sebesar 2861,661 dapat kita katakan relatif kecil sehingga model yang kita peroleh sepertinya sudah cukup bagus.

Python

from statsmodels.tools.eval_measures import  rmse

rmse = rmse(data_test["nilai_jt_usd"], preds)

print(f"RMSE : {rmse}")



data.plot(figsize=(12, 6), alpha=0.5, marker="s")
preds.plot(linewidth=2, marker="o", legend=True)

Output

RMSE : 2861.6614013795875

Dalam evalausi model time series, terdapat teknik lainnya yang dapat digunakan untuk validasi dan evaluasi yaitu walk forward validation. Teknik ini mirip seperti K-Fold CV atau LOOVC namun tentunya dengan mempertimbangkan adanya keterkaitan antar waktu pada data time series. Topik ini akan dibahas lebih lanjut pada kesempatan lainnya.