Data Source:

https://www.kaggle.com/datasets/robertritz/ub-meat-prices?resource=download

Time Series¶

시간 데이터 : 특정 간격을 갖는 시간이 필요합니다.(Lag: 시간의 간격)

Lag데이터 에서 Data의 경향성을 찾아야 합니다.

일변량 정상 시계열 - ARIMA
시계열 성분 특징

① Trend (추세)
② Seasonality (계절성)
③ Cycle (주기)
④ Noise (잡음) → White Noise 평균과 분산이 변하지 않아 분석가능한 데이터

시계열 분석은 외부적인 변동 요인에 취약하기에 모델이 불확실할 수 있다.

In [1]:

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt

mpl.rc('font', family = 'Malgun Gothic')

In [2]:

df1 = pd.read_csv('meat_prices_20180103_20211027.csv')
df1.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 796 entries, 0 to 795
Data columns (total 3 columns):
 #   Column     Non-Null Count  Dtype  
---  ------     --------------  -----  
 0   meat_type  796 non-null    object 
 1   date       796 non-null    object 
 2   price      796 non-null    float64
dtypes: float64(1), object(2)
memory usage: 18.8+ KB

In [3]:

df1.head()

Out[3]:

	meat_type	date	price
0	Mutton, kg	2018-01-03	5090.0
1	Beef, kg	2018-01-03	6470.0
2	Horse meat, kg	2018-01-03	4900.0
3	Goat meat, kg	2018-01-03	3435.0
4	Mutton, kg	2018-01-10	5100.0

In [4]:

df1.tail()

Out[4]:

	meat_type	date	price
791	Goat meat, kg	2021-10-20	9000.0
792	Mutton, kg	2021-10-27	10974.0
793	Beef, kg	2021-10-27	13322.0
794	Horse meat, kg	2021-10-27	9937.0
795	Goat meat, kg	2021-10-27	8986.0

In [5]:

# 날짜 데이터 날짜 타입으로 변경
df1['Datetime'] = pd.to_datetime(df1['date'])

In [6]:

df1['Datetime'].describe(datetime_is_numeric=True)

Out[6]:

count                              796
mean     2019-12-02 00:07:14.170854144
min                2018-01-03 00:00:00
25%                2018-12-12 00:00:00
50%                2019-12-04 00:00:00
75%                2020-11-18 00:00:00
max                2021-10-27 00:00:00
Name: Datetime, dtype: object

In [7]:

df1['Year'] = df1['Datetime'].dt.year
df1['Month'] = df1['Datetime'].dt.month
df1['Week'] = df1['Datetime'].dt.isocalendar().week
df1['Day_of_week'] = df1['Datetime'].dt.day_name()

In [8]:

df1.head()

Out[8]:

	meat_type	date	price	Datetime	Year	Month	Week	Day_of_week
0	Mutton, kg	2018-01-03	5090.0	2018-01-03	2018	1	1	Wednesday
1	Beef, kg	2018-01-03	6470.0	2018-01-03	2018	1	1	Wednesday
2	Horse meat, kg	2018-01-03	4900.0	2018-01-03	2018	1	1	Wednesday
3	Goat meat, kg	2018-01-03	3435.0	2018-01-03	2018	1	1	Wednesday
4	Mutton, kg	2018-01-10	5100.0	2018-01-10	2018	1	2	Wednesday

추세 확인 - lineplot or pointplot

합을 계산해서 확인해봅니다.

In [9]:

sns.lineplot(data = df1, x = 'Datetime', y = 'price', estimator = sum)

Out[9]:

<AxesSubplot:xlabel='Datetime', ylabel='price'>

위 방식은 데이터가 방대해지면 데이터를 비교하기에 시간이 오래걸립니다.

In [10]:

pivot_df1 = pd.pivot_table(data = df1, index = "Datetime",
               values = 'price', aggfunc = 'sum').reset_index()

In [11]:

sns.lineplot(data = pivot_df1, x= 'Datetime', y='price')

Out[11]:

<AxesSubplot:xlabel='Datetime', ylabel='price'>

In [12]:

pivot_df2 = pd.pivot_table(data = df1, index = ["Datetime","meat_type"],
               values = 'price', aggfunc = 'sum').reset_index()
sns.lineplot(data = pivot_df2, x= 'Datetime', y='price', hue= 'meat_type')

Out[12]:

