Python 시계열 데이터 분석 - AR, MA, ARIMA, SARIMAX, 분해, 지수평활

시계열 데이터 분석이라는 다소 넓은 범위의 주제이지만, 그 중에서도 전통적인 시계열 분석에서 언급되는 주요 분석 기법에 대해서 간략히 알아보고, 파이썬(Python) 코드로 구현하는 방법에 대해 알아보겠습니다.

필요한 모듈 import

시계열 데이터 분석에 필요한 패키지를 설치합니다.

# !pip install pmdarima -q
# !pip install statsmodels -q

import os
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
import warnings

# SEED 설정
SEED = 123

#####################    Pandas   #####################

# DataFrame 행과 열의 출력 개수 지정. None은 전체 출력.
pd.set_option('display.max_columns', None)

# 1,000 단위 표기
pd.options.display.float_format = '{:,.3f}'.format

##################### Matplotlib #####################

# Unicode warning 제거 (폰트 관련 경고메시지)
plt.rcParams['axes.unicode_minus']=False

### 한글 폰트 설정 ###
# 나눔고딕 폰트
plt.rcParams['font.family'] = "NanumGothic"
# 맑은고딕 폰트
# plt.rcParams['font.family'] = "Malgun Gothic"
# 애플고딕 (Mac OS)
# plt.rcParams['font.family'] = "AppleGothic"

# 그래프 출력 사이즈 설정
plt.rcParams["figure.figsize"] = (10, 6)

# Seaborn
# sns.set_style("white")

COLORS = {
    'green': '#278F45',
    'red': '#DB4A25', 
    'blue': '#23D3DB',
    'black': '#000000', 
    'yellow': '#DBBE2E',
}

%matplotlib inline
######################## Others  #########################

# 노트북 출력에 나타나는 경고(warning) 메시지 무시
warnings.filterwarnings('ignore')

데이터 로드

• 금융 데이터 로드를 위한 Finance DataReader 모듈 설치

# Finance DataReader 설치
# !pip install -U finance-datareader

import FinanceDataReader as fdr

df = fdr.DataReader('BTC/USD')
df.head()

	Open	High	Low	Close	Adj Close	Volume
Date
2014-09-17	465.864	468.174	452.422	457.334	457.334	21056800
2014-09-18	456.860	456.860	413.104	424.440	424.440	34483200
2014-09-19	424.103	427.835	384.532	394.796	394.796	37919700
2014-09-20	394.673	423.296	389.883	408.904	408.904	36863600
2014-09-21	408.085	412.426	393.181	398.821	398.821	26580100

freq에 원하는 시계열 데이터의 주기를 입력할 수 있습니다.

• 참고링크: https://pandas.pydata.org/docs/user_guide/timeseries.html#timeseries-offset-aliases

# 주간(Weekly) 데이터만 추출하기 위하여 date_range()를 주 단위로 생성합니다.
target_index = pd.date_range('20130107', '20230701', freq='W-MON')
target_index

DatetimeIndex(['2013-01-07', '2013-01-14', '2013-01-21', '2013-01-28',
               '2013-02-04', '2013-02-11', '2013-02-18', '2013-02-25',
               '2013-03-04', '2013-03-11',
               ...
               '2023-04-24', '2023-05-01', '2023-05-08', '2023-05-15',
               '2023-05-22', '2023-05-29', '2023-06-05', '2023-06-12',
               '2023-06-19', '2023-06-26'],
              dtype='datetime64[ns]', length=547, freq='W-MON')

# 추출한 시계열 인덱스 데이터로 필터링
df.loc[df.index.isin(target_index)].head(10)

	Open	High	Low	Close	Adj Close	Volume
Date
2014-09-22	399.100	406.916	397.130	402.152	402.152	24127600
2014-09-29	376.928	385.211	372.240	375.467	375.467	32497700
2014-10-06	320.389	345.134	302.560	330.079	330.079	79011800
2014-10-13	377.921	397.226	368.897	390.414	390.414	35221400
2014-10-20	389.231	390.084	378.252	382.845	382.845	16419000
2014-10-27	354.777	358.632	349.809	352.989	352.989	13033000
2014-11-03	325.569	334.002	325.481	327.554	327.554	12948500
2014-11-10	362.265	374.816	357.561	366.924	366.924	30450100
2014-11-17	388.349	410.199	377.502	387.408	387.408	41518800
2014-11-24	366.948	387.209	366.669	376.901	376.901	30930100

