AR模型全称为Autoregressive Models，即自回归模型，用于刻画因变量能由它的多个滞后项表示。p阶自回归模型可以写成：
下面模拟一个AR(1)模型。
import pandas as pdimport numpy as npimport statsmodels.tsa.api as smt #tsa为Time Series analysis缩写import statsmodels.api as smimport scipy.stats as scsfrom arch import arch_model#画图import matplotlib.pyplot as pltimport matplotlib as mpl
%matplotlib inline#正常显示画图时出现的中文和负号from pylab import mpl
mpl.rcParams['font.sans-serif']=['SimHei']
mpl.rcParams['axes.unicode_minus']=False
#先定义一个画图函数，后面都会用到
def ts_plot(data, lags=None,title=''):
if not isinstance(data, pd.Series):
data = pd.Series(data)#matplotlib官方提供了五种不同的图形风格，#包括bmh、ggplot、dark_background、fivethirtyeight和grayscalewith plt.style.context('ggplot'):
fig = plt.figure(figsize=(10, 8))
layout = (3, 2)
ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
acf_ax = plt.subplot2grid(layout, (1, 0))
pacf_ax = plt.subplot2grid(layout, (1, 1))
qq_ax = plt.subplot2grid(layout, (2, 0))
pp_ax = plt.subplot2grid(layout, (2, 1))
data.plot(ax=ts_ax)
ts_ax.set_title(title+'时序图')
smt.graphics.plot_acf(data, lags=lags, ax=acf_ax, alpha=0.5)
acf_ax.set_title('自相关系数')
smt.graphics.plot_pacf(data, lags=lags, ax=pacf_ax, alpha=0.5)
pacf_ax.set_title('偏自相关系数')
sm.qqplot(data, line='s', ax=qq_ax)
qq_ax.set_title('QQ 图')
scs.probplot(data, sparams=(data.mean(),
data.std()), plot=pp_ax)
pp_ax.set_title('PP 图')
plt.tight_layout()return
# 模拟AR(1) 过程 #设置随机种子(括号里数字无意义)
np.random.seed(1)#模拟次数
n=5000#AR模型的参数
a = 0.8#扰动项为正态分布
x = w = np.random.normal(size=n)for t in range(1,n):
x[t] = a*x[t-1] + w[t]#画图
ts_plot(x, lags=30)
模拟的AR(1)模型是正态的。自相关系数图(ACF)显示滞后值之间存在显著的序列相关性，偏自相关系数图(PACF)则显示在滞后1期时截尾(迅速降为0)。下面使用statsmodels构建AR(p)模型，先用AR模型拟合上述模拟的数据，并返回估计的系数参数)，然后选择最佳滞后阶数，最后与原模型设置对比看是否选择了正确的滞后项。假如AR模型是正确的，那估计的系数参数将很接近真实的系数0.8，选择的阶数也会等于1。#估计数据的AR模型参数和滞后阶数
def simu_ar(data,a,maxlag=30,true_order = 1):
'''data:要拟合的数据；a为参数,可以为列表；maxlag:最大滞后阶数'''
# 拟合AR(p)模型
result = smt.AR(data).fit(maxlag=maxlag, ic='aic', trend='nc')
#选择滞后阶数
est_order = smt.AR(data).select_order(maxlag=maxlag,
ic='aic', trend='nc')
#参数选择标准ic : 有四个选择 {‘aic’,’bic’,’hqic’,’t-stat’}
#趋势项：trend：c是指包含常数项，nc为不含常数项
#打印结果
print(f'参数估计值：{result.params.round(2)}，
估计的滞后阶数：{est_order}')
print(f'真实参数值：{a}，真实滞后阶数 {true_order}')
simu_ar(x,a=0.8)
参数估计值：[0.8]，估计的滞后阶数：1
真实参数值：0.8，真实滞后阶数 1
看下如何用AR(p)模型来拟合沪深300的对数收益
# Select best lag order for hs300 returns
import tushare as ts
token='输入token'
pro=ts.pro_api(token)
df=pro.index_daily(ts_code='000300.SH')
df.index=pd.to_datetime(df.trade_date)
del df.index.name
df=df.sort_index()
df['ret']=np.log(df.close/df.close.shift(1))
max_lag = 30
Y=df.ret.dropna()
ts_plot(Y,lags=max_lag,title='沪深300')
result = smt.AR(Y.values).fit(maxlag=max_lag, ic='aic', trend='nc')
est_order = smt.AR(Y.values).select_order(maxlag=max_lag,
ic='aic', trend='nc')
print(f'沪深300拟合AR模型的参数：{result.params.round(2)}')
print(f'沪深300拟合AR模型的最佳滞后阶数 {est_order}')
沪深300拟合AR模型的参数：[ 0.03 -0.03 ...]
沪深300拟合AR模型的最佳滞后阶数 15
最好的阶数选择是15或者说有15个参数！任何模型有这么多参数在实际中不可能有用。显然有比这个模型更好的模型可以解释沪深300收益率走势。03MA模型MA(q)模型与AR(p)模型非常相似。不同之处在于，MA(q)模型是对过去的白噪声误差项的线性组合，而不是过去观测的线性组合。MA模型的动机是我们可以直接通过拟合误差项的模型来观察误差过程中的“冲击”。在一个AR(p)模型中，通过在一系列过去的观察中使用ACF间接观察到这些冲击。MA(q)模型的公式是：
下面使用Python模拟MA(1) 过程。
#这里使用arma模型进行模拟，设定ar阶数为0，即得到ma模型
alphas = np.array([0.])
betas = np.array([0.6])
ar = np.r_[1, -alphas]
ma = np.r_[1, betas]
#模拟MA的样本数据
ma_sample = smt.arma_generate_sample(ar=ar, ma=ma, nsample=1000)
ts_plot(ma_sample, lags=30,title='MA(1)模型')
ACF函数显示滞后1阶系数显著异于0，表明MA(1)模型适合拟合的数据。
# 对上述模拟数据进行ARMA模型拟合
max_lag = 30
result = smt.ARMA(ma1, order=(0, 1)).fit(maxlag=max_lag,
method='mle', trend='nc')
print(result.summary())
模型估计d 滞后系数为0.6277，与真实值0.6比较接近。注意到，95％置信区间确实包含该真实值。
下面尝试用MA(3)模型去拟合沪深300股价的对数收益，但这次并不知道真实的参数值。结果显示，拟合的残差自相关系数和偏自相关系数比较符合白噪声过程，但由于存在厚尾，MA模型并不是预测沪深300未来回报的最佳模型。
max_lag = 30
result=smt.ARMA(Y.values,order(0,3)).fit(maxlag=max_lag,
method='mle', trend='nc')
print(result.summary())
resid=pd.Series(result.resid,index=Y.index)
ts_plot(resid, lags=max_lag,title='沪深300指数MA拟合残差')
ARMA模型全称为自回归移动平均模型Autoregressive Moving Average Models - ARMA(p, q)，是AR(p)和MA(q)模型之间的结合，从金融的角度理解，AR和MA模型的理论意义在于：AR(p)模型试图捕捉(解释)交易市场中经常观察到的动量和均值回复效应。MA(q)模型尝试捕捉(解释)在白噪声条件下观察到的冲击效应。这些冲击效应可以被认为是影响观察过程的意外事件。ARMA模型的弱点在于忽视了大多数金融时间序列中的波动聚集效应。模型的公式可以表示为:
print(result.summary())
# 下面使用ARMA(2, 2) 模型进行模拟分析
max_lag = 30
n = 5000
burn = int(n/10)
alphas = np.array([0.5, -0.25])
betas = np.array([0.5, -0.3])#注意ar模型1代表0阶(自身)，然后在其他系数前加负号
ar = np.r_[1, -alphas]
ma = np.r_[1, betas]
arma22 = smt.arma_generate_sample(ar=ar, ma=ma, nsample=n, burnin=burn)
_ = ts_plot(arma22, lags=max_lag)
result = smt.ARMA(arma22, order=(2, 2)).fit(maxlag=max_lag,
method='mle', trend='nc', burnin=burn)
结果显示模型估计的参数与真实参数基本上吻合。下面使用ARMA模型来拟合沪深300的收益数据。ACF和PACF没有显示出明显的自相关性。QQ和概率图显示残差大致为正态分布但厚尾。总体而言，这个模型的残差看起来不像白噪声，说明模型还是没有很好的拟合其波动性特性。#不事先确定滞后阶数，而是通过信息准则选择最佳的滞后阶数#先将初始值设置为无穷大
best_aic = np.inf
best_order = None
best_mdl = None
rng = range(5)for i in rng:for j in rng:try:
tmp_mdl = smt.ARMA(Y.values, order=(i,j))
.fit(method='mle', trend='nc')
tmp_aic = tmp_mdl.aicif tmp_aic best_aic = tmp_aic
best_order = (i, j)
best_mdl = tmp_mdlexcept: continue
print(f'最佳滞后阶数:{best_order}')
print(best_mdl.summary())
resid=pd.Series(best_mdl.resid,index=Y.index)
ts_plot(resid, lags=30,title='沪深300指数ARMA拟合残差')
最佳滞后阶数:(4, 4)
ARIMA模型全称是差分移动自回归模型(Autoregressive Integrated Moving Average Models)，是ARMA模型的拓展。由于现实中很多时间序列不是平稳的，但可以通过差分来实现平稳，即通过一阶差分可以将非平稳机游走其转化为平稳的白噪声。由于前三个模型都有时间序列平稳的假设在，如果时间序列存在明显的上升或者下降趋势，模型预测的效果大大折扣。对于有明显下降或者上升趋势的数据集，可以使用差分的方式将其转化为平稳序列，然后使用ARMA模型进行拟合。假设模型经过d次差分通过了时间序列平稳的检验，ARMA的系数为p，q，ARIMA模型为ARIMA(p，d，q)。
下面通过迭代(p，d，q)的不同组合，找到拟合沪深300收益率数据的最佳ARIMA模型。通过AIC信息准则来评估每个模型，最后选取AIC最小的。
#原理与拟合ARMA模型类似
best_aic = np.inf
best_order = None
best_mdl = None#假定最多滞后4阶
pq_rng = range(5) #假定最多差分一次
d_rng = range(2) for i in pq_rng:for d in d_rng:for j in pq_rng:try:
tmp_mdl = smt.ARIMA(Y.values, order=(i,d,j))
.fit(method='mle', trend='nc')
tmp_aic = tmp_mdl.aicif tmp_aic best_aic = tmp_aic
best_order = (i, d, j)
best_mdl = tmp_mdlexcept: continue
print(f'ARIMA模型最佳阶数选择：{best_order}')# 对拟合残差进行可视化
print(best_mdl.summary())
resid=pd.Series(best_mdl.resid,index=Y.index)
_ = ts_plot(resid, lags=30,title='沪深300指数ARIMA残差')
ARIMA模型最佳阶数选择：(4, 0, 4)
最好的模型是差分为0，因为我们使用的是收益率数据，相对于已经采用了第一次对数差分来计算股票收益率。模型残差图结果与上面使用的ARMA模型基本相同。显然，ARIMA模型同样无法解释时间序列中的条件波动性。到这一步，时间序列的基本模型和建模步骤基本上大家已经熟知，下面利用模型的forecast()方法进行预测。
# 对沪深300收益率未来20天进行预测
n_steps = 20#分别设置95%和99%的置信度
f, err95, ci95 = best_mdl.forecast(steps=n_steps)
_, err99, ci99 = best_mdl.forecast(steps=n_steps, alpha=0.01)
date=(df.index[-1]).strftime('%Y%m%d')
cal=pro.trade_cal(exchange='', start_date=date)
idx = cal[cal.is_open==1][:20]['cal_date'].values
fc_95 = pd.DataFrame(np.column_stack([f, ci95]),
index=idx,columns=['forecast', 'lower_ci_95', 'upper_ci_95'])
fc_99 = pd.DataFrame(np.column_stack([ci99]),
index=idx, columns=['lower_ci_99', 'upper_ci_99'])
fc_all = fc_95.combine_first(fc_99)
#fc_all.head()
# 对预测的20日收益率数据进行可视化
plt.style.use('ggplot')
fig = plt.figure(figsize=(12,7))
ax = plt.gca()
ts = df['ret'][-500:].copy()
ts.plot(ax=ax, label='HS300收益率')# 样本内预测
pred=best_mdl.predict(np.arange(len(ts)) [0], np.arange(len(ts))[-1])
pf=pd.Series(pred,index=ts.index)
pf.plot(ax=ax, style='r-', label='样本内预测')
fc_all.index=pd.to_datetime(fc_all.index)
fc_all.plot(ax=ax)
plt.fill_between(fc_all.index, fc_all.lower_ci_95,
fc_all.upper_ci_95, color='gray', alpha=0.7)
plt.fill_between(fc_all.index, fc_all.lower_ci_99,
fc_all.upper_ci_99, color='gray', alpha=0.2)
plt.title('{} 天HS300收益率预测\nARIMA{}'.format(n_steps,
best_order))
plt.legend(loc='best', fontsize=10)
plt.show()
参考资料：
1. statsmodels官方文档
2. Time Series Analysis (TSA) in Python - Linear Models to GARCH
关于Python金融量化专注于分享Python在金融量化领域的应用。加入知识星球，可以免费获取量化投资视频资料、量化金融相关PDF资料、公众号文章Python完整源码、量化投资前沿分析框架，与博主直接交流、结识圈内朋友等。
