请　[注册]　或　[登录]　 | 返回主站

量化交易吧 /  量化平台 帖子：3354526 新帖：23

【QLS-47】自回归与AR模型

你在说什么呢发表于：5 月 9 日 20：24回复(1)

金融市场上，时间序列是我们天天打交道的数据形式。QLS课程中，有专门一连串课程来介绍时间序列的各种模型，本文就是该系列的第一篇。让我们一起来学习最简单的时间序列模型：AR模型。

AR模型¶

AR模型，或者说AR$(p)$模型，是对一个时间序列在它的历史数据，即滞后项上做回归。最简单的自回归模型是AR$(1)$模型，只使用了一阶滞后项。AR$(1)$模型的数学形式是：

$$ x_t = b_0 + b_1 x_{t - 1} + \epsilon_t $$

$x_{t - 1}$代表$(t - 1)$时刻序列的取值， $\epsilon_t$是残差项。我们可以将这个数学表达式推广，得到 AR$(p)$模型：

$$ x_t = b_0 + b_1 x_{t-1} + b_2 x_{t - 2} \ldots + b_p x_{t - p} + \epsilon_t $$

为了让AR模型得到较为可信的结果，我们必须检验试剂序列是否是平稳的。平稳意味着需要满足以下三个条件：

1.序列的期望在任何时刻都是有限常数，即无论t取什么值， $E[y_t] = \mu$ and $\mu < \infty$ ；

2.序列的方差在任何时间段里都是有限常数；

3.在任意相同时间段上，序列的协方差是有限常数，即 $$ COV(y_t, y_{t - s}) = \lambda, \ |\lambda| < \infty, \text{ $\lambda$ 是常数}, \ t = 1, 2, \ \ldots, T; \ s = 0, \pm 1, \pm 2, \ldots, \pm T $$

如果上面的条件没有满足，我们估计的结果就会失去现实意义：估测的参数会有偏差，任何关于模型的检验都不可用。不幸的是，在金融市场里想要找到一个平稳时间序列是非常困难的。例如，我们考察苹果公司的股价，我们可以很清楚地看到它的上升趋势。股价的平均值随着时间的推移而增长。然而，我们有办法将非平稳序列平稳化。当我们要完成这种转化后，得到的新序列就能符合以上三个准则。

模拟数据¶

下面我们从一个AR$(3)$过程中提取数据，绘制时间序列曲线。

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import statsmodels.tsa as tsa
import matplotlib.pyplot as plt

# 设置随机数种子，确保每次随机数相同
np.random.seed(100)

# 下面这个函数用于模拟AR过程
# b1 ... bk是自相关系数, 再加上服从正态分布的残差项
def AR(b, X, mu, sigma):
    l = min(len(b)-1, len(X))
    b0 = b[0]
    
    return b0 + np.dot(b[1:l+1], X[-l:]) + np.random.normal(mu, sigma)

b = np.array([0, 0.8, 0.1, 0.05])
X = np.array([1])

mu = 0
sigma = 1

for i in range(10000):
    X = np.append(X, AR(b, X, mu, sigma))
plt.figure(figsize=(16,8))    
plt.plot(X)
plt.xlabel('Time')
plt.ylabel('AR Series Value');

注意到上面的曲线围绕着某个中间值上下波动，这个值就是我们时间序列的平均值。如果我们有一个序列，在任何时间里它都在均值两侧波动（波动范围在一定距离内），那么我们就可以假设这个序列是平稳的。当然，我们还希望能更严格地对这点进行检验。不过在这个例子里，我们的确已经由AR(3)模型构造了平稳的序列，所以不用担心需要对它进行平稳性检验。

厚尾风险¶

自相关过程通常会比正态分布具有更多的极端值。这是因为每一个时间点的值都会被最近时间的值所影响。这种现象被称作自相关过程的肥尾分布，也称厚尾风险。

厚尾风险的讨论集中在，“尾部事件确实会发生，但由于它的低概率性，人们很难去建模”这一点上。如果我们有理由去推测一个过程是自相关的，那么我们就应该考虑它的厚尾风险并做相应调整。

AR模型只是厚尾风险来源的一种，所以也不能因为某个过程不是AR风险，就判定它不具有厚尾风险。

我们下面来验证一下厚尾分布极端事件发生的概率是否确实更高。