# 2020년 ~ 2022년 3년치의 시세 데이터를 조회합니다. (Close)는 종가를 의미
data = df.loc[(df.index >= '2020-01-01') & (df.index <= '2022-12-31'), 'Close']
data

Date
2020-01-01    7,200.174
2020-01-02    6,985.470
2020-01-03    7,344.884
2020-01-04    7,410.657
2020-01-05    7,411.317
                ...    
2022-12-27   16,717.174
2022-12-28   16,552.572
2022-12-29   16,642.342
2022-12-30   16,602.586
2022-12-31   16,547.496
Name: Close, Length: 1096, dtype: float64

시계열 데이터 결측치 처리

Simple Imputer

• 단순하게 평균, 중앙값 등으로 채웁니다.

• 하지만, 시계열 데이터의 특성상 앞,뒤 데이터의 추세를 반영하지 못합니다.

# 결측치 확인
data.isnull().sum()

0

# 평균 값으로 채우는 경우
data.fillna(data.mean())

Date
2020-01-01    7,200.174
2020-01-02    6,985.470
2020-01-03    7,344.884
2020-01-04    7,410.657
2020-01-05    7,411.317
                ...    
2022-12-27   16,717.174
2022-12-28   16,552.572
2022-12-29   16,642.342
2022-12-30   16,602.586
2022-12-31   16,547.496
Name: Close, Length: 1096, dtype: float64

Interpolate

• 앞, 뒤 데이터의 추세를 보고 interpolation 값을 계산하여 채웁니다.

# interpolation 방식으로 결측치를 채웁니다.
data = data.interpolate()
data

Date
2020-01-01    7,200.174
2020-01-02    6,985.470
2020-01-03    7,344.884
2020-01-04    7,410.657
2020-01-05    7,411.317
                ...    
2022-12-27   16,717.174
2022-12-28   16,552.572
2022-12-29   16,642.342
2022-12-30   16,602.586
2022-12-31   16,547.496
Name: Close, Length: 1096, dtype: float64

# 결측치 확인
data.isnull().sum()

0

지수평활 (Exponential Smoothing)

지수평활(Exponential Smoothing) 은 시계열 데이터 분석에서 일반적으로 사용되는 기법 중 하나입니다. 이는 데이터의 패턴을 예측하는 데 도움이 되는, 데이터의 노이즈를 줄이고 추세를 더 잘 파악 할 수 있는 방법입니다.

지수평활의 기본적인 아이디어는 최근의 관찰값에 더 많은 가중치 를 주고, 과거의 관찰값에는 적은 가중치 를 주는 것입니다. ‘지수’라는 용어는 각 관찰값에 곱해지는 가중치가 지수적으로 감소하기 때문에 사용됩니다.

이 방법은 다음과 같은 이유로 사용됩니다:

• 간단함: 지수평활은 계산이 간단하고 이해하기 쉬운 방법입니다. 이는 데이터가 복잡한 패턴이나 추세를 보이지 않을 때 특히 유용합니다.

• 노이즈 감소: 지수평활은 원래의 시계열 데이터에서 노이즈를 줄이고 주요한 패턴을 강조함으로써 데이터를 더 잘 이해하는 데 도움이 됩니다.

• 예측: 지수평활은 미래 값의 예측에 사용할 수 있습니다. 이 방법은 과거의 데이터를 사용하여 미래의 데이터를 예측하는 방법 중 하나로, 데이터의 패턴이 시간이 지남에 따라 일정하게 유지될 것으로 예상되는 경우 특히 유용합니다.

• 빠른 업데이트: 지수평활은 새로운 데이터가 추가될 때마다 쉽게 업데이트될 수 있습니다. 이는 특히 실시간 분석이 필요한 경우에 유용합니다.

지수평활에는 다양한 변형이 있습니다.

• 가장 간단한 형태의 지수평활은 단순 지수평활(Simple Exponential Smoothing, SES)이며, 이는 데이터에 선형적인 패턴이 없을 때 사용됩니다.
• 더 복잡한 형태의 지수평활은 홀트의 선형 지수평활(Holt's Linear Exponential Smoothing) 이나 홀트-윈터스의 계절적 지수평활(Holt-Winters' Seasonal Exponential Smoothing) 과 같이 트렌드와 계절성을 포함하는 데이터를 다룰 수 있습니다.

fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 6)
fig.set_dpi(300)
ax.plot(data, label='Close', color='black', linestyle='solid', linewidth=0.5)
ax.plot(data.ewm(alpha=0.1).mean(), label='alpha=0.1', color=COLORS['red'], linestyle='solid', linewidth=0.8)
ax.plot(data.ewm(alpha=0.5).mean(), label='alpha=0.5', color=COLORS['blue'], linestyle='solid', linewidth=0.7)
ax.plot(data.ewm(alpha=0.8).mean(), label='alpha=0.8', color=COLORS['green'], linestyle='solid', linewidth=0.5)

plt.xticks(rotation=30)
plt.legend()
plt.show()

이동평균(Moving Average)

이동 평균은 주어진 시점에서 이전 n 개의 데이터 포인트의 평균을 계산 하는 것을 의미합니다.

이는 시계열 데이터의 잡음을 줄이고 데이터의 기본적인 패턴을 확인하는 데 도움이 됩니다.

이동평균(moving average) 를 구하는 방식은 다음과 같습니다.

(예시)

• 10일 이동평균을 계산

data.rolling(window=10).mean()

fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 6)
fig.set_dpi(300)
ax.plot(data, label='Close', color='black', linestyle='solid', linewidth=0.5)
ax.plot(data.rolling(window=10).mean().bfill(), label='ma10', color=COLORS['red'], linestyle='solid', linewidth=0.8)
ax.plot(data.rolling(window=20).mean().bfill(), label='ma20', color=COLORS['blue'], linestyle='solid', linewidth=0.7)
ax.plot(data.rolling(window=60).mean().bfill(), label='ma60', color=COLORS['green'], linestyle='solid', linewidth=0.5)

plt.xticks(rotation=30)
plt.legend()
plt.show()

정상성(stationary) / 비정상성(non-stationary)

이미지 출처: https://t1.daumcdn.net/cfile/tistory/99E7CC3E5D13859D30

정상 시계열 (stationary)

정상 시계열은 시간의 흐름에 따라 그 통계적 특성이 변하지 않는 시계열 을 의미합니다.

다시 말해, 정상 시계열의 평균과 분산(또는 표준편차), 그리고 공분산은 시간이 변해도 일정 합니다.

이러한 특성 때문에 정상 시계열은 분석이 용이하며, 많은 통계적 시계열 모델들(AR, MA, ARMA 등)은 데이터가 정상성을 가정 하고 있습니다. (여기서, 정상성을 가정한다는 점이 매우 중요합니다!!)

(예시)

• 시간에 따른 일정한 주기를 가진 신호(일정한 주기로 진동하는 심장 박동)

• 백색 잡음(랜덤한 값들의 연속) 등이 있습니다.

비정상 시계열 (non-stationary)

비정상 시계열은 시간의 흐름에 따라 그 통계적 특성이 변하는 시계열 을 의미합니다. 즉, 평균, 분산, 공분산 등이 시간에 따라 변화합니다.

비정상 시계열 데이터는 시간의 흐름에 따라 통계적 속성이 변하는 데이터를 의미합니다. 이는 추세(trend) 또는 계절성(seasonality) 때문에 발생할 수 있습니다.

• 추세(trend): 시계열 데이터가 시간에 따라 증가하거나 감소하는 경향 을 보일 때, 이를 추세 라고 합니다.

• 계절성(seasonality): 시계열 데이터가 일정한 주기로 반복되는 패턴 을 보일 때, 이를 계절성 이라고 합니다.

대부분의 실제 시계열 데이터는 비정상적인 성질을 가지고 있습니다. 이 경우에는 데이터를 정상적으로 만드는 변환(예: 차분, 로그 변환 등) 이 필요할 수 있습니다.

(예시)

• 시간에 따라 트렌드가 있는 주가 데이터

• 인구 증가 데이터 등

이러한 데이터는 시간에 따라 평균과 분산이 변하기 때문에 비정상적(non-stationary) 입니다.

정상성(stationary) / 비정상성(non-stationary) 확인

시각화를 통한 확인

ACF(Auto Correlation Function) Plot - 시각화

• 링크: https://www.statsmodels.org/stable/generated/statsmodels.graphics.tsaplots.plot_acf.html

분석 tip 1. 통계적 유의성 확인

• 자기 상관이 통계적으로 유의미한지 여부를 나타내는 “신뢰구간”을 회색 음영 영역(shading area) 으로 표시

• 기본적으로 95% 신뢰구간 을 사용하여 shading area를 그립니다.

