Seasonal ARIMA dengan Python (Step by Step)

SARIMA Model dengan Python

Pada sebuah unggahan beberapa waktu sebelumnya, kita telah membahas bersama bagaimana praktik memodelkan Seasonal Autoregressive Integrated and Moving Average (SARIMA) dengan R. Pada unggahan kali ini, kita akan mencoba melakukan praktik pemodelan SARIMA menggunakan Python dengan GUI-nya Jupyter Notebook.

SARIMA merupakan salah satu model runtun waktu (time series) yang sering digunakan dalam melakukan estimasi atau peramalan baik di masa lampau maupun masa depan dengan berbekal data masa lalu. Prinsip dasar dari SARIMA adalah melakukan pemodelan berdasarkan pola musiman dengan memerhatikan rata-rata dan varians data pada periode sebelumnya.

Tipe data runtun waktu yang mengandung pola musiman ditengarai dengan adanya efek-efek musim yang terlihat dalam data. Periode musiman bisa terjadi secara bulanan, triwulanan, semesteran, atau tahunan. Sebagaimana ketika kita amati produksi buah semangka. Sebagai buah musiman, produksi semangka tidak selalu banyak di setiap bulannya. Dalam setahun, buah ini biasanya dipanen sebanyak 3-4 kali.

Itu artinya, pada periode-periode tertentu, produksi semangka dapat mencapai puncaknya. Sedangkan di lain waktu, semangka justru sangat langka dan tidak mudah ditemui di pasaran.

Dalam praktikum pemodelan SARIMA kali ini, data yang kita pakai bukanlah data produksi semangka. Melainkan data harga emas yang dapat diperoleh pada tautan berikut. Setelah datanya telah diunduh dan siap dimodelkan, berikut langkah-langkah pemodelan SARIMA dengan Python sebagai berikut:

In [54]:

#Import dan aktivasi package
import statsmodels as sm
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.stattools import adfuller
import matplotlib.pyplot as plt
import numpy as np

In [55]:

#Import pandas untuk memanggil data
import pandas as pd

In [56]:

#Import data harga emas
data = pd.read_excel("C:\\Users\\56848\\Documents\\emasku.xlsx")
data

Out[56]:

	tgl	hrg
0	2018-03-09	18204469.0
1	2018-03-16	18014528.0
2	2018-03-23	18558165.9
3	2018-03-30	18226102.9
4	2018-04-06	18337944.9
...	...	...
257	2023-02-09	28364076.4
258	2023-02-16	27720479.9
259	2023-02-23	27737700.2
260	2023-03-02	28057135.3
261	2023-03-03	28169593.5

262 rows × 2 columns

In [57]:

#Visualisasi data harga emas
plt.figure(figsize=[15, 7.5]); #Menentukan dimensi visualisasi
plt.plot(data['tgl'], data['hrg'])
plt.title('Harga Emas per Oz') #judul grafik
plt.ylabel('Harga (US$)') #judul ordinat
plt.xlabel('tanggal') #judul absis
plt.xticks(rotation=90)
plt.grid(True)
plt.show()

In [58]:

#Plot ACF dan PACF
#Plot ACF untuk menentukan order MA
#Plot PACF untuk menentukan order AR
plot_pacf(data['hrg']);
plot_acf(data['hrg']);

C:\Users\56848\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning: The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
  warnings.warn(

In [59]:

#Uji ADF untuk menentukan order difference (d)
adf0 = adfuller(data['hrg'])
print(f'ADF Statistic: {adf0[0]}')
print(f'p-value: {adf0[1]}')

ADF Statistic: -1.1962536051436001
p-value: 0.675205236451486

In [60]:

#Karena tidak stasioner pada level d = 0, maka uji ADF pada diff 1 atau d = 1
adf1 = adfuller(data['hrg'].diff().dropna())
print(f'ADF Statistic: {adf1[0]}')
print(f'p-value: {adf1[1]}')

ADF Statistic: -6.627815218860708
p-value: 5.816158480071784e-09

In [61]:

def optimize_SARIMA(parameters_list, d, D, s, exog): #D adalah difference musiman, s adalah periode musiman
    
    results = []
    
    for param in tqdm_notebook(parameters_list):
        try: 
            model = SARIMAX(exog, order=(param[0], d, param[1]), seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
        except:
            continue
            
        aic = model.aic
        results.append([param, aic])
        
    result_df = pd.DataFrame(results)
    result_df.columns = ['(p,q)x(P,Q)', 'AIC']
    #Melakukan sorting berdasarkan nilai AIC terkecil
    result_df = result_df.sort_values(by='AIC', ascending=True).reset_index(drop=True)
    
    return result_df

In [62]:

#Aktivasi kombinasi order SARIMA
from itertools import product

In [63]:

#Menentukan sejumlah kombinasi order SARIMA(p, d, q)(P, D, Q)
p = range(0, 4, 1)
d = 1
q = range(0, 4, 1)
P = range(0, 4, 1)
D = 1
Q = range(0, 4, 1)
s = 4
parameter = product(p, q, P, Q)
daftar_parameter = list(parameter)

#Banyak kombinasi order
print(len(daftar_parameter))

256

In [64]:

from tqdm import tqdm_notebook

In [65]:

#Model tentatif berdasarkan nilai AIC
hasil_ten = optimize_SARIMA(daftar_parameter, 1, 1, 4, data['hrg'])
hasil_ten

C:\Users\56848\AppData\Local\Temp\ipykernel_6784\1988398620.py:5: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0
Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`
  for param in tqdm_notebook(parameters_list):

  0%|          | 0/256 [00:00<?, ?it/s]

C:\Users\56848\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to 

Out[65]:

	(p,q)x(P,Q)	AIC
0	(0, 0, 0, 3)	7602.104865
1	(0, 0, 3, 3)	7603.533947
2	(0, 0, 1, 1)	7605.428019
3	(0, 0, 2, 1)	7607.059058
4	(0, 0, 1, 2)	7608.304828
...	...	...
251	(0, 2, 0, 0)	7669.222716
252	(2, 0, 0, 0)	7669.397991
253	(0, 3, 0, 0)	7670.220241
254	(1, 2, 0, 0)	7670.998407
255	(3, 0, 0, 0)	7671.116292

256 rows × 2 columns

In [66]:

#Mengurutkan hasil order tentatif menurut nilai AIC terkecil ke terbesar
tabulasi_urut = result_df.sort_values(by='AIC',ascending = True).reset_index(drop = True)
tabulasi_urut

Out[66]:

	(p,q)x(P,Q)	AIC
0	(0, 0, 0, 3)	7602.104865
1	(0, 0, 3, 3)	7603.533947
2	(0, 0, 1, 1)	7605.428019
3	(0, 0, 2, 1)	7607.059058
4	(0, 0, 1, 2)	7608.304828
...	...	...
251	(0, 2, 0, 0)	7669.222716
252	(2, 0, 0, 0)	7669.397991
253	(0, 3, 0, 0)	7670.220241
254	(1, 2, 0, 0)	7670.998407
255	(3, 0, 0, 0)	7671.116292

256 rows × 2 columns

In [88]:

#Memodelkan berdasarkan order dengan AIC terendah
mod_terbaik = SARIMAX(data['hrg'], order=(0, 1, 2), seasonal_order=(0, 0, 3, 4)).fit(dis=-1)
print(mod_terbaik.summary())

                                         SARIMAX Results                                         
=================================================================================================
Dep. Variable:                                       hrg   No. Observations:                  262
Model:             SARIMAX(0, 1, 2)x(0, 0, [1, 2, 3], 4)   Log Likelihood               -3788.114
Date:                                   Mon, 20 Mar 2023   AIC                           7588.227
Time:                                           15:40:17   BIC                           7609.614
Sample:                                                0   HQIC                          7596.824
                                                   - 262                                         
Covariance Type:                                     opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ma.L1         -0.0120      0.050     -0.238      0.812      -0.110       0.086
ma.L2          0.0388      0.046      0.841      0.400      -0.052       0.129
ma.S.L4       -0.0091      0.037     -0.244      0.808      -0.082       0.064
ma.S.L8       -0.0057      0.058     -0.099      0.921      -0.119       0.108
ma.S.L12      -0.0174      0.053     -0.329      0.742      -0.121       0.086
sigma2      2.428e+11   6.77e-14   3.58e+24      0.000    2.43e+11    2.43e+11
===================================================================================
Ljung-Box (L1) (Q):                   0.16   Jarque-Bera (JB):               125.45
Prob(Q):                              0.69   Prob(JB):                         0.00
Heteroskedasticity (H):               2.80   Skew:                             0.09
Prob(H) (two-sided):                  0.00   Kurtosis:                         6.39
===================================================================================

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

C:\Users\56848\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\base\optimizer.py:17: FutureWarning: Keyword arguments have been passed to the optimizer that have no effect. The list of allowed keyword arguments for method lbfgs is: m, pgtol, factr, maxfun, epsilon, approx_grad, bounds, loglike_and_score. The list of unsupported keyword arguments passed include: dis. After release 0.14, this will raise.
  warnings.warn(

In [89]:

#Pengecekan residual model SARIMA
mod_terbaik.plot_diagnostics(figsize=(15,12));

In [92]:

#Visualisasi hasil model SARIMA terbaik, data aktual, dan hasil forecast
data['arima_model'] = mod_terbaik.fittedvalues
data['arima_model'][:4+1] = np.NaN
fc = mod_terbaik.predict(start=data.shape[0], end=data.shape[0] + 8)
fc = data['arima_model'].append(fc)
plt.figure(figsize=(15, 7.5))
plt.plot(fc, color='r', label='hasil model')
plt.axvspan(data.index[-1], fc.index[-1], alpha=0.5, color='lightgrey')
plt.plot(data['hrg'], label=' harga aktual')
plt.legend()
plt.show()

C:\Users\56848\AppData\Local\Temp\ipykernel_6784\430267613.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['arima_model'][:4+1] = np.NaN
C:\Users\56848\AppData\Local\Temp\ipykernel_6784\430267613.py:5: FutureWarning: The series.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  fc = data['arima_model'].append(fc)

Hasil pemodelan SARIMA pada kasus ini dari sisi akurasi memang terbilang cukup baik. Terlihat dari perbandingan antara data hasil prediksi model dengan data aktual harga emas berhimpitan. Kendati demikian, uji normalitas residual model belum normal. Terlihat nilai p-value uji Jarque Bera < alpha atau 0,05 sehingga residual model tidak mengikuti distribusi normal sehingga bila digunakan untuk peramalan jangka panjang justru biasnya semakin besar.

Demikian sedikit sharing kali ini. Semoga bermanfaat. Jangan lupa untuk stay tune di blog sederhana ini. Selamat memahami dan mempraktikkan!