基于梯度提升的时序预测: Skforecast、XGBoost、LightGBM、Scikit-learn 和 CatBoost

# 数据处理
# ==============================================================================
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from feature_engine.timeseries.forecasting import WindowFeatures

# 绘图
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from skforecast.plot import plot_residuals, calculate_lag_autocorrelation
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
poff.init_notebook_mode(connected=True)
plt.style.use('seaborn-v0_8-darkgrid')

# 建模与预测
# ==============================================================================
import xgboost
import lightgbm
import catboost
import sklearn
import shap
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from catboost import CatBoostRegressor
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import (
    OneHotEncoder,
    OrdinalEncoder,
    FunctionTransformer,
    PolynomialFeatures,
)
from sklearn.feature_selection import RFECV
from sklearn.pipeline import make_pipeline
from sklearn.compose import make_column_transformer, make_column_selector

import skforecast
from skforecast.recursive import ForecasterEquivalentDate, ForecasterRecursive
from skforecast.model_selection import (
    TimeSeriesFold,
    OneStepAheadFold,
    bayesian_search_forecaster,
    backtesting_forecaster,
)
from skforecast.preprocessing import RollingFeatures
from skforecast.feature_selection import select_features
from skforecast.metrics import calculate_coverage

# 警告配置
# ==============================================================================
import warnings
warnings.filterwarnings('once')

color = '\033[1m\033[38;5;208m' 
print(f"{color}版本 skforecast: {skforecast.__version__}")
print(f"{color}版本 scikit-learn: {sklearn.__version__}")
print(f"{color}版本 lightgbm: {lightgbm.__version__}")
print(f"{color}版本 xgboost: {xgboost.__version__}")
print(f"{color}版本 catboost: {catboost.__version__}")
print(f"{color}版本 pandas: {pd.__version__}")
print(f"{color}版本 numpy: {np.__version__}")

# 数据下载
# ==============================================================================
data = fetch_dataset('bike_sharing', raw=True)

# 数据预处理（设置索引与频率）
# ==============================================================================
data = data[['date_time', 'users', 'holiday', 'weather', 'temp', 'atemp', 'hum', 'windspeed']]
data['date_time'] = pd.to_datetime(data['date_time'], format='%Y-%m-%d %H:%M:%S')
data = data.set_index('date_time')
if pd.__version__ < '2.2':
    data = data.asfreq('H')
else:
    data = data.asfreq('h')
data = data.sort_index()
data.head()

# 训练-验证-测试 划分
# ==============================================================================
end_train = '2012-04-30 23:59:00'
end_validation = '2012-08-31 23:59:00'
data_train = data.loc[: end_train, :]
data_val   = data.loc[end_train:end_validation, :]
data_test  = data.loc[end_validation:, :]

print(f"训练集日期 : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"验证集日期 : {data_val.index.min()} --- {data_val.index.max()}  (n={len(data_val)})")
print(f"测试集日期 : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")

# 基准模型：使用“前一天同一小时”的值
# ==============================================================================
forecaster = ForecasterEquivalentDate(
                 offset    = pd.DateOffset(days=1),
                 n_offsets = 1
             )

# 训练预测器
# ==============================================================================
forecaster.fit(y=data.loc[:end_validation, 'users'])
forecaster

# 回测
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(data.loc[:end_validation]))
metric_baseline, predictions = backtesting_forecaster(
                                   forecaster = forecaster,
                                   y          = data['users'],
                                   cv         = cv,
                                   metric     = 'mean_absolute_error'
                               )
metric_baseline

# 创建预测器
# ==============================================================================
window_features = RollingFeatures(stats=["mean"], window_sizes=24 * 3)
forecaster = ForecasterRecursive(
                estimator       = LGBMRegressor(random_state=15926, verbose=-1),
                lags            = 24,
                window_features = window_features
             )

# 训练预测器
# ==============================================================================
forecaster.fit(y=data.loc[:end_validation, 'users'])
forecaster

# 预测
# ==============================================================================
forecaster.predict(steps=10)