• 각 Lag 에서 Shading 영역을 벗어나 자기상관계수(Autocorrelation) 추정치는 95% 신뢰구간을 벗어났기 때문에, 통계적 유의성을 찾을 수 없다 는 의미로 해석할 수 있습니다.

분석 tip 2. 정상성 확인

• 보통은 시차가 커질수록 상관계수는 점차 0에 가까워지며, 시차가 커질 수록 자기상관성이 떨어지는 양상을 보여주는 것이 일반적입니다. 하지만, Seasonal 한 데이터의 경우, 값이 다시 튀어오르는 현상은 있을 수 있습니다.

• ACF는 시계열 데이터의 정상성을 판단할 때 유용 하게 활용될 수 있습니다.

  • 빠르게 0으로 수렴하는 경우 -> 정상(stationary) 시계열, 천천히 감소합니다.

  • 큰 상관계수를 가지는 경우 비정상(non-stationary) 시계열 데이터로 평가할 수 있습니다.

수식

• ACF는 시계열 데이터의 lagged 복사본과 원본 사이의 상관 관계를 측정합니다.

$\rho_k = \frac{Cov(X_t, X_{t+k})}{\sqrt{Var[X_t]*Var[X_{t+k}]}}$

PACF(Partial Autocorrelation Function) Plot - 시각화

• 시간차(lag)가 늘어나면서 관측치 간의 상관 관계가 어떻게 변하는지 설명 합니다.

• 원리적으로, PACF는 선형 회귀분석을 통해 계산됩니다. 각 시차의 관측치를 이전 시차들의 관측치로 선형 회귀분석 하고, 이때의 회귀 계수를 PACF로 간주합니다. 이렇게 하면, PACF는 이전 시차들의 영향을 배제하고 특정 시차만의 자기상관을 측정할 수 있습니다.

분석 tip 1. AR 모델의 차수 결정

• PACF를 사용하면 AR 모델의 차수 를 결정하는 데 도움이 될 수 있습니다.

수식

• PACF는 다른 모든 시점의 영향을 배제하고 두 시점 사이의 상관 관계를 측정합니다.

$\phi_{kk} = Cor r(X_t, X_{t+k} | X_{t+1}, …, X_{t+k-1})$

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf


fig, axes = plt.subplots(1, 2)
fig.set_size_inches(12, 4)

# ACF Plot
plot_acf(data, lags=30, ax=axes[0])
# PACF Plot
plot_pacf(data, lags=10, zero=False, ax=axes[1])

for ax in axes:
    ax.set_ylim(-0.5, 1.25)
plt.show()

정상성 검증

kpss test

• 귀무가설: 해당 시계열은 정상(stationary) 시계열이다. <-> 대립가설: 해당 시계열은 비정상(non-stationary) 시계열이다.

  • p-value <= 0.05 : 귀무가설 기각/ 대립가설 채택 -> 비정상(non-stationary) 시계열.

  • p-value > 0.05 : 귀무가설 채택/ 대립가설 기각 -> 정상(stationary) 시계열.

from statsmodels.tsa.stattools import kpss

def kpss_test(series, **kw):    
    stats, p_value, nlags, critical_values = kpss(series, **kw)
    
    print(f'KPSS Stat: {stats:.5f}')
    print(f'p-value: {p_value:.2f}')
    print(f'Lags: {nlags}')
    
    print(f'검증결과: {"비정상(non-stationary)" if p_value <= 0.05 else "정상(stationary)"} 시계열 데이터입니다.')

kpss_test(data, nlags=30)

KPSS Stat: 1.18454
p-value: 0.01
Lags: 30
검증결과: 비정상(non-stationary) 시계열 데이터입니다.

# 1차 차분에 대한 결과
diff_1 = data.diff(periods=1).iloc[1:]

kpss_test(diff_1, nlags=30)

KPSS Stat: 0.29060
p-value: 0.10
Lags: 30
검증결과: 정상(stationary) 시계열 데이터입니다.

# 1차 차분 결과에 대한 시각화
fig, axes = plt.subplots(1, 1)
fig.set_size_inches(8, 4)

# ACF Plot
plot_acf(diff_1, lags=30, ax=axes)

axes.set_ylim(-0.5, 1.25)
plt.show()

Adfuller 테스트

• kpss test 와 귀무가설이 반대입니다.

