# 数据处理
# ==============================================================================
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset

# 绘图
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from skforecast.plot import plot_residuals
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')
plt.rcParams.update({'font.size': 8})

# 建模与预测
# ==============================================================================
import skforecast
import lightgbm
import sklearn
from lightgbm import LGBMRegressor
from sklearn.preprocessing import PolynomialFeatures
from sklearn.feature_selection import RFECV
from feature_engine.timeseries.forecasting import WindowFeatures
from skforecast.preprocessing import CalendarFeatures, RollingFeatures
from skforecast.recursive import ForecasterEquivalentDate, ForecasterRecursive
from skforecast.direct import ForecasterDirect
from skforecast.model_selection import (
    TimeSeriesFold,
    bayesian_search_forecaster,
    backtesting_forecaster
)
from skforecast.feature_selection import select_features
from skforecast.stats import calculate_lag_autocorrelation
from skforecast.metrics import calculate_coverage
import shap

# 警告配置
# ==============================================================================
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 pandas: {pd.__version__}')
print(f'{color}Version numpy: {np.__version__}')

# 数据下载
# ==============================================================================
data = fetch_dataset(name='vic_electricity', raw=True)
data.info()

# 数据预处理
# ==============================================================================
data['Time'] = pd.to_datetime(data['Time'], utc=True)
data = data.set_index('Time')
data = data.tz_convert('Australia/Melbourne')
data = data.sort_index()
data = data.asfreq('30min')
data.head(2)

# 验证时间索引是否完整
# ==============================================================================
start_date = data.index.min()
end_date = data.index.max()
complete_date_range = pd.date_range(
    start=start_date, end=end_date, freq=data.index.freq
)
is_index_complete = data.index.equals(complete_date_range)
print(f'Index complete: {is_index_complete}')
print(f'Number of rows with missing values: {data.isnull().any(axis=1).sum()}')

# 填补时间索引中的缺口
# ==============================================================================
# data.asfreq(freq='30min', fill_value=np.nan)

# 以 1 小时为间隔进行聚合
# ==============================================================================
# 删除 Date 列，避免聚合时报错。
data = data.drop(columns='Date')
data = (
    data
    .resample(rule='h', closed='left', label='right')
    .agg({
        'Demand': 'mean',
        'Temperature': 'mean',
        'Holiday': 'mean',
    })
)
data

# 划分数据为 train-val-test
# ==============================================================================
data = data.loc[:'2014-12-31 23:00:00', :].copy()
end_train = '2013-12-31 23:59:00'
end_validation = '2014-09-30 23:59:00'
data_train = data.loc[: end_train, :].copy()
data_val   = data.loc[end_train:end_validation, :].copy()
data_test  = data.loc[end_validation:, :].copy()

print(
    f'Train dates      : {data_train.index.min()} --- {data_train.index.max()}  '
    f'(n={len(data_train)})'
)
print(
    f'Validation dates : {data_val.index.min()} --- {data_val.index.max()}  '
    f'(n={len(data_val)})'
)
print(
    f'Test dates       : {data_test.index.min()} --- {data_test.index.max()}  '
    f'(n={len(data_test)})'
)

# 时间序列交互式绘图
# ==============================================================================
fig = go.Figure()
for partition, name in zip(
    [data_train, data_val, data_test], ['Train', 'Validation', 'Test']
):
    fig.add_trace(
        go.Scatter(x=partition.index, y=partition['Demand'], mode='lines', name=name)
    )
fig.update_layout(
    title='Hourly energy demand',
    xaxis_title='Time',
    yaxis_title='Demand',
    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 = ('2013-05-01 14:00:00','2013-06-01 14:00:00')
fig, axs = plt.subplots(2, 1, figsize=(8, 4), gridspec_kw={'height_ratios': [1, 2]})
data['Demand'].plot(ax=axs[0], color='black', alpha=0.5)
axs[0].axvspan(zoom[0], zoom[1], color='blue', alpha=0.7)
axs[0].set_title('Electricity demand')
axs[0].set_xlabel('')
data.loc[zoom[0] : zoom[1], 'Demand'].plot(ax=axs[1], color='blue')
axs[1].set_title(f'Zoom: {zoom[0]} to {zoom[1]}', fontsize=10)
plt.tight_layout()
plt.show()

# 年度、每周和每日季节性
# ==============================================================================
fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex=False, sharey=True)
axs = axs.ravel()
flierprops = {'markersize': 3, 'alpha': 0.3}

