使用梯度提升进行时间序列预测 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

# 图表绘制
# ==============================================================================
from skforecast.plot import set_dark_theme
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, PolynomialFeatures
from sklearn.feature_selection import RFECV
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}Version skforecast: {skforecast.__version__}")
print(f"{color}Version scikit-learn: {sklearn.__version__}")
print(f"{color}Version lightgbm: {lightgbm.__version__}")
print(f"{color}Version xgboost: {xgboost.__version__}")
print(f"{color}Version catboost: {catboost.__version__}")
print(f"{color}Version pandas: {pd.__version__}")
print(f"{color}Version 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"Dates train      : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Dates validacion : {data_val.index.min()} --- {data_val.index.max()}  (n={len(data_val)})")
print(f"Dates test       : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")

# 时间序列交互图
# ==============================================================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=data_train.index, y=data_train['users'], mode='lines', name='Train'))
fig.add_trace(go.Scatter(x=data_val.index, y=data_val['users'], mode='lines', name='Validation'))
fig.add_trace(go.Scatter(x=data_test.index, y=data_test['users'], mode='lines', name='Test'))
fig.update_layout(
    title  = 'Number of users',
    xaxis_title="Time",
    yaxis_title="Users",
    legend_title="Partition:",
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    legend=dict(orientation="h", yanchor="top", y=1, xanchor="left", x=0.001)
)
#fig.update_xaxes(rangeslider_visible=True)
fig.show()

# 带缩放功能的时间序列静态图
# ==============================================================================
zoom = ('2011-08-01 00:00:00', '2011-08-15 00:00:00')
fig, axs = plt.subplots(2, 1, figsize=(8, 4), gridspec_kw={"height_ratios": [1, 2]})
data_train['users'].plot(ax=axs[0], label='train', alpha=0.5)
data_val['users'].plot(ax=axs[0], label='validation', alpha=0.5)
data_test['users'].plot(ax=axs[0], label='test', alpha=0.5)
axs[0].axvspan(zoom[0], zoom[1], color="blue", alpha=0.7)
axs[0].set_title("Number of users")
axs[0].set_xlabel("")
# 局部放大图
data.loc[zoom[0] : zoom[1], "users"].plot(ax=axs[1], color="blue")
axs[1].set_title(f"Zoom: {zoom[0]} to {zoom[1]}", fontsize=10)
plt.tight_layout()
plt.show()

# 年度、周度与日度季节性
# ==============================================================================
set_dark_theme()
fig, axs = plt.subplots(2, 2, figsize=(8.5, 5.5), sharex=False, sharey=True)
axs = axs.ravel()
flierprops=dict(marker='o', markerfacecolor='white', markeredgecolor='black', markersize=4)
    
# Users distribution by month
data['month'] = data.index.month
data.boxplot(column='users', by='month', ax=axs[0], flierprops=flierprops)
data.groupby('month')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[0])
axs[0].set_ylabel('Users')
axs[0].set_title('Users distribution by month', fontsize=10)