• 귀무가설: 해당 시계열은 비정상(non-stationary) 시계열이다. <-> 대립가설: 해당 시계열은 정상(stationary) 시계열이다.

  • p-value <= 0.05 : 귀무가설 기각/ 대립가설 채택 -> 정상(stationary) 시계열.

  • p-value > 0.05 : 귀무가설 채택/ 대립가설 기각 -> 비정상(non-stationary) 시계열.

# 데이터 정상 상태 수치 확인 
from statsmodels.tsa.stattools import adfuller #ADF Test를 위한 함수 호출 

st_result = adfuller(data)

def adfuller_test(series, **kw):    
    adf, p_value, nlags, number_of_observations, critical_values, _ = adfuller(series, **kw)
    
    print(f'ADF: {adf:.5f}')
    print(f'p-value: {p_value:.2f}')
    print(f'Lags: {nlags}')
    print(f'Number of Observations: {number_of_observations}')
    
    print(f'검증결과: {"비정상(non-stationary)" if p_value > 0.05 else "정상(stationary)"} 시계열 데이터입니다.')

# Adfuller 테스트
adfuller_test(data, maxlag=30)

ADF: -1.44857
p-value: 0.56
Lags: 0
Number of Observations: 1095
검증결과: 비정상(non-stationary) 시계열 데이터입니다.

# 1차 차분에 대한 Adfuller 테스트
adfuller_test(diff_1, maxlag=30)

ADF: -10.34401
p-value: 0.00
Lags: 8
Number of Observations: 1086
검증결과: 정상(stationary) 시계열 데이터입니다.

ARIMA Model

시계열 데이터 분석에서는 AR(AutoRegressive), MA(Moving Average), ARMA(AutoRegressive Moving Average), ARIMA(AutoRegressive Integrated Moving Average) 등의 모델이 있습니다.

• AR (AutoRegressive): 이전 관측치의 값들이 현재 관측치에 영향을 주는 모델입니다. 예를 들어, AR(1) 모델은 현재 관측치가 바로 이전 관측치에 의해 영향을 받는다는 것을 나타냅니다. AR 모델은 파라미터 p로 몇 개의 이전 관측치를 고려할지를 결정합니다.

• MA (Moving Average): 이전 오차항들이 현재 관측치에 영향을 주는 모델입니다. MA(1) 모델은 현재 관측치가 바로 이전의 오차항에 의해 영향을 받는다는 것을 나타냅니다. MA 모델은 파라미터 q로 몇 개의 이전 오차항을 고려할지를 결정합니다.

AR & MA 모델의 차이

AR(AutoRegressive) 모델과 MA(Moving Average) 모델은 둘 다 시계열 분석에서 널리 사용되는 모델이지만, 그 방식이 다릅니다.

AR 모델

• AR 모델은 ‘자기회귀 모델’ 로, 과거의 자신의 값에 의존하여 현재 값을 예측하는 방법을 사용합니다. 즉, 이전 시점의 관측값이 현재 시점의 관측값에 영향을 줍니다.
• 예를 들어, 주식 시장에서 어제 주가가 상승했다면 오늘 주가도 상승할 가능성이 있다는 것 이 AR 모델의 관점입니다.

AR 모델은 파라미터 p로 이전 몇 개의 관측치를 고려할지 결정 하게 됩니다. AR(1) 모델은 바로 이전의 한 단계 관측치를 고려하고, AR(2) 모델은 바로 이전의 두 단계 관측치를 고려하게 됩니다.

(예시)

예를 들어, 전력 소비를 예측하는 문제를 생각해 보겠습니다. 사람들은 특정 시간대에 일관되게 전력을 사용하는 경향이 있습니다. 아침에 사람들이 일어나서 활동을 시작하면 전력 사용량이 증가하고, 밤에 사람들이 잠에 들면 전력 사용량이 감소합니다. 이러한 경우, AR 모델은 이전 시간대의 전력 사용 패턴을 학습하여 다음 시간대의 전력 사용량을 예측하는 데 사용 될 수 있습니다.

MA 모델

• MA 모델은 ‘이동 평균 모델’ 로, 과거의 오차항에 의존하여 현재 값을 예측합니다.
• 예를 들어, 우리가 날씨를 예측한다고 가정해 봅시다. 오늘의 기온 예측치가 실제 기온보다 더 높았다면, 이는 오차를 의미합니다. MA 모델에 따르면, 이 오차는 내일의 기온 예측에 반영되어야 합니다. 즉, 오늘 예측이 너무 높았다면, 내일의 예측치는 그 오차를 보정하기 위해 조금 더 낮아질 것입니다.

