分布变化下的自适应共形推断

自适应共形推理(ACI)[Gibbs 等人,2021] 是一种旨在在数据分布随时间变化时(例如时间序列预测)在顺序预测框架中纠正共形预测方法覆盖范围的方法。ACI 可以应用于任何共形预测方法之上,适用于回归和分类。在本笔记本中,我们展示了它与共形化分位数回归(CQR)[Romano 等人,2019] 结合用于股票价格趋势估计任务。

下载并准备数据

我们下载了一些大型科技公司过去一年的每日调整后收盘股票价格。我们将每个时间序列数据依次分为训练、验证和测试数据集。

[1]:
import yfinance as yf

df = yf.download(
    tickers="AMZN AAPL GOOGL NFLX MSFT",
    period="1y",
    interval="1d",
    ignore_tz=True,
    prepost=False,
)
tickers = df["Adj Close"].columns
dates = df.index
y_all = df["Adj Close"].to_numpy().T
y_all /= y_all.max(1)[:, None]
[*********************100%%**********************]  5 of 5 completed
[2]:
import matplotlib.pyplot as plt

n_tickers = len(tickers)
plt.figure(figsize=(12, 1.5 * n_tickers))
for i in range(n_tickers):
    plt.subplot(n_tickers, 1, i + 1)
    plt.plot(df.index, y_all[i])
    plt.title(tickers[i], fontsize=12)
plt.tight_layout()
../_images/examples_adaptive_conformal_inference_5_0.svg
[3]:
import numpy as np

period = 60
y_train = y_all[:, : -2 * period, None]
y_val = y_all[:, -2 * period : -period, None]
y_test = y_all[:, -period:, None]
n_data = y_all.shape[1]
x_all = np.broadcast_to(np.linspace(0, 1, n_data), shape=y_all.shape)
x_train = x_all[:, : -2 * period, None]
x_val = x_all[:, -2 * period : -period, None]
x_test = x_all[:, -period:, None]

建模和训练

为了捕捉数据中的总体趋势,我们为每个股票价格的时间序列在时间步长上构建了一个简单的二次模型。

[4]:
import flax.linen as nn


class QuadraticModel(nn.Module):
    output_dim: int

    @nn.compact
    def __call__(self, x, **kwargs):
        params = self.param("coeffs", nn.initializers.zeros, (3,))
        return params[0] + params[1] * x + params[2] * x**2

然后我们使用Fortuna从二次模型开始构建概率模型,为每个股票价格时间序列训练模型,进行预测并估计可信区间。

[5]:
from fortuna.data import DataLoader
from fortuna.prob_model import ProbRegressor
from fortuna.prob_model import MAPPosteriorApproximator
from fortuna.prob_model import FitConfig, FitOptimizer, FitMonitor
from fortuna.model import ScalarConstantModel
import optax

target_error = 0.05

n_train = y_train.shape[1]
val_cred_intervals = np.zeros((n_tickers, period, 2))
test_cred_intervals = np.zeros((n_tickers, period, 2))
test_conformal_intervals = np.zeros((n_tickers, period, 2))
train_preds, val_preds, test_preds = (
    np.zeros((n_tickers, n_train, 1)),
    np.zeros((n_tickers, period, 1)),
    np.zeros((n_tickers, period, 1)),
)
statuses = []

for i in range(n_tickers):
    print(f"Ticker: {tickers[i]}.")
    train_data_loader = DataLoader.from_array_data(
        (x_train[i], y_train[i]), batch_size=128, prefetch=True
    )
    val_data_loader = DataLoader.from_array_data(
        (x_val[i], y_val[i]), batch_size=128, prefetch=True
    )
    test_data_loader = DataLoader.from_array_data(
        (x_test[i], y_test[i]), batch_size=128, prefetch=True
    )

    prob_model = ProbRegressor(
        model=QuadraticModel(1),
        likelihood_log_variance_model=ScalarConstantModel(output_dim=1),
        posterior_approximator=MAPPosteriorApproximator(),
    )

    statuses.append(
        prob_model.train(
            train_data_loader=train_data_loader,
            fit_config=FitConfig(
                optimizer=FitOptimizer(method=optax.adam(1e-1), n_epochs=1000),
                monitor=FitMonitor(early_stopping_patience=2),
            ),
        )
    )

    train_preds[i] = prob_model.predictive.mean(
        inputs_loader=train_data_loader.to_inputs_loader()
    )
    val_preds[i] = prob_model.predictive.mean(
        inputs_loader=val_data_loader.to_inputs_loader()
    )
    test_preds[i] = prob_model.predictive.mean(
        inputs_loader=test_data_loader.to_inputs_loader()
    )
    val_cred_intervals[i] = prob_model.predictive.credible_interval(
        inputs_loader=val_data_loader.to_inputs_loader(), error=target_error
    )
    test_cred_intervals[i] = prob_model.predictive.credible_interval(
        inputs_loader=test_data_loader.to_inputs_loader(), error=target_error
    )
