Analisis Time Series ARIMA dengan Python

Melanjutkan pembahasan kita mengenai pemanfaatan Python di Era Data Science dan Big Data, kali ini kita akan belajar bersama mengenai satu alat analisis dari sekian alat analisis yang sangat sering digunakan dalam riset, terutama mengenai topik-topik peramalan (forecasting).

Alat analisis tersebut adalah Autoregressive Integrated and Moving Average atau biasa disingkat dengan ARIMA. Pada bagian terdahulu, kita telah membahas bersama bagaimana pemodelan ARIMA dan SARIMAX dengan menggunakan paket program R (bisa dicek pada tautan [1] dan [2]).

Adapun data yang akan kita gunakan kali ini adalah data inflasi bulanan (month to month) dari situs resmi Badan Pusat Statistik Provinsi Jawa Timur (jatim.bps.go.id). Inflasi ini sangat sering digunakan dalam analisis perekonomian suatu wilayah karena sifatnya yang series dan selalu tersedia datanya karena BPS selalu merilis angkanya pada setiap awal bulan. Dengan data inflasi, kita bisa melakukan baik nowcasting atau forecasting dengan alat-alat analisis yang ada, terutama menggunakan ARIMA.

Pada bagian unggahan ini, data yang diperlukan dapat diunduh di tautan berikut. Setelah datanya siap, langkah-langkah pemodelan ARIMA menggunakan Python (di sini menggunakan GUI jupyter notebook) adalah sebagai berikut:

Untitled

In [87]:

#Import data inflasi month to month BPS Provinsi Jawa Timur (jatim.bps.go.id)
import pandas as pd
df = pd.read_excel("E:\\R\\inflasi.xlsx")
df

Out[87]:

	bulan	inf
0	2015-01-01	0.20
1	2015-02-01	-0.52
2	2015-03-01	0.31
3	2015-04-01	0.39
4	2015-05-01	0.41
...	...	...
93	2022-10-01	0.42
94	2022-11-01	0.32
95	2022-12-01	0.60
96	2023-01-01	0.36
97	2023-02-01	0.10

98 rows × 2 columns

In [88]:

#Visualisasi data inflasi
df.inf.plot(figsize=(12,6))

Out[88]:

<AxesSubplot:>

In [89]:

#Uji Stasioneritas data inflasi
#Jika p-value < alpha (0,05) maka data stasioner
from statsmodels.tsa.stattools import adfuller
hasil = adfuller(df.inf)
print('ADF Statistic: %f' % hasil[0])
print('p-value: %f' % hasil[1])

ADF Statistic: -1.019581
p-value: 0.746056

In [90]:

#Uji Stasioneritas data difference 1 inflasi
#Jika p-value < alpha (0,05) maka data stasioner
from statsmodels.tsa.stattools import adfuller
hasil = adfuller(df.inf.diff().dropna())
print('ADF Statistic: %f' % hasil[0])
print('p-value: %f' % hasil[1])

ADF Statistic: -5.331544
p-value: 0.000005

In [91]:

#Melihat plot ACF untuk menentukan order MA
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(df.inf)

Out[91]:

In [92]:

#Melihat plot PACF untuk menentukan order AR
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(df.inf, lags = 15)