# 按月份统计需求分布
data['month'] = data.index.month
data.boxplot(column='Demand', by='month', ax=axs[0], flierprops=flierprops)
data.groupby('month')['Demand'].median().plot(style='o-', linewidth=0.8, ax=axs[0])
axs[0].set_ylabel('Demand')
axs[0].set_title('Demand distribution by month', fontsize=9)

# 按星期几统计需求分布（1 = 周一）
data['week_day'] = data.index.day_of_week + 1
data.boxplot(column='Demand', by='week_day', ax=axs[1], flierprops=flierprops)
data.groupby('week_day')['Demand'].median().plot(style='o-', linewidth=0.8, ax=axs[1])
axs[1].set_ylabel('Demand')
axs[1].set_title('Demand distribution by week day', fontsize=9)

# 按一天中的小时统计需求分布（0 至 23，当地时间）
data['hour_day'] = data.index.hour
data.boxplot(column='Demand', by='hour_day', ax=axs[2], flierprops=flierprops)
# 箱线图绘制在位置 1 至 24，因此中位数也绘制在相同位置
median_hour = data.groupby('hour_day')['Demand'].median()
axs[2].plot(range(1, 25), median_hour.to_numpy(), 'o-', linewidth=0.8)
axs[2].set_ylabel('Demand')
axs[2].set_title('Demand distribution by the hour of the day', fontsize=9)

# 按星期几和小时统计需求分布
mean_day_hour = data.groupby(['week_day', 'hour_day'])['Demand'].mean()
mean_day_hour.plot(ax=axs[3])
axs[3].set(
    title       = 'Mean Demand during week',
    xticks      = [i * 24 for i in range(7)],
    xticklabels = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun'],
    xlabel      = 'Day and hour',
    ylabel      = 'Mean Demand'
)
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['Demand'], ax=ax, lags=60, fft=True)
plt.show()

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

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

# 创建基线模型：前一天同一小时的值
# ==============================================================================
# 偏移量以步数表示（24 小时）。对于带时区信息的索引，
# pd.DateOffset(days=1) 这类日历偏移量在夏令时切换时会失败，
# 因为某些当地时间并不存在。
forecaster = ForecasterEquivalentDate(
                 offset    = 24,
                 n_offsets = 1
             )

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

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

# 创建预测器
# ==============================================================================
# 滞后：过去 24 小时的需求
lags = 24
window_features = RollingFeatures(stats=['mean'], window_sizes=24 * 3)
forecaster = ForecasterRecursive(
                 estimator       = LGBMRegressor(random_state=15926, verbose=-1),
                 lags            = lags,
                 window_features = window_features
             )

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

# 回测
# ==============================================================================
metric, predictions = backtesting_forecaster(
                          forecaster    = forecaster,
                          y             = data['Demand'],
                          cv            = cv,
                          metric        = 'mean_absolute_error',
                          verbose       = True,  # 设为 False 可避免打印
                      )
metric_recursive_no_exog = metric

# 绘制预测值与真实值对比图
# ======================================================================================
fig = go.Figure()
fig.add_trace(
    go.Scatter(x=data_test.index, y=data_test['Demand'], name='test', mode='lines')
)
fig.add_trace(
    go.Scatter(
        x=predictions.index, y=predictions['pred'], name='prediction', mode='lines'
    )
)
fig.update_layout(
    title='Real value vs predicted in test data',
    xaxis_title='Date time',
    yaxis_title='Demand',
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    legend=dict(orientation='h', yanchor='top', y=1.01, xanchor='left', x=0)
)
fig.show()

