量化交易全流程（三）

輕栀

已于 2023-10-02 00:17:43 修改

阅读量637

点赞数

文章标签： python

于 2023-09-17 19:15:33 首次发布

本文链接：https://blog.csdn.net/mnwl12_0/article/details/132949766

版权

本节目录

收益率、波动率

夏普比率

索提诺比率

阿尔法和贝塔

最大回撤

利率的计量

债券定价

债券的剥息定价

远期利率

久期

期权

平稳性

自相关系数

混成检验

AR模型（及p的确定）

信息准则（拟合优度）

MA模型（ARMA）

ARCH和GARCH模型

-----------------------------------------------金融基础概念-------------------------------------------------------

主要介绍关于金融量化分析的一些基本概念，比如对数收益率、年化收益、波动率等。在严格的量化分析中，这些基本概念是必须要掌握并熟练使用的。
市面上，很多策略回测笼统地使用所谓的"胜率"来确定策略的好坏，这是一种非常业余而且不可靠的分析方法。在衡量基金业绩好坏的时候，大部分人也只是看基金的年化收益，忽略了基金的风险情况（波动率）。
市场中充斥着大量类似的业余的分析方法，这些方法导致很多所谓的回测看起来很美好，其实在统计上根本站不住脚，即所谓的"统计骗局"。
在对基本的概念掌握得足够熟练之后，我们很容易就可以识破很多统计骗局。经典的统计骗局有如下几种：只告诉你基金的年化收益，而不提基金的波动率或者最大回撤；使用周收益率来计算夏普比率，而不是使用日收益率来计算。
熟练掌握基本的概念，可以让我们对基金业绩的评估有更为全面理性的认识。

收益率
在学术界，建模一般不直接使用资产价格，而是使用资产收益率（Returns)。因为收益率比价格具有更好的统计特性，更便于建模。下面我们来看一下经典的收益率算法。
假设Pt表示在时刻 t 时一种资产的价格，在没有利息的情况下，从时刻t-1到时刻t这一持有阶段的收益率为：

Rt=(Pt-Pt-1)/Pt-1

其中，分子表示资产在持有期内的收人或利润，如果该值为负，则表示亏损。分母Pt-1表示持有资产初期的原始投资。

对数收益率（Log returns)，用rt表示,定义如下：

$r_{t}=ln(1+R_{t})=ln(\frac{P_{t}}{P_{t-1}})$

其中In(x）表示自然对数，即以e为底的对数。对数收益率比简单的收益率更为常见，因为对数收益率具有三个良好的数学性质，具体如下。

口当x比较小的时候（比如小于10％时）,In(x）和x的值是很接近的。
口使用对数收益率，可以简化多阶段收益。k阶段总的对数收益就是k阶段的对数收益之和。
口将对数收益绘制成图表，在直观上更接近真实的表现。比如股票价格从1元涨到10元，相当于翻了10倍，再从10元涨到100元，也是翻了10倍。如果单纯绘制股票价格，那么从10元涨到100元的这一段明显会"看起来涨了更多"。但是如果换算成对数价格，那么就不会存在这种直观偏差了。

年化收益
年化收益（Annualized Returns）表示资产平均每年能有多少收益。我们在对比资产的收益的时候，需要有一个统一的标准。计算方式是：

（最终价值／初始价值 - 1)／交易日数量 X 252

其中，252代表每年有252个交易日，这个数字每年都不一样，但业界为了方便，一般都将其固定为252，即：

(Pt-Pt-1)/days X 252

年化收益的一个直观的理解是，假设按照某种盈利能力，换算成一年的收益大概能有多少。这个概念常常会存在误导性，比如，这个月股票赚了5%，在随机波动的市场中，这是很正常的现象。如果据此号称年化收益为5%x12个月＝60%，这就显得不太可信了，实际上每个月的收益不可能都这么稳定。所以在听到有人说年化收益的时候，需要特别留意一下具体的情况，否则很容易被误导。

波动率
波动率（Volatility）和风险，可以算是一对同义词，都是用来衡量收益率的不确定性的。波动率可以定义为收益率的标准差，即：

δ= Std(r)

假设不同时间段的收益率没有相关性（称为没有自相关性），那么可以证明的是，收益率的方差Var(r）具有时间累加性。时间累加性的意思是，不同时间段t1，... ,tn的方差，加总即可得到t1，... ,tn这段时间的方差。

换句话说，随着时间的增加，方差将会成正比增加，波动率（标准差）将会按时间开根号的比例增加。举个例子，假设股票收益率的日波动率为δ，那么股票每年的波动率就为√252δ（假设年有252个交易日）。
这种不同周期间的波动率换算，在投资计算中非常常见。最常用的波动率是年化波动率，我们经常需要将日波动率、月波动率换算成年化波动率。

夏普比率

在研发策略的时候，经常会接触到各种各样的指标，这些指标代表了策略（资产、基金）的表现，比如下面将要介绍的夏普比率（Sharpe Ratio)。
关于夏普比率（Sharpe ratio)，具有投资常识的人都明白，投资光看收益是不够的，还要看承受的风险，也就是收益风险比。夏普比率描述的正是这个概念，即每承受一单位的总风险，会产生多少超额的报酬。用数学公式描述就是：

$SharpeRatio=\frac{E(R_{p})-R_{f}}{\delta_{ p}}$

其中各参数说明如下。
E(Rp)：表示投资组合预期收益率。
Rf：表示无风险利率。
δp：表示投资组合的波动率（亦即投资组合的风险）。
上面三个值一般是指年化后的值，比如预期收益率是指预期年化收益率。