Ticker: AAPL.
Epoch: 1000 | loss: -329.26602: 100%|██████████| 1000/1000 [00:04<00:00, 224.31it/s]
Ticker: AMZN.
Epoch: 1000 | loss: -356.22971: 100%|██████████| 1000/1000 [00:04<00:00, 231.92it/s]
Ticker: GOOGL.
Epoch: 1000 | loss: -328.9819: 100%|██████████| 1000/1000 [00:04<00:00, 230.54it/s]
Ticker: MSFT.
Epoch: 1000 | loss: -230.51251: 100%|██████████| 1000/1000 [00:04<00:00, 233.62it/s]
Ticker: NFLX.
Epoch: 1000 | loss: -292.70657: 100%|██████████| 1000/1000 [00:04<00:00, 234.52it/s]
[6]:
plt.figure(figsize=(10, 2))
for i in range(n_tickers):
    plt.plot(statuses[i]["fit_status"]["loss"])
plt.legend(tickers)
plt.title("Loss decays")
[6]:
Text(0.5, 1.0, 'Loss decays')
../_images/examples_adaptive_conformal_inference_12_1.svg

顺序符合区间

随着测试数据点被顺序观察,我们希望为下一个目标变量构建保形区间。为此,我们使用最初用于验证的所有数据的可信区间估计作为验证数据,再加上已经观察到的测试数据点的可信区间估计。我们通过两种方式估计保形区间:使用Romano等人,2019年的保形化分位数回归(CQR),以及使用Gibbs等人,2021年的自适应保形推断(ACI)提供的额外覆盖误差更新功能的CQR。

[7]:
from fortuna.conformal import (
    QuantileConformalRegressor,
    AdaptiveConformalRegressor,
)

errors = target_error * np.ones((n_tickers, period + 1))
is_ins = np.zeros((n_tickers, period))
conformal_intervals = np.zeros((n_tickers, period, 2))
adaptive_conformal_intervals = np.zeros((n_tickers, period, 2))
qcr = QuantileConformalRegressor()
aqcr = AdaptiveConformalRegressor(conformal_regressor=qcr)

for i in range(n_tickers):
    for j in range(period):
        if j == 0:
            val_lower_bounds, val_upper_bounds = val_cred_intervals[i].T
            val_targets = np.copy(y_val[i])
        else:
            val_lower_bounds = np.concatenate(
                (val_cred_intervals[i, :, 0], test_cred_intervals[i, :j, 0])
            )
            val_upper_bounds = np.concatenate(
                (val_cred_intervals[i, :, 1], test_cred_intervals[i, :j, 1])
            )
            val_targets = np.concatenate((y_val[i, :], y_test[i, :j]))

        conformal_intervals[i, j] = qcr.conformal_interval(
            val_lower_bounds=val_lower_bounds,
            val_upper_bounds=val_upper_bounds,
            test_lower_bounds=test_cred_intervals[None, i, j, 0],
            test_upper_bounds=test_cred_intervals[None, i, j, 1],
            val_targets=val_targets,
            error=target_error,
        )[0]

        adaptive_conformal_intervals[i, j] = aqcr.conformal_interval(
            val_lower_bounds=val_lower_bounds,
            val_upper_bounds=val_upper_bounds,
            test_lower_bounds=test_cred_intervals[None, i, j, 0],
            test_upper_bounds=test_cred_intervals[None, i, j, 1],
            val_targets=val_targets,
            error=errors[i, j],
        )[0]

        errors[i, j + 1] = aqcr.update_error(
            conformal_interval=adaptive_conformal_intervals[i, j],
            error=errors[i, j],
            target=y_test[i, j],
            target_error=target_error,
        )

