(dacon) 농산물 가격 예측 (시계열 분석)

 

 

 

dacon 농산물 가격 예측 튜토리얼

https://dacon.io/edu/21025

농산물 가격 예측 프로젝트

 

개요

프로젝트의 초점은 2021년도에 열렸던 농산물 가격 예측 대회의 데이터를 기반으로 시계열 데이터 분석 및 시계열 예측의 기본기를 다질 수 있도록 구성되어 있습니다.

목표

1. 시계열 데이터 이해와 분석: 시계열 데이터의 기본 구조와 특성, 그리고 주요 분석 기법(e.g., 이동평균, 지수이동평균, 계절성 분해 등)을 이해하고 적용합니다.

2. 데이터 전처리 및 정상성: 누락된 데이터, 이상치, 0 값 등을 처리하고, 시계열 데이터의 정상성을 검정합니다.

3. 기초부터 고급 예측 모델까지: Naive Forecasting으로 시작하여, ARIMA와 LightGBM 같은 고급 예측 모델까지 이해하고 적용합니다.

4. 파라미터 튜닝과 성능 평가: ACF, PACF 그래프 및 Optuna를 활용하여 모델의 성능을 최적화하고, NMAE 같은 성능 지표로 평가합니다.

5. Feature Engineering과 모델 고도화: 다양한 Feature Engineering 기법을 활용하여 모델의 성능을 향상시킵니다.

6. 실전 적용과 결과 해석: 학습한 모델을 실제 데이터에 적용하여 예측 결과를 생성하고, 이를 시각적으로 분석 및 해석합니다.

이러한 학습 목표를 통해 참가자들은 시계열 데이터 분석의 전반적인 과정을 체계적으로 배우고, 실제 문제에 적용하는 능력을 키울 수 있습니다.

설명

본 프로젝트의 목표는 농산물 가격 예측을 위한 시계열 데이터 분석과 모델링을 주제로 합니다. 총 6개의 스테이지로 구성되어 있으며, 각 스테이지는 다양한 학습 목표와 단계를 가지고 있습니다. 이 교재를 통해 학습자는 시계열 데이터의 특성을 이해하고, 이를 분석, 전처리, 모델링하는 전반적인 과정을 체계적으로 배울 수 있습니다.

 

 

 

데이터 검토

데이터 기간

• train: 2016/1/1 ~ 2021/9/27
• test: 2021/9/28 ~ 2021/11/4 (훈련 데이터 직후 기간)

컬럼

• date
• 배추_거래량(kg)
• 배추_가격(원/kg)
• 무_거래량(kg)
• 무_가격(원/kg)
• 깻잎_거래량(kg)
• 깻잎_가격(원/kg)
• 시금치_거래량(kg)
• 시금치_가격(원/kg)
• 토마토_거래량(kg)
• 토마토_가격(원/kg)
• 청상추_거래량(kg)
• 청상추_가격(원/kg)
• 캠벨얼리_거래량(kg)
• 캠벨얼리_가격(원/kg)
• 샤인마스캇_거래량(kg)
• 샤인마스캇_가격(원/kg)

 

Naive Forecasting

단순히 가장 최근 데이터로 미래 값을 예측하는 방법 ;;;

2021-09-28 데이터를 +1w, +2w, +4w 예측에 그대로 사용

 

누락 데이터 채우기

시계열 데이터의 누락 데이터는 이전 데이터를 복사해서 채우기도 한다.

DataFrame.replace에 ffill method 인자를 사용하면 된다.

튜토리얼에서는 아래와 같이 0으로 입력된 데이터를 replace하고 있다.

train_prep.replace(0, method="ffill", inplace=True)

 

그런데, 매뉴얼에 들어가보니 2.1.0 버전부터 method 인자가 deprecated되었다.

replace 대신 ffill 함수를 사용하면 되는데, ffill은 None만 치환할 수 있고, 치환할 value를 직접 설정할 수 없었다.