需要注意的是，虽然公式看起来很简单，但是计算起来其实并不容易。原因就是预期收益率E(R）和波动率δp，其实是无法准确得知的。我们只能用统计方法来估计这两个值，然而估计方法也有很多种。估计E(R）和δp最简单的方法就是计算历史年化收益率及其标准差。然而，即使是同一种方法，针对不同周期计算出来的结果也可能存在很大的差别，从而产生误导。

知乎网站上有一个很经典的问题："夏普比率和最大回撤到底怎么计算"?
如题，这里列举一个极端例子来进行说明，三个交易日。
我9月1日一大早进场，带了100万，9月1日亏了999999元。
9月2日﹣大早我就只剩1元钱了，9月2日赚了一元钱，变成了2元钱。
9月3日我又赚了一元钱，变成了3元钱。
9月3日收盘时我决定退市。
上面所列举的这个极端情况中，最大回撤和夏普比率怎么计算？

提问者理解了公式的意思，但在具体计算的时候，却产生了迷惑。提问者之所以会产生迷惑，是因为他将收益率和波动率的周期弄混了。

同一项资产，用不同周期频率的收益率，计算出来的夏普值，根本就不是一回事！比如，利用每日的收益率计算夏普值，与利用每年的收益率计算夏普，就不是一回事。而且在计算的时候，收益率和波动率周期是要一致的，你不能用日线数据计算收益率，然后用周线数据计算波动率。

在上面的问题里，提问者其实是将周期的概念弄混了，如果你要按三天一个样本来计算收益，那么必须也要按三天一个样本的频率来计算波动率，然而，这里按三天一个周期来算的话，就只有一个样本，是无法计算波动率的。所以只能按每日的收益来进行计算。

事实上，自己手动算过一遍后，就不容易被误导了，也很容易发现别人计算中存在的问题。

举个例子来说明一下，我们先生成一组收益率数据，代码如下：

import pandas as pd
import numpy as np
year_list=[]
month_list=[]
rtn_list=[]

# 生成对数收益率，以半年为周期

for year in range(2006,2017):
    for month in [6,12]:
        year_list.append(year)
        month_list.append(month) # 除法运算，将月份除以6再除以10。
        rtn=round((-1)**(month/6) * (month/6/10),3)+(np.random.random()-0.5)*0.1
        rtn_list.append(rtn)

# 生成半年为周期的收益率df 22个数据点，
df=pd.DataFrame ()
df['year']=year_list
df['month']=month_list
df['rtn'] = rtn_list

# 计算SP
# 年化的值，收益率乘以2，波动率乘以sqrt(2)
round(df['rtn'].mean()/df['rtn'].std() *np.sqrt (2),3) # 假设无风险收益率为0

现在我们将数据变换成以年为频率的收益率。使用groupby方法，代码如下：

# 生成每年的收益数据dt_year（对数收益率可以直接相加）
df_year = df.groupby(['year']).sum()
del df_year['month']

# 计算其夏普比率
round(df_year['rtn'].mean()/df_year['rtn'].std(),3)

可以看到，同样的收益率数据，使用不同的周期，计算出来的结果之间的差距非常大。一般来说，周期频率越小，越难以保持收益稳定，每天都盈利比每年都盈利困难太多了。我们可以想象一种极端情况，10年中，每年的收益都是10%，夏普值就是无穷大，因为收益完全稳定，没有任何波动，然而每月的收益又不完全相同，所以从每月的收益率来看，夏普值并不是无穷大的。

所以在看夏普值的时候，一定要留意这个夏普值的计算方式，否则很容易产生误判。

自行计算的话，并没有强行的标准，只是需要注意如下两点。

一是要结合自己的实际，比如，高频策略当然得用日收益率，每周调仓的策略则可以使用周收益率。

二是对比策略优劣的时候，周期要一致，比如，对比每日调仓的策略和每月调仓的策略，一定要换算到同一个周期上，才有可比性。

索提诺比率

索提诺比率与夏普比率相似，不一样的是，索提诺比率是使用下行风险来衡量波动率的。在夏普比率中，资产大涨与资产大跌都可视为波动风险。

实际上，有时候大涨并不算风险，大跌才是风险，比如基金净值，所以索提诺比率只考虑大跌的风险，这也可以看作是对夏普比率的一种修正方式。

但是在某些品种中，大涨大跌都可能是风险，比如可以做多做空的期货。美股也可以做多或者做空，这种情况下，涨或跌都是风险。索提诺比率之所以没有流行，大概是因为美股可以做空。

阿尔法和贝塔

关注量化的读者，可能经常会听到有人谈论阿尔法（Alpha）策略。所谓的阿尔法策略，其实是来源于资本资产定价模型（CAPM)。这个模型将股票的收益分为了两个部分，一部分是由大盘涨跌带来的，另一部分则是由股票自身的特性带来的。大盘的那部分影响就是贝塔（Beta）值，剔除大盘的影响，剩下的股票自身就是Alpha值。
所以在谈论Alpha策略的时候，其实就是在谈论股票与大盘无关的那部分收益。如果Alpha策略做得好，对冲掉大盘风险后，可以取得相当稳定的收益。

最大回撤

顾名思义，是指投资一项资产，可能产生的最大亏损，即所谓的"买在最高点，抛在最低点"。
计算公式如下：

max(1﹣当日净值／当日之前最高净值）

这个max需要对每个交易日进行循环。

下面给出一段计算最大回撤的代码，算法复杂度为O(n)。这里的df是常见的K线数据（包含open、high、low、close字段）。

第一步是循环计算出截至每个交易日之前的最高值，第二步是计算出每个交易日的最大回撤，第三步是在所有交易日中选出极值。示例代码如下：

# 记录截至每个时点的前最高值

for idx,row in df.iterrows():
    max_c=max(row['high'].max_c)
    df.loc[idx,'max_c']=max_c

# 计算每个交易日的最大回撤

df['max_dd']=df['low']/df['max_c']-1

# 得到所有交易日中的最大回撤

max_dd=abs(df['max_dd'].min())

--------------------------------------------------------资产定价------------------------------------------------------

利率
利率是决定金融产品价格的一个基本因素。讨论利率的基本分析方式，首先，定义利率所用的复利频率和连续复利利率的含义，然后介绍零息利率、平价收益率、收益率曲线以及债券定价分析的内容。

利率定义了资金借入方承诺支付给资金借出方的资金数量。利率包含了很多种类，比如存款利率、贷款利率、国债利率等。利率与信用风险有关，信用风险越高，利率就越高。国债收益率是投资者投资国债所得到的收益率。国债是政府借入以本国货币为计量单位而发行的产品。通常我们认为政府不会对其债务违约，因此国债利率可以视为无风险利率。

LIBOR 是伦敦同业银行的拆借利率（London Interbank Offered Rate)，这是一个参照利率，每天的拆借利率由英国银行家协会提供。该利率是指伦敦一流银行之间短期资金借贷的利率，是国际金融市场中大多数浮动利率的基础利率。SHIBOR是我国仿照LIBOR的模式，建立的上海同业拆借利率。SHIBOR报价银行现由18家商业银行组成。全国银行间同业拆借中心授权SHIBOR的报价计算和信息发布。每个交易日根据各报价行的报价，剔除最高、最低报价各4家，然后对其余报价进行算术平均计算之后，得出每一期限品种的SHIBOR，并于每天早上的9:30对外发布。

无风险利率广泛应用于金融产品的定价。无风险利率的选取并没有一个完全客观的标准，在美国，有交易员用一年期的国债利率作为无风险利率，但也有交易员认为政府债券隐含的利率是偏低的，所以很多金融机构将LIBOR作为无风险利率。2007年发生信用危机导致LIBOR利率激增，许多市场参与者开始使用隔夜指数互换利率（Overnight Indexed Swap,OIS）作为对无风险利率的近似取值。在国内，SHIBOR和10年期国债利率都可以作为无风险利率。

利率的计量

当我们谈论年利率的时候，根据复利频率的不同，会有不同的表现形式。比如，假设年利率是10%，那么可以是一年按10％复利一次，也可以是半年按5％复利两次，还可以是每3个月按2.5％复利4次。

推广以上结果，假设将A资金投资n年。如果利率是按年复利的，那么投资的终值为：

$A(1+R)^{^{n}}$

如果利率是对应于一年复利m次，那么投资终值为：

$A(1+\frac{R}{m})^{^{n}}$

下面用Python程序来对复利m次的投资终值进行一个模拟操作。假设初始投资为100元，年利率为10%，投资1年，复利次数从1到15。程序代码具体如下：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data = pd.DataFrame()
data['frequency'] = np.arange(1,15,1)
data['final_value'] = 100*((1+0.1/data['frequency'])**data['frequency'])
data.plot(figsize = (10,6),x = 'frequency', y='final_value',style = 'o-')

由上图可以看到，随着复利频率的增加，投资终值趋近于一个极限。可以从数学上进行证明，当复利频率趋于无穷大的时候，投资的终值为：

$Ae^{^{Rn}}$

其中e=2.71828。这种情况下对应的利率被称为连续复利（continuous compounding）利率。在实际情况下，由于每天计算复利的频率非常高（每年复利365次），非常接近于连续复利，因此可以认为普通复利的计算方法与连续复利的计算方法是等价的。

对一笔资金，以利率R连续复利n年，相当于乘以 $e^{^{Rn}}$ ：反过来我们可以计算贴现值。"贴现值"是一个较为复杂的概念，为了便于理解，下面举例说明。

比如，我给你一张支票，这张支票你必须要等一年之后才能拿到1万元钱，你愿意花多少钱买这张支票？这个钱肯定比1万元要少，因为当前的1万元现金放在银行是有利息的，放在银行一年后可能是10400 元。那么这张支票值多少钱呢？未来到手的资产在当前值多少钱，就是某资产的贴现值。一笔n年后的资金如果要计算贴现值，则以利率为R按连续复利进行贴现，相当于乘以 $e^{^{-Rn}}$ 。

连续复利在衍生品产品定价中的使用非常广泛，所以除非特别指明，利率均将按照连续复利来计算。

假设Rc是连续复利利率，Rm是与之等价的每年m次复利利率。可以得到：

用这个公式可以对m次复利利率和连续复利利率进行相互转换，公式如下，即：

零息利率
N年的零息利率（zero rate）是指今天投入的资金在N年后所得的收益率。所有的利息以及本金在N年末支付给投资者，在N年满期之前，投资不支付任何利息收益。零息利率也称作即期利率（spot rate)。
假如一个10年期连续复利的零息利率是每年5%，那么100元投资在10年后将增长到：

100 х $e^{0.05X10}$ =164.87

市场上很多都不是零息利率，大多会在利息周期内支付部分利息。

债券定价

任何金融产品的理论价格都等于未来现金流的折现值，大多数债券提供周期性的券息，发行人在到期前会把债权本金偿还给投资者，有时，投资者使用同样的贴现率对债券的现金流进行贴现，但更精确的办法是对不同的现金流采用不同的零息贴现率。