# 回测误差
# ==============================================================================
metric_recursive_no_exog

# 日历特征
# ==============================================================================
calendar_transformer = CalendarFeatures(
    features = ['month', 'week', 'day_of_week', 'hour'],
    encoding = 'cyclical',
    keep_original_columns = False,
)
calendar_features = calendar_transformer.fit_transform(data)
calendar_features.head(2)

# 日照特征
# ==============================================================================
location = LocationInfo(
    latitude  = -37.8,
    longitude = 144.95,
    timezone  = 'Australia/Melbourne'
)
# 日出和日落只取决于日期，因此每天只计算一次
dates = pd.Series(data.index.date, index=data.index)
sun_by_date = {
    date: sun(location.observer, date=date, tzinfo=location.timezone)
    for date in dates.unique()
}
sunrise_hour = dates.map(lambda date: sun_by_date[date]['sunrise'])
sunset_hour = dates.map(lambda date: sun_by_date[date]['sunset'])
sunrise_hour = sunrise_hour.dt.round('h').dt.hour
sunset_hour = sunset_hour.dt.round('h').dt.hour
sunrise_hour_sin = np.sin(2 * np.pi * sunrise_hour / 24)
sunrise_hour_cos = np.cos(2 * np.pi * sunrise_hour / 24)
sunset_hour_sin = np.sin(2 * np.pi * sunset_hour / 24)
sunset_hour_cos = np.cos(2 * np.pi * sunset_hour / 24)

daylight_hours = sunset_hour - sunrise_hour
# 索引以及日出和日落的小时均为墨尔本当地时间
is_daylight = np.where(
    (data.index.hour >= sunrise_hour) & (data.index.hour < sunset_hour), 1, 0,
)

sun_light_features = pd.DataFrame({
                        'sunrise_hour_sin': sunrise_hour_sin,
                        'sunrise_hour_cos': sunrise_hour_cos,
                        'sunset_hour_sin': sunset_hour_sin,
                        'sunset_hour_cos': sunset_hour_cos,
                        'daylight_hours': daylight_hours,
                        'is_daylight': is_daylight
                     })
sun_light_features.head(2)

# 节假日特征
# ==============================================================================
# 公共假日是提前已知的，因此将次日的值用作特征并不构成数据泄露。
# 平移量以步数表示（24 小时），因此在每年夏令时切换的两天会相差一个小时。
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)
holiday_features.head(2)

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

# 合并所有外生变量
# ==============================================================================
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)
exogenous_features = pd.concat([
                         calendar_features,
                         sun_light_features,
                         temp_features,
                         holiday_features
                     ], axis=1)

# 由于创建了滑动平均，序列开头会有缺失值。
# 而 holiday_next_day 会导致序列结尾也有缺失值。
exogenous_features = exogenous_features.iloc[7 * 24:, :]
exogenous_features = exogenous_features.iloc[:-24, :]
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',
    'Temperature_window_1D_mean',
    'Temperature_window_1D_min',
    'Temperature_window_1D_max',
    'Temperature_window_7D_mean',
    'Temperature_window_7D_min',
    'Temperature_window_7D_max',
    'Temperature',
    'Holiday'
]
poly_features = transformer_poly.fit_transform(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 == exogenous_features.index)
exogenous_features = pd.concat([exogenous_features, poly_features], axis=1)
exogenous_features.head(3)

# 选择要纳入模型的外生变量
# ==============================================================================
exog_features = []
# 选择以 _sin 或 _cos 结尾的列（周期性特征及其交互项）
exog_features.extend(
    exogenous_features.filter(regex='_sin$|_cos$').columns.tolist()
)
# 选择以 Temperature_ 开头的列
exog_features.extend(
    exogenous_features.filter(regex='^Temperature_.*').columns.tolist()
)
# 选择以 holiday_ 开头的列
exog_features.extend(
    exogenous_features.filter(regex='^holiday_.*').columns.tolist()
)
# 包含原始特征
exog_features.extend(['Temperature', 'Holiday', 'daylight_hours', 'is_daylight'])