<AxesSubplot:xlabel='Datetime', ylabel='price'>

시계열 분석 목적

이전 시간에 따라서 나온 데이터가 다음 시간에 어떤 영향을 줄까를 파악하는 분석

외부적인 작용이 없는 시간에 있어 독립적인 데이터만 모아서 분석해야합니다.

위 데이터에서는 눈에 띄는 데이터는 없어보입니다.

In [13]:

# 특정 데이터만 뽑아서 보기
cond1 = (df1['meat_type'] == 'Beef, kg')
beef_df = df1.loc[cond1]

ARIMA (Auto Regressive Integrated Moving Average Model)¶

AR (Auto Regressive model) :

자기회귀모델 - 특정 시점(p) 전의 자료가 현재 시점의 데이터에 영향을 주는 자기 회귀 모델

특정 간격 내 불특정한 변동된 자료에 대한 설명력이 떨어집니다.

ACF (자기상관함수) auto correlation fuction / 특정 구간 내 데이터 간 상관관계 표현 (외부적인 요인으로 두 개의 독립적인 변수가 상관관계를 보임)

PACF (부분자기상관함수) partial auto correlation fuction / 연속적인 연관성을 배제하고 부분적인 데이터간의 관계를 표현

MA (Moving Average model) :

이동평균모델 - 일정한 구간 데이터의 평균을 계산하여 미래를 예측하는 모델

특정 구간 내 데이터의 평균을 이용하므로 불규칙적인 변동을 제거해줍니다.

Difference 차분 :

기존에는 시간의 변화에 있어서 평균이나 분산이 일정한 데이터를 다루었습니다.

실무에서는 외부적인 요인이 많기에 정상성을 만족하지 않은 경우가 많습니다.

차분을 이용해 모델이 이해가능하도록 변화 시켜줍니다.

ARIMA model

단기 예측에 적합합니다.
계절적 변동요인 (주기적 변동 not random)이 있는 경우에 사용합니다.
Sample > 50 이상이여야 사용가능하다
정성적 자료 (일정한 분산)

not 정성적 → 차분

In [14]:

import statsmodels.tsa.api as tsa # 통계 라이브러리

In [15]:

beef_df

Out[15]:

	meat_type	date	price	Datetime	Year	Month	Week	Day_of_week
1	Beef, kg	2018-01-03	6470.0	2018-01-03	2018	1	1	Wednesday
5	Beef, kg	2018-01-10	6475.0	2018-01-10	2018	1	2	Wednesday
9	Beef, kg	2018-01-17	6525.0	2018-01-17	2018	1	3	Wednesday
13	Beef, kg	2018-01-24	6525.0	2018-01-24	2018	1	4	Wednesday
17	Beef, kg	2018-01-31	7075.0	2018-01-31	2018	1	5	Wednesday
...	...	...	...	...	...	...	...	...
777	Beef, kg	2021-09-29	13367.0	2021-09-29	2021	9	39	Wednesday
781	Beef, kg	2021-10-06	13304.0	2021-10-06	2021	10	40	Wednesday
785	Beef, kg	2021-10-13	13405.0	2021-10-13	2021	10	41	Wednesday
789	Beef, kg	2021-10-20	13081.0	2021-10-20	2021	10	42	Wednesday
793	Beef, kg	2021-10-27	13322.0	2021-10-27	2021	10	43	Wednesday

199 rows × 8 columns

In [16]:

# 시계열 - 인덱스 (시간), 하나의 feature (분석할 대상)
# 만약, 인덱스가 시간이 아닌 경우 set_index('Datetime') 해주어야 합니다.
df_time = pd.pivot_table(data = beef_df, index = 'Datetime', values = 'price', aggfunc = 'sum')
df_time.head()

Out[16]:

	price
Datetime
2018-01-03	6470.0
2018-01-10	6475.0
2018-01-17	6525.0
2018-01-24	6525.0
2018-01-31	7075.0

In [17]:

# 하루 단위
y = df_time['price'].resample('1w').mean()

In [18]:

y.isnull().sum()

Out[18]:

In [19]:

y1 = y.fillna(method= 'ffill')

시계열 분해

Trend (추세요인), seasonality(계절요인), Residual (불규칙 or 순환요인)

In [20]:

from pylab import rcParams
rcParams['figure.figsize'] = 15, 10