假设下表给出的是零息利率，其中的利率是连续复利。假设一个两年期债券面值为100元，券息为6%，每半年付息一次。为了计算第一个3元券息的现值，我们用5％的6个月贴现率进行贴现操作；为了计算第二个3元券息的现值，我们用5.8％的1年贴现率进行贴现操作，以此类推。因此债券理论价格为：
3exp(-0.05x05)+3exp(-0.058x1)+ 3exp(-0.064x1.5)+ 103(-0.068x2)=98.39

国债零息利率
期限（年）零息利率（连续复利，%)
0.5 5
1.0    5.8
1.5    6.4
2.0   6.8

债券收益率
债券收益率又称为到期收益率（Yield To Maturity,YTM)。债券收益率等于对所有现金流贴现并使债券的价格与市场价格相等的贴现率。假定上面考虑的债券理论价格也等于市场价格，即98.39元，如果用y表示按连续复利的债券收益率，那么我们可以得到：

3exp(-y*0.5)+ 3exp(-y*1)+3exp(-y*1.5)+103exp(-y*2)=98.39

下面使用Python 的Scipy 库来解这个方程，可以得到收益率y=6.76%，程序代码具体如下：

# 用 `fsolve` 函数，来找到满足方程 `3*e**(-y*0.5)+3*e**(-y*1)+3*e**(-y*1.5)+103*e**(-y*2)-98.39 = 0` 的 `y` 的值
# `0.1` 是作为 `fsolve` 函数的初始猜测值（initial guess）传递给方程的参数 `y`
# `fsolve` 函数通过迭代的方式寻找方程的解，并根据初始猜测值逐步逼近解的准确值
# 初始猜测值的选择可以影响算法的收敛速度和结果的准确性,`0.1` 是作为 `y` 的初始猜测值，用于启动解的搜索过程
from scipy.optimize import fsolve
import math
def func(y):
    e=math.e
    return 3*e**(-y*0.5)+3*e**(-y*1)+3*e**(-y*1.5)+103*e**(-y*2)-98.39
y=fsolve(func, 0.1)
print(y)

平价收益率
债券平价收益率（par yield）是使债券价格等于面值（par value）的券息率。假定债券每年的券息为c．每半年支付一次即c/2。采用之前的零息利率，可以得到以下方程：

解出c=6.87%，即平价收益率为6.87%。

国债零息利率确定
确定零息收益率曲线的一种常用方法是从国债价格入手来进行计算。最流行的方法就是所谓的票息剥离法（bootstrap method)。为了说明这种方法，下面首先来看一下下表中关于5个债券价格的数据。

票息剥离法原始数据，用如下代码创建dataframe，得到的结果如下：

import pandas as pd
import numpy as np
data=pd.DataFrame()
data['maturity']=[0.25,0.5,1,1.5,2]
data['coupon']=[0,0,0,8,12]
data['principal']=[100]*5
data['price']=[97.5,94.9,90,96,101.6]
data['zero_rate']=[np.nan]*5
print(data)

前三个债券不支付券息，对应的零息利率可以直接进行计算。比如，对于第一个债券（期限为0.25年），满足公式：

$100=97.5e^{^{R*0.25}}$
得到：R=10.127%

类似的0.5年和1年计算得到： $100=94.9e^{^{R*0.5}}$ ， $100=90e^{^{R*1}}$ 。

用代码计算具体如下：

# 由价格和现金流计算零息时的利率:face_value=price*exp(R*ΔT)➡R=ln(face_value/price)/ΔT
df=data
df.loc[0:2,'zero_rate']=df.loc[0:2,'principal']/df.loc[0:2,'price']
df.loc[0:2, 'zero_rate']= df['zero_rate'].apply(math.log)
df.loc[0:2,'zero_rate']=df.loc[0:2, 'zero_rate']/df.loc[0:2,'maturity']
print(df.loc[0:2,'zero_rate'])

得到的结果如下：

第四个债券，假设半年支付一次券息。那么债券带来的现金流如下。
6个月：4元
1年：4元
1.5年：104元
债券价格必须等于债券未来现金收入的现值。假设1.5年所对应的零息利率为R，那么
可以得到
4exp(-0.104693*0.5)+ 4exp(-0.10536*1)+4exp(-R*1.5)=96
由此得出R=0.10681，因此，1.5年所对应的零息利率为10.681%。
对应Python代码具体如下：

#计算半年的券息，半年付息
c=df.loc[3,'coupon']/2
# 计算券息的现值
coupon_pv=0
e=math.e

#计算每半年的付息现值 = ∑coupon*exp(-r*t)
for i in [1,2]:
    coupon_pv+=c*e**(-df.loc[i,'zero_rate']*df.loc[i,'maturity'])

# 计算第三年的零息利率
df.loc[3,'zero_rate']=-math.log((df.loc[3, 'price']-coupon_pv)/(100+c))/df.loc[3,'maturity']

print(df.loc[3,'zero_rate'])

同样地，对于2年期零息利率，假设该利率为R，用同样的方法，可以得到R=0.10808，即10.808%。

#计算半年的券息，半年付息
c=df.loc[4,'coupon']/2
print(c)
# 计算券息的现值
coupon_pv=0
e=math.e
print(coupon_pv)
#计算每半年的付息现值 = ∑coupon*exp(-r*t)
for i in [1,2,3]:  # 这里不是范围,是循环这里面的数
    coupon_pv+=c*e**(-df.loc[i,'zero_rate']*df.loc[i,'maturity'])
print(coupon_pv)
# 计算第三年的零息利率
df.loc[4,'zero_rate']=-math.log((df.loc[4, 'price']-coupon_pv)/(100+c))/df.loc[4,'maturity']

print(df.loc[4,'zero_rate'])

下表总结了计算结果。我们可以以期限为x轴，零息利率为y轴，绘制曲线图，该图形称为零息利率曲线（如图9-2所示）。由票息剥离法计算只能计算几个标准期限对应的零息利率，对于其他非标准期限，比如1.23年，则需要使用其他方法。

这里假设数值节点之间为线性关系，进而绘制零息利率曲线图。

print(df)
df.plot(figsize=(10,6),x='maturity',y='zero_rate',style='o-')

远期利率

远期利率（forward interest rate）是由当前零息利率所蕴含的将来一定期限的利率。下
面举例说明远期利率的计算方法，假设一组零息利率如下表所示。

期限（年）零息利率（每年，%) 远期利率（每年，%)
1 3

2 4 5

3 4.6 5.8

4 5 6.2

5 5.3 6.5
如果这些利率为连续复利，那么1年期的3％利率即意味着今天投资100元，一年后投资者收到 $100*e^{^{0.03x1}}$ =103.05元。2年期4％利率意味着今天投资100元，两年后投资收到 $100*e^{^{0.04x2}}$ =108.33元，以此类推。
上表中第2年的远期利率为每年5%。这个利率的意思是，在未来的第1年年末投资，到第2年年末的隐含利率。这个利率可以由1年期和2年期的零息利率计算得出。为了说明远期利率为5%，假定投资100元，第1年的利率是3%，第2年的利率是5%，那么在第2年年末的收益为
$100*e^{^{0.03x1}}*e^{^{0.05x1}}$ =108.33元
如果按照4％的零息利率计算，那么连续投资两年的收益为：
$100*e^{^{0.04x2}}$ =108.33元
可以看到，这两种投资最终的收益是相等的，所以1年期末到2年期末的远期利率为5%。进一步可得到一般结论：当利率按连续复利表示时，将相互衔接的时间段的利率结合在一起，整个时段的等价利率即为各个时段利率的平均值。比如，在上面的例子中，第一年的利率3％和第二年的利率5％求平均值将会得到两年的平均利率为4%。对于非连续复利的利率，这一结论只是近似成立。
第3年的远期利率可以由2年的零息利率4％与3年的零息利率4.6％隐含得出，结果为5.8%。其他远期利率也可以用类似的方法进行计算，结果如上表中的第3列所示。

$100*e^{^{0.04x2}}*e^{^{0.058x1}}$ =114.797

$100*e^{^{0.046x3}}$ =114.797

一般来说，如果 R和R2分别对应期限为 $T_{1}$ 和 $T_{2}$ 的零息利率， $R_{f}$ 为 $T_{1}$ 和 $T_{2}$ 之间的远期利率，那么有：

$R_{f}=\frac{R_{2}*T_{2}-R_{1}*T_{1}}{T_{2}-T_{1}}$