C:\Users\Joko Ade\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(

Out[92]:

In [93]:

#Membangun model ARIMA dengan order p=1, d=1 dan q=1, sesuai uji ADF dan plot ACF PACF
from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(df.inf, order=(1, 1, 1))
model_fit = model.fit()

#Ringkasan model ARIMA(1, 1, 1)
model_fit.summary()

C:\Users\Joko Ade\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'

Out[93]:

SARIMAX Results
Dep. Variable:	inf	No. Observations:	98
Model:	ARIMA(1, 1, 1)	Log Likelihood	-26.306
Date:	Sat, 11 Mar 2023	AIC	58.612
Time:	04:57:32	BIC	66.336
Sample:	0	HQIC	61.735
	- 98
Covariance Type:	opg

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	0.2709	0.106	2.554	0.011	0.063	0.479
ma.L1	-0.9998	4.456	-0.224	0.822	-9.733	7.734
sigma2	0.0966	0.424	0.228	0.820	-0.735	0.928

Ljung-Box (L1) (Q):	0.11	Jarque-Bera (JB):	42.53
Prob(Q):	0.74	Prob(JB):	0.00
Heteroskedasticity (H):	0.75	Skew:	1.00
Prob(H) (two-sided):	0.43	Kurtosis:	5.55

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

In [94]:

#Melakukan forecast 5 bulan ke depan
model_fit.forecast(5)

Out[94]:

98     0.219970
99     0.252469
100    0.261273
101    0.263657
102    0.264303
Name: predicted_mean, dtype: float64

In [95]:

#Membuat iterasi order ARIMA dengan p maksimal 1, d maksimal 1, dan q maksimal 1
from itertools import product

#List order p, d, dan q
p = range(0, 2)  # 0-1
d = range(0, 2)  # 0-1
q = range(0, 2)  # 0-1

#Hasil kombinasi p, d, dan q
pdq = list(product(p, d, q))

print(pdq)

[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)]

In [96]:

#Membuat sejumlah model ARIMA dengan order kombinasi dan mendapatkan model terbaik berdasarkan nilai AIC terkecil
from statsmodels.tsa.arima.model import ARIMA
#List penyimpanan nilai AIC
aic_scores = []
#Membentuk ARIMA dengan kombinasi p, d, q yang optimal
for param in pdq:
    #Model fit ARIMA
    model = ARIMA(df.inf, order=param)
    model_fit = model.fit()
    #Menambahkan variabel AIC untuk setiap model ARIMA yang terbentuk
    aic_scores.append({'par': param, 'aic': model_fit.aic})
    
#Model ARIMA dengan AIC terkecil
best_aic = min(aic_scores, key=lambda x: x['aic'])
best_aic

C:\Users\Joko Ade\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'

Out[96]:

{'par': (0, 0, 1), 'aic': 52.271878948035464}

In [97]:

#Membuat model ARIMA dengan AIC terkecil
model = ARIMA(df.inf, order=best_aic['par']) # order=(3, 0, 2)
model_fit = model.fit()

#Ringkasan model
model_fit.summary()

Out[97]:

SARIMAX Results
Dep. Variable:	inf	No. Observations:	98
Model:	ARIMA(0, 0, 1)	Log Likelihood	-23.136
Date:	Sat, 11 Mar 2023	AIC	52.272
Time:	04:58:01	BIC	60.027
Sample:	0	HQIC	55.409
	- 98
Covariance Type:	opg

	coef	std err	z	P>\|z\|	[0.025	0.975]
const	0.2655	0.046	5.740	0.000	0.175	0.356
ma.L1	0.3163	0.104	3.038	0.002	0.112	0.520
sigma2	0.0938	0.009	10.137	0.000	0.076	0.112

Ljung-Box (L1) (Q):	0.01	Jarque-Bera (JB):	66.12
Prob(Q):	0.92	Prob(JB):	0.00
Heteroskedasticity (H):	0.79	Skew:	1.04
Prob(H) (two-sided):	0.51	Kurtosis:	6.45

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

In [98]:

#Melakukan forecast 5 bulan ke depan
model_fit.forecast(5)

Out[98]:

98     0.212974
99     0.265470
100    0.265470
101    0.265470
102    0.265470
Name: predicted_mean, dtype: float64

In [99]:

#Memecah data menjadi 2 bagian, data train dan test
data_train = df[:-5]
data_test = df[-5:]

#Melihat data test
data_test

Out[99]:

	bulan	inf
93	2022-10-01	0.42
94	2022-11-01	0.32
95	2022-12-01	0.60
96	2023-01-01	0.36
97	2023-02-01	0.10

In [100]:

