시계열 예측

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')
 
# 재현 가능한 결과를 위한 시드 설정
np.random.seed(42)
 
# 2년치 일별 매출 데이터 생성
date_range = pd.date_range(start='2022-01-01', end='2023-12-31', freq='D')
n_days = len(date_range)
 
# 기본 매출 (추세 + 노이즈)
base_sales = 10000 + np.linspace(0, 2000, n_days)  # 상승 추세
 
# 주간 계절성 (주말 매출 증가)
weekly_pattern = np.array([1.0, 0.95, 0.9, 0.95, 1.1, 1.3, 1.2])
weekly_seasonality = np.tile(weekly_pattern, n_days // 7 + 1)[:n_days]
 
# 월간 계절성 (월말 매출 증가)
monthly_seasonality = 1 + 0.1 * np.sin(2 * np.pi * np.arange(n_days) / 30)
 
# 연간 계절성 (11-12월 연말 시즌)
yearly_seasonality = np.ones(n_days)
for i, date in enumerate(date_range):
    if date.month in [11, 12]:
        yearly_seasonality[i] = 1.3
    elif date.month in [1, 2]:
        yearly_seasonality[i] = 0.8
 
# 최종 매출
daily_sales = pd.DataFrame({
    'ds': date_range,
    'y': (base_sales * weekly_seasonality * monthly_seasonality *
          yearly_seasonality * (1 + np.random.normal(0, 0.05, n_days)))
})
 
print(f"데이터 기간: {daily_sales['ds'].min().date()} ~ {daily_sales['ds'].max().date()}")
print(f"총 일수: {len(daily_sales)}")
print(f"평균 일별 매출: ${daily_sales['y'].mean():,.0f}")
print(daily_sales.head())

# 전체 시계열 시각화
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
 
# 전체 매출 추이
axes[0].plot(daily_sales['ds'], daily_sales['y'], linewidth=0.8)
axes[0].set_title('일별 매출 추이', fontsize=14, fontweight='bold')
axes[0].set_xlabel('날짜')
axes[0].set_ylabel('매출 ($)')
axes[0].grid(True, alpha=0.3)
 
# 월별 집계
monthly_sales = daily_sales.set_index('ds').resample('M').sum()
axes[1].bar(monthly_sales.index, monthly_sales['y'], width=20, color='steelblue')
axes[1].set_title('월별 매출', fontsize=14, fontweight='bold')
axes[1].set_xlabel('월')
axes[1].set_ylabel('매출 ($)')
axes[1].grid(True, alpha=0.3)
 
# 요일별 평균
daily_sales['dayofweek'] = daily_sales['ds'].dt.day_name()
dow_order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
dow_sales = daily_sales.groupby('dayofweek')['y'].mean().reindex(dow_order)
axes[2].bar(dow_sales.index, dow_sales.values, color='coral')
axes[2].set_title('요일별 평균 매출', fontsize=14, fontweight='bold')
axes[2].set_xlabel('요일')
axes[2].set_ylabel('평균 매출 ($)')
axes[2].tick_params(axis='x', rotation=45)
 
plt.tight_layout()
plt.show()

from statsmodels.tsa.seasonal import seasonal_decompose
 
# 시계열 분해 (곱셈 모델)
ts_data = daily_sales.set_index('ds')['y']
decomposition = seasonal_decompose(ts_data, model='multiplicative', period=7)
 
fig, axes = plt.subplots(4, 1, figsize=(14, 12))
 
# 원본 데이터
axes[0].plot(decomposition.observed)
axes[0].set_title('원본 데이터 (Observed)', fontsize=12)
axes[0].set_ylabel('매출')
 
# 추세
axes[1].plot(decomposition.trend, color='orange')
axes[1].set_title('추세 (Trend)', fontsize=12)
axes[1].set_ylabel('매출')
 
# 계절성
axes[2].plot(decomposition.seasonal, color='green')
axes[2].set_title('계절성 (Seasonal)', fontsize=12)
axes[2].set_ylabel('계절 지수')
 
# 잔차
axes[3].plot(decomposition.resid, color='red')
axes[3].set_title('잔차 (Residual)', fontsize=12)
axes[3].set_ylabel('잔차')
 
plt.tight_layout()
plt.show()
 
print("=== 시계열 분해 요약 ===")
print(f"추세 범위: ${decomposition.trend.min():,.0f} ~ ${decomposition.trend.max():,.0f}")
print(f"계절성 범위: {decomposition.seasonal.min():.2f} ~ {decomposition.seasonal.max():.2f}")

# 이동 평균 계산
daily_sales['ma_7'] = daily_sales['y'].rolling(window=7).mean()
daily_sales['ma_30'] = daily_sales['y'].rolling(window=30).mean()
 
# 시각화
plt.figure(figsize=(14, 6))
plt.plot(daily_sales['ds'][-180:], daily_sales['y'][-180:],
         alpha=0.5, label='실제 매출')
plt.plot(daily_sales['ds'][-180:], daily_sales['ma_7'][-180:],
         linewidth=2, label='7일 이동평균')
plt.plot(daily_sales['ds'][-180:], daily_sales['ma_30'][-180:],
         linewidth=2, label='30일 이동평균')
plt.title('이동 평균 (최근 180일)', fontsize=14, fontweight='bold')
plt.xlabel('날짜')
plt.ylabel('매출 ($)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

from statsmodels.tsa.holtwinters import SimpleExpSmoothing, ExponentialSmoothing
 
# 학습/테스트 분리
train_size = int(len(daily_sales) * 0.8)
train = daily_sales['y'][:train_size]
test = daily_sales['y'][train_size:]
 
print(f"학습 데이터: {train_size}일")
print(f"테스트 데이터: {len(test)}일")
 
# 단순 지수 평활
ses_model = SimpleExpSmoothing(train).fit(smoothing_level=0.2, optimized=False)
ses_forecast = ses_model.forecast(len(test))
 
# 평가
from sklearn.metrics import mean_absolute_error, mean_squared_error
 
mae_ses = mean_absolute_error(test, ses_forecast)
rmse_ses = np.sqrt(mean_squared_error(test, ses_forecast))
print(f"\n=== 단순 지수 평활 결과 ===")
print(f"MAE: ${mae_ses:,.2f}")
print(f"RMSE: ${rmse_ses:,.2f}")

# Holt-Winters 모델 (추세 + 계절성)
hw_model = ExponentialSmoothing(
    train,
    trend='add',           # 가법적 추세
    seasonal='mul',        # 승법적 계절성
    seasonal_periods=7     # 주간 계절성
).fit()
 
hw_forecast = hw_model.forecast(len(test))
 
# 평가
mae_hw = mean_absolute_error(test, hw_forecast)
rmse_hw = np.sqrt(mean_squared_error(test, hw_forecast))
print("=== Holt-Winters 결과 ===")
print(f"MAE: ${mae_hw:,.2f}")
print(f"RMSE: ${rmse_hw:,.2f}")
 
# 예측 시각화
plt.figure(figsize=(14, 6))
plt.plot(daily_sales['ds'][:train_size], train, label='학습 데이터', alpha=0.7)
plt.plot(daily_sales['ds'][train_size:], test, label='실제 값', alpha=0.7)
plt.plot(daily_sales['ds'][train_size:], hw_forecast,
         label='Holt-Winters 예측', linewidth=2, color='red')
plt.axvline(x=daily_sales['ds'].iloc[train_size], color='gray',
            linestyle='--', label='예측 시작')
plt.title('Holt-Winters 예측 결과', fontsize=14, fontweight='bold')
plt.xlabel('날짜')
plt.ylabel('매출 ($)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
 
# 정상성 검정 (ADF 테스트)
adf_result = adfuller(daily_sales['y'])
print("=== ADF 검정 (정상성 테스트) ===")
print(f"ADF 통계량: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
print(f"정상성 여부: {'정상' if adf_result[1] < 0.05 else '비정상 (차분 필요)'}")
 
# ACF, PACF 시각화
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(daily_sales['y'], lags=30, ax=axes[0])
axes[0].set_title('자기상관함수 (ACF)')
plot_pacf(daily_sales['y'], lags=30, ax=axes[1])
axes[1].set_title('편자기상관함수 (PACF)')
plt.tight_layout()
plt.show()

from statsmodels.tsa.arima.model import ARIMA
 
# ARIMA 모델 (p=1, d=1, q=1)
arima_model = ARIMA(train, order=(1, 1, 1)).fit()
 
print("=== ARIMA 모델 요약 ===")
print(f"AIC: {arima_model.aic:.2f}")
print(f"BIC: {arima_model.bic:.2f}")
 
# 예측
arima_forecast = arima_model.forecast(len(test))
 
# 평가
mae_arima = mean_absolute_error(test, arima_forecast)
rmse_arima = np.sqrt(mean_squared_error(test, arima_forecast))
print(f"\nMAE: ${mae_arima:,.2f}")
print(f"RMSE: ${rmse_arima:,.2f}")

from statsmodels.tsa.statespace.sarimax import SARIMAX
 
# SARIMA 모델
sarima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 7)  # 주간 계절성
).fit(disp=False)
 
sarima_forecast = sarima_model.forecast(len(test))
 
# 평가
mae_sarima = mean_absolute_error(test, sarima_forecast)
rmse_sarima = np.sqrt(mean_squared_error(test, sarima_forecast))
print("=== SARIMA 결과 ===")
print(f"MAE: ${mae_sarima:,.2f}")
print(f"RMSE: ${rmse_sarima:,.2f}")

# 모델별 성능 비교
models_comparison = pd.DataFrame({
    '모델': ['단순 지수평활', 'Holt-Winters', 'ARIMA', 'SARIMA'],
    'MAE': [mae_ses, mae_hw, mae_arima, mae_sarima],
    'RMSE': [rmse_ses, rmse_hw, rmse_arima, rmse_sarima]
}).round(2)
 
models_comparison['MAPE'] = (models_comparison['MAE'] / test.mean() * 100).round(2)
 
print("=== 시계열 모델 성능 비교 ===")
print(models_comparison.to_string(index=False))
 
# 시각화
fig, ax = plt.subplots(figsize=(10, 6))
x = np.arange(len(models_comparison))
width = 0.35
 
bars1 = ax.bar(x - width/2, models_comparison['MAE'], width, label='MAE', color='steelblue')
bars2 = ax.bar(x + width/2, models_comparison['RMSE'], width, label='RMSE', color='coral')
 
ax.set_xlabel('모델')
ax.set_ylabel('오차 ($)')
ax.set_title('시계열 모델 성능 비교', fontsize=14, fontweight='bold')
ax.set_xticks(x)
ax.set_xticklabels(models_comparison['모델'])
ax.legend()
ax.grid(True, alpha=0.3)
 
plt.tight_layout()
plt.show()

# 전체 데이터로 최종 모델 학습
final_model = ExponentialSmoothing(
    daily_sales['y'],
    trend='add',
    seasonal='mul',
    seasonal_periods=7
).fit()
 
# 90일 예측
forecast_days = 90
future_forecast = final_model.forecast(forecast_days)
future_dates = pd.date_range(
    start=daily_sales['ds'].max() + timedelta(days=1),
    periods=forecast_days
)
 
# 신뢰구간 (단순 추정)
forecast_std = daily_sales['y'].std() * 0.1
upper_bound = future_forecast + 1.96 * forecast_std
lower_bound = future_forecast - 1.96 * forecast_std
 
# 시각화
plt.figure(figsize=(14, 6))
plt.plot(daily_sales['ds'][-90:], daily_sales['y'][-90:],
         label='실제 데이터', alpha=0.7)
plt.plot(future_dates, future_forecast,
         label='예측', color='red', linewidth=2)
plt.fill_between(future_dates, lower_bound, upper_bound,
                 alpha=0.3, color='red', label='95% 신뢰구간')
plt.axvline(x=daily_sales['ds'].max(), color='gray',
            linestyle='--', label='예측 시작')
plt.title('90일 매출 예측', fontsize=14, fontweight='bold')
plt.xlabel('날짜')
plt.ylabel('매출 ($)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 월별 예측 합계
forecast_df = pd.DataFrame({
    'ds': future_dates,
    'yhat': future_forecast,
    'yhat_lower': lower_bound,
    'yhat_upper': upper_bound
})
 
forecast_df['month'] = forecast_df['ds'].dt.to_period('M')
monthly_forecast = forecast_df.groupby('month').agg(
    예측_매출=('yhat', 'sum'),
    하한=('yhat_lower', 'sum'),
    상한=('yhat_upper', 'sum')
).round(0)
 
print("=== 월별 매출 예측 ===")
print(monthly_forecast)
 
# 요일별 평균 예측
forecast_df['dayofweek'] = forecast_df['ds'].dt.day_name()
daily_pattern = forecast_df.groupby('dayofweek')['yhat'].mean().reindex(dow_order)
print("\n=== 요일별 평균 예측 매출 ===")
for day, value in daily_pattern.items():
    print(f"{day}: ${value:,.0f}")

# 주요 휴일 정의
holidays = pd.DataFrame({
    'holiday': ['New Year', 'Christmas', 'Black Friday', 'Thanksgiving'],
    'ds': pd.to_datetime(['2023-01-01', '2023-12-25', '2023-11-24', '2023-11-23'])
})
 
# 휴일 주변 매출 분석
def analyze_holiday_effect(data, holiday_date, window=7):
    mask = (data['ds'] >= holiday_date - timedelta(days=window)) & \
           (data['ds'] <= holiday_date + timedelta(days=window))
    holiday_sales = data.loc[mask, 'y'].mean()
    normal_sales = data['y'].mean()
    effect = (holiday_sales - normal_sales) / normal_sales * 100
    return effect
 
print("=== 휴일 효과 분석 ===")
for _, row in holidays.iterrows():
    effect = analyze_holiday_effect(daily_sales, row['ds'])
    print(f"{row['holiday']}: {effect:+.1f}%")