# 合并目标变量和外生变量到同一 DataFrame
# ==============================================================================
data = data[['Demand']].merge(
           exogenous_features[exog_features],
           left_index  = True,
           right_index = True,
           how         = 'inner'  # 仅使用所有变量均有值的日期
       )
data = data.astype('float32')

# 划分数据为 train-val-test
data_train = data.loc[: end_train, :].copy()
data_val   = data.loc[end_train:end_validation, :].copy()
data_test  = data.loc[end_validation:, :].copy()

# 回测模型
# ==============================================================================
# 由于 `data` 的长度发生了变化，需要重新创建折
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_validation]),
        refit              = False
)

metric, predictions = backtesting_forecaster(
                          forecaster = forecaster,
                          y          = data['Demand'],
                          exog       = data[exog_features],
                          cv         = cv,
                          metric     = 'mean_absolute_error'
                      )
metric_recursive_exog = metric
display(metric_recursive_exog)
predictions.head()

# 超参数搜索
# ==============================================================================
forecaster = ForecasterRecursive(
                 estimator       = LGBMRegressor(random_state=15926, verbose=-1),
                 lags            = 24,  # 该值会在搜索中被替换
                 window_features = window_features
             )

# 用作特征的滞后
lags_grid = [
    # 前一天
    24,
    # 前两天
    48,
    # 前一天以及前一周的同一小时（+-1）
    list(range(1, 25)) + [167, 168, 169],
    # 最近几小时以及前两天的同一小时（+-1）
    [1, 2, 3, 23, 24, 25, 47, 48, 49],
    # 同上，再加上前一周的同一小时（+-1）
    [1, 2, 3, 23, 24, 25, 47, 48, 49, 167, 168, 169],
]

# 估计器超参数搜索空间
def search_space(trial):
    params = {
        '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, 0.5),
        '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 params

# 训练和验证折
cv_search = TimeSeriesFold(
                steps              = 24,
                initial_train_size = len(data.loc[:end_train]),
                refit              = False,
            )

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster   = forecaster,
    y            = data.loc[:end_validation, 'Demand'],
    exog         = data.loc[:end_validation, exog_features],
    cv           = cv_search,
    metric       = 'mean_absolute_error',
    search_space = search_space,
    n_trials     = 30,  # 增大以获得更全面搜索
    return_best  = True
)

# 搜索结果
# ==============================================================================
best_params = results_search.at[0, 'params']
best_params = best_params | {'random_state': 15926, 'verbose': -1}
best_lags = results_search.at[0, 'lags']
results_search.head(3)

# 最优模型
# ==============================================================================
forecaster

# 在测试集上回测最终模型
# ==============================================================================
metric, predictions = backtesting_forecaster(
                          forecaster = forecaster,
                          y          = data['Demand'],
                          exog       = data[exog_features],
                          cv         = cv,
                          metric     = 'mean_absolute_error'
                      )
metric_recursive_exog_tuned = metric
display(metric_recursive_exog_tuned)
predictions.head()

# 创建预测器
# ==============================================================================
estimator = LGBMRegressor(
                n_estimators = 100,
                max_depth    = 4,
                random_state = 15926,
                verbose      = -1
            )

forecaster = ForecasterRecursive(
                 estimator       = estimator,
                 lags            = best_lags,
                 window_features = window_features
             )

# 递归特征消除与交叉验证
# ==============================================================================
warnings.filterwarnings('ignore', message='X does not have valid feature names.*')
selector = RFECV(
    estimator = estimator,
    step      = 1,
    cv        = 3,
)
lags_select, window_features_select, exog_select, _ = select_features(
    forecaster      = forecaster,
    selector        = selector,
    y               = data_train['Demand'],
    exog            = data_train[exog_features],
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.5,  # 子采样以加快计算
    random_state    = 123,
    verbose         = True,
)

