시계열 : 서울시 평균기온 예측 프로젝트 3(ARIMA 모델을 사용한 단변량 예측)
1. 데이터 로드
csv 파일을 읽어올 때, parse_dates 매개변수에 일시 칼럼을 입력하여 일시 칼럼의 자료형을 datetime으로 변경합니다.
import pandas as pd
train = pd.read_csv('train.csv', parse_dates=['일시'])
display(train.head(2))
# 데이터 타입 확인
print(f'일시 칼럼의 데이터타입은 {train["일시"].dtypes} 입니다.')
2. 시계열 데이터 인덱싱 및 학습/검증 세트 분할
# '일시'를 인덱스로 설정하고 빈도를 'D'로 설정
train = train.set_index('일시')
# 학습/ 검증 데이터 분할
train_data = train[train.index < '2022-01-01']
validation_data = train[train.index >= '2022-01-01']
train_data.index.freq = 'D'
validation_data.index.freq = 'D'
※ 코드 해석
예측할 기간이 '2023-01-01'부터 '2023-12-24'까지 약 1년인 것을 고려할 때, 검증 데이터 세트를 동일한 기간 동안 설정하는 것이 중요합니다.
이렇게 하면, 모델이 예측해야 할 실제 기간과 유사한 패턴과 변동성을 학습할 수 있으며, 이는 최종 예측 결과의 정확성 향상시킬 수 있습니다.
시계열이 아닌 프로젝트에서는 일반적으로 사이킷런의 train_test_split() 함수를 사용했습니다.
그런데 시계열 데이터의 경우, 데이터 포인트들 간의 시간적 순서가 중요합니다.
train_test_split() 함수는 데이터를 무작위로 분할합니다. 이는 시계열 데이터의 시간적 연속성을 무시합니다.
시간적 연속성을 무시할 경우 시계열 데이터의 패턴, 특성과 구조를 학습하는 데 부적합합니다.
그러므로, 시계열 분석에서는 특정 시점을 기준으로 데이터를 나눈 후 분할함으로써 모델이 시간에 따른 패턴을 인지하고 예측하는 데 필요한 학습을 지속할 수 있게 해야 합니다.
이는 모델이 과거의 데이터를 기반으로 미래를 예측하는 시계열 분석의 본질에 더 적합합니다.
마지막으로, train_data와 validation_data의 인덱스는 빈도를 일별인 'D'로 설정합니다.
ARIMA 모델의 입력 데이터의 정확한 시간 구조를 제공하여 필요할 경고 메시지를 방지합니다.
3. ARIMA 모델 초기 설정 및 학습
※ ARIMA 모델
ARIMA는 자기회귀 누적 이동 평균(AutoRegressive Integrated Moving Average)의 약자로, 시계열 데이터의 과거 값과 오류를 사용하여 미래 값을 예측하는 모델입니다.
이 모델은 자기회귀(AR), 차분(I), 이동평균(MA) 세 가지의 구성 요소로 이루어져 있습니다.
- 자기회귀(AR) : '자기회귀'는 과거의 값들이 미래 값에 어떤 영향을 미치는지를 설명합니다.예를 들어, 지난주의 판매량이 이번 주의 판매량에 영향을 미치는 경우, 이는 자기회귀 관계에 해당합니다.AR부분에서는 'p'라는 파라미터를 사용하는데, 이는 과거 데이터 포인트 중 얼마만큼 고려할 것인지를 결정합니다.
- 차분(I) : '차분'은 비정상적인 시계열 데이터를 정상적인 상태로 만들기 위해 사용됩니다.비정상 시계열은 평균, 분산 등이 시간에 따라 변하는 특성을 가집니다.이러한 데이터를 정상화하기 위해 차분(d)을 수행하여 일정한 수준을 유지하도록 합니다.'d'라는 파라미터는 몇 번 차분을 수행할 것인지를 나타냅니다.
- 이동평균(MA) : '이동평균' 부분은 과거 예측 오차가 미래 값에 어떻게 영향을 미치는지를 나타냅니다. MA에서는 'q'라는 파라미터를 사용하여 과거 예측 오차 중 몇 개를 고려할 것인지를 결정합니다.
model = ARIMA(train_data['평균기온'], order=(2, 1, 3))
위 코드는 ARIMA 모델을 초기화합니다.
여기서 '평균기온'은 모델에 사용할 시계열 데이터입니다.
order=(2, 1, 3)은 ARIMA 모델의 세 가지 주요 파라미터를 설정합니다.
- 2 → 자기회귀(AR) 부분의 차수(p)
- 1 → 차분(I)의 차수(d)
- 3 → 이동평균(MA) 부분의 차수(q)
summary() 메서드에서 확인할 수 있는 주요 정보
- Log Likelihood
- 모델이 데이터를 얼마나 잘 설명하는지를 나타내는 지표
- 값이 클수록 모델이 데이터를 더 잘 설명한다는 의미
- AIC (Akaike Information Criterion) / BIC (Bayesian Information Criterion)
- 모델의 복잡도를 평가하여 모델을 비교하는 지표
- 값이 작을수록 더 적합한 모델
- 자기회귀 계수 (ar.L1, ar.L2 등)
- ar.L1 계수가 1.534라면, 양수이므로 첫번째 지연값은 현재 값에 긍정적인 영향을 미칩니다.
- P-value (유의확률)
- 각 계수의 통계적 유의성을 평가
- p-value가 0.05 이하라면 해당 계수가 유의미하다고 판단 가능합니다.
- Ljung-Box / Jarque-Bera / Heteroskedasticity(이분산성) 테스트
- 잔차(오차)들이 정상성을 띄는지 평가
- 분산이 일정하면 정상적인 시계열로 볼 수 있음
- Prob(Q, JB, H)
- p-값이 낮을수록 (0.05 또는 0.01 이하), 잔차가 독립적이고 정규분포를 따르고 시차에 따라 일정한 분산을 가진다고 볼 수 있습니다.
이와 같은 정보를 활용하여 ARIMA 모델의 적합성을 평가하고, 시계열 예측의 신뢰도를 높일 수 있습니다.
from statsmodels.tsa.arima.model import ARIMA
# ARIMA 모델 설정 (p=2, d=1, q=3)
model = ARIMA(train_data['평균기온'], order=(2, 1, 3))
# 모델 학습
fitted_model = model.fit()
# 모델 요약 정보 출력
print(fitted_model.summary())
SARIMAX Results
==============================================================================
Dep. Variable: 평균기온 No. Observations: 22646
Model: ARIMA(2, 1, 3) Log Likelihood -49847.449
Date: Tue, 04 Feb 2025 AIC 99706.898
Time: 05:17:34 BIC 99755.064
Sample: 01-01-1960 HQIC 99722.562
- 12-31-2021
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 1.5315 0.009 174.175 0.000 1.514 1.549
ar.L2 -0.5380 0.009 -61.186 0.000 -0.555 -0.521
ma.L1 -1.5370 0.009 -166.594 0.000 -1.555 -1.519
ma.L2 0.1953 0.016 12.543 0.000 0.165 0.226
ma.L3 0.3487 0.007 50.592 0.000 0.335 0.362
sigma2 4.7803 0.033 143.299 0.000 4.715 4.846
===================================================================================
Ljung-Box (L1) (Q): 0.21 Jarque-Bera (JB): 5264.71
Prob(Q): 0.65 Prob(JB): 0.00
Heteroskedasticity (H): 0.96 Skew: -0.67
Prob(H) (two-sided): 0.08 Kurtosis: 4.94
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
※ 결과 해석
summary() 메서드를 사용하여 출력된 결과는 ARIMA 모델에 대한 통계적 정보를 나타냅니다.
Log Likelihood: 값이 약 -49847 이므로, 모델의 적합도가 최적이 아닐 수 있음을 시사합니다.
AIC/BIC: AIC는 약 99706.897, BIC는 약 99755.064로 나타나 있습니다. 이 값들은 비교적 높게 나타나 있어, 다른 파라미터 설정이 필요할 수 있음을 암시합니다.
coef: `ar.L1`, `ar.L2`, `ma.L1`, `ma.L2`, `ma.L3`의 계수들은 각각의 파라미터가 모델에 미치는 영향의 크기와 방향을 나타냅니다. 이 계수들은 통계적으로 유의미한 값을 보여주고 있습니다.
P>|z|: 각 계수의 p-value 값은 모든 파라미터에서 0.000으로 나타나 있어, 통계적으로 매우 유의미하다는 것을 나타냅니다.
Prob(Q, JB, H): Ljung-Box Test와 Heteroskedasticity Test의 p-값은 모델의 잔차가 무작위 노이즈를 따르지 않고, 시차에 따라 일정한 분산을 가진다는 귀무가설을 기각하지 않습니다. 반면, Jarque-Bera Test는 잔차가 정규 분포를 따르지 않는다는 귀무가설을 기각하고 있습니다.
따라서 각 항목을 고려하였을 때 모델이 가정을 충족하고 있으나, 전반적인 모델의 성능과 적합도에 대한 추가적인 검토가 필요함을 나타냅니다.
4. ARIMA 모델 검증 및 예측 정확도 평가
from sklearn.metrics import mean_absolute_error
# 검증 데이터 기간 동안 예측
validation_predictions = fitted_model.predict(start=validation_data.index[0], end=validation_data.index[-1])
# 실제 값과 예측 값의 차이 계산
mae = mean_absolute_error(validation_data['평균기온'], validation_predictions)
print(f'Mean Absolute Error (MAE): {mae}')
Mean Absolute Error (MAE): 32.58036858617591
5. 시계열 예측 결과의 시각적 비교
import matplotlib.pyplot as plt
# 축 객체 생성
fig, ax = plt.subplots(figsize=(12, 6))
# 예측 결과 시각화
ax.plot(validation_data.index, validation_data['평균기온'], label='실제', color='blue')
ax.plot(validation_data.index, validation_predictions, label='예측', color='red')
ax.set_xlabel('일시')
ax.set_ylabel('평균기온')
ax.legend()
plt.show()
※ 결과 해석
결과에서의 빨간색 선(예측)은 실제 파란색 선(실제)과 상당히 다른 패턴을 보이고 있습니다.
빨간색 선은 시각 부분에서는 실제 데이터가 하락하는 패턴을 보이는 반면 예측값은 그렇지 않습니다.
또한, 빨간색 선은 시간이 지남에 따라 점차 파란색 선과 멀어지며, 특히 그 차이는 시간이 끝날수록 더욱 커지는 것을 볼 수 있습니다.
이는 예측 모델이 실제 데이터의 변동성을 잘 포착하지 못했고, 특히 추세나 계절성과 같은 시간에 따른 변화를 반영하지 못했음을 나타냅니다.
이는 모델이 실제 데이터의 중요한 패턴을 학습하지 못했음을 의미합니다.
결론적으로, 예측 모델의 성능이 실제 데이터의 변동을 정확히 캐치하지 못하고 있으며, 모델의 개선이 필요함을 시사합니다.
이를 위해서는 모델 매개변수의 재조정, 데이터의 추가적인 전처리, 또는 다른 모델링 기법의 적용 등이 고려되어야 합니다.
6. 기간 설정 및 학습/검증 세트 재분할
# 원본 데이터를 복사하고 2015년 1월 1일 이후의 데이터를 사용합니다.
copy_train = train.copy()
slice_train = copy_train[copy_train.index >= '2015-01-01']
# 학습 데이터와 검증 데이터 분할
train_data = slice_train[slice_train.index < '2022-01-01']
validation_data = slice_train[slice_train.index >= '2022-01-01']
train_data.index.freq = 'D'
validation_data.index.freq = 'D'
※ 결과 해석
이 코드에서는 2015년 1월 1일 이후의 데이터를 선택하여 사용합니다.
과거의 모든 데이터를 사용하는 것보다 최신 데이터에 초점을 맞추면 과적합의 위험을 줄일 수 있습니다.
과거 데이터 중 일부는 현재 상황과 더 이상 관련이 없거나 노이즈를 포함할 수 있기 때문에, 이를 제외함으로써 일반화 능력을 향상시킬 수 있습니다.
데이터 양의 축소는 계산 효율성을 증가시킵니다.
적은 양의 데이터로 모델을 훈련시키면 학습 시간이 단축되고, 모델의 하이퍼파라미터 튜닝과 교차 검증 과정이 더 빠르게 이루어질 수 있습니다.
특히 많은 시계열 모델을 실험해야 하는 상황에서 유리합니다.
7. 시계열 데이터의 정상성 검증
adfuller() 함수는 시계열 데이터의 정상성을 검증하기 위한 Augmented Dickey-Fuller 테스트를 수행합니다.
정상성이란 시계열 데이터의 통계적 속성이 시간에 따라 변하지 않음을 의미하는데, 데이터의 평균과 분산이 일정함을 의미합니다.
ADF 테스트는 '단위근'이 존재하는 특정 유형의 비정상성을 확인합니다.
이 테스트에서 계산되는 ADF 통계치는 데이터에 단위근이 없다는 대립 가설에 대한 증거를 제공합니다.
p-value는 통계치가 얼마나 유의미한지를 나타냅니다.
일반적으로 p-value가 0.05 이하일 경우, 우리는 데이터가 정상성을 가진다고 간주하고 귀무 가설을 기각할 수 있습니다.
테스트 결과 ADF 통계량과 p-value가 출력됩니다.
p-value가 0.05 이하일 경우, 데이터는 정상성을 가진다고 볼 수 있습니다.
from statsmodels.tsa.stattools import adfuller
# ADF 테스트
adf_result = adfuller(train_data['평균기온'])
print(f'ADF 통계값: {adf_result[0]}')
print(f'p-value: {adf_result[1]}')
ADF 통계값: -2.9047450206736
p-value: 0.04479899198599979
※ 결과 해석
코드 실행 결과, ADF 통계값은 -2.9040로 나타났고, p-value는 0.044로 나타났습니다.
일반적으로 p-value가 0.05 이하일 경우에는 귀무 가설을 기각하고 데이터가 정상성을 가진다고 결론 내릴 수 있습니다.
그러나 이 경우 p-value가 0.050에 근접하여, 데이터가 정상 시계열이라고 결론짓기에는 약간의 불확실성이 남습니다.
이 경우, 차분을 수행하여 정상성을 확인하거나 모델의 성능을 검증하며 평가할 수 있습니다.
대신, 가능한 범위 내에서 여러 ARIMA 모델을 피팅하고 각 모델의 AIC(Akaike Information Criterion) 값을 비교함으로써 최적의 파라미터 조합을 찾는 방법을 사용할 것입니다.
AIC는 모델의 적합도와 복잡도를 동시에 고려하는 척도로, 낮은 AIC 값을 가지는 모델이 주어진 데이터에 대해 좋은 예측 성능을 제공할 가능성이 높습니다.
8. 자기상관 분석
ACF(자기상관 함수)는 시계열 데이터의 현재 값과 과거 값 사이의 상관 관계를 시간 지연에 따라 보여줍니다.
PACF(부분 자기상관 함수 )는 시간 지연 간의 상관 관계를 보여줍니다. 이는 다른 지연의 영향을 배제한 후의 상관 관계를 나타냅니다.
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# ACF와 PACF 플롯
fig, axes = plt.subplots(1, 2, figsize=(15,4))
plot_acf(train_data['평균기온'], ax=axes[0])
plot_pacf(train_data['평균기온'], ax=axes[1])
plt.show()
※ 결과 해석
ACF 그래프는 시차(lag)가 증가함에 따라 상관계수가 천천히 감소하는 패턴을 보였습니다.
일반적으로 이런 패턴은 데이터에 추세성이나 계절성 같은 비정상적인 성분이 포함될 수 있음을 나타냅니다.
그러나 ADF 테스트를 통해 데이터가 정상 시계열임을 확인하였기 때문에, 이는 단순히 높은 시차 간의 자기상관 패턴으로 볼 수 있으며, 데이터를 정상적인 시계열로 간주하고 ARIMA 모델에 사용할 것입니다.
이로 인해 데이터에 추가적인 차분이 필요 없으며, 이것이 모델링을 신뢰하고 해당 데이터로 모델링을 진행할 것임을 의미합니다.
PACF 그래프에서 신뢰구간 안으로 처음으로 들어가는 시차(lag)의 인덱스는 9번입니다.
해당 지점까지의 부분 자기상관계수가 통계적으로 유의미하다고 해석할 수 있습니다.
이는 ARIMA 모델에서 p(AR 차수)를 결정할 때 참고할 수 있는 정보를 제공합니다.
첫 번째 라그 이후 상관계수가 점차 감소하는 것은 이전 시점의 데이터가 현재 값에 영향을 미칠 수 있음을 나타내지만, 9번째 라그 이후에는 그러한 영향이 통계적으로 유의하지 않게 되어 추가적인 AR구성 요소가 필요하지 않을 수 있음을 의미합니다.
이정보를 기반으로 ARIMA 모델링시 p를 최대 8로 설정할 수 있으며, 이는 시계열 데이터의 현재 값이 과거 8개의 시점에 영향을 받을 수 있다는 가정 하에 모델을 구축하는 것입니다.
이후 모델 선택과정에서 AIC와 같은 정보 기준을 통해 더 낮은 p값이 더 적합한 모델을 제공하는지 평가할 수 있습니다.
9. 잔차 분석을 기반으로 한 ARIMA 모델 설정 및 학습
※ ARIMA 모델의 변형
AR 모델 (ARIMA(p, 0, 0))
이 모델은 ARIMA 모델의 간소화된 형태로, 차분(d)과 이동평균(MA) 부분이 없습니다.
이 모델은 과거의 값들이 미래 값에 어떻게 영향을 미치는지만을 고려합니다.
MA 모델 (ARIMA(0, 0, q))
이 모델은 AR 모델의 또 다른 변형으로, 자기회귀 부분이 없으며 과거의 예측 오차가 현재 값에 어떻게 영향을 미치는지를 모델링합니다.
ARMA 모델 (ARIMA(p, 0, q))
이 모델은 ARIMA 모델에서 차분 부분이 없는 형태로, 과거 값과 과거의 예측 오차가 현재 값에 어떻게 영향을 미치는지를 동시에 고려합니다.
AR 모델은 데이터에 명확한 자기상관이 있는 경우에 적합합니다.
MA 모델은 시계열 데이터의 오차 항에 관심이 있을 때 유용합니다.
ARMA 모델은 시계열 데이터가 정상적이고 AR과 MA 두 특성을 보일 때 적합합니다.
ARIMA 모델은 비정상적인 데이터에 적용되며, 가장 일반적인 형태의 시계열 모델링 방법 중 하나입니다.
- 이전 스텝에서 확인했던 내용을 바탕으로 p와 q, d의 범위를 설정합니다.
- itertools 모듈의 product() 함수를 사용하여 가능한 모든 p,d,q 조합을 생성하고 이를 pdq_combinations 변수에 저장합니다.
- for 반복문을 통해 생성된 각 조합에 대해 ARIMA 모델을 학습시키고,
AIC 값이 가장 낮은 모델을 best_model_results에 저장합니다. - 최종적으로, 가장 우수한 성능을 보인 최적의 모델을 요약 정보와 함께 출력하여 확인합니다.
import itertools
# p와 q의 최댓값을 8로 설정할 경우 약 5분 정도가 소요됩니다.
# 시간적 여유가 있을 경우에 p와 q의 범위를 range(0,9) 과 같이 설정해주세요.
p = range(0,4)
q = range(0,4)
d = [0] # d는 차분을 하지 않기 때문에 0으로 고정합니다.
# 가능한 모든 p,d,q의 조합 생성
pdq_combinations = list(itertools.product(p, d, q))
results_df = pd.DataFrame(columns=['p', 'd', 'q', 'AIC'])
best_aic = float('inf')
best_order = None
best_model_results = None
for order in pdq_combinations:
model = ARIMA(train_data['평균기온'], order=order)
fitted_model = model.fit()
aic = fitted_model.aic
results_df = results_df.append({'p': order[0], 'd': order[1], 'q': order[2], 'AIC': aic}, ignore_index=True)
if aic < best_aic:
best_aic = aic
best_order = order
best_model_results = fitted_model
results_df = results_df.sort_values(by='AIC')
if best_model_results:
print(best_model_results.summary())
SARIMAX Results
==============================================================================
Dep. Variable: 평균기온 No. Observations: 2557
Model: ARIMA(3, 0, 3) Log Likelihood -5661.670
Date: Tue, 04 Feb 2025 AIC 11339.340
Time: 07:22:13 BIC 11386.113
Sample: 01-01-2015 HQIC 11356.302
- 12-31-2021
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 8.6867 1.467 5.920 0.000 5.811 11.563
ar.L1 2.5156 0.026 97.568 0.000 2.465 2.566
ar.L2 -2.0337 0.052 -39.426 0.000 -2.135 -1.933
ar.L3 0.5179 0.026 20.061 0.000 0.467 0.568
ma.L1 -1.5255 0.026 -58.000 0.000 -1.577 -1.474
ma.L2 0.1450 0.046 3.175 0.001 0.055 0.234
ma.L3 0.3847 0.021 18.588 0.000 0.344 0.425
sigma2 4.8833 0.105 46.632 0.000 4.678 5.089
===================================================================================
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 894.52
Prob(Q): 0.88 Prob(JB): 0.00
Heteroskedasticity (H): 1.11 Skew: -0.83
Prob(H) (two-sided): 0.13 Kurtosis: 5.38
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
※ 결과 해석
제공된 모델 요약 결과에 따르면, 선택된 ARIMA(3, 0, 3) 모델은 자기회귀(AR) 부분에 3개의 래그, 이동평균(MA) 부분에도 3개의 래그를 사용하며, 차분은 적용되지 않은 모델입니다.
모델의 매개변수들 중 다수가 통계적으로 유의미한 것으로 나타나, 모델이 데이터의 패턴을 잘 포착하고 있음을 시사합니다.
ARIMA모델의 p,d,q 를 (3,0,3)으로 설정했을 때 가장 우수한 성능의 모델임을 확인했습니다.
10. 시계열 예측 성능 평가 및 결과 시각화
# ARIMA 모델 정의 및 피팅
model = ARIMA(train_data['평균기온'], order=(3, 0, 3))
fitted_model = model.fit()
# 동적 예측 수행
dynamic_forecast = fitted_model.get_prediction(start=validation_data.index[0], end=validation_data.index[-1], dynamic=False)
# 예측 결과와 신뢰 구간 추출
forecast = dynamic_forecast.predicted_mean
# 결과 시각화
plt.figure(figsize=(8, 6))
plt.plot(validation_data['평균기온'], label='실제 기온', color='blue')
plt.plot(forecast, label='동적 예측 기온', color='red', linestyle='--')
plt.title('실제 기온 vs 동적 예측 기온')
plt.xlabel('Date')
plt.ylabel('평균기온')
plt.legend()
plt.show()
※ 결과 해석
이전 스텝에서는 예측값이 실제값의 추세를 따라가지 못하는 경향이 있었으나,
최근의 접근 방식을 통해 예측값이 실제 기온의 장향성을 어느 정도 반영하고 있음을 확인할 수 있습니다.
그럼에도 불구하고, 모델은 여전히 데이터의 변동성, 특히 극단적인 고점과 저점에서의 움직임을 완전 포함하지 못하고 있습니다.
이는 모델이 데이터 내 주요 패턴을 일부 반영할 수 있다는 긍정적인 신호이지만, 데이터의 모든 변동 특성을 잡아내는 데에는 제한이 있음을 의미합니다.
이러한 예측의 한계를 차분을 수행하지 않은 비정상 시계열의 특성에서 비롯될 수 있습니다.
비록 ADF 테스트 등을 통해 데이터가 정상 시계열로 판단되더라도, 추가적인 차분이나 다른 비정상 시계열 처리 방법을 탐색함으로써 모델의 예측 성능은 더욱 개선될 수 있을 것입니다.