猜涨跌预测(三):置信区间

之前 猜涨跌预测(二):简化版 中给出的是上证指数上涨概率的点估计,如果能以区间的形式给出上涨概率的估计,在这个场景中会更加的实用,因为不同的人对风险的偏好是不一致的(保守的投资者希望获得稳定可靠的收益,激进的投资者希望获得丰厚的回报)。

置信区间的定义

置信区间(Confidence interval)是一个十分古典的概念(频率学派)。95% 置信区间的古典定义为(按照我的理解):上帝可以知道统计模型里参数的真值(true value),根据此模型随机生成 NN 组样本数据,通过一些估计方法算出关于某个参数的 NN 个不同的置信区间,然后有 95% 的置信区间包含该参数真值(真值只有上帝知道)。换种说法,收集一组样本数据后,通过统计推断得到的 95% 置信区间包含了该参数真值的概率为 95%。

实际应用

计算置信区间对我们日常的决策很有帮助。举个例子,如果统计模型给出今天猜涨跌上证指数上涨的概率为 60%,95% 置信区间为 [50%, 70%],以及猜涨的赔率为 100 支付宝积分赔付 100 支付宝积分。那么一个保守的投资者可能做出如下决策:通过 95% 的置信区间构造今天的收益区间:[200×0.5100,200×0.7100][200\times 0.5-100, 200\times 0.7-100] = [0,40][0, 40],因为该投资者的策略为保守型,因此他可能做出今天不参与猜涨跌的决策(极端厌恶损失)。但是一个激进的投资者可能做出不同的决策:今天的期望收益为:200×0.6100=20200\times 0.6-100=20,而且 95% 的置信区间的下界都大于等于 50%,(Loosely speaking)因此 95% 概率下,今天猜上涨能获得超额利润,果断压上 100 支付宝积分赌今日上证指数会上涨。 可见计算上涨概率的置信区间相比于仅仅计算上涨概率(点估计),我们能获得更多的信息,从而帮助我们做出相应的决策。

在逻辑回归模型下计算上涨概率的置信区间

简单回顾一下知识点,逻辑回归(Logistic Regression)模型为

P(y=1x)=11+e(β0+β1x)P(y=1|x) = \frac{1}{1+e^{-(\beta_0 + \beta_1 x)}}

对于数据 {x1,y1;x2,y2;...;xN,yN}\{x_1, y_1;x_2, y_2;...;x_N, y_N\}(其中 yi{0,1}y_i\in\{0, 1\},模型通过极大似然估计(Maximum likelihood estimation, MLE)求解参数 β0^\hat{\beta_0}β1^\hat{\beta_1},使得在这组参数下,生成该组数据的概率最大,严格地说

β0^,β1^=maxβ0,β1log(i=1Npiyi(1pi)1yi)\hat{\beta_0}, \hat{\beta_1} = \max_{\beta_0, \beta_1} \log ( \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i} )

其中 pi=11+e(β0+β1xi)p_i = \frac{1}{1+e^{-(\beta_0 + \beta_1 x_i)}}
熟悉统计推断的同学应该立马意识到 MLE 具有优良的性质(under mild assumption):

  1. 一致性(Consistency):θ^mlePθtrue\hat{\theta}_{mle} \overset{P}{\rightarrow} \theta_{true}
  2. 渐进正态性(Asymptotic normality):存在 σmle2\sigma^2_{mle},使得 θ^mleθtruedN(0,σmle2)\hat{\theta}_{mle} - \theta_{true} \overset{d}{\rightarrow} N(0, \sigma^2_{mle})

更进一步地,σmle2=I(θtrue)1\sigma^2_{mle} = I(\theta_{true})^{-1},其中 I(θture)=2θ2loglik(Yθtrue)I(\theta_{ture}) = \frac{\partial^2}{\partial \theta^2} \log \text{lik}(Y|\theta_{true}) 是 Fisher information matrix。

