统计估计不稳定性

作者：chen_h
微信号 & QQ：862251340
微信公众号：coderpai

---

参数

参数是模型用于约束其预测的任何参数。通常，参数是有助于描述数据集或者分布的数据。例如，正太分布的均值是一个参数；实际上，我们说正太分布时通过其均值和方差来参数化的。如果我们采用从正太分布中抽取的一组样本的均值，我们得到分布均值的估计。类似的，一组观察值的平均值是对基础分布参数的估计（通常认为是正太的）。其他参数包括中值，与另一个系列的相关系数，标准偏差以及数据集的每个其他测量值。

你永远不知道，你只是估计

当你取数据集的平均值时，你不知道真正的平均值。你已从所拥有的数据中尽可能的估算出均值。估计可能会被取消。对于你估计的任何参数都是如此。要真正理解发生了什么，你需要通过查看其稳定性、标准误差、置信区间来确定你的估计有多好。

估计不稳定

每当我们考虑一组观察时，我们对参数的计算只能是估计。它会随着我们进行更能多测量或者随着时间的推移而改变，我们会得到新的观察结果。我们可以通过查看参数在参数数据的不同子集时如何变化来量化我们估计中的不确定性。利润，标准偏差描述了一组的平均值和每个观察的平均值的差异，即每个观察本身的平均值。在金融应用程序中，数据通常按时间顺序排列。在这种情况下，我们可以估计不同时间点的参数；比如，前 30 天。通过观察这个移动估计在我们改变时间窗口时波动的程度，我们可以计算估计参数的不稳定性。

# We'll be doing some examples, so let's import the libraries we'll need
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

例子：均值和标准差

首先，让我们看一下正太分布的一些样本。我们知道分布的均值为 0 ，标准差为 1；但是如果我们从观察中测量参数，我们将只得到大约 0 和 大约 1；我们可以看到这些估计如何随着我们采取越来越多的样本而改变；

# Set a seed so we can play with the data without generating new random numbers every time
np.random.seed(123)

normal = np.random.randn(500)
print(np.mean(normal[:10]))
print(np.mean(normal[:100]))
print(np.mean(normal[:250]))
print(np.mean(normal))

# Plot a stacked histogram of the data
plt.hist([normal[:10], normal[10:100], normal[100:250], normal], normed=1, histtype='bar', stacked=True);
plt.ylabel('Frequency')
plt.xlabel('Value');

plt.show()

-0.26951611032632805
0.027109073490359778
-0.020616059111720507
-0.038643973513210604

print(np.std(normal[:10]))
print(np.std(normal[:100]))
print(np.std(normal[:250]))
print(np.std(normal))

1.236304801499023
1.128240470477961
1.0174604368340197
1.0032028561568238

请注意，尽管平均值和标准偏差越接近 0 和 1 的概率随着样本数量的增加而增加，但我们并不总是通过获取更多数据点来获得更好的估计值。无论我们的期望是什么，我们总能获得不同的结果，我们的目标通常是计算结果与预期显着不同的概率。

对于时间序列数据，我们通常只关心数据的连续子集。移动平均线（也称为运行或者滚动）将前 n 个数据点的平均值分配给每个时间点。下面，我们计算股票价格的 90 天移动平均线，并绘制它以查看它如何变化。开始时没有结果，因为我们首先必须累积至少 90 天的数据。

例子：非正态基础分布

如果基础数据不正太分布会发生什么呢？一个含义是这是非常具有欺骗性的。因此，测试数据的正常性非常重要。我们将以 Jarque-Bera 测试为例。

#Generate some data from a bi-modal distribution
def bimodal(n):
    X = np.zeros((n))
    for i in range(n):
        if np.random.binomial(1, 0.5) == 0:
            X[i] = np.random.normal(-5, 1)
        else:
            X[i] =  np.random.normal(5, 1)
    return X
            
X = bimodal(1000)