# 在测试集上进行回测
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(data.loc[:end_validation]))
metric, predictions = backtesting_forecaster(
                          forecaster = forecaster,
                          y          = data['users'],
                          cv         = cv,
                          metric     = 'mean_absolute_error'
                      )
predictions.head()

# 回测误差
# ==============================================================================
metric

# 超参数搜索
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator        = LGBMRegressor(random_state=15926, verbose=-1),
                lags             = 72,
                window_features  = window_features,
                transformer_exog = transformer_exog,
                fit_kwargs       = {"categorical_feature": "auto"}
             )

# 滞后项网格
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

# 回归器超参数搜索空间
def search_space(trial):
    search_space  = {
        'n_estimators'    : trial.suggest_int('n_estimators', 300, 1000, step=100),
        'max_depth'       : trial.suggest_int('max_depth', 3, 10),
        'min_data_in_leaf': trial.suggest_int('min_data_in_leaf', 25, 500),
        'learning_rate'   : trial.suggest_float('learning_rate', 0.01, 0.5),
        'feature_fraction': trial.suggest_float('feature_fraction', 0.5, 1),
        'max_bin'         : trial.suggest_int('max_bin', 50, 250),
        'reg_alpha'       : trial.suggest_float('reg_alpha', 0, 1),
        'reg_lambda'      : trial.suggest_float('reg_lambda', 0, 1),
        'lags'            : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

# 训练/验证折叠
cv_search = TimeSeriesFold(steps = 36, initial_train_size = len(data_train))

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'], # 不使用测试集数据
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20, # 若需更全面的搜索可增大该值
    return_best   = True
)

best_params = results_search['params'].iat[0]
best_params = best_params | {'random_state': 15926, 'verbose': -1}
best_lags   = results_search['lags'].iat[0]

# 搜索结果
# ==============================================================================
results_search.head(3)

# 回测模型
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(data.loc[:end_validation]))
metric, predictions = backtesting_forecaster(
                        forecaster = forecaster,
                        y          = data['users'],
                        exog       = data[exog_features],
                        cv         = cv,
                        metric     = 'mean_absolute_error'
                     )
metric

# 预测值 vs 真实值（绘图）
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="测试集", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="预测", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="测试集：真实值 vs 预测值",
    xaxis_title="日期时间",
    yaxis_title="用户数",
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    legend=dict(orientation="h", yanchor="top", y=1.1, xanchor="left", x=0.001)
)
fig.show()

# 创建预测器
# ==============================================================================
estimator = LGBMRegressor(
    n_estimators = 100,
    max_depth    = 5,
    random_state = 15926,
    verbose      = -1
)
forecaster = ForecasterRecursive(
    estimator        = estimator,
    lags             = [1, 2, 3, 23, 24, 25, 167, 168, 169],
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"}
)

# 带交叉验证的递归特征消除（RFECV）
# ==============================================================================
warnings.filterwarnings("ignore", message="X does not have valid feature names.*")
selector = RFECV(
    estimator = estimator,
    step      = 1,
    cv        = 3,
)
selected_lags, selected_window_features, selected_exog = select_features(
    forecaster      = forecaster,
    selector        = selector,
    y               = data_train['users'],  
    exog            = data_train[exog_features],
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.5,
    random_state    = 123,
    verbose         = True,
)

# 使用已选特征创建预测器
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator        = LGBMRegressor(**best_params),
    lags             = selected_lags,
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"}
)

# 带外生变量的回测（测试集）
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(data.loc[:end_validation]))
metric_lgbm, predictions = backtesting_forecaster(
    forecaster = forecaster,
    y          = data['users'],
    exog       = data[selected_exog],
    cv         = cv,
    metric     = 'mean_absolute_error'
)
metric_lgbm

# 更新外生变量集合
# ==============================================================================
exog_features = selected_exog

# 创建并训练预测器
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator        = LGBMRegressor(**best_params),
    lags             = [1, 2, 3, 23, 24, 25, 167, 168, 169],
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"},
    binner_kwargs    = {"n_bins": 5}
)
forecaster.fit(
    y    = data.loc[:end_train, 'users'],
    exog = data.loc[:end_train, exog_features],
    store_in_sample_residuals = True
)