因为我们的参数 β0^,β1^\hat{\beta_0}, \hat{\beta_1} 是从 MLE 估计出的,由渐近正态性可以认为 logit = β0^+β1^x\hat{\beta_0} + \hat{\beta_1} x 近似服从于正态分布,从而构造 logit 的 95% 置信区间 [logitlower,logitupper][\text{logit}_{lower}, \text{logit}_{upper}],从而给出支付宝猜涨跌上涨概率 pp 的 95% 置信区间

[11+exp(logitlower),11+exp(logitupper)][\frac{1}{1+\exp(-\text{logit}_{lower})}, \frac{1}{1+\exp(-\text{logit}_{upper})}]

数学推导 Fisher information matrix

定义 log-likelihood

l(β0,β1)=log(i=1Npiyi(1pi)1yi)=i=1Nyilogpi1pi+log(1pi)l(\beta_0, \beta_1) = \log ( \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i} ) = \sum_{i=1}^N y_i \log \frac{p_i}{1-p_i} + \log(1-p_i)

其中 pi=11+exp(β0β1xi)p_i = \frac{1}{1+\exp(-\beta_0-\beta_1 x_i)},代入上式,得到

l(β0,β1)=i=1Nyi(β0+β1xi)i=1Nlog(1+exp(β0+β1xi))l(\beta_0, \beta_1) = \sum_{i=1}^N y_i(\beta_0+\beta_1 x_i) - \sum_{i=1}^N \log(1+\exp(\beta_0+\beta_1 x_i))

l(β0,β1)l(\beta_0, \beta_1) 求偏导

lβ0=i=1N(yiexp(β0+β1xi)1+exp(β0+β1xi))lβ1=i=1N(yiexp(β0+β1xi)1+exp(β0+β1xi))xi\begin{aligned} \frac{\partial l}{\partial \beta_0} & = \sum_{i=1}^N ( y_i - \frac{\exp(\beta_0+\beta_1 x_i)}{1+\exp(\beta_0+\beta_1 x_i)}) \\ \frac{\partial l}{\partial \beta_1} & = \sum_{i=1}^N (y_i - \frac{\exp(\beta_0+\beta_1 x_i)}{1+\exp(\beta_0+\beta_1 x_i)}) x_i \end{aligned}

l(β0,β1)l(\beta_0, \beta_1) 求二阶偏导

2lβ02=i=1Nexp(β0+β1xi)(1+exp(β0+β1xi))22lβ0β1=i=1Nexp(β0+β1xi)(1+exp(β0+β1xi))2xi2lβ12=i=1Nexp(β0+β1xi)(1+exp(β0+β1xi))2xi2\begin{aligned} \frac{\partial^2 l}{\partial \beta_0^2} & = -\sum_{i=1}^N \frac{\exp(\beta_0+\beta_1 x_i)}{(1+\exp(\beta_0+\beta_1 x_i))^2} \\ \frac{\partial^2 l}{\partial \beta_0 \beta_1} & = -\sum_{i=1}^N \frac{\exp(\beta_0+\beta_1 x_i)}{(1+\exp(\beta_0+\beta_1 x_i))^2} x_i \\ \frac{\partial^2 l}{\partial \beta_1^2} & = -\sum_{i=1}^N \frac{\exp(\beta_0+\beta_1 x_i)}{(1+\exp(\beta_0+\beta_1 x_i))^2} x_i^2 \end{aligned}

W(β0,β1)=(2lβ022lβ0β12lβ0β12lβ12)W(\beta_0, \beta_1) = \begin{pmatrix} \frac{\partial^2 l}{\partial \beta_0^2} & \frac{\partial^2 l}{\partial \beta_0 \beta_1} \\ \frac{\partial^2 l}{\partial \beta_0 \beta_1} & \frac{\partial^2 l}{\partial \beta_1^2} \end{pmatrix}

由 MLE 的性质,我们知道

(β0^,β1^)N((β0,β1),W(β0,β1)1)(\hat{\beta_0}, \hat{\beta_1}) \sim N((\beta_0, \beta_1), W(\beta_0, \beta_1)^{-1})