#Let's see how it looks
plt.hist(X, bins=50)
plt.ylabel('Frequency')
plt.xlabel('Value')
plt.show()
print('mean:', np.mean(X))
print('standard deviation:', np.std(X))

mean: 0.009847581282146528
standard deviation: 5.060708740105227

果然，平均值对数据中发生的事情的信息量非常大。我们已将所有数据合并为单一估算，并且丢失了大量信息。如果我们假设它是正太分布的，那么这分布应该是正确的。

mu = np.mean(X)
sigma = np.std(X)

N = np.random.normal(mu, sigma, 1000)

plt.hist(N, bins=50)
plt.ylabel('Frequency')
plt.xlabel('Value');

plt.show()

我们将使用 Jarque-Bera 测试来测试我们的数据，看它是否正太分布。显著的 p 值表示非正太。

from statsmodels.stats.stattools import jarque_bera

jarque_bera(X)

(142.12550136207705,
 1.3735343038981241e-31,
 -0.007644415681800414,
 1.1531707484649847)

果然这个值 < 0.05，我们说 X 是非正太的。这使得我们免于意外制作的可怕预测。

例子：sharpe ratio

一个通常用于描述资产和投资组合表现的统计数据是夏普比率，它衡量投资组合所获得的每单位额外风险的额外回报，相对于无风险回报来说。

$R = \frac{E[r_{a} - r_{b}]}{\sqrt{Var(r_{a} - r_{b})}}$R=Var(ra​−rb​)​E[ra​−rb​]​

其中，$r_{a}$ra​ 是我们资产的回报，$r_{b}$rb​ 是无风险回报，与均值和标准差一样，我们可以计算滚动夏普比率，以了解我们的估计如何随着时间变化。

def sharpe_ratio(asset, riskfree):
    return np.mean(asset - riskfree)/np.std(asset - riskfree)

start = '2012-01-01'
end = '2015-01-01'
# Use an ETF that tracks 3-month T-bills as our risk-free rate of return
treasury_ret = get_pricing('BIL', fields='price', start_date=start, end_date=end).pct_change()[1:]
pricing = get_pricing('AMZN', fields='price', start_date=start, end_date=end)
returns = pricing.pct_change()[1:] # Get the returns on the asset

# Compute the running Sharpe ratio
running_sharpe = [sharpe_ratio(returns[i-90:i], treasury_ret[i-90:i]) for i in range(90, len(returns))]

# Plot running Sharpe ratio up to 100 days before the end of the data set
_, ax1 = plt.subplots()
ax1.plot(range(90, len(returns)-100), running_sharpe[:-100]);
ticks = ax1.get_xticks()
ax1.set_xticklabels([pricing.index[i].date() for i in ticks[:-1]]) # Label x-axis with dates
plt.xlabel('Date')
plt.ylabel('Sharpe Ratio');

夏普比率看起来相当不稳定，很明显，仅将其报告为单一值对于预测未来价值并不是很有帮助。相反，我们可以计算上述数据的均值和标准差，然后看看它是否有助于我们预测接下来 100 天的夏普比率。

# Compute the mean and std of the running Sharpe ratios up to 100 days before the end
mean_rs = np.mean(running_sharpe[:-100])
std_rs = np.std(running_sharpe[:-100])

# Plot running Sharpe ratio
_, ax2 = plt.subplots()
ax2.set_xticklabels([pricing.index[i].date() for i in ticks[:-1]]) # Label x-axis with dates
ax2.plot(range(90, len(returns)), running_sharpe)

# Plot its mean and the +/- 1 standard deviation lines
ax2.axhline(mean_rs)
ax2.axhline(mean_rs + std_rs, linestyle='--')
ax2.axhline(mean_rs - std_rs, linestyle='--')

# Indicate where we computed the mean and standard deviations
# Everything after this is 'out of sample' which we are comparing with the estimated mean and std
ax2.axvline(len(returns) - 100, color='pink');
plt.xlabel('Date')
plt.ylabel('Sharpe Ratio')
plt.legend(['Sharpe Ratio', 'Mean', '+/- 1 Standard Deviation'])