为了说明上式，考虑表中第4年的远期利率：T=3,T2=4,R1=0.046,R2=0.05，由公式可以得出R=0.062。

上式可以写成：

$R_{f}=R_{2}+\frac{(R_{2}-R_{1})*T_{1}}{T_{2}-T_{1}}$

上式说明如果零息利率曲线在 $T_{1}$ 和 $T_{2}$ 之间向上倾斜，即 $R_{2}>R_{1}$ ，那么 $R_{f}>R_{2}$ 。类似地，如果零息利率曲线向下倾斜，即 $R_{2}<R_{1}$ ，那么 $R_{f}<R_{2}$ 。令 $T_{2}$ 趋近于 $T_{1}$ ，并将共同值记为T，可以得到

$R_{f}=R_{2}+T*\frac{\partial R}{\partial T}$

上式R为期限为T的零息利率。Rf称为期限为T的瞬时远期利率（instantaneous forward rate)。这是用于在T开始的一段很短时间内的远期利率。

久期

债券的久期（duration）是指投资者收到所有的现金流所要等待的平均时间。一个n年期的零息债券久期为n年，而一个n年带息债券的久期小于n年，因为持有人在n年之前就已经收到部分现金付款了。假定一个债券在 $t_{i}$ 时刻给债券持有人提供的现金流为 $c_{i}$ ,1≤i≤n。债券价格B与连续复利收益率关系式为：

$B=\sum_{_{i=1}}^{n}c_{i}*e^{-yt_{i}}$

$D=\frac{\sum_{_{i=1}}^{n}t_{i}*c_{i}*e^{-yt_{i}}}{B}$

期权

期权是一项在一段有限时间内按一定价格买进或者卖出某种资产的权利。这里所说的某种资产称为标的资产（underlying asset)。
期权包含两种基本类型：看涨期权（call option）赋予了期权持有者在将来某一特定时刻以一定的价格买入某项资产的权利，看跌期权（put option）赋予了期权持有者在将来某一特定时刻以一定的价格卖出某项资产的权利。看涨期权也称为认购期权。看跌期权也称为认沽期权。目前上海证券交易所官方使用的就是认购期权和认沽期权。

期权合约中注明的日期称为到期日（expiration date）或者满期日（maturity date)，期权合约中所注明的价格称为执行价格（exercise price）或者敲定价格（strike price)。
期权可以是美式期权（American option）或者是欧式期权（European option)。美式期权可以在到期日之前的任何时刻行使，而欧式期权只能在到期日才能行使。

每个期权合约均包含如下4个基本属性。
类型（看涨期权或者看跌期权）。
标的资产。
到期日。
执行价格。
举个例子，国内的50ETF期权，其中一个合约代码是10000691.SH，简称"50ETF购12月2.348A"，这个合约包含如下4个属性。
口类型：看涨期权（认购期权）。
口标的资产：上证50交易型开放式指数证券投资基金（"50ETF")。
口到期日：2016年12月28日。
口执行价格：2.348元。
这个期权的意思是，如果投资者买入一手期权，那么可以在2016年12月28日，以2.348元的价格买人相应份额的50ETF。

看涨期权和看跌期权

考虑某投资者以5元的价格买人执行价格为100元的看涨期权。如果在到期日，股票价格小于 100元，那么投资者明显不会行使期权，此时投资损失的就是他买人期权的费用5元。如果股票价格大于100元，比如是110元，那么投资者将会行使期权，可以获取110-100=10元的收益，再减去5元的期权费用，最后获得5元的利润。
看涨期权持有人希望股票价格上涨，而看跌期权持有人则希望股票价格下跌。考虑某投资者以5元的价格购买了执行价格为100元的看跌期权。如果在到期日，股票价格大于100元，那么投资者不会行使期权。如果股票价格低于100元，比如是90元，那么投资者将会行使期权，可以获取100-90=10元的收益，再减去5元的期权费用，最后获得5元的利润。

期权价格与股票价格的关系

期权可以分为实值期权（in-the-money option)、平值期权（at-the-money option）及虚值期权（out-of-the-money option)。如果股票价格低于一手看涨期权的执行价格，那么这手看涨期权称为虚值期权；如果股票价格等于看涨期权的执行价格，那么看涨期权称为平值期权；如果股票价格大于看涨期权的执行价格，那么看涨期权称为实值期权。

假设S为股票价格，K为执行价格，对于看涨期权，当S<K时为虚值期权，S=K时为平值期权，S>K 时为实值期权。相反，对于看跌期权，S<K时为实值期权，S=K时为平值期权，S>K时为虚值期权。只有期权为实值期权时才会被行使。

一个期权的内涵价值（intrinsic value）是期权被行使时能产生的额外价值。如果一个期权是虚值的，那么内涵价值就是0。一个看涨期权的内涵价值为max(S-K,0)；一个看跌期权的内涵价值为max(K-S,0)。

影响期权价格的因素

有6种因素会影响股票期权的价格，具体列举如下。
当前股票价格 $S_{0}$ 。
执行价格K。
期权期限T。
股票价格的波动率δ。
无风险利率r（例如1年期国债利率）。
期权期限内预期发放的股息。
前4项是决定一个期权价格的主要因素，后两项一般没那么重要。不过，对红利（股息）高的股票来说，股息率可能会有相当大的影响。
为期权定价，最常用的公式就是布莱克—斯科尔斯—默顿（Black-Scholes-Merton）定价公式：

函数N(x）为标准正态分布的累积概率分布函数。c与p分别为欧式看涨期权与看跌期权的价格， $S_{0}$ 为股票在0时刻的价格，K为执行价格，r为连续复利的无风险利率，δ为股票价格的波动率，T为期权期限。
下面我们用Python来实现期权定价的函数。option_pricer 函数共有6个参数，每个参数在注释里面都有说明。此外，这里用到了scipy.stats 的 norm 模块，程序代码具体如下：

from scipy.stats import *
import numpy as np
def option_pricer(call_or_put, spot, strike, maturity, r, vol):
    """
    参数:call_or_put:1代表计算看涨期权价格,0代表计算看跌期权价格
    spot:标的资产价格    strike:执行价格      maturity:到期时间(年化)
    r:无风险利率         vol:波动率     返回：相应期权的价格
    """
    d1=(np.log(spot/strike) + (r + 0.5*vol*vol)*maturity)/vol/np.sqrt(maturity)
    d2=d1-vol*np.sqrt (maturity)

    if call_or_put==1:
        call_price=spot*norm.cdf(d1)-strike*np.exp(-r*maturity)*norm.cdf(d2)
        return call_price
    
    elif call_or_put==0:
        put_price=strike*np.exp(-r*maturity)*norm.cdf(-d2) -spot*norm.cdf(-d1)
        return put_price

现在考虑一个6个月的期权。股票的当前价格为50元，期权执行价格为40元。无风险利率为3%，年化波动率为20%。也就是说 $S_{0}$ =50,K=40, r=0.03, δ=0.2,T=6/12=0.5,调用以上函数。
假设该期为看涨期权，其价格为：

option_pricer(1, 50, 40, 0.5, 0.03, 0.2)

输出结果如下：

假设该期权为看跌期权，其价格为：

option_pricer(0, 50, 40, 0.5, 0.03, 0.2)

输出结果如下：

本部分介绍了利率、债券和期权的一些基本概念，并且使用Python实现了部分概念的计算。

---------------------------------------------金融时间序列分析---------------------------------------------

金融时间序列分析是一门比较成熟，同时也比较难掌握的学科。时间序列分析对投资者的数学功底有一定的要求，特别是关于理论方面的知识。为了突出Python在时间序列分析中的应用，将会尝试尽量少地涉及数学性的证明和推导。但也不可能完全不讨论理论问题，因为有些理论如果不理解清楚，就只是稀里糊涂地当黑盒使用，也很容易出错。虽然条条大路通罗马，量化不止这一条路，但多掌握些分析工具总是好的。

为什么用收益率而不是价格
我们知道，业界非常流行的技术分析是基于资产的价格数据来进行分析，比如求MACD、KDJ等，都是基于价格本身来进行的。然而在学术界，金融时间序列分析，通常是基于收益率数据的，而不是基于价格数据。
主要原因是可以假设收益率数据是平稳的（stationary)，而价格数据是不平稳的。平稳性是金融时间序列分析中一个非常重要的概念。平稳的时间序列有着更好的统计特性，更便于建模分析。平稳性，直观地理解，就是时间序列的统计特性不随着时间的变化而变化。后文将对平稳性给出严格的定义。这里举一个直观的例子，比如股票的收益率，每天的收益率都在﹣10％到10％之间波动，可以预计的是将来收益率也不会超出这个区间太远。然而股票的价格是完全不可预期的，很可能10年就涨了100倍。
所以对于习惯于使用技术分析的投资者来说，如果想要使用更严谨的金融时间序列分析方法，首先需要转变思维的习惯，即从分析价格转变为分析收益率。