# 仅保留选定的窗口特征
# ==============================================================================
# `select_features` 返回选定窗口特征的名称。
# 仅用选定的统计量重新创建 RollingFeatures 对象。
# 如果没有选中任何窗口特征，则返回 None（不含窗口特征的预测器）。
def filter_rolling_features(window_features, selected_names):
    keep = [
        i for i, name in enumerate(window_features.features_names)
        if name in selected_names
    ]
    if not keep:
        return None

    return RollingFeatures(
        stats          = [window_features.stats[i] for i in keep],
        window_sizes   = [window_features.window_sizes[i] for i in keep],
        min_periods    = [window_features.min_periods[i] for i in keep],
        features_names = [window_features.features_names[i] for i in keep],
        fillna         = window_features.fillna,
        kwargs_stats   = window_features.kwargs_stats,
    )

window_features_select = filter_rolling_features(
    window_features, window_features_select
)
print(f'Selected window features: {window_features_select}')

# 用选定特征创建预测器
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator       = LGBMRegressor(**best_params),
                lags            = lags_select,
                window_features = window_features_select
             )
# 在测试数据上回测包含外生变量的模型
# ==============================================================================
metric, predictions = backtesting_forecaster(
                            forecaster = forecaster,
                            y          = data['Demand'],
                            exog       = data[exog_select],
                            cv         = cv,
                            metric     = 'mean_absolute_error'
                      )
metric_recursive_exog_selection = metric
display(metric_recursive_exog_selection)
predictions.head()

# 创建并训练预测器
# ==============================================================================
forecaster = ForecasterRecursive(
                 estimator       = LGBMRegressor(**best_params),
                 lags            = lags_select,
                 window_features = window_features_select,
                 binner_kwargs   = {'n_bins': 5}
             )
forecaster.fit(
    y    = data.loc[:end_train, 'Demand'],
    exog = data.loc[:end_train, exog_select],
    store_in_sample_residuals = True
)

# 预测区间
# ==============================================================================
# 由于模型是用外生变量训练的，预测时必须提供这些变量。
predictions = forecaster.predict_interval(
                  exog     = data.loc[end_train:, exog_select],
                  steps    = 24,
                  interval = [0.05, 0.95],
                  method   = 'conformal',
              )
predictions.head()

# 在验证数据上进行回测以获得样本外残差
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_train]),
        refit              = False,
)
metric_val, predictions_val = backtesting_forecaster(
                                  forecaster = forecaster,
                                  y          = data.loc[:end_validation, 'Demand'],
                                  exog       = data.loc[:end_validation, exog_select],
                                  cv         = cv,
                                  metric     = 'mean_absolute_error'
                              )
metric_val

# 样本外残差分布
# ==============================================================================
residuals = data.loc[predictions_val.index, 'Demand'] - predictions_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())
_ = plot_residuals(residuals=residuals, figsize=(7, 4))

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

# 在测试数据中使用样本外残差进行区间预测的回测
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_validation]),
        refit              = False,
)
metric, predictions = backtesting_forecaster(
    forecaster              = forecaster,
    y                       = data['Demand'],
    exog                    = data[exog_select],
    cv                      = cv,
    metric                  = 'mean_absolute_error',
    interval                = [0.05, 0.95],
    interval_method         = 'conformal',
    use_in_sample_residuals = False,  # 使用样本外残差
    use_binned_residuals    = True,   # 区间根据预测值的范围进行调整
)
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['Demand'], 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='Demand',
    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)
)
fig.show()

# 预测区间覆盖率（测试数据）
# ==============================================================================
coverage = calculate_coverage(
              y_true       = data.loc[predictions.index, 'Demand'],
              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            = lags_select,
                 window_features = window_features_select
             )
forecaster.fit(
    y    = data.loc[:end_validation, 'Demand'],
    exog = data.loc[:end_validation, exog_select]
)