#차트 기본 크기 설정
mpl.rcParams['axes.labelsize'] = 15
mpl.rcParams['xtick.labelsize'] = 12
mpl.rcParams['ytick.labelsize'] = 12
mpl.rcParams['text.color'] = 'k'

In [21]:

# 시계열 모델 생성
model_series = tsa.seasonal_decompose(y1, model = 'additive')
# 모델 시각화
fig = model_series.plot()
plt.show()

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\backends\backend_agg.py:240: RuntimeWarning: Glyph 8722 missing from current font.
  font.set_text(s, 0.0, flags=flags)
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\backends\backend_agg.py:203: RuntimeWarning: Glyph 8722 missing from current font.
  font.set_text(s, 0, flags=flags)

In [22]:

import itertools # 반복수를 만드는 라이브러리

In [23]:

p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

In [24]:

seasonal_pdq

Out[24]:

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

In [25]:

param_list = []
param_seasonal_list = []
results_AIC_list = []

In [34]:

for param in pdq:
     for param_seasonal in seasonal_pdq:
        try:
            mod = tsa.statespace.SARIMAX(y1,order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False)
            results = mod.fit()
            param_list.append(param)
            param_seasonal_list.append(param_seasonal)
            results_AIC_list.append(results.aic)
        except:
             continue

C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\base\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

In [27]:

ARIMA_list = pd.DataFrame({'Parameter':param_list,'Seasonal':param_seasonal_list,'AIC':results_AIC_list})
ARIMA_list.to_excel('arima_model_list.xlsx')

In [35]:

ARIMA_list.sort_values(by='AIC').head() # 낮을수록 좋은 모델

Out[35]:

	Parameter	Seasonal	AIC
27	(0, 1, 1)	(0, 1, 1, 12)	2510.739906
59	(1, 1, 1)	(0, 1, 1, 12)	2511.992587
31	(0, 1, 1)	(1, 1, 1, 12)	2512.299133
63	(1, 1, 1)	(1, 1, 1, 12)	2513.453048
19	(0, 1, 0)	(0, 1, 1, 12)	2523.460945

Likelihood (최대 우도) - 특정 데이터가 모수로부터 추출되었을 가능도

즉, 특정 값이 집단에서의 분포에서 속할 확률 추정 (연속 확률 밀도 함수 pdf의 y값)

AIC (Akaike Information Criterion) : 데이터에 대한 모델의 상대적 품질 수준 수치화
- AIC = -2 ln(L최대우도) + 2k
- 값이 낮을 수록 모형 적합도가 높습니다.
BIC (Bayes Infomation Criterion)
- BIC = =2 ln(L) + log(n)p
- 변수가 더 많을 때, AIC에 더 많은 패널티를 부여해 계산합니다.
- 값이 낮을 수록 모형 적합도가 높습니다.
HQIC (Hannan Quinn Information Criterion)
- HQIC = -2 ln(L) + 2kln(ln(n))
- 값이 낮을 수록 모형 적합도가 높습니다.

In [36]:

model = tsa.statespace.SARIMAX(y1,order=(0, 1, 1),seasonal_order=(0, 1, 1, 12),
                              enforce_stationarity=False, enforce_invertibility=False)
results = model.fit()
print(results.summary())

                                     SARIMAX Results                                      
==========================================================================================
Dep. Variable:                              price   No. Observations:                  200
Model:             SARIMAX(0, 1, 1)x(0, 1, 1, 12)   Log Likelihood               -1252.370
Date:                            Mon, 30 Jan 2023   AIC                           2510.740
Time:                                    00:24:03   BIC                           2520.200
Sample:                                01-07-2018   HQIC                          2514.578
                                     - 10-31-2021                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ma.L1          0.0824      0.082      1.001      0.317      -0.079       0.244
ma.S.L12      -0.9998      0.043    -23.042      0.000      -1.085      -0.915
sigma2      9.786e+04   4.43e-07   2.21e+11      0.000    9.79e+04    9.79e+04
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):               901.47
Prob(Q):                              0.97   Prob(JB):                         0.00
Heteroskedasticity (H):               0.69   Skew:                             1.66
Prob(H) (two-sided):                  0.16   Kurtosis:                        13.68
===================================================================================

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