def compare_tails_to_normal(X):
    A = np.zeros((2,4))    
    for k in range(4):             
        #厚尾分布 vs 正态分布
        A[0, k] = len(X[X > (k + 1)]) / float(len(X))      
        A[1, k] = 1 - stats.norm.cdf(k + 1) 
    print ('Frequency of std events in X \n1: %s\t2: %s\t3: %s\t4: %s' % tuple(A[0]))
    print ('Frequency of std events in a normal process \n1: %s\t2: %s\t3: %s\t4: %s' % tuple(A[1]))
    return A

compare_tails_to_normal(X);

Frequency of std events in X 
1: 0.32756724327567244	2: 0.18548145185481452	3: 0.09099090090990901	4: 0.0357964203579642
Frequency of std events in a normal process 
1: 0.15865525393145707	2: 0.02275013194817921	3: 0.0013498980316301035	4: 3.167124183311998e-05

显然，厚尾分布确实更容易发生极端事件，这在金融领域更是不容小觑。

错误的方差估计¶

由于AR过程是厚尾且非正态的，对AR过程的方差估计会出现偏差。由于方差在很多统计量的计算中都要使用（特别是置信区间和p值），所以这一点非常致命。对于AR过程，我们无法使用平常通用的区间估计和p值检验。

让我们来确认这一点。首先先定义一个用来计算输入序列95%水平置信区间的函数。

def compute_unadjusted_interval(X):
    T = len(X)
    # 计算 mu and sigma 
    mu = np.mean(X)
    sigma = np.std(X)
    # 95%的置信水平对应着1.96个标准差
    lower = mu - 1.96 * (sigma/np.sqrt(T))
    upper = mu + 1.96 * (sigma/np.sqrt(T))
    return lower, upper

# 下面这个函数的作用是：如果计算值是0，则返回True
def check_unadjusted_coverage(X):
    l, u = compute_unadjusted_interval(X)
    # 用于确认 l <= 0 <= u
    if l <= 0 and u >= 0:
        return True
    else:
        return False
    
def simululate_AR_process(b, T):
    X = np.array([1])

    mu = 0
    sigma = 1

    for i in range(T):
        X = np.append(X, AR(b, X, mu, sigma))
        
    return X

现在，我们用AR模型来模拟1000次AR过程，并估算它的平均值和95%水平的置信区间。然后我们看看置信区间中是否包含0（长期情况下序列的真实平均值）。

trials = 1000
outcomes = np.zeros((trials, 1))

for i in range(trials):
    #note these are the same values we used to generate the initial AR array
    Z = simululate_AR_process(np.array([0, 0.8, 0.1, 0.05]), 100)
    if check_unadjusted_coverage(Z):
        # The internal contains 0, the true value
        outcomes[i] = 1
    else:
        outcomes[i] = 0

让我们来确认一下，如果我们的检验是正确的，那么应该大约有95%的实验包含0。但实际上，正如下面结果所示，包含0的实验不过五分之一。

np.sum(outcomes) / trials

0.193

看样子肯定是哪里出了什么问题。事实上，AR序列移动得比想象得多：方差估计假设序列平稳，完全没有考虑进这些移动。结果就是置信区间比实际的会更小，无法将所有本来的正确值都囊括。平稳性检验通常能够敏锐地捕捉到AR序列的种种行为，并让我们知晓方差估计是否出错。

实践上，我们一般很难去精准地估测AR序列的方差。一种方法是Newey-West估计，关于它的详细内容可参见维基百科“Newey–West estimator”。

AR模型的检验¶

为了确定AR模型滞后项的阶数$p$，我们关注实践序列的自相关性。自相关就是序列当前的值与它过去的值的相关性， $k$阶自相关的数学表达式：