# Users distribution by week day
data['week_day'] = data.index.day_of_week + 1
data.boxplot(column='users', by='week_day', ax=axs[1], flierprops=flierprops)
data.groupby('week_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[1])
axs[1].set_ylabel('Users')
axs[1].set_title('Users distribution by week day', fontsize=10)

# Users distribution by the hour of the day
data['hour_day'] = data.index.hour + 1
data.boxplot(column='users', by='hour_day', ax=axs[2], flierprops=flierprops)
data.groupby('hour_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[2])
axs[2].set_ylabel('Users')
axs[2].set_title('Users distribution by the hour of the day', fontsize=10)

# Users distribution by week day and hour of the day
mean_day_hour = data.groupby(["week_day", "hour_day"])["users"].mean()
mean_day_hour.plot(ax=axs[3])
axs[3].set(
    title       = "Mean users during week",
    xticks      = [i * 24 for i in range(7)],
    xticklabels = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
    xlabel      = "Day and hour",
    ylabel      = "Number of users"
)
axs[3].title.set_size(10)

fig.suptitle("Seasonality plots", fontsize=12)
fig.tight_layout()

# 自相关图
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_acf(data['users'], ax=ax, lags=72)
plt.show()

# 偏自相关图
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(data['users'], ax=ax, lags=72, method='ywm')
plt.show()

# 偏自相关绝对值最大的前 10 个滞后阶数
# ==============================================================================
calculate_lag_autocorrelation(
    data    = data['users'],
    n_lags  = 60,
    sort_by = "partial_autocorrelation_abs"
).head(10)

# 创建基线：前一天同一小时的值
# ==============================================================================
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

# 日历特征
# ==============================================================================
features_to_extract = [
    'month',
    'week',
    'day_of_week',
    'hour'
]
calendar_transformer = DatetimeFeatures(
    variables           ='index',
    features_to_extract = features_to_extract,
    drop_original       = True,
)
calendar_features = calendar_transformer.fit_transform(data)[features_to_extract]

# 日照特征
# ==============================================================================
location = LocationInfo(
    name      = 'Washington DC',
    region    = 'USA',
    timezone  = 'US/Eastern',
    latitude  = 40.516666666666666,
    longitude = -77.03333333333333
)
sunrise_hour = [
    sun(location.observer, date=date, tzinfo=location.timezone)['sunrise'].hour
    for date in data.index
]
sunset_hour = [
    sun(location.observer, date=date, tzinfo=location.timezone)['sunset'].hour
    for date in data.index
]
sun_light_features = pd.DataFrame({
                        'sunrise_hour': sunrise_hour,
                        'sunset_hour': sunset_hour}, 
                        index = data.index
                     )
sun_light_features['daylight_hours'] = (
    sun_light_features['sunset_hour'] - sun_light_features['sunrise_hour']
)
sun_light_features["is_daylight"] = np.where(
    (data.index.hour >= sun_light_features["sunrise_hour"])
    & (data.index.hour < sun_light_features["sunset_hour"]),
    1,
    0,
)

# 节假日特征
# ==============================================================================
holiday_features = data[['holiday']].astype(int)
holiday_features['holiday_previous_day'] = holiday_features['holiday'].shift(24)
holiday_features['holiday_next_day'] = holiday_features['holiday'].shift(-24)

# 温度的滚动窗口
# ==============================================================================
wf_transformer = WindowFeatures(
    variables   = ["temp"],
    window      = ["1D", "7D"],
    functions   = ["mean", "max", "min"],
    freq        = "h",
)
temp_features = wf_transformer.fit_transform(data[['temp']])


# 合并所有外生变量
# ==============================================================================
assert all(calendar_features.index == sun_light_features.index)
assert all(calendar_features.index == temp_features.index)
assert all(calendar_features.index == holiday_features.index)
df_exogenous_features = pd.concat([
    calendar_features,
    sun_light_features,
    temp_features,
    holiday_features
], axis=1)

# 由于创建了移动平均，序列开头存在缺失值；
# 由于 holiday_next_day，序列末尾也存在缺失值。
df_exogenous_features = df_exogenous_features.iloc[7 * 24:, :]
df_exogenous_features = df_exogenous_features.iloc[:-24, :]
df_exogenous_features.head(3)

# 日历与日照特征的周期性编码
# ==============================================================================
features_to_encode = [
    "month",
    "week",
    "day_of_week",
    "hour",
    "sunrise_hour",
    "sunset_hour",
]
max_values = {
    "month": 12,
    "week": 52,
    "day_of_week": 7,
    "hour": 24,
    "sunrise_hour": 24,
    "sunset_hour": 24,
}
cyclical_encoder = CyclicalFeatures(
    variables     = features_to_encode,
    max_values    = max_values,
    drop_original = False
)

df_exogenous_features = cyclical_encoder.fit_transform(df_exogenous_features)
df_exogenous_features.head(3)

# 外生变量之间的交互
# ==============================================================================
transformer_poly = PolynomialFeatures(
                        degree           = 2,
                        interaction_only = True,
                        include_bias     = False
                    ).set_output(transform="pandas")
poly_cols = [
    'month_sin', 
    'month_cos',
    'week_sin',
    'week_cos',
    'day_of_week_sin',
    'day_of_week_cos',
    'hour_sin',
    'hour_cos',
    'sunrise_hour_sin',
    'sunrise_hour_cos',
    'sunset_hour_sin',
    'sunset_hour_cos',
    'daylight_hours',
    'is_daylight',
    'holiday_previous_day',
    'holiday_next_day',
    'temp_window_1D_mean',
    'temp_window_1D_min',
    'temp_window_1D_max',
    'temp_window_7D_mean',
    'temp_window_7D_min',
    'temp_window_7D_max',
    'temp',
    'holiday'
]
poly_features = transformer_poly.fit_transform(df_exogenous_features[poly_cols])
poly_features = poly_features.drop(columns=poly_cols)
poly_features.columns = [f"poly_{col}" for col in poly_features.columns]
poly_features.columns = poly_features.columns.str.replace(" ", "__")
assert all(poly_features.index == df_exogenous_features.index)
df_exogenous_features = pd.concat([df_exogenous_features, poly_features], axis=1)
df_exogenous_features.head(3)

# 将类别变量存储为 category 类型
# ==============================================================================
data["weather"] = data["weather"].astype("category")

# 独热编码转换器
# ==============================================================================
one_hot_encoder = make_column_transformer(
    (
        OneHotEncoder(sparse_output=False, drop='if_binary'),
        make_column_selector(dtype_include=['category', 'object']),
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
).set_output(transform="pandas")  

# 创建带有外生特征转换器的预测器
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator        = LGBMRegressor(random_state=15926, verbose=-1),
                lags             = 72,
                window_features  = window_features,
                transformer_exog = one_hot_encoder
             )

# 查看训练矩阵
# ==============================================================================
exog_features = ['weather']         
X_train, y_train = forecaster.create_train_X_y(
                        y    = data.loc[:end_validation, 'users'],
                        exog = data.loc[:end_validation, exog_features]
                   )
X_train.head(3)

# 创建带有自动类别检测功能的预测器
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator            = LGBMRegressor(random_state=15926, verbose=-1),
                lags                 = 72,
                categorical_features = 'auto'
             )

# 选择要纳入模型的外生变量
# ==============================================================================
exog_features = []
# 选择以 _sin 或 _cos 结尾的列
exog_features.extend(df_exogenous_features.filter(regex='_sin$|_cos$').columns.tolist())
# 选择以 temp_ 开头的列
exog_features.extend(df_exogenous_features.filter(regex='^temp_.*').columns.tolist())
# 选择以 holiday_ 开头的列
exog_features.extend(df_exogenous_features.filter(regex='^holiday_.*').columns.tolist())
exog_features.extend(['temp', 'holiday', 'weather'])

df_exogenous_features = df_exogenous_features.filter(exog_features, axis=1)

# 将目标变量与外生变量合并到同一数据框
# ==============================================================================
data = data[['users', 'weather']].merge(
            df_exogenous_features,
            left_index=True,
            right_index=True,
            how='inner' # To use only dates for which we have all the variables
        )
data = data.astype({col: np.float32 for col in data.select_dtypes("number").columns})
data_train = data.loc[: end_train, :].copy()
data_val   = data.loc[end_train:end_validation, :].copy()
data_test  = data.loc[end_validation:, :].copy()

# 在测试集上对含外生变量的模型进行回测
# ==============================================================================
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',
                        n_jobs        = 'auto',
                        verbose       = False,
                        show_progress = True
                     )