SARIMA는 Order (p,d,q) * Seasonal Order (P,D,Q,M)
- p : AR에서의 p값 (p 구간 내 데이터 사이의 상관관계)
- d : 차분
- q : MA모델 PACF 편상관계수 q 값
- P : 계절성 주기에서의 패턴
- D : 계절성 여부
- Q : 주기의 패턴
- M : 계절성 주기
Ljung - Box Test : 일정 기간동안 관측치가 랜덤이고, 독립적인지 여부를 검정
- 귀무 : 데이터가 상관관계를 나타내지 않는다.
- 대립 : 데이터가 상관관계를 나타낸다.
- P.value(귀무가설이 참일 확률) < 0.05 (유의수준) 일 경우 대립가설 참
- 위 결과 데이터가 상관관계를 나타내지 않는다.
Jarque Bera Test : 왜도와 첨도가 정규분포와 일치하는지 가설검정
- SARIMAX : 잔차의 분포가 정규분포 인가
- 귀무 가설 : 해당 잔차(residual)는 정규분포와 일치한다.
- 대립 가설 : 해당 잔차(residual)는 정규분포와 일치하지 않는다.
- P.value < 0.05 , 해당 잔차(residual)는 정규분포와 일치하지 않는다.
Heteroskedasticity : 분산의 이분산 검정
- 잔차의 분산이 실제의 분산과 같으냐 다르냐
Warnings :
- 모수의 추정치가 불안전하다 (오차의 변동이 심함)

예측

In [30]:

results.plot_diagnostics(figsize=(16, 8))
plt.show()

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\backends\backend_agg.py:240: RuntimeWarning: Glyph 8722 missing from current font.
  font.set_text(s, 0.0, flags=flags)
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\backends\backend_agg.py:203: RuntimeWarning: Glyph 8722 missing from current font.
  font.set_text(s, 0, flags=flags)

In [37]:

results.get_prediction()

Out[37]:

<statsmodels.tsa.statespace.mlemodel.PredictionResultsWrapper at 0x23e6a746b20>

In [47]:

y1.head()

Out[47]:

Datetime
2018-01-07    6470.0
2018-01-14    6475.0
2018-01-21    6525.0
2018-01-28    6525.0
2018-02-04    7075.0
Freq: W-SUN, Name: price, dtype: float64

In [46]:

pred = results.get_prediction(start=pd.to_datetime('2018-07-29'), dynamic=False)
pred_ci = pred.conf_int()

ax = y1.plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7, figsize=(14, 7))
ax.fill_between(pred_ci.index,pred_ci.iloc[:, 0],pred_ci.iloc[:, 1], color='k', alpha=.2)

plt.legend()
plt.show()

In [49]:

pred = results.get_prediction(start=pd.to_datetime('2018-07-29'), dynamic=False)
pd.DataFrame(pred.predicted_mean).reset_index()

Out[49]:

	Datetime	predicted_mean
0	2018-07-29	9359.863615
1	2018-08-05	8634.595815
2	2018-08-12	8699.159155
3	2018-08-19	8667.909590
4	2018-08-26	7955.054619
...	...	...
166	2021-10-03	13512.735065
167	2021-10-10	13371.699994
168	2021-10-17	13317.223676
169	2021-10-24	13442.343605
170	2021-10-31	13102.780673

171 rows × 2 columns

In [50]:

y_forecasted = pred.predicted_mean
y_truth = y['2018-01-07':]
mse = ((y_forecasted - y_truth) ** 2).mean()
print('MSE {}'.format(round(mse, 2)))

MSE 115737.47

In [51]:

pred_uc = results.get_forecast(steps=50)

In [52]:

pred_uc.predicted_mean.head()

Out[52]:

2021-11-07    13312.710236
2021-11-14    13227.484620
2021-11-21    13280.164685
2021-11-28    13446.366205
2021-12-05    13514.876586
Freq: W-SUN, Name: predicted_mean, dtype: float64

In [54]:

pred_uc = results.get_forecast(steps=50)
pred_ci = pred_uc.conf_int() #추정된 계수의 신뢰구간 계산
ax = y1.plot(label='observed', figsize=(14, 7))
pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
ax.fill_between(pred_ci.index,
                 pred_ci.iloc[:, 0],
                 pred_ci.iloc[:, 1], color='k', alpha=.25)

plt.legend()
plt.show()

Kang's Note

[일변량 정상 시계열 분석]

Time Series¶

ARIMA (Auto Regressive Integrated Moving Average Model)¶

댓글

티스토리툴바