金融时间序列定义
前面给出了收益率和对数收益率的定义，在进行实际分析的时候，这两种收益率都可以使用。这里针对对数收益率进行分析，下面先使用akShare获取沪深300指数2022到2023年的数据：

import akshare as ak
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime


data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()

print(df.head())

df['date'] = pd.to_datetime(df['date'])
plt.plot(df['rtn'])
plt.xlabel('date')
indices = np.linspace(0, len(df['date'])-1, 8, dtype=int)
dates = df['date'].tolist()
plt.xticks(indices, [dates[i].strftime('%Y-%m-%d') for i in indices], rotation=45, ha='right')
plt.show()

相对于价格而言，收益率序列要更加平稳一些。下面就来看一下平稳性的严格定义。

平稳性

通常我们会看到一个时间序列的波动表现为随机性，但是在一段时间内其统计特性却是保持不变的。对于一个随机过程，如果随着时间的变化，其表现的各个统计特性都不变，则称这个随机过程具有强平稳性。
在数学上，强平稳性是这样定义的，对于时间序列（rt)，r1, r2,…,rn,的分布与r1+m,r2+m,…,rn+m的分布相同。也就是说，这n个观测序列的概率分布不取决于它们的初始时间。强平稳性是一个很强的假设，因为它需要所有的统计特性在时间上保持不变。一般来说，我们会假设金融序列是弱平稳性的。
弱平稳性并不要求所有的统计特性一直保持不变，只要求均值、方差以及协方差随着时间变化不发生改变即可。这个条件用数学公式表达具体如下：