MA 모델은 파라미터 q로 이전 몇 개의 오차항을 고려할지 결정하게 됩니다. MA(1) 모델은 바로 이전의 한 단계 오차항 을 고려하고, MA(2) 모델은 바로 이전의 두 단계 오차항 을 고려하게 됩니다.

즉, AR 모델은 자신의 과거 값을, MA 모델은 과거의 예측 오차를 사용하여 현재 값을 예측하는 방법론입니다. 어떤 모델을 사용할지는 주로 데이터의 특성과 목표에 따라 결정됩니다.

(예시)

매장의 판매량 예측을 예로 들어 보겠습니다. 이전의 판매 예측이 과소평가되었다면 (즉, 실제 판매량이 예측보다 높았다면), 그 다음날의 판매량 예측을 증가시키는 것이 합리적일 수 있습니다. 반대로, 이전의 판매 예측이 과대평가 되었다면 (즉, 실제 판매량이 예측보다 낮았다면), 그 다음날의 판매량 예측을 감소시키는 것이 합리적 일 수 있습니다. 이런 상황에서는 MA 모델이 유용하게 사용될 수 있습니다.

ARMA, ARIMA

• ARMA (AutoRegressive Moving Average): AR 모델과 MA 모델의 결합으로, 이전 관측치와 이전 오차항 모두를 고려하는 모델입니다. ARMA 모델은 두 파라미터 p와 q를 모두 사용합니다.

• ARIMA (AutoRegressive Integrated Moving Average): ARMA 모델에 비정상 시계열을 정상 시계열로 변환하는 과정을 포함한 모델입니다. 이 변환 과정은 차분(differencing)이라고 하며, 몇 차 차분을 사용할지를 결정하는 파라미터가 d입니다. 따라서 ARIMA 모델은 세 파라미터 p, d, q를 모두 사용합니다.

ARIMA 적용시 고려사항

• AR(p) = ARIMA(p, 0, 0)

• MA(q) = ARIMA(0, 0, q)

• ARMA(p, q) = ARIMA(p, 0, q)

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import itertools
from tqdm import tqdm
import warnings

warnings.filterwarnings('ignore')

p = range(0, 4)
d = range(0, 3)
q = range(0, 5)

pdq = list(itertools.product(p,d,q))

aic = []
params = []

with tqdm(total=len(pdq)) as pg:
    for i in pdq:
        pg.update(1)
        try:
            model = ARIMA(data, order=(i)).fit()
            aic.append(round(model.aic, 2))
            params.append((i))
        except:
            continue 

100%|███████████████████████████████████████████| 60/60 [00:22<00:00,  2.66it/s]

가장 Optimal 한 (p, d, q) 하이퍼 파라미터를 찾습니다.

optimal = [(params[i],j) for i,j in enumerate(aic) if j == min(aic)]
arima = ARIMA(data, order = optimal[0][0], freq='D').fit()
arima.summary()

SARIMAX Results

Dep. Variable:	Close	No. Observations:	1096
Model:	ARIMA(2, 2, 3)	Log Likelihood	-9369.370
Date:	Mon, 03 Jul 2023	AIC	18750.740
Time:	01:46:18	BIC	18780.726
Sample:	01-01-2020	HQIC	18762.087
	- 12-31-2022
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
ar.L1	-1.9073	0.014	-137.751	0.000	-1.934	-1.880
ar.L2	-0.9732	0.014	-71.464	0.000	-1.000	-0.946
ma.L1	0.8959	0.030	29.584	0.000	0.837	0.955
ma.L2	-0.9402	0.046	-20.326	0.000	-1.031	-0.850
ma.L3	-0.9557	0.028	-34.322	0.000	-1.010	-0.901
sigma2	1.637e+06	8.1e-08	2.02e+13	0.000	1.64e+06	1.64e+06

Ljung-Box (L1) (Q):	0.03	Jarque-Bera (JB):	1219.10
Prob(Q):	0.86	Prob(JB):	0.00
Heteroskedasticity (H):	5.39	Skew:	-0.19
Prob(H) (two-sided):	0.00	Kurtosis:	8.16

arima_model = ARIMA(data, order=optimal[0][0], freq='D').fit()
forecast = arima_model.forecast(steps=10)

fig, ax = plt.subplots(1, 1)
fig.set_size_inches(10, 5)
fig.