我们用 W^=W(β0^,β1^)\hat{W}=W(\hat{\beta_0}, \hat{\beta_1}) 来估计 W(β0,β1)1)W(\beta_0, \beta_1)^{-1}),这就估计出了 β0^,β1^\hat{\beta_0}, \hat{\beta_1} 的协方差矩阵。

数值计算

首先使用 Logistic Regression 计算 β0^,β1^\hat{\beta_0}, \hat{\beta_1}

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
from sklearn.linear_model import LogisticRegression

with open('stock.csv', 'r') as f:
stock_dataset = np.loadtxt(f, delimiter=',', skiprows=1, usecols=[1,2,3])
x = stock_dataset[:, 1] - stock_dataset[:, 0]
y = stock_dataset[:, 2]
N = len(y)
x = np.reshape(x, (N,1))

logr = LogisticRegression()
logr.fit(x, y)
print('beta_0: {0:.5f}, beta_1: {1:.5f}'.format(logr.intercept_[0], logr.coef_[0][0]))

计算结果为

1
beta_0: -0.08161, beta_1: 0.09779

接着计算 Fisher information matrix 的估计 W(β0^,β1^)W(\hat{\beta_0}, \hat{\beta_1})

1
2
3
4
5
6
b0 = -0.08161
b1 = 0.09779
W_11 = np.sum(np.exp(b0 + b1*x) / np.square(1 + np.exp(b0 + b1*x)))
W_12 = np.sum(x*np.exp(b0 + b1*x) / np.square(1 + np.exp(b0 + b1*x)))
W_22 = np.sum(x*x*np.exp(b0 + b1*x) / np.square(1 + np.exp(b0 + b1*x)))
print('[[{0:.2f}, {1:.2f}]\n [{1:.2f}, {2:.2f}]]'.format(W_11, W_12, W_22))

计算结果为

1
2
[[27.37, -6.219]
[-6.219, 4553.]]

计算 W(β0^,β1^)W(\hat{\beta_0}, \hat{\beta_1}) 的矩阵逆

1
2
W = [[27.37, -6.219], [-6.219, 4553]]
print(np.linalg.inv(W))

得到

1
2
[[0.036540994 0.000049912]
[0.000049912 0.00021967 ]]

因此 Var(logit)=Var(β0^+β1^x)=Var(β0^)+2xCov(β0^,β1^)+x2Var(β1^)\mathrm{Var}(\text{logit}) = \mathrm{Var}(\hat{\beta_0} + \hat{\beta_1} x) = \mathrm{Var}(\hat{\beta_0}) + 2x\mathrm{Cov}(\hat{\beta_0}, \hat{\beta_1}) + x^2 \mathrm{Var}(\hat{\beta_1}),开根号,得到 std(logit)=0.036540994+0.000049912x+0.00021967x2\text{std}(\text{logit}) = \sqrt{0.036540994 + 0.000049912 x + 0.00021967 x^2}
至此,我们可以根据 [logit1.96std(logit),logit+1.96std(logit)][\text{logit} - 1.96 * \text{std}(\text{logit}), \text{logit} + 1.96 * \text{std}(\text{logit})] 构造 logit 的 95% 置信区间,于是支付宝猜涨跌上涨概率 pp 的 95% 置信区间为

[11+exp(logit+1.96std(logit)),11+exp(logit1.96std(logit))]\left[ \frac{1}{1+\exp(-\text{logit} + 1.96 * \text{std}(\text{logit}))}, \quad \frac{1}{1+\exp(-\text{logit} - 1.96 * \text{std}(\text{logit}))} \right ]

在线计算

将我们构造的 95% 置信区间计算过程在线部署后,大家可以访问 猜涨跌预测 页面,输入今天的上证指数开盘价和中午收盘价,来估计下午收盘后,今天的指数相比于昨天收盘价的涨跌的概率以及相对应的 95% 置信区间,欢迎大家踊跃测试,撒花~