$$ \rho_k = \frac{COV(x_t, x_{t - k})}{\sigma_x^2} = \frac{E[(x_t - \mu)(x_{t - k} - \mu)}{\sigma_x^2} $$

$k$代表滞后项的阶数，我们无法直观测到自相关系数，因此我们用下面的形式来估测它：

$$ \hat{\rho}_k = \frac{\sum_{t = k + 1}^T[(x_t - \bar{x})(x_{t - k} - \bar{x})]}{\sum_{t = 1}^T (x_t - \bar{x})^2} $$

为了达到我们的目的，我们可以使用ACF和PACF这对工具来决定我们模型的阶数。这两个函数在很多统计包中都有，我们可以很方便地使用它们来确定$p$的合适取值。我们用下面的例子来演示这两个函数：

from statsmodels.tsa.stattools import acf, pacf

X = simululate_AR_process(np.array([0, 0.8, 0.1, 0.05]), 1000)

nlags = 40
X_acf = acf(X, nlags=nlags)
print ('Autocorrelations:\n' + str(X_acf) + '\n')
X_pacf = pacf(X, nlags=nlags)
print ('Partial Autocorrelations:\n' + str(X_pacf))

Autocorrelations:
[1.         0.57070553 0.58001725 0.88553976 0.55669882 0.57640436
 0.78307405 0.54759415 0.55254723 0.70158722 0.53082572 0.52673077
 0.62850967 0.50091633 0.48541475 0.56544198 0.46391848 0.45390828
 0.51491108 0.43375521 0.43938017 0.46783793 0.40843913 0.42270738
 0.42566802 0.38298857 0.38379558 0.40029717 0.36548704 0.36088258
 0.37743521 0.34745858 0.34065687 0.35239327 0.33194099 0.31758452
 0.33113111 0.31688646 0.2818213  0.31237461 0.31252031]

Partial Autocorrelations:
[ 1.          0.57127624  0.37827414  0.80672256  0.00162478  0.00490887
 -0.04497956  0.04013524 -0.05301909  0.00706455 -0.00528104  0.0051476
 -0.03421094 -0.03637203 -0.06841025  0.00396849 -0.02433695  0.03906542
  0.02344142  0.02958713  0.07306916 -0.01201434  0.00279546  0.01328444
 -0.00679875 -0.01547705 -0.09186641  0.05847768  0.03298422  0.04395628
 -0.01605778  0.00286143 -0.00734411 -0.01496441  0.01152254 -0.01185183
  0.01895763  0.01901218 -0.09532658  0.00791091  0.06173332]

让我们用图来分别展示ACF和PACF的值。

plt.figure(figsize=(16,8))
plt.plot(X_acf, 'ro')
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title("ACF")

Text(0.5, 1.0, 'ACF')
plt.figure(figsize=(16,8))
plt.plot(X_pacf, 'ro')
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.title("PACF");

统计检验¶

仅仅看上面的散点图当然会一头雾水，我们需要使用统计方法进行严格的检验。acf和pacf函数会返回自相关系数的置信区间，我们可以通过它们来确认置信区间是否包含0.如果它们确实如此，我们就说0在正确参数值的置信区间集内，并且忽略不显著的滞后项。

X_acf, X_acf_confs = acf(X, nlags=nlags, alpha=0.05)
X_pacf, X_pacf_confs = pacf(X, nlags=nlags, alpha=0.05)

def plot_acf(X_acf, X_acf_confs, title='ACF'):
    errorbars = np.ndarray((2, len(X_acf)))
    errorbars[0, :] = X_acf - X_acf_confs[:,0]
    errorbars[1, :] = X_acf_confs[:,1] - X_acf
    plt.figure(figsize=(16,8))
    plt.plot(X_acf, 'ro')
    plt.errorbar(range(len(X_acf)), X_acf, yerr=errorbars, fmt='none', ecolor='gray', capthick=2)
    plt.xlabel('Lag')
    plt.ylabel('Autocorrelation')
    plt.title(title);
plot_acf(X_acf, X_acf_confs)

plot_acf(X_pacf, X_pacf_confs, title='PACF')

这个例子里，ACF图不太好判断。但我们注意到PACF图中，只有前三阶的值是显著大于0的，于是我们可以认为AR模型的阶数为3。不过，如果你生成了另一个随机序列，图中的结果或许未必就是3了。

对于现实世界中的时间序列，我们也用以上的方法来决定我们模型的形式。我们会试图用极大似然函数的方法来拟合模型。

拟合模型¶

我们对上面的X序列进行建模， model 这个变量类型包含很多有用的信息，可使用? 或者help来了解关于它的更多信息。我们主要关注model.params模型中的估测参数和model.bse估测参数的标准差。

model = tsa.api.AR(X)
model = model.fit()

print ('Parameters')
print (model.params)
print ('Standard Error')
print (model.bse)

Parameters
[ 0.01001762  0.04631952  0.06669125  0.84444655 -0.02662476  0.07031018
 -0.03883539  0.05362627 -0.06257934  0.03183551  0.02255835  0.06366225
 -0.02663838 -0.00824025 -0.09032174 -0.02259851 -0.05894006 -0.03575162
  0.01899755  0.0265497   0.08158745 -0.00494997]
Standard Error
[0.03210283 0.03228516 0.03221285 0.03226938 0.04224386 0.04226218
 0.04228528 0.04223737 0.04222499 0.04230633 0.04232241 0.04227582
 0.04230175 0.04232492 0.04221253 0.04221065 0.04218858 0.04214434
 0.04212185 0.03217742 0.03220652 0.03227438]

model_confs = np.asarray((model.params - model.bse, model.params + model.bse)).T
plot_acf(model.params, model_confs, title='Model Estimated Parameters')

确定滞后项阶数¶

我们可以看到，上面的模型产生了不少估测参数，我们知道它肯定太多了，因为我们是用一个AR(3)过程来模拟得到的数据。用AR模型建模会产生多于实际的滞后项，是因为间接的依赖性。如果$X_t$ 依赖于 $X_{t-1}$， 则会对$X_{t-2}$ 有一个更小程度的依赖。生成数据时如果不止一个滞后项，我们建模时往往就会得到更为复杂的模型结构。这些间接依赖性可以用简单的方法来估测。

总体而言，你很难从有很多参数的模型中获得有用的信息，因为过拟合会导致适得其反的结果。在这种情况下，我们想要找到合适的阶数，使得我们相信它能够解释什么正在发生，并且避免过拟合的风险。

观察ACF和PACF，我们知道只有前3阶的系数是可能有用得多。然而，当我们将阶数扩大到10来再次确认我们的初始数据。我们将使用信息标准，具体地说用Akaike信息标准（AIC）和Bayes信息标注（BIC）来决定合适的系数。

计算我们考察范围内所有阶数模型的AIC和BIC，并将最小值分别记为$AIC_{min}$ 和 $BIC_{min}$。它们是模型能够在两种考察标准下最小的信息损失。对每一种IC我们可以计算每一个模型的相对似然（relative likelihood）：

$$l = e^{(IC_{min} - IC_{i})/2}$$

我们可以用过$l$来在AIC模型上最小化信息损失,听起来可能比较难理解，让我们来看看下面的实例。

N = 10
AIC = np.zeros((N, 1))

for i in range(N):
    model = tsa.api.AR(X)
    model = model.fit(maxlag=(i+1))
    AIC[i] = model.aic
    
AIC_min = np.min(AIC)
model_min = np.argmin(AIC)

print ('Relative Likelihoods')
print (np.exp((AIC_min-AIC) / 2))
print ('Number of parameters in minimum AIC model %s' % (model_min+1))

Relative Likelihoods
[[0.54739169]
 [0.59052713]
 [0.99861675]
 [0.997749  ]
 [0.99628437]
 [0.99619915]
 [0.9960884 ]
 [1.        ]
 [0.99850924]
 [0.99809322]]
Number of parameters in minimum AIC model 8

N = 10
BIC = np.zeros((N, 1))

for i in range(N):
    model = tsa.api.AR(X)
    model = model.fit(maxlag=(i+1))
    BIC[i] = model.bic
    
BIC_min = np.min(BIC)
model_min = np.argmin(BIC)

print ('Relative Likelihoods')
print (np.exp((BIC_min-BIC) / 2))
print ('Number of parameters in minimum BIC model %s' % (model_min+1))

Relative Likelihoods
[[0.5508575 ]
 [0.59280494]
 [1.        ]
 [0.99666667]
 [0.99274501]
 [0.99020386]
 [0.98763995]
 [0.98905709]
 [0.9851272 ]
 [0.98226451]]
Number of parameters in minimum BIC model 3

综合AIC和BIC来看，我们得出结论：AR模型的阶数应该取3。

检验残差¶

最后一步是对残差项进行检验。AIC和BIC其实已经进行了这项工作的一部分，因为它们是测量模型损失了多少信息（残差中的信息）。此处我们再对残差的正态性进行检验

model = tsa.api.AR(X)
model = model.fit(maxlag=3)

from statsmodels.stats.stattools import jarque_bera

score, pvalue, _, _ = jarque_bera(model.resid)

if pvalue < 0.10:
    print ('We have reason to suspect the residuals are not normally distributed.')
else:
    print ('The residuals seem normally distributed.')

The residuals seem normally distributed.

如上，残差项服从正态分布，这也与AR模型的数学表达式一致。

       

全部回复

0/140

量化课程

    移动端课程