以下图表显示了针对每家公司股票价格的训练、验证和测试数据集的股票趋势预测。在左侧列中,我们可视化了通过CQR获得的95%置信区间;在右侧列中,我们看到了使用CQR与ACI结合计算的针对95%覆盖率的置信区间。

[8]:
plt.figure(figsize=(12, 2 * n_tickers))
for i in range(n_tickers):
    plt.subplot(n_tickers, 2, 2 * i + 1)
    plt.plot(df.index, y_all[i], label="data")
    plt.plot(
        df.index[:n_train],
        train_preds[i],
        color="C1",
        label="training prediction",
    )
    plt.plot(
        df.index[n_train : n_train + period],
        val_preds[i],
        color="C2",
        label="validation prediction",
    )
    plt.plot(
        df.index[n_train + period :],
        test_preds[i],
        color="C3",
        label="test prediction",
    )
    plt.fill_between(
        df.index[n_train + period :],
        *conformal_intervals[i].T,
        alpha=0.3,
        color="C0",
        label="conformal interval",
    )
    plt.title(tickers[i], fontsize=12)
    if i == 0:
        plt.legend(loc="upper left")

    plt.subplot(n_tickers, 2, 2 * i + 2)
    plt.plot(df.index, y_all[i], label="data")
    plt.plot(
        df.index[:n_train],
        train_preds[i],
        color="C1",
        label="training prediction",
    )
    plt.plot(
        df.index[n_train : n_train + period],
        val_preds[i],
        color="C2",
        label="validation prediction",
    )
    plt.plot(
        df.index[n_train + period :],
        test_preds[i],
        color="C3",
        label="test prediction",
    )
    plt.fill_between(
        df.index[n_train + period :],
        *adaptive_conformal_intervals[i].T,
        alpha=0.3,
        color="C0",
        label="adaptive conformal interval",
    )
    plt.title(tickers[i], fontsize=12)
    if i == 0:
        plt.legend(loc="upper left")
plt.tight_layout()
../_images/examples_adaptive_conformal_inference_17_0.svg

下表显示了每个股票价格时间序列的两种保形区间估计的覆盖率估计。由于股票价格分布随时间变化,使用CQR估计的保形区间的有效性无法保证,因此得出的覆盖率估计往往远低于期望的95%。ACI对此进行了改进,使其更接近期望的覆盖率。

[9]:
from fortuna.metric.regression import prediction_interval_coverage_probability
from tabulate import tabulate

table = [
    [
        "",
        "estimated coverage of conformal intervals",
        "estimated coverage of adaptive conformal intervals",
    ]
]
for i in range(n_tickers):
    table.append(
        [
            tickers[i],
            float(
                prediction_interval_coverage_probability(
                    *conformal_intervals[i].T, y_test[i]
                )
            ),
            float(
                prediction_interval_coverage_probability(
                    *adaptive_conformal_intervals[i].T, y_test[i]
                )
            ),
        ]
    )
print(tabulate(table, headers="firstrow", tablefmt="fancy_grid"))
╒═══════╤═════════════════════════════════════════════╤══════════════════════════════════════════════════════╕
│       │   estimated coverage of conformal intervals │   estimated coverage of adaptive conformal intervals │
╞═══════╪═════════════════════════════════════════════╪══════════════════════════════════════════════════════╡
│ AAPL  │                                    0.133333 │                                             0.5      │
├───────┼─────────────────────────────────────────────┼──────────────────────────────────────────────────────┤
│ AMZN  │                                    0.283333 │                                             0.65     │
├───────┼─────────────────────────────────────────────┼──────────────────────────────────────────────────────┤
│ GOOGL │                                    0.383333 │                                             0.683333 │
├───────┼─────────────────────────────────────────────┼──────────────────────────────────────────────────────┤
│ MSFT  │                                    0.7      │                                             0.816667 │
├───────┼─────────────────────────────────────────────┼──────────────────────────────────────────────────────┤
│ NFLX  │                                    0.316667 │                                             0.683333 │
╘═══════╧═════════════════════════════════════════════╧══════════════════════════════════════════════════════╛
[ ]: