 |
| 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:
#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
#Import pandas untuk memanggil data
import pandas as pd
#Import data harga emas
data = pd.read_excel("C:\\Users\\56848\\Documents\\emasku.xlsx")
data
| 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
#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()
#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(
#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
#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
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
#Aktivasi kombinasi order SARIMA
from itertools import product
#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
from tqdm import tqdm_notebook
#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
| (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
#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
| (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
#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(
#Pengecekan residual model SARIMA
mod_terbaik.plot_diagnostics(figsize=(15,12));
#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!