$E(r_{i})=\mu$ (μ是常数）
$Corr(r_{i},r_{i-l})=\gamma _{l}$ ( $\gamma _{l}$ 只依赖于l)

从而，均值和方差不会随着时间的变化而变化，并且两个观测量之间的联系只取决于它们之间的间隔。
平稳性是很多时间序列模型的基础。也就是说，在使用这些模型之前，一定要检验数据是否平稳。否则，建模的结果就是不可信的，比如，可能会出现"伪回归问题"。

即使假定一个时间序列是平稳的，建模也依然非常困难，因为存在着无穷多的参数需要估计。所以我们需要进一步简化模型，将其简化成一个只有有限个参数的模型，这样才能进行建模和估计。ARIMA 模型就是这样一类模型。最简单的ARIMA是AR(p）模型。当p=1的时候，就是最简单的AR(1）模型，下面就来介绍AR模型。在介绍AR模型之前，需要引入一个特殊的随机过程，即白噪声序列。另外，我们还需要了解自相关系数的概念，以下就来对这两个概念进行简单介绍。

白噪声序列
时间序列 { ${r_{i}}$ } 如果具有有限均值和有限方差，并且是独立同分布随机变量序列，那么就称 { ${r_{i}}$ } 为白噪声序列。直观地理解，每一个观测值，都像是重新摇了一次骰子而生成的。
白噪声的一个特例是高斯白噪声，即 { ${r_{i}}$ } 服从均值为0，方差为 $\delta ^{2}$ 的正态分布。数学表达式为：

$E(r_{i})=\mu$ (对任意i）

$Var(r_{i})=\delta ^{2}$ (对任意i）

$Corr(r_{i},r_{j})=0$ (对任意i≠j）

白噪声序列是平稳的，虽然才噪声本身无趣，但对于其他序列来说，他是一个基础。

自相关系数：衡量两个随机变量的线性相关性和相关程度，计算公式为：

$\rho _{x,y}=\frac{COV(X, Y)}{\delta _{x}\delta _{y}}$

相关系数很容易理解，重要的是自相关函数(ACF):

考虑一个弱平稳收益率序列 $r_{t}$ 。当r与它的过去值 $r_{t-l}$ 是线性相关关系时，就可以将相关系数的概念推广到自相关系数。 $r_{t}$ 与 $r_{t-l}$ 的相关系数则称为 $r_{t}$ 的间隔为l的自相关系数，通常记为 $p_{l}$ ，在弱平稳性的假设下，它只是l的函数。具体地说，定义：

$\rho _{l}=\frac{COV(r_{t}r_{t-l})}{\delta _{r_{t}}\delta _{r_{t-l}}}=\frac{\gamma _{l}}{\gamma _{0}}$

对于一个平稳时间序列{ $r_{t}$ },1≤t≤T，间隔为l的自相关系数的估计为：

我们称函数

$\widetilde{\rho _{1}};\widetilde{\rho _{2}};...\widetilde{\rho _{l}}$

为 $r_{t}$ 的样本自相关函数（ACF)。
可以通过如下代码绘制出收益率的ACF图：

import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline
fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot(111)
fig=sm.graphics.tsa.plot_acf(df['rtn'], ax=ax1, lags=50)

得到的结果如图所示：

可以看到，在10阶之前自相关度还比较高。

混成检验：

很多时候，我们倾向于认为通过历史收益率是很难预测将来的收益率的，也就是说，收益率的时间序列自相关系数为0。为了验证这个结论，需要使用相应的统计检验手段。混成检验（Portmanteau Test）经常用来检验 $r_{t}$ 的几个相关系数是否同时为零。

原假设H0： ${\rho _{1}}={\rho _{2}}=...={\rho _{l}}=0$
备择假设H1：对某i∈{1，...，m}，存在 ${\rho _{i}}\neq 0$

混成检验统计量：

$Q^{*}(m)=T(T+2)\sum_{l-1}^{m}\frac{\rho ^{2}_{l}}{T-l}$

决策规则是：

Q(m)> $\chi ^{2}_{\alpha }$ ，拒绝原假设H0

我们可以用Python 计算出Q(m）的 p-value，当 p-value 小于等于显著性水平α时拒绝H0，也就是说序列自相关系数为0不成立，即序列是具有自相关性的。

以下是对沪深300指数进行混成检验的代码：acf是导入import statsmodels.api as sm

import akshare as ak
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import statsmodels.api as sm

data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()

# 检验10个自相关系数
m = 10
acf,q,p = sm.tsa.acf(df['rtn'],nlags=m,qstat=True)
out = np.c_[range(1,11), acf[1:], q,p]
output = pd.DataFrame(out, columns=['lag','AC','Q','P-value'])
output = output.set_index('lag')
output

得到的结果如图所示。

取显著性水平为0.05，可以看到，所有的p-value 都大于0.05。所以我们没有理由拒绝 H0假设，认为此段收益率是序列不相关的。

AR(p）模型简介

最简单的具有相关性的平稳过程是自回归过程，在该过程中，观测值 $r_{t}$ 的定义是过去观测值的加权平均值与白噪声残差之和。以下用数据公式来进行说明。
假设当前的观测值，与前一个观测值之间存在某种相关关系。那么存在如下公式：

$r_{t}=\varphi _{0}+\varphi _{1}r_{t-1}+\alpha_{t}$

其中， $\alpha_{t}$ 是均值为0，方差为 $\delta ^{2}_{\alpha }$ 的白噪声序列，即服从分布White Noise(0, $\delta ^{2}_{\alpha }$ )。我们称 $r_{t}$ 为AR(1）的过程，其中AR代表Autoregressive，所以这个模型又称为一阶自回归模型。直观的理解是，当前观测值是前一个观测值线性的"记忆"函数加上了一个白噪声的随机扰动。
将AR(1）推广，我们可以得到AR(p）模型：

$r_{t}=\varphi _{0}+\varphi _{1}r_{t-1}+...+\varphi _{p}r_{t-p}+\alpha_{t}$

当这个模型表示给定过去的数据时，过去的p个值 $r_{t-i}$ (i=1,…,p）联合决定了 $r_{t}$ ，即当前观测值的条件期望。

AR(p）平稳性检验

要想使用AR(p）模型进行拟合和预测，首先要检测时间序列是否平稳。这里先介绍一下如何判断AR(p）模型是否平稳。我们先假定序列是弱平稳性的，则有如下公式：

$E(r_{i})=\mu$ (对任意i）

$Var(r)=\gamma _{0}$

$Corr(r_{t},r_{t-l})=\gamma _{l}$ ( $\gamma _{l}$ 只依赖于l)

由于 $\alpha_{t}$ 是白噪声序列，因此有：

$E(a_{t})=0$ (对任意i）

$Var(a_{t})=\gamma _{0}$ (对任意i）

$Corr(r_{t},r_{t-l})=\gamma _{l}$ (对任意i≠j）

因此由AR(p）两边同时求期望可以得到：

$E(r_{t})=\varphi _{0}+\varphi _{1}E(r_{t-1})+...+\varphi _{p}E(r_{t-p})$

根据平稳性性质有： $E(r_{t})=E(r_{t-1})=...=\mu$ ，从而有

$u=\varphi _{0}+\varphi _{1}u+...+\varphi _{p}u$

对应的方程为： $r_{t}=\varphi _{0}+\varphi _{1}r_{t-1}+...+\varphi _{p}r_{t-p}+\alpha_{t}$ ，即：

$r_{t}=\varphi _{0}+\varphi _{1}r_{t}L^{1}+\varphi _{2}r_{t}L^{2}+...+\varphi _{p}r_{t}L^{p}$

$1-\varphi _{1}x-\varphi _{2}x^{2}-...-\varphi _{p}x^{p}=0$

我们称这个方程的解的倒数为该模型的特征根。如果所有特征根的模都小于1，则该序列是平稳的。在Python 中，statsmodel 提供了AR模型的函数，用于拟合AR 模型。这里使用前面谈到的收益率数据来拟合AR 模型。示例代码如下：AR模型是导入import statsmodels as sm

import akshare as ak
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import statsmodels as sm

data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()

# 检验10个自相关系数
# m = 10
# acf,q,p = sm.tsa.acf(df['rtn'],nlags=m,qstat=True)
# out = np.c_[range(1,11), acf[1:], q,p]
# output = pd.DataFrame(out, columns=['lag','AC','Q','P-value'])
# output = output.set_index('lag')
# output

print(df)

# tsa库的AR模型参数只能使用array类型
rtn=np.array(df['rtn'])

# 生成模型，要指定阶数
model=sm.tsa.ar_model.AutoReg(rtn,15)

# 拟合
fit_AR =model.fit()

# 绘图，蓝色的是收益率，红色的是AR模型拟合的曲线
plt.figure(figsize=(10,4))
plt.plot (rtn, 'b', label='return')
plt.plot(fit_AR.fittedvalues, 'r', label='AR fit')
plt.legend()

# 查看滞后阶数,对应特征返程的根
len(fit_AR.roots)

得到的结果如图所示：

这里调用的模型要指定滞后阶数p的值。查看p的值示例代码如下：

len(fit_AR. roots)

也就是说这个模型是15 阶的，即p= 15。

下面画出模型的特征根，以检验平稳性，示例代码如下：

# 画单位圆
pi,sin,cos = np.pi,np.sin,np.cos
r1=1
theta=np.linspace(0,2*pi,360)
x1=r1*cos(theta)
y1=r1*sin(theta)

plt.figure(figsize=(6,6))
plt.plot (x1,y1,'k')

# 这里fit_AR.roots是计算的特征方程的解，特征根应该取倒数
roots=1/fit_AR.roots

# 画特征根
for i in range(len(roots)):
    plt.plot(roots[i].real,roots[i].imag,'.r',markersize=8)

plt.show()

得到的结果如图所示：

可以看到，15个特征根都在单位圆内，说明模都小于1，序列是平稳的。

AR(p) 如何确定参数 p

在实际应用中，AR(p）模型中的p是未知的，必须根据实际数据来决定。求解阶p的过程称为AR模型的定阶。解决这个问题一般有两种方法：第一种是利用偏相关函数（Partial Auto Correlation Function, PACF)；第二种方法利用某个信息准则函数。

偏自相关函数
考虑如下一连串的AR模型：

$r_{t}=\varphi _{0,1}+\varphi _{1,1}r_{t-1}+e_{1t}$

$r_{t}=\varphi _{0,2}+\varphi _{1,2}r_{t-1}+\varphi _{2,2}r_{t-2}+e_{2t}$

$r_{t}=\varphi _{0,3}+\varphi _{1,3}r_{t-1}+\varphi _{2,3}r_{t-2}+\varphi _{3,3}r_{t-3}+e_{3t}$

...

其中 $\varphi _{0,j}$ 是常数项， $\varphi _{i,j}$ 是 $r_{t-i}$ 的系数， $e_{jt}$ 是AR(j）模型的误差项。

第一个式子中的估计 $\widetilde{\varphi _{1,1}}$ 称为 $r_{t}$ 的间隔为1的样本偏自相关函数；第二个式子中的估计 $\widetilde{\varphi _{2,2}}$ 称为 $r_{t}$ 的间隔为2的样本偏自相关函数；第三个式子中的估计 $\widetilde{\varphi _{3,3}}$ 称为 $r_{t}$ 的间隔为3的样本偏自相关函数.

在具体确定AR(p）模型的p值时，往往会采用绘制PACF图的方式。示例代码如下：

fig=plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(111)
fig=sm.graphics.tsaplots.plot_pacf(df['rtn'],ax=ax1,lags=80)


# fig=plt.figure(figsize=(10,5))
# ax1 = fig.add_subplot(111)
# fig=sm.graphics.tsaplots.plot_acf(df['rtn'],ax=ax1,lags=80)

得到的结果如图所示：

在观察PACF图的时候，其实是存在主观判断的。比如，在上图中，我们既可以认为p=20（因为20之后的数据都在范围之内），也可以认为p=20是偶然情况，因为从8到20之间的数据都在范围之内。所以也可以认为p=8。实践中往往会结合其他的方法，综合考虑判断p的值。

信息准则：

有几种信息准则可用来决定AR过程的阶p，它们都是基于似然函数的。例如，常用的Akaike 信息准则（AIC)(Akaike,1973）的定义如下：

AIC= $\frac{-2}{T}$ ln(似然函数的最大值)+ $\frac{2}{T}$ (参数的个数)

其中，T是样本容量。对高斯AR(I）模型，AIC可以简化为：

AIC(l)= $ln(\widetilde{\delta ^{2}_{l}})+\frac{2l}{T}$

在实际应用中，我们首先要计算AIC(l)，其中，l=0,1,2, …, P。然后选择阶k，使AIC达到最小值。类似的还有BIC准则，这里不再详细介绍。
在Python中，我们可以使用ARMA(p,0）模型进行拟合，ARMA模型在后文中会详细介绍，此处了解即可。ARMA模型有两个参数，p代表了AR(p）模型，0代表没有MA部分。所以ARMA(p, 0) 模型实际就是AR(p) 模型。下面依然使用沪深300的收益率数据进行实验。示例代码如下：

import akshare as ak
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
import statsmodels as sm  # 拟合AR时是sm没有api
from statsmodels.tsa.ar_model import AR
from statsmodels.tsa.arima.model import ARIMA

data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()

rtn=np.array(df['rtn'])

aicList=[]
bicList=[]
# 从1到15，对每一阶都进行一个AR(p）模型的拟合，并得到相应的AIC数值
for i in range(1,16):
# 这里使用了ARMA模型，order 代表了模型的 (p,q）值，我们令 q始终为O，就只考虑了AR情况。
    order=(i,0)
    # tempModel=sm.ARMA(rtn,order).fit()
    tempModel = sm.tsa.ar_model.AutoReg(rtn,i).fit()
    aicList.append(tempModel.aic)
    bicList.append(tempModel.bic)
plt.figure(figsize=(10,4))
plt.plot (aicList,'r',label='AIC value')
plt.plot (bicList, 'b',label='BIC value')
plt.legend(loc=0)

可以看到，当p=1时，AIC和BIC同时都最小，所以可以认为p=1是合理的。

拟合优度：

根据AR(p)模型来看，模型拟合较好的一个必要条件是残差序列应该是白噪声，也就是说，模型拟合出来的数据减去收益率序列本身，得到的差值，应该是白噪声。
下面先求出残差序列，并绘图，程序代码如下：

# 查看滞后阶数,对应特征返程的根
len(fit_AR.roots)

# 计算残差
delta=fit_AR.fittedvalues-rtn[15:]
plt.figure(figsize= (10,6))
plt.plot (delta,'r',label=' residual error')
plt.legend()

绘图结果如图所示：

然后计算其acf和对应的p值：

#计算自相关系数及p-value
lags=10
acf,q,p = sm.tsa.api.acf (delta, nlags=lags,qstat=True)
out =np.c_[range(1,lags+1),acf[1:],q,p]
df=pd.DataFrame(out, columns=['lag','AC', 'O','p-value'])
df

得到的结果如图所示：

从p-value的值可以知道，该序列没有相关性，即序列接近白噪声。
衡量平稳模型拟合是否足够好的一个常用统计量是拟合优度，拟合优度的公式如下：

$R^{2}=1-\frac{\delta ^{2}_{\xi }}{\delta ^{2}_{rt}}$

它的值在0~1之间，越接近1，拟合效果越好。现在我们来计算AR模型的拟合优度，示例代码如下：

value=1-delta.var()/rtn[15:].var()
print (value)

得到的结果为:

可以看到，该模型的拟合优度不算好。

预测

我们可以将样本分为训练集和测试集。用训练集来拟合模型，用测试集来测试模型的准确度。测试集一般选取最后的t_len（这里的t_len=10）个样本。示例代码如下：

import statsmodels as sm
# 测试集长度
t_len=10
# 测试集
train=rtn[:-t_len]
# 训练集
test=rtn[-t_len:]
# 拟合模型
fit=sm.tsa.ar_model.AutoReg(train,15).fit()
# 获取对应的预测值
predicts=fit.predict(len(rtn)-t_len,len(rtn)-1,dynamic=True)
# 生成对比数据
df=pd.DataFrame()
df['real']=rtn[-t_len:]
df['predict']=predicts
df

结果如图：

可以看到，预测效果并不是很好，这很正常，找到一个好的收益率模型是非常困难的。

ARMA模型：ARMA是AR和MA模型的一种结合形式，因此先讨论MA模型：

MA模型：MA(q）模型的表达公式如下：

$r_{t}=c_{0}+\alpha _{t}-\theta _{1}\alpha _{t-1}--\theta _{2}\alpha _{t-2}-...--\theta _{q}\alpha _{t-q}$

其中， $c_{0}$ 是一个常数， $\alpha _{0}$ 是一个白噪声序列。
Python并未为MA模型专门提供模型，如果需要拟合MA模型，只需要使用ARMA 模型，将p置为0即可。这里演示一个简单的建模过程。
先使用ACF函数判断阶次q，程序代码如下：

import akshare as ak
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()
df=df.set_index('date')

fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot(111)
fig=sm.graphics.tsa.plot_acf(df['rtn'],ax=ax1,lags=60)

得到的结果如图所示：

可以看到，在8之前存在大量的高自相关系数（超过范围的圆点都是高自相关系数）。虽然21点也超过了，但10到21之间的点都没有超过，故可以认为21是意外情况，8才是临界值，所以可以假定阶次q为8。

使用MA模型拟合，并且查看预测效果，代码如下：

import akshare as ak
import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
import matplotlib.pyplot as plt

data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()
df=df.set_index('date')

fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot(111)
fig=sm.graphics.tsa.plot_acf(df['rtn'],ax=ax1,lags=60)
#  MA模型
rtn=df.rtn.values
# 测试集长度
t_len=20
# 测试集
train = rtn[:-t_len]
# 训练集
test=rtn[-t_len:]
# 拟合模型
order = (0, 0, 10)  # Example order parameter for ARMA model
arima_model = ARIMA(train, order=order).fit()

# Get corresponding predictions
predicts = arima_model.predict(start=len(train), end=len(train) + len(test) - 1, dynamic=True)

comp = pd.DataFrame()
comp['original'] = test
comp['predict'] = predicts
comp.plot(figsize=(10, 5))

得到的结果如图所示:

可以看到，预测并不准确。这是很容易理解的，MA模型过于简单，一般不会产生有效的拟合和预测。

ARMA模型公式
在某些应用中，我们需要通过高阶的AR或者MA模型才能充分地描述数据的动态结构。这时会有很多参数需要估计，问题也就变得烦琐了。为了克服这个困难，人们提出了自回归滑动平均（ARMA）模型。该模型的基本思想是将AR和MA模型的想法结合在一个紧凑的形式中，使所用参数的个数保持很小。对于金融中的收益率序列，直接使用ARMA 模型的机会较少。然而，ARMA模型的概念与波动率建模存在密切关系。事实上，广义的自回归条件异方差（GARCH）模型就可以认为是 ${\alpha ^{2}_{t}}$ 的ARMA 模型。

ARMA(p,q）模型的表达形式如下：

$r_{t}=\varphi_{0} +\sum_{i=1}^{p}\varphi_{i}r_{t-i}+{\alpha_{t}}-\sum_{i=1}^{q}\theta _{i}\alpha_{t-i}$

其中， $\alpha_{t}$ 是白噪声序列，p和q都是非负整数。AR和MA模型是ARMA(p,q）的特殊情形。当q=0的时候，ARMA(p,0）就是AR(p）模型。当p=0的时候，ARMA(0,q）就是MA(q）模型。

ARMA 模型阶次判定

一般来说，判断AR（p）p的阶数要使用PACF函数，判断MA（q）q的阶数要使用ACF函数，这里可以同时绘制PACF和ACF图形，示例代码如下：

import akshare as ak
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()
df = df.set_index('date')
fig = plt.figure(figsize=(10,8))
ax1 = fig.add_subplot (211)
fig = sm.graphics.tsa.plot_acf(df.rtn, lags=30, ax=ax1)
ax2 = fig.add_subplot (212)
fig = sm.graphics.tsa.plot_pacf(df.rtn, lags=30,ax=ax2)

结果如图：

观察上图可知，阶次为（8,6）是合理的。不过一般来说，阶次太高的话，拟合的计算量会比较大。下面我们再用信息准则判断一下，示例代码如下：

用指数来验证效果不太好，于是换成了600519贵州茅台，结果会好一点：

import akshare as ak
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

# data = ak.index_zh_a_hist(symbol="000001", period="daily", start_date="20220901", end_date="20230901")
data = ak.stock_zh_a_hist('600519', period="daily", start_date = "20220901", end_date = "20230901", adjust="qfq")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()
df = df.set_index('date')
fig = plt.figure(figsize=(10,8))
ax1 = fig.add_subplot (211)
fig = sm.graphics.tsa.plot_acf(df.rtn, lags=30, ax=ax1)
ax2 = fig.add_subplot (212)
fig = sm.graphics.tsa.plot_pacf(df.rtn, lags=30,ax=ax2)

sm.tsa.arma_order_select_ic(df.rtn.values,max_ar=10, max_ma=10,ic='aic') ['aic_min_order']

计算出来的结果是（4,10)，这是用信息准则判断出来的阶次。

建立ARMA 模型
下面用上述计算出来的阶次（4,10）来建立ARMA模型，并利用拟合的模型进行预测：

import akshare as ak
import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
import matplotlib.pyplot as plt
rtn = df.rtn.values
# 测试集长度
t_len = 20
# 测试集
train=rtn[:-t_len]
# 训练集
test = rtn[-t_len:]
# 拟合模型
order = (4,0,10)
arima_model = ARIMA(train, order=order).fit()

# 获取对应的预测值
predicts=arima_model.predict (len(rtn)-t_len, len(rtn)-1, dynamic=True)
comp=pd.DataFrame()
comp['original']=test
comp['predict'] = predicts
comp.plot(figsize = (10,5))

下图是预测值和实际值的对比图，可以看到，预测效果仍然不好。
效果仍然有差距。

ARCH和GARCH模型

之前我们讨论了AR模型和ARMA模型，这两个模型是对收益率数据的拟合，并且进行了预测。实际上，时间序列模型对波动率的应用更为广泛。
股票波动的一个特殊性是它不能被直接观测。我们只能通过收益率数据对波动率进行估计，常用方法是计算一段时间内收益率的标准差。

波动率的特征
对于金融时间序列，波动率往往具有以下特征。
1、存在波动率聚集现象（volatility cluster)，也就是说，波动率在一段时间内都比较高，在另一段时间内都比较低。
2、波动率以连续方式随时间发生变化，即波动率很少发生跳跃。
3、波动率不会发散到无穷，即波动率只在固定的范围内发生变化。从统计学角度来说，这意味着波动率往往是平稳的。
4、波动率对价格大幅上升和价格大幅下降的反应不同，这种现象称为杠杆效应（leverage effect)。

波动率模型框架

在对波动率进行建模的时候，我们需要同时考虑收益率的条件均值和条件方差。给定t-1时刻已知的信息集的条件均值和条件方差为：

$\mu _{t}=E(r_{t}|F_{t-1}) , \delta ^{2}_{t}=Var(r_{t}|F_{t-1})$

我们假设 $r_{t}$ 服从一个简单的时间序列模型，得到的平稳ARMA(p,q）模型如下：

$r_{t}=\mu_{t}+\alpha _{t} , \delta ^{2}_{t}=Var(r_{t}|F_{t-1})$

其中， $\alpha _{t}$ 为资产收益率在t时刻的"扰动"或"新息",ut为 rt的均值方程， $\delta ^{2}_{t}$ 的模型称为rt的波动率方程。可以看到，条件异方差性建模就是对时间序列增加一个动态方程，用来刻画资产收益率的条件方差随时间的演变规律。
对资产收益率序列建立一个波动率模型（ARCH）需要如下4个步骤。
1）通过检验数据的序列相关性建立一个均值方程，如有必要，对收益率序列建立一个计量经济模型（如ARMA模型）来消除任何的线性依赖。
2）对均值方程的残差进行ARCH效应检验。
3）如果ARCH效应在统计上是显著的，则指定一个波动率模型，并对均值方程和波动率方程进行联合估计。
4）仔细地检验所拟合的模型，如有必要则对其进行改进。

ARCH模型
在较早的经济学模型中，干扰项的方差被假设为常数。但是实际上，有很多经济时间序列呈现出的是波动聚集性，这种情况下假设方差为常数是不恰当的。
为波动率建模提供一个系统框架的第一个模型是Engle(1982年）提出的ARCH模型。ARCH采用自回归形式来刻画方差的变异。其基本思想具体如下。

1）资产收益率的扰动a，是序列不相关的，但不是独立的。
2)at的不独立性可以用其延迟值的简单二次函数来描述。
具体地说，一个ARCH(m）模型假定：

$a_{t}=\delta _{t}\xi _{t} , \delta^{2}_{t}=\alpha _{0}+\alpha _{1}a^{2}_{t-1}+...+\alpha _{m}a^{2}_{t-m}$

其中， $\xi _{t}$ 是均值为0、方差为1的独立同分布随机变量序列，a0>0，对于i>0，有ai≥0。系数a必须满足一些正则性条件以保证at的无条件方差是有限的。实际中，通常假定 $\xi _{t}$ 服从标准正态分布。
从模型的结构上看，大的"扰动"会倾向于紧接着出现另一个大的"扰动"。这一点与在资产收益率中所观察到的"波动率聚集"现象相似。
现在，我们用Python来进行ARCH建模。
ARCH和GARCH模型在另外一个模块中，Anaconda没有包含这两个模型，需要自己重新安装。安装方法为：

pip install arch

先导入相关的库，获取数据，计算收益率，代码如下：找了好几只股票效果不好，最终选了300002：

import akshare as ak
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import arch #条件异方差模型相关的库

# data = ak.index_zh_a_hist(symbol="000002", period="daily", start_date="20220901", end_date="20230901")
data = ak.stock_zh_a_hist('300002', period="daily", start_date = "20220901", end_date = "20230901", adjust="qfq")
df = data[['日期','开盘','最高','最低','收盘','成交量']]                  # 选取日期、高开低收价格、成交量数据
df = df.rename(columns={'日期': 'date','开盘': 'open','最高': 'high','最低': 'low','收盘': 'close','成交量':'volume'})
df['rtn'] = np.log(df['close'])-np.log(df['close'].shift(1))
df = df.dropna()
df = df.set_index('date')
rtn=df.rtn.values

# 然后检验收益率序列是否平稳，代码如下：
t=sm.tsa.stattools.adfuller(df.rtn) # ADF
print ("p-value: ",t[1])

# 假设该序列的均值方程是AR(p）模型
fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot (111)
fig=sm.graphics.tsa.plot_pacf (df.rtn, lags = 25,ax=ax1)

得到的结果为p-value:1.0249703008901893e-23。p-value小于0.05，说明序列是平稳的。

为了简化问题，我们假设该序列的均值方程是AR(p）模型。先判定阶次p，示例代码如下：

# 假设该序列的均值方程是AR(p）模型
fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot (111)
fig=sm.graphics.tsa.plot_pacf (df.rtn, lags = 25,ax=ax1)

通过上图可以假定阶次 p 为 1。现在来生成均值方程AR(1）模型，并绘制残差 $a_{t} , a^{2}_{t}$ 的图形，示例代码如下：

import statsmodels as sm

mean_model = sm.tsa.ar_model.AutoReg(rtn,1).fit()

rtn =rtn[:len(mean_model.fittedvalues)] #可能数据长度不一样
at = rtn - mean_model.fittedvalues
at2=np.square(at)
plt.figure(figsize=(10,6))
plt.subplot (211)
plt.plot (at,label='at')
plt.legend()
plt.subplot (212)
plt.plot (at2,label='at^2')
plt.legend(loc=0)

绘制结果如图所示:

下面来检验是否具有相关性，代码如下：

import statsmodels as sm
m=10 # 我们检验10个自相关系数
acf,q,p = sm.tsa.api.acf(at2,nlags=m,qstat=True) # 计算自相关系数及p-value
out = np.c_[range(1,m+1),acf[1:], q, p]
output = pd.DataFrame(out, columns=['lag', 'AC','Q','p-value'])
output=output.set_index('lag')
output

得到的结果如图所示:

从上图可以看出，p-value小于显著性水平0.05，所以认为序列具有相关性，即具有ARCH效应。现在我们来判定ARCH模型的阶次，代码如下：

fig=plt.figure(figsize=(10,5))
axl=fig.add_subplot (111)
sm.graphics.tsaplots.plot_pacf(at2,lags = 30,ax=ax1)

结果如图所示：

根据上图，我们可以将阶次定为5阶（从第6个点开始，都没有超出范围）。所以
最终的选择是均值方程为AR(5）模型，波动率模型为ARCH(5）模型。
下面来拟合模型，示例代码如下：通过以下命令显示拟合模型的结果：

am=arch.arch_model (rtn, mean='AR', lags=5, vol='ARCH',p=5)
res=am.fit()
res.summary()

得到的结果如图所示：

虽然上面演示了 ARCH模型的基本语法。但在实际应用中，我们通常更多地是使用GARCH模型来进行波动率预测。但是两个系数并不显著，因此拟合效果并不好。

GARCH模型
如果说ARCH模型对应着AR模型，那么GARCH模型对应的就是ARMA模型。GARCH模型是ARCH模型的推广形式，是由Bollerslev(1986年）提出来的。
对于收益率序列 $r_{t}$ ，令 $a_{t}=r_{t}-\mu _{t}$ ，若 $a_{t}$ 满足下式，则称 $a_{t}$ 服从GARCH(m,s)模型：

$a_{t}=\delta _{t}\xi _{t} , \delta^{2}_{t}=\alpha _{0}+\sum_{i=1}^{m}\alpha _{i}a^{2}_{t-i}+\sum_{j=1}^{s}\beta _{j}\delta^{2}_{t-j}$

其中， $\xi _{t}$ 是均值为0、方差为1的独立同分布随机变量序列，a0>0,ai≥0,βj≥0，
若s=0,GARCH(m, s）就简化成一个ARCH(m）模型。ai，βj分别称为ARCH参数和GARCH参数。
在实际应用中，GARCH模型的阶不太容易确定，所以往往直接使用低阶的GARCH 模型，比如GARCH(1,1）模型。
下面的示例仍然使用演示ARCH 模型时所用的数据，那么均值方程也仍然可以选择AR(5）模型。拟合模型的示例代码如下：

am=arch.arch_model(rtn,mean='AR',lags=5, vol='GARCH')

res=am.fit()
# 可通过如下命令来查看拟合的结果：
res.summary()

得到的结果如图所示:

效果并不是很好，画图命令如下：

res.plot()

绘制的图形如图所示：

輕栀

关注

0
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