metric

# 超参数搜索
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator            = LGBMRegressor(random_state=15926, verbose=-1),
                lags                 = 72,
                window_features      = window_features,
                categorical_features = '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'], # Test data not used
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20, # Increase this value for a more exhaustive search
    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

# 绘制预测值与真实值对比图
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Real value vs predicted in test data",
    xaxis_title="Date time",
    yaxis_title="Users",
    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                 = best_lags,
    window_features      = window_features,
    categorical_features = 'auto'
)

# 递归特征消除与交叉验证
# ==============================================================================
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,
    categorical_features = '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                 = best_lags,
    window_features      = window_features,
    categorical_features = '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% prediction interval
   interval_method         = 'conformal',
   use_in_sample_residuals = False,  # Use out-sample residuals
   use_binned_residuals    = True,   # Use residuals conditioned on predicted values
)
predictions.head(5)

# 绘制预测区间与真实值对比图
# ==============================================================================
fig = go.Figure([
    go.Scatter(name='Prediction', x=predictions.index, y=predictions['pred'], mode='lines'),
    go.Scatter(
        name='Real value', x=data_test.index, y=data_test['users'], mode='lines',
    ),
    go.Scatter(
        name='Upper Bound', x=predictions.index, y=predictions['upper_bound'], mode='lines',
        marker=dict(color="#444"), line=dict(width=0), showlegend=False
    ),
    go.Scatter(
        name='Lower Bound', 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="Real value vs predicted in test data",
    xaxis_title="Date time",
    yaxis_title="users",
    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),
    # Initial zoom on x axis betwee 1 oct to 10 oct
    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"Total area of the interval: {round(area, 2)}")
print(f"Predicted interval coverage: {round(100 * coverage, 2)} %")

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

# 提取特征重要性
# ==============================================================================
importance = forecaster.get_feature_importances()
importance.head(10)

# 预测器用于训练内部估计器的训练矩阵
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
                       y    = data.loc[:end_validation, 'users'],
                       exog = data.loc[:end_validation, exog_features]
                   )
display(X_train.head(3))
display(y_train.head(3))