튜토리얼처럼 0을 치환하고 싶을 때에는 아래와 같이 0을 None으로 바꿔주는 과정을 거쳐야 한다.

train_prep.replace(0, None, inplace=True)  # 0을 None으로 치환
train_prep.ffill(inplace=True)  # None을 이전 데이터로 채우기
train_prep.fillna(0, inplace=True)  # 시작일 데이터가 None으로 남아있으면 다시 0으로 치환

 

원본 데이터 (2016-01-03 데이터 누락)

2016-01-03 데이터를 기존 값으로 채움

 

이동평균(MA), 지수이동평균(EMA)

DataFrame에서 이동 평균(MA)과 지수 이동 평균(EMA)은 각각 rolling().mean()과 ewm().mean()으로 구할 수 있다.

rolling()의 window 인자는 윈도우 크기, 즉, 평균을 구할 기간이다.

ewm()의 span 인자는 지수이동평균 계산에 사용할 기간 인자이다.

# 다양한 윈도우 크기로 MA와 EMA 계산
window_size = 7

train_selected[f'{item_selected}_가격(원/kg)_MA'] = train_selected[f'{item_selected}_가격(원/kg)'].rolling(window=window_size).mean()
train_selected[f'{item_selected}_가격(원/kg)_EMA'] = train_selected[f'{item_selected}_가격(원/kg)'].ewm(span=window_size, adjust=False).mean()

# ... 차트 표시

 

ewm을 매뉴얼(DataFrame.ewm, Exponentially weighted window)에서 찾아보니 "span은 지수이동평균 관련 기간이다"라는 설명은 부족한 것 같다.
지수이동평균은 해당일까지의 기간 전체 데이터를 사용하는데, 얼마나 최근 데이터에 가중치를 줄 것인지에 대한 파라미터라고 생각하고 일단 넘어간다.

 

계절성 분해 (Seasonal Decomposition)

seasonal_decompose 함수를 이용해서 시계열 데이터를 추세(trend), 계절성(seasonal), 잔차(residual)로 분해해 볼 수 있다.

from statsmodels.tsa.seasonal import seasonal_decompose

item_selected = '배추'
start_date_train = '2021-01'

# 배추 가격 추출
train_selected = train_prep[train_prep['date'] >= start_date_train].copy()
train_selected_item = train_selected[['date', f'{item_selected}_가격(원/kg)']].copy()

# 계절성 분해 수행 (연별 계절성 가정)
decomposition = seasonal_decompose(train_selected_item[f'{item_selected}_가격(원/kg)'], model='additive', period=30)

# ... 차트 표시

 

 

period 인자에 30을 사용해서 월간 seasonality가 계산되었다. 배추 가격은 연 단위로 봐야하겠지만 데이터가 부족해서 예제에서는 30을 사용한 듯하다. (period=365로 바꿔봤더니 데이터가 cycle의 두 배인 730개 이상이어야 하는데 부족하다고 오류 발생)

 