#Membuat sejumlah model ARIMA dengan order kombinasi dan mendapatkan model terbaik berdasarkan nilai AIC terkecil
from statsmodels.tsa.arima.model import ARIMA
#List penyimpanan nilai AIC
aic_scores = []
#Membentuk ARIMA dengan kombinasi p, d, q yang optimal
for param in pdq:
    #Model fit ARIMA
    model = ARIMA(data_train.inf, order=param)
    model_fit = model.fit()
    #Menambahkan variabel AIC untuk setiap model ARIMA yang terbentuk
    aic_scores.append({'par': param, 'aic': model_fit.aic})
    
#Model ARIMA dengan AIC terkecil
best_aic = min(aic_scores, key=lambda x: x['aic'])
best_aic

C:\Users\Joko Ade\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'

Out[100]:

{'par': (0, 0, 1), 'aic': 52.79079707247058}

In [101]:

#Membuat model ARIMA dengan AIC terkecil
model = ARIMA(data_train.inf, order=best_aic['par'])
model_fit = model.fit()

#Rangkuman model
model_fit.summary()

Out[101]:

SARIMAX Results
Dep. Variable:	inf	No. Observations:	93
Model:	ARIMA(0, 0, 1)	Log Likelihood	-23.395
Date:	Sat, 11 Mar 2023	AIC	52.791
Time:	04:58:19	BIC	60.389
Sample:	0	HQIC	55.859
	- 93
Covariance Type:	opg

	coef	std err	z	P>\|z\|	[0.025	0.975]
const	0.2658	0.049	5.373	0.000	0.169	0.363
ma.L1	0.3516	0.113	3.114	0.002	0.130	0.573
sigma2	0.0967	0.010	9.721	0.000	0.077	0.116

Ljung-Box (L1) (Q):	0.01	Jarque-Bera (JB):	62.20
Prob(Q):	0.94	Prob(JB):	0.00
Heteroskedasticity (H):	0.83	Skew:	1.04
Prob(H) (two-sided):	0.60	Kurtosis:	6.43

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

In [103]:

#Melakukan forecast 5 bulan ke depan
model_fit.forecast(5)

Out[103]:

93    0.701042
94    0.265799
95    0.265799
96    0.265799
97    0.265799
Name: predicted_mean, dtype: float64

In [104]:

#Melihat RMSE untuk hasil forecast 5 bulan ke depan
from statsmodels.tools.eval_measures import  rmse
rmse = rmse(data_test["inf"], preds)
print(rmse)

0.21446480374898652

In [105]:

#Visualisasi data hasil forecast dan data aktual
df.inf.plot(figsize=(12, 6), alpha=0.5, marker="s")
preds.plot(linewidth=2, marker="o", legend=True)

Out[105]:

<AxesSubplot:>

Dari hasil pemodelan ARIMA(0, 0, 1) di atas, terlihat bahwa model yang terbentuk dan digunakan dalam forecasting cukup baik (nilai p value Ljung Box > alpha 0,05 demikian pula uji homoskedastisitas residual) meski residual modelnya melanggar asumsi kenormalan. Terlihat pada nilai p value uji Jarque Bera kurang dari 0,05 yang menyatakan gagal terima Ho (residual mengikuti distribusi normal).

Meski demikian, kalau dilihat berdasarkan nilai Root Mean Square Error (RMSE) relatif kecil kendati masih di atas 10 persen. Oleh karena itu, direkomendasikan periode forecasting tidak terlalu panjang karena berisiko bias akan terus membesar sebagai efek ketidaknormalan residual model.

Demikian sedikit sharing kita pada kesempatan kali ini. Jangan lupa untuk terus mengikuti setiap unggahan baru, menarik, dan unik dalam blog sederhana ini. Semoga sedikit atau banyak memberi manfaat bagi seluruh pembaca. Selamat memahami dan mempraktikkan!

Blog belajar dan berbagi tentang bahasa pemrograman R, Python, dan visualisasi data Maret 11, 2023

Data Data Mining Modelling Python Tutorial