# 创建 SHAP 解释器（适用于树模型）
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.estimator)

# 采样 50% 的数据以加快计算速度
rng = np.random.default_rng(seed=785412)
sample = rng.choice(X_train.index, size=int(len(X_train)*0.5), replace=False)
X_train_sample = X_train.loc[sample, :]
shap_values = explainer.shap_values(X_train_sample)

# SHAP 汇总图（前 10 个特征）
# ==============================================================================
shap.initjs()
shap.summary_plot(shap_values, X_train_sample, max_display=10, show=False)
fig, ax = plt.gcf(), plt.gca()
ax.set_title("SHAP Summary plot")
ax.tick_params(labelsize=8, colors='white')
fig.set_size_inches(10, 4.5)

# 回测，并返回预测所用的特征值
# ==============================================================================
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',
                        return_predictors       = True,
                     )

predictions.head(3)

# 针对单次具体预测的瀑布图
# ==============================================================================
predictions = predictions.astype(data[exog_features].dtypes) # Ensure that the types are the same
iloc_predicted_date = predictions.index.get_loc('2012-10-06 12:00:00')
shap_values_single = explainer(predictions.iloc[:, 2:])
shap.plots.waterfall(shap_values_single[iloc_predicted_date], show=False)
fig, ax = plt.gcf(), plt.gca()
fig.set_size_inches(8, 3.5)
ax_list = fig.axes
ax = ax_list[0]
ax.tick_params(labelsize=8)
ax.set
plt.show()

# 针对单次具体预测的力图
# ==============================================================================
shap.force_plot(
    base_value  = shap_values_single.base_values[iloc_predicted_date],
    shap_values = shap_values_single.values[iloc_predicted_date],
    features    = predictions.iloc[iloc_predicted_date, 2:],
)

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

# 创建预测器
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = XGBRegressor(tree_method='hist', random_state=123),
    lags                 = 24,
    window_features      = window_features,
    categorical_features = '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),
        'learning_rate'   : trial.suggest_float('learning_rate', 0.01, 1),
        'subsample'       : trial.suggest_float('subsample', 0.1, 1),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.1, 1),
        'gamma'           : trial.suggest_float('gamma', 0, 1),
        '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

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
)

# 在测试集上对含外生变量的模型进行回测
# ==============================================================================
metric_xgboost, predictions = backtesting_forecaster(
    forecaster = forecaster,
    y          = data['users'],
    exog       = data[exog_features],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error'
)
metric_xgboost

# 创建预测器
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = HistGradientBoostingRegressor(random_state=123),
    lags                 = 24,
    window_features      = window_features,
    categorical_features = 'auto'
)

# 超参数搜索
# ==============================================================================
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

def search_space(trial):
    search_space  = {
        'max_iter'          : trial.suggest_int('max_iter', 300, 1000, step=100),
        'max_depth'         : trial.suggest_int('max_depth', 3, 10),
        'learning_rate'     : trial.suggest_float('learning_rate', 0.01, 1),
        'min_samples_leaf'  : trial.suggest_int('min_samples_leaf', 1, 20),
        'l2_regularization' : trial.suggest_float('l2_regularization', 0, 1),
        'lags'              : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

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
)

# 在测试集上对含外生变量的模型进行回测
# ==============================================================================
metric_histgb, predictions = backtesting_forecaster(
    forecaster = forecaster,
    y          = data['users'],
    exog       = data[exog_features],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error'
)
metric_histgb

# 创建预测器
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = CatBoostRegressor(
                                random_state=123, 
                                silent=True, 
                                allow_writing_files=False,
                                boosting_type = 'Plain',         # 更快的训练速度
                                leaf_estimation_iterations = 3,  # 更快的训练速度
                            ),
    lags                 = 24,
    window_features      = window_features,
    categorical_features = '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', 100, 1000, step=100),
        'max_depth'     : trial.suggest_int('max_depth', 3, 10),
        'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
        'lags'          : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

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
)

# 在测试集上进行回测
# ==============================================================================
metric_catboost, predictions = backtesting_forecaster(
    forecaster = forecaster,
    y          = data['users'],
    exog       = data[exog_features],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error'
)
metric_catboost

metrics = pd.concat(
    [metric_baseline, metric_lgbm, metric_xgboost, metric_histgb, metric_catboost]
)
metrics.index = [
    "Baseline",
    "LGBMRegressor",
    "XGBRegressor",
    "HistGradientBoostingRegressor",
    "CatBoostRegressor",
]
metrics.round(2).sort_values(by="mean_absolute_error")

import session_info
session_info.show(html=False)