print('Mean of running Sharpe ratio:', mean_rs)
print('std of running Sharpe ratio:', std_rs)

Mean of running Sharpe ratio: 0.0646215053325
std of running Sharpe ratio: 0.0778015776531

在这种情况下，标准偏差约为范围的四分之一，因此该数据极不稳定。虽然我们仍然观察到查过一个标准偏差的数据，但是在向前看时考虑到这一点可以提供比仅适用均值更好的预测。我们还可以计算夏普比率的滚动均值来尝试和跟踪趋势；但在这种情况下，我们也应该记住标准偏差。

例子：移动平均线

假设你使用回滚窗口取平均值；你如何确定该估计的标准误差？让我们从一个显示 90 天移动平均线的例子开始。

# Load time series of prices
start = '2012-01-01'
end = '2015-01-01'
pricing = get_pricing('AMZN', fields='price', start_date=start, end_date=end)

# Compute the rolling mean for each day
mu = pd.rolling_mean(pricing, window=90)

# Plot pricing data
_, ax1 = plt.subplots()
ax1.plot(pricing) 
ticks = ax1.get_xticks()
ax1.set_xticklabels([pricing.index[i].date() for i in ticks[:-1]]) # Label x-axis with dates
plt.ylabel('Price')
plt.xlabel('Date')

# Plot rolling mean
ax1.plot(mu);
plt.legend(['Price','Rolling Average']);

这让我们可以看到均值的不稳定性 / 标准误差，并有助于预测数据的未来变化。我们可以通过计算滚动平均值的均值和标准差来量化这种可变性。

print 'Mean of rolling mean:', np.mean(mu)
print 'std of rolling mean:', np.std(mu)

Mean of rolling mean: 288.399003348
std of rolling mean: 51.1188097398

事实上，我们用来量化变异性的标准差本身就是一个变量。下面我们绘制滚动标准偏差（对于 90 天窗口），并计算器平均值和标准偏差。

# Compute rolling standard deviation
std = pd.rolling_std(pricing, window=90)

# Plot rolling std
_, ax2 = plt.subplots()
ax2.plot(std)
ax2.set_xticklabels([pricing.index[i].date() for i in ticks[:-1]]) # Label x-axis with dates
plt.ylabel('Standard Deviation of Moving Average')
plt.xlabel('Date')

print 'Mean of rolling std:', np.mean(std)
print 'std of rolling std:', np.std(std)

Mean of rolling std: 17.3969897999
std of rolling std: 7.54619079684

为了查看这个变化的标准偏差对我们的数据意味着什么，让我们再次绘制数据以及布林带：滚动平均值，一个滚动标准偏差（数据）高于平均值，一个标准偏差低于平均值。

请注意，尽管标准偏差为我们提供了有关数据传播的更多信息，但我们无法在不假设基础流程的特定分布的情况下为我们对未来观察的预期分配精确概率。

# Plot original data
_, ax3 = plt.subplots()
ax3.plot(pricing)
ax3.set_xticklabels([pricing.index[i].date() for i in ticks[:-1]]) # Label x-axis with dates

# Plot Bollinger bands
ax3.plot(mu)
ax3.plot(mu + std)
ax3.plot(mu - std);
plt.ylabel('Price')
plt.xlabel('Date')
plt.legend(['Price', 'Moving Average', 'Moving Average +1 Std', 'Moving Average -1 Std'])

结论

每当我们计算数据集的参数时，我们也应该计算其波动率。否则，我们不知道是否应该期望新数据点与此参数对齐。计算波动率的一种好方法是将数据划分为子集并从每个子集中估计参数，然后找出结果之间的可变性。在我们的样本期后可能仍然存在外部因素，而这些因素是我们无法预测的。然而，标准误差的不稳定性分析和测试对于高速我们我们应该多少不信任我们的估计仍然非常有用。