# 模型特定的特征重要性
# ==============================================================================
feature_importances = forecaster.get_feature_importances()
feature_importances.head(10)

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

# 创建 SHAP 解释器
# ==============================================================================
shap.initjs()
explainer = shap.TreeExplainer(forecaster.estimator)

# 抽取 50% 的数据以加快计算
X_train_sample = X_train.sample(frac=0.5, random_state=785412)
shap_values = explainer.shap_values(X_train_sample)

# SHAP 摘要图（前 10 个特征）
# ==============================================================================
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)
fig.set_size_inches(6, 4.5)

# 回测并返回特征
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_validation]),
        refit              = False,
)
_, predictions = backtesting_forecaster(
                        forecaster        = forecaster,
                        y                 = data['Demand'],
                        exog              = data[exog_select],
                        cv                = cv,
                        metric            = 'mean_absolute_error',
                        return_predictors = True,
                )

# 预测值和特征
# ==============================================================================
predictions.head(3)

# 回测期间生成的单个预测的瀑布图
# ==============================================================================
predictors = predictions.drop(columns=['fold', 'pred'])
# 确保类型与训练矩阵中的一致
predictors = predictors.astype(data[exog_select].dtypes)
iloc_predicted_date = predictions.index.get_loc('2014-12-16 12:00:00')
shap_values_single = explainer(predictors)
shap.plots.waterfall(shap_values_single[iloc_predicted_date], show=False)
fig = plt.gcf()
fig.set_size_inches(8, 3.5)
fig.axes[0].tick_params(labelsize=8)
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    = predictors.iloc[iloc_predicted_date, :],
)

# 直接法预测器
# ==============================================================================
forecaster = ForecasterDirect(
                 estimator       = LGBMRegressor(**best_params),
                 steps           = 24,
                 lags            = lags_select,
                 window_features = window_features_select
             )

# 回测模型
# ==============================================================================
metric, predictions = backtesting_forecaster(
                          forecaster = forecaster,
                          y          = data['Demand'],
                          exog       = data[exog_select],
                          cv         = cv,
                          metric     = 'mean_absolute_error'
                      )
metric_direct_exog_selection = metric
display(metric_direct_exog_selection)
predictions.head()

# initial_train_size + 12 位置的末尾
# ==============================================================================
data.iloc[:len(data.loc[:end_validation]) + 12].tail(2)

# 预测器
# ==============================================================================
forecaster = ForecasterRecursive(
                 estimator       = LGBMRegressor(**best_params),
                 lags            = lags_select,
                 window_features = window_features_select
             )

# 带间隔的回测
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_validation]) + 12,
        refit              = False,
        gap                = 12,
)

metric, predictions = backtesting_forecaster(
                          forecaster = forecaster,
                          y          = data['Demand'],
                          exog       = data[exog_select],
                          cv         = cv,
                          metric     = 'mean_absolute_error'
                      )
display(metric)
predictions.head(5)

# 结果
# ======================================================================================
forecaster_type = [
    'ForecasterEquivalentDate (Baseline)', 'ForecasterRecursive',
    'ForecasterRecursive', 'ForecasterRecursive', 'ForecasterRecursive',
    'ForecasterDirect'
]
exog_included = [
    'False', 'False', 'True', 'True', 'True (Feature Selection)',
    'True (Feature Selection)'
]
tuned = ['False', 'False', 'False', 'True', 'True', 'True']
metrics = pd.concat(
    [
        metric_baseline, metric_recursive_no_exog, metric_recursive_exog,
        metric_recursive_exog_tuned, metric_recursive_exog_selection,
        metric_direct_exog_selection
    ],
    axis=0,
)
metrics.insert(0, 'Forecaster', forecaster_type)
metrics.insert(1, 'Exogenous Variables', exog_included)
metrics.insert(2, 'Tuned Hyperparameters', tuned)
metrics = (
    metrics.reset_index(drop=True).round(2).sort_values(by='mean_absolute_error')
)
metrics

import session_info
session_info.show(html=False)