# 预测区间
# ==============================================================================
# 由于模型在训练时使用了外生变量，进行预测时也必须提供相应的外生变量。
predictions = forecaster.predict_interval(
    exog     = data.loc[end_train:, exog_features],
    steps    = 24,
    interval = [5, 95],
    method   = 'conformal',
)
predictions.head()

# 在验证集上回测以获得样本外残差
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(data.loc[:end_train]))
_, predictions_val = backtesting_forecaster(
    forecaster = forecaster,
    y          = data.loc[:end_validation, 'users'],
    exog       = data.loc[:end_validation, exog_features],
    cv         = cv,
    metric     = 'mean_absolute_error'
)

# 样本外残差的分布
# ==============================================================================
residuals = data.loc[predictions_val.index, 'users'] - predictions_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(
        y_true = data.loc[predictions_val.index, 'users'],
        y_pred = predictions_val['pred'],
        figsize=(7, 4)
    )

# 在预测器中存储样本外残差
# ==============================================================================
forecaster.set_out_sample_residuals(
    y_true = data.loc[predictions_val.index, 'users'],
    y_pred = predictions_val['pred']
)

# 使用样本外残差在测试集上进行带预测区间的回测
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(data.loc[:end_validation]))
metric, predictions = backtesting_forecaster(
   forecaster              = forecaster,
   y                       = data['users'],
   exog                    = data[exog_features],
   cv                      = cv,
   metric                  = 'mean_absolute_error',
   interval                = [5, 95],  # 90% 预测区间
   interval_method         = 'conformal',
   use_in_sample_residuals = False,  # 使用样本外残差
   use_binned_residuals    = True,   # 根据预测值范围选择（分箱）残差
)
predictions.head(5)

# 预测区间 vs 真实值（绘图）
# ==============================================================================
fig = go.Figure([
    go.Scatter(name='预测值', x=predictions.index, y=predictions['pred'], mode='lines'),
    go.Scatter(
        name='真实值', x=data_test.index, y=data_test['users'], mode='lines',
    ),
    go.Scatter(
        name='上界', x=predictions.index, y=predictions['upper_bound'], mode='lines',
        marker=dict(color="#444"), line=dict(width=0), showlegend=False
    ),
    go.Scatter(
        name='下界', x=predictions.index, y=predictions['lower_bound'], marker=dict(color="#444"),
        line=dict(width=0), mode='lines', fillcolor='rgba(68, 68, 68, 0.3)', fill='tonexty', showlegend=False
    )
])
fig.update_layout(
    title="测试集：真实值 vs 预测值（含预测区间）",
    xaxis_title="日期时间",
    yaxis_title="用户数",
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    hovermode="x",
    legend=dict(orientation="h", yanchor="top", y=1.1, xanchor="left", x=0.001),
    # 初始 x 轴缩放范围：10 月 1 日至 10 日
    xaxis=dict(range=['2012-10-01', '2012-10-10'])
)
fig.show()

# 预测区间覆盖率（测试集）
# ==============================================================================
coverage = calculate_coverage(
              y_true       = data.loc[end_validation:, "users"],
              lower_bound  = predictions["lower_bound"],
              upper_bound  = predictions["upper_bound"]
            )
area = (predictions['upper_bound'] - predictions['lower_bound']).sum()
print(f"区间总面积: {round(area, 2)}")
print(f"预测区间覆盖率: {round(100 * coverage, 2)} %")

# 创建并训练预测器
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator        = LGBMRegressor(**best_params),
    lags             = [1, 2, 3, 23, 24, 25, 167, 168, 169],
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"}
)
forecaster.fit(
    y    = data.loc[:end_validation, 'users'],
    exog = data.loc[:end_validation, exog_features]
)

# 用于超参数搜索与回测的折叠（folds）
# ==============================================================================
cv_search = OneStepAheadFold(initial_train_size = len(data_train))
cv_backtesting = TimeSeriesFold(steps = 36, initial_train_size = len(data[:end_validation]))