seasonal_decompose가 실제로도 많이 쓰이는지 궁금해서 검색해보았다.
아직 쓰는 곳이 많지만, 오래된 방법이므로 이후에 나온 X11, SEATS, STL 등을 쓰는 것이 좋다고 한다.
(https://yoongaemii.github.io/seasonal_decomposition/)

 

ARIMA(AutoRegressive Integrated Moving Average) 모델

시계열 데이터 예측 모델 중 하나인 ARIMA 모델에 대해 알아본다.

 

- ARIMA 모델 참고: https://siddp6.medium.com/arima-model-424a1b0f0970

ARIMA Model

 

AR(AutoRegressive) 모델

과거 데이터로 미래를 예측하는 시계열 모델이다.

예를 들어, 오늘 날씨가 맑았다면, 내일도 맑을 것이라고 예측하는 방식이다.

AR(p) 모델의 p 파라미터는 PACF(Partial AutoCorrelation Function, 부분 자기상관 함수)로 구하는데, PACF 값이 급격히 0으로 떨어지는 시차(lag)를 찾아서 직전 값을 사용한다.

PACF는 시점 t와 t-n 사이의 관계를 중간의 t-1 ~ t-(n-1) 데이터를 무시하고 상관관계를 구하는 함수이다.

 

Integrated

ARIMA 모델은 stationary한 데이터에서 동작한다. stationary하다는 것은 평균과 분포(표준 편차)가 시간이 지나도 일정하다는 뜻이다.

그런데, 가격 데이터는 점차 증가하는 경향이 있으므로 stationary하지 않다.

이러한 데이터를 stationary하게 만들기 위해 데이터를 차분하는 방법을 사용한다.

 

MA(Moving Average, 이동 평균)

이동 평균 MA는 단기 추세를 확인하는데 도움을 준다.

MA(q) 모델의 q 파라미터는 ACF(AutoCorrelation Function, 자기상관 함수)로 구하며, ACF 값이 급격히 0으로 떨어지는 시차(lag)를 찾아서 직전 값을 사용한다.

ACF는 한 시점의 값이 이전 여러 시점들의 값과 어떤 상관관계를 가지는지 측정하는 함수이다.

 

AR(p), MA(q)의 p, q 파라미터를 찾기 위해 우선 ACF, PACF 그래프를 그려 본다.

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

n_diff = 1
features = train_selected.columns[1:]
train_selected_diff = train_selected[features].diff(n_diff).dropna()

fig, ax = plt.subplots(1,2 ,figsize = (12, 4))

# ACF (AutoCorrelation Function) 그래프 그리기
plot_acf(train_selected_diff[f'{item_selected}_가격'], lags=50, title=f'ACF ({item_selected}_가격)', ax=ax[0])

# PACF (Partial AutoCorrelation Function) 그래프 그리기
plot_pacf(train_selected_diff[f'{item_selected}_가격'], lags=50, title=f'PACF ({item_selected}_가격)', ax=ax[1])

plt.tight_layout()

 

그리고, (95%) 신뢰구간을 벗어나는 ACF, PACF 값을 확인해본다.

from statsmodels.tsa.stattools import acf, pacf

acf_values = acf(train_selected_diff[f'{item_selected}_가격'], nlags=50)
pacf_values = pacf(train_selected_diff[f'{item_selected}_가격'], nlags=50)

# 신뢰 구간 계산 (95% 신뢰 구간)
confidence_interval = 1.96 / np.sqrt(len(train_selected_diff))

q_param = []
for lag, value in enumerate(acf_values):
    if np.abs(value) > confidence_interval:
        q_param.append(lag)

p_param = []
for lag, value in enumerate(pacf_values):
    if np.abs(value) > confidence_interval:
        p_param.append(lag)


display(f"Q parameters of MA model is {q_param},  P of AR model is {p_param}")

'Q parameters of MA model is [0, 1, 7, 14, 25, 29], P of AR model is [0, 1, 6, 25, 29, 31, 44]'

 

p, q 파라미터를 각각 7, 6으로 선정하여 ARIMA 모델을 학습하는 코드 및 summary 내용은 다음과 같다.

from statsmodels.tsa.arima.model import ARIMA

start_date_valid = '2021-05-01'

_train_item = train_1week[ train_1week['date'] < start_date_valid].reset_index(drop=True)
# _valid_item = train_1week[ train_1week['date'] >= start_date_valid].reset_index(drop=True)

train_item_arima = _train_item.copy()
train_item_arima.set_index('date', inplace = True)

model_arima = ARIMA(train_item_arima[f"{item_selected}_가격"], order = (7,1,6))

fit_result = model_arima.fit()
display(fit_result.summary())

 

위에서 AIC와 BIC가 모델 복잡성 및 적합도를 동시에 고려하는 지표이며, 낮을수록 모델이 데이터에 잘 적합되었다고 볼 수 있다.

 

p, q 파라미터는 다음과 같이 grid search 방식으로 찾을 수도 있다.

import itertools
import warnings
warnings.filterwarnings('ignore')

# p, d, q 범위 설정
p = range(0, 7)
d = range(1, 2)
q = range(0, 3)

# 가능한 모든 p, d, q 조합 생성
pdq = list(itertools.product(p, d, q))

# AIC 값을 저장할 리스트 초기화
aic = []

train_item_arima = _train_item.copy()
train_item_arima.set_index('date', inplace = True)


# 각 조합에 대해 ARIMA 모델 적합 및 AIC 계산
for i in pdq:
    try:
        # 모델 생성 및 적합
        model = ARIMA(train_item_arima[f"{item_selected}_가격"], order=i)  # 빈도를 'B' (영업일)로 설정
        model_fit = model.fit()

        # AIC 값 출력 및 저장
        print(f'ARIMA : {i} >> AIC : {round(model_fit.aic, 2)}')
        aic.append(round(model_fit.aic, 2))
    except:
        continue

optimal_pdq_info = [ (pdq[i], j) for i, j in enumerate(aic) if j == min(aic)]
optimal_pdq = optimal_pdq_info[0][0]
display(f"optimal pdq parameter : {optimal_pdq_info}")

ARIMA : (0, 1, 0) >> AIC : 1396.25
ARIMA : (0, 1, 1) >> AIC : 1397.45
ARIMA : (0, 1, 2) >> AIC : 1378.85
ARIMA : (1, 1, 0) >> AIC : 1397.98
ARIMA : (1, 1, 1) >> AIC : 1383.4
ARIMA : (1, 1, 2) >> AIC : 1378.5
ARIMA : (2, 1, 0) >> AIC : 1387.26
ARIMA : (2, 1, 1) >> AIC : 1378.86
ARIMA : (2, 1, 2) >> AIC : 1380.14
ARIMA : (3, 1, 0) >> AIC : 1388.66
ARIMA : (3, 1, 1) >> AIC : 1380.42
ARIMA : (3, 1, 2) >> AIC : 1382.1
ARIMA : (4, 1, 0) >> AIC : 1389.76
ARIMA : (4, 1, 1) >> AIC : 1379.7
ARIMA : (4, 1, 2) >> AIC : 1383.49
ARIMA : (5, 1, 0) >> AIC : 1360.47
ARIMA : (5, 1, 1) >> AIC : 1361.47
ARIMA : (5, 1, 2) >> AIC : 1356.71
ARIMA : (6, 1, 0) >> AIC : 1360.49
ARIMA : (6, 1, 1) >> AIC : 1359.88
ARIMA : (6, 1, 2) >> AIC : 1355.34
'optimal pdq parameter : [((6, 1, 2), 1355.34)]'

 

ARIMA 모델로 예측한 1, 2, 4주 후 가격과 실제 가격 비교 코드 및 그래프는 다음과 같다.

import matplotlib.pyplot as plt
import pandas as pd


def plot_predict_behavior(df_all, df_valid_week1, df_valid_week2, df_valid_week4, pred_1week, pred_2week, pred_4week, item_selected):

    predicted_1week_df = pd.DataFrame({'date_1week': list(df_valid_week1['date_1week']), 'pred_1week': pred_1week}).dropna()
    predicted_2week_df = pd.DataFrame({'date_2week': list(df_valid_week2['date_2week']), 'pred_2week': pred_2week}).dropna()
    predicted_4week_df = pd.DataFrame({'date_4week': list(df_valid_week4['date_4week']), 'pred_4week': pred_4week}).dropna()

    start_date_valid = df_valid_week1['date'].min()
    end_date_valid = df_valid_week1['date'].max()

    plt.figure(figsize=(12, 6))

    # 전체 기간 그래프
    plt.subplot(2, 1, 1)
    plt.plot(df_all['date'], df_all[f'{item_selected}_가격'], label='Real', color='blue')
    plt.plot(predicted_1week_df['date_1week'], predicted_1week_df['pred_1week'], label='Predicted 1week', color='orange', linestyle='--')
    plt.plot(predicted_2week_df['date_2week'], predicted_2week_df['pred_2week'], label='Predicted 2week', color='green', linestyle='--')
    plt.plot(predicted_4week_df['date_4week'], predicted_4week_df['pred_4week'], label='Predicted 4week', color='black', linestyle='--')
    plt.axvline(x=start_date_valid, color='r', linestyle='--', label='Start of Validation')
    plt.axvline(x=end_date_valid, color='r', linestyle='--', label='End of Validation')
    plt.title('Real vs Predicted Prices (From 2021)')
    plt.legend()

    y_min_zoomed_graph = min(df_valid_week1[f'{item_selected}_가격'].min(),
                             predicted_1week_df['pred_1week'].min(),
                             predicted_2week_df['pred_2week'].min(),
                             predicted_4week_df['pred_4week'].min())

    y_max_zoomed_graph = max(df_valid_week1[f'{item_selected}_가격'].max(),
                             predicted_1week_df['pred_1week'].max(),
                             predicted_2week_df['pred_2week'].max(),
                             predicted_4week_df['pred_4week'].max())

    # 검증 기간 확대 그래프
    plt.subplot(2, 1, 2)
    plt.plot(df_all['date'], df_all[f'{item_selected}_가격'], label='Real', color='blue')
    plt.plot(predicted_1week_df['date_1week'], predicted_1week_df['pred_1week'], label='Predicted 1week', color='orange', linestyle='--')
    plt.plot(predicted_2week_df['date_2week'], predicted_2week_df['pred_2week'], label='Predicted 2week', color='green', linestyle='--')
    plt.plot(predicted_4week_df['date_4week'], predicted_4week_df['pred_4week'], label='Predicted 4week', color='black', linestyle='--')
    plt.axvline(x=start_date_valid, color='r', linestyle='--', label='Start of Validation')
    plt.axvline(x=end_date_valid, color='r', linestyle='--', label='End of Validation')
    plt.xlim(left=start_date_valid - pd.DateOffset(days=14))
    plt.ylim(top=y_max_zoomed_graph, bottom=y_min_zoomed_graph)
    plt.title('Real vs Predicted Prices (Validation Period)')
    plt.legend()

    plt.tight_layout()
    plt.show()

plot_predict_behavior(train_selected, _valid_data_1week, _valid_data_2week, _valid_data_4week, pred_1week, pred_2week, pred_4week, item_selected)

 

 

LightGBM을 사용한 예측 모델

통계적 접근법이 아닌 decision tree 방식의 예측 모델 예제.

모델을 학습한 후 1, 2, 4주 후의 가격을 예측한 그래프는 아래와 같다.

단순 예제이기는 하지만 여기서는 주목할 만한 정확도는 보여주지 않는 듯하다.

 

LightGBM 시계열 예측 정확도 향상을 위한 피쳐 엔지니어링

1. 날짜로부터 월, 주차, 요일 피쳐 추가

# 날짜 관련 feature 추가 하기

train_selected['month'] = train_selected['date'].dt.month
train_selected['week'] = train_selected['date'].dt.isocalendar().week.astype(np.int32)
train_selected['weekday']  = train_selected['date'].dt.weekday

# 날자 관련 피처를 저장해 둔다.
features_date = ['month', 'week', 'weekday']

display(train_selected.columns)

 

2. 상관관계가 높은 다른 품목의 과거 1주일 변동량 피쳐 추가

상관관계가 높은 다른 품목(피쳐)의 과거 1주일 변동량 데이터가 모델 예측 정확도를 높이는데 도움이 될 수 있다.

 

3. EMA(Exponential Moving Average)를 이용한 피쳐 추가

7, 14, 28일 window로 지수 이동평균(EMA)을 계산해서, 실제 가격과의 오차율을 피쳐로 추가한다.