跳转至




Index

节前迎来揪心一幕!谁来告诉我,A股现在有没有低估?

节前迎来揪心一幕,主要指数均创出今年最低周收盘。很自然,我们也想知道,现在处于什么状态,存在着低估机会吗?这篇文章,我们从市盈利的角度来探讨是存在机会,还是要警惕陷阱。

我们通过akshare来获取沪指导市盈率。实际上,akshare中的这个数据又来自乐咕乐股网站。

1
2
3
4
5
6
7
import akshare as ak

pe = ak.stock_market_pe_lg(symbol="上证")
pe.set_index("日期", inplace=True)
pe.index.name = "date"
pe.rename(columns={"平均市盈率": "pe", "指数": "price"}, inplace=True)
pe.tail(15)

我们将得到从1999年以来的所有数据。它的格式如下:

表1 市盈率与指数

低估的八月?

我们可以通过quantile函数,找出PE的25%, 50%和75%分位数:

1
2
3
4
5
percentiles = []
for i in range(1, 4):
    percentiles.append(pe["pe"].quantile(i/4))

percentiles

我们发现,从1999年以来,25%, 50%和75%分位数分别是 13.9, 17.4和 33.5。

那么,到上个月底(2024年8月),沪指的市盈率处于什么位置呢?

我们可以用以下方法来计算出这个位置:

1
2
3
rank = pe.rank().loc[datetime.date(2024,8,30), "pe"]
percentile = rank / len(pe)
percentile

结果显示,当前的市盈率处在10.6%分位数左右,也就是属于1999年以来较低的位置。

如果仅根据统计数据来看,显然沪指在2024年8月底,是被低估了,换句话说,这里出现了买入机会。

但是,分位数是静态的,它不能反映数据的运行趋势

趋势中的PE

如果你打开东财的客户端,找到沪指的资料页,就可以看到PE的历史走势。近十年的走势图如下:

不过这似乎也给出不了更多信息。从这个图来看,沪指目前也是在低估中。

那么,我们把历年的价格走势叠加到PE走势上,看看PE的峰谷是否真的会对应价格的峰谷。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

fig, ax1 = plt.subplots(figsize=(60,6))

color = "tab:blue"
ax1.plot(pe["price"], label="Index", color=color)
ax1.set_xlabel("Year")
ax1.set_ylabel("Index", color="tab:blue")
ax1.xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2, 5, 8, 11]))
ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) 

color = "tab:red"
ax2 = ax1.twinx()
ax2.plot(pe["pe"], label="PE", color=color)
ax2.set_ylabel("PE", color=color)

for i in range(1, 4):
    quantile = pe["pe"].quantile(i/4)
    ax2.axhline(quantile, color='gray', linestyle='--', label=f"{i/4:02.0%}")

plt.title("Index vs PE")

fig.tight_layout()
for date in pe.index[::3]:
    ax1.axvline(date, color='gray', linestyle=':', linewidth=0.5)
plt.gcf().autofmt_xdate()

plt.legend(loc="upper left")
plt.show()

我们会得到如下所示的绘图:

1999年-2024年

从图中可以看出,从1999年到2024年间,总体上存在沪指上涨,PE下降的趋势。沪指上涨应该是反应了GDP增长的长期趋势;而PE下降,则主要是由于资产供应数量的增加,导致资产泡沫不断被挤出的趋势。

为了能研究细节部分,我们将上图切割成为5个子图:

1999年1月-2004年5月

2004年6月-2009年5月

2009年6月-2014年7月

2014年6月-2019年10月

2019年8月-2024年8月

思考与结论

  • 思考: 在哪些阶段,PE走势与沪指完全一致?在这种走势中,谁是因变量,谁是自变量?它反应什么样的现象?

  • 回答:在1999年到2001年,2005年8月到2008年5月,以及2014年5月到1015年5月等时间段PE的走势与沪指完全一致。在这种走势中,由于市场炒作十分激烈,公司盈利的变化相对于股价的变化而言,完全可以忽略不计,PE的走势完全由股价决定,从而导致PE曲线与股价走势极为相似。一旦出现这种现象,就说明市场已经过度投机了。

  • 思考:大约从2023年8月起,似乎可以看出,沪指的下降速度快于PE的下降速度。比如,2024年2月,沪指出现阶段低点,PE也出现阶段低点。但2024年9月,沪指创新低之后,PE并没有创新低。这种背离反映了什么现象?能否用一个常见的术语来表述它?

  • 回答: 在2024年9月,沪指创新低之后,PE并没有创新低,出现了背离。这种情况其实从峰值分析也可以看出来。比如,2024年6月(参见表我1),PE值创2023年8月以来的新高,但沪指却低于2023年8月后任何一个高点。这种背离可以用低市盈率陷阱来描述。这个指标提醒我们,当前指数虽然在下跌,但上市公司的盈利能力也可能在下降,从而出现指数小涨,PE大涨;指数跌得多,PE跌得少的情况。

因此,如果仅从分位数统计来看,当下的A股是低估的。但如果考虑到市盈率总体上一直在下降的趋势,以及最近一年来PE与指数涨跌的背离情况,判断A股是否已经低估还存有疑问,应该纳入更多维度进行判断。

文章来源于《因子分析与机器学习策略》第2课的习题。

[0915] QuanTide Weekly

本周要闻

  • 主要指数创今年最低周收盘,也是5年最低周收盘
  • 8月有增有降,广义货币增长6.5%,狭义货币下降7.3%。
  • 一意孤行!美提高部分对华301关税 中方:强烈不满 坚决反对
  • 茅台业绩说明会之后,本周白酒指数再跌3.21%

下周看点

  • 周三(晚),美联储公布议息决议。
  • 周五央行公布最新LPR利率。
  • A股将有2只新股发行,创业板1只,北交所1支

本周精选

  • 连载!量化人必会的 Numpy 编程 (3)

本周要闻

  • 本周大盘继续调整,主要指数均创出今年最低周收盘,也是5年最低周收盘,其中,上证指数2700点大关岌岌可危,深证成指跌破8000点大关。本周留下一个9点的跳空缺口,这是今年第一次周K线留下跳空缺口新民晚报
  • 9月13日,央行发布8月社会融资规模统计数据,央行有关部门负责人进行了解读,指出央行认真贯彻落实党中央、国务院决策部署,稳健的货币政策灵活适度、精准有效,强化逆周期调节,为经济社会发展营造良好的货币金融环境。红网.财富频道
  • 9月13日,美国贸易代表办公室宣布将提高部分中国商品的301关税。中方对此强烈不满,坚决反对。自9月27日起,中国制造的电动汽车的关税将上调至100%,太阳能电池上调至50%,电动汽车电池、关键矿产、钢铁、铝、口罩和岸边集装箱起重机将上调至25%,而包括半导体芯片在内的其他产品的关税上调也将在未来两年内生效。东方财富

量化人必会的Numpy编程(3)

1. 处理包含 np.nan 的数据

在量化分析中,我们常常会遇到数据为 np.nan 情况。比如,某公司上年利润为负数,今年利润实现正增长,请问要如何表示公司的 YoY 的利润增长呢?

Info

np.nan 是 numpy 中的一个特殊值,表示“Not a Number”,即“不是数字”。注意,在 Numpy 中,尽管 np.nan 不是一个数字,但它确实数字类型。确切地说,它是 float 类型。此外,在 float 类型中,还存在 np.inf(正无穷大)和负无穷大 (np.NINF,或者-np.inf)。

又比如,在计算个股的 RSI 或者移动平均线时,最初的几期数据是无法计算出来的(在回测框架 backtrader 中,它把这种现象称之为技术指标的冷启动)。如果不要求返回的技术指标的取值与输入数据长度一致,则会返回短一些、但全部由有效数据组成的数组;否则,此时我们常常使用 np.NaN 或者 None 来进行填充,以确保返回的数据长度与输入数据长度一致。

但是,如果我们要对返回的数组进行统计,比如求均值、最大值、排序,对包含 np.nan 或者 None 的数组,应该如何处理?


1.1. 包含 np.nan 和 np.inf 的数组运算

在 numpy 中,提供了对带 np.nan 的数组进行运算的支持。比如有以下数组:

1
2
3
4
import numpy as np

x = np.array([1, 2, 3, np.nan, 4, 5])
print(x.mean())

我们将得到一个 nan。实际上,多数情况下,我们希望忽略掉 nan,只对有效数据进行运算,此时得到的结果,我们往往仍然认为是有意义的。

因此,Numpy 提供了许多能处理包含 nan 数据的数组输入的运算函数。下面是一个完整的列表:

在这里,我们以输入 np.array([1, 2, 3, np.nan, np.inf, 4, 5]) 为例

函数 nan 处理 inf 处理 输出
nanmin 忽略 inf 1.0
nanmax 忽略 inf inf
nanmean 忽略 inf inf
nanmedian 忽略 inf 3.5
nanstd 传递 inf nan
nanvar 传递 inf nan
nansum 忽略 inf inf
nanquantile 忽略 inf 2.25
nancumsum 忽略 inf inf
nancumprod 忽略 inf inf

对 np.nan 的处理中,主要是三类,一类是传递,其结果导致最终结果也是 nan,比如,在计算方差和标准差时;一类是忽略,比如在找最小值时,忽略掉 np.nan,在余下的元素中进行运算;但在计算 cumsum 和 cumprod 时,"忽略"意味着在该元素的位置上,使用前值来填充。我们看一个不包含 np.inf 的示例:

1
2
3
x = np.array([1, 2, 3, np.nan, 4, 5])
np.nancumprod(x)
np.nancumsum(x)
输出结果是:

1
2
3
array([  1.,   2.,   6.,   6.,  24., 120.])

array([ 1.,  3.,  6.,  6., 10., 15.])
结果中的第 4 个元素都是由第 3 个元素复制而来的。

如果一个数组中包含 inf,则在任何涉及到排序的操作(比如 max, median, quantile)中,这些元素总是被置于数组的最右侧;如果是代数运算,则结果会被传导为 inf。这些地方,Numpy 的处理方式与我们的直觉是一致的。

除了上述函数,np.isnan 和 np.isinf 函数,也能处理包含 np.nan/np.inf 元素的数组。它们的作用是判断数组中的元素是否为 nan/inf,返回值是一个 bool 数组。

1.2. 包含 None 的数组运算

在上一节中,我们介绍的函数能够处理包含 np.nan 和 np.inf 的数组。但是,在 Python 中,None 是任何类型的一个特殊值,如果一个数组包含 None 元素,我们常常仍然会期望能对它进行 sum, mean, max 等运算。但是,Numpy 并没有专门为此准备对应的函数。


但是,我们可以通过 astype 将数组转换为 float 类型,在此过程中,所有的 None 元素都转换为 np.nan,然后就可以进行运算了。

1
2
x = np.array([3,4,None,55])
x.astype(np.float64)

输出为:array([3., 4., nan, 55.])

1.3. 性能提升

当我们调用 np.nan *函数时,它的性能会比普通的函数慢很多。因此,如果性能是我们关注的问题,我们可以使用 bottleneck 这个库中的同名函数。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
from bottleneck import nanstd
import numpy as np
import random

x = np.random.normal(size = 1_000_000)
pos = random.sample(np.arange(1_000_000).tolist(), 5)
x[pos] = np.nan

%timeit nanstd(x)
%timeit np.nanstd(x)

我们担心数组中 np.nan 元素个数会影响到性能,所以,在上面的示例中,在随机生成数组时,我们只生成了 5 个元素。在随后的一次测试中,我们把 nan 元素的个数增加了 10 倍。实验证明,nan 元素的个数对性能没有什么影响。在所有的测试中,bottlenect 的性能比 Numpy 都要快一倍。

Info

根据 bottleneck 的文档,它的许多函数,要比 Numpy 中的同名函数快 10 倍左右。


2. 随机数和采样

随机数和采样是量化中的高频使用的操作。在造数据方面非常好用。我们在前面的示例中,已经使用过了 normal() 函数,它是来自 numpy.random 模块下的一个重要函数。借由这个函数,我们就能生成随机波动、但总体上来看又是上涨、下跌或者震荡的价格序列。

Tip

我们会在何时需要造价格序列?除了前面讲过的例子外,这里再举一例:我们想知道夏普为\(S\)的资产,它的价格走势是怎么样的?价格走势与夏普的关系如何?要回答这个问题,我们只能使用“蒙”特卡洛方法,造出若干模拟数据,然后计算其夏普并绘图。此时我们一般选造一个符合正态分布的收益率数组,然后对它进行加权(此时即可算出夏普),最后通过 np.cumprod() 函数计算出价格走势,进行绘图。

我们通过一个例子来说明夏普与股价走势之间的关系:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import numpy as np
from empyrical import sharpe_ratio
import matplotlib.pyplot as plt

returns_ = np.random.normal(0, 0.02, size=100)
legend = []
for alpha in (-0.01, 0, 0.01):
    returns = returns_ + alpha
    prices = np.cumprod(returns + 1)
    sharpe = sharpe_ratio(returns)
    _ = plt.plot(prices)
    legend.append(f"{sharpe:.1f}")

lines = plt.gca().lines
plt.legend(lines, legend)

从绘制的图形可以看出,当 alpha 为 1%时,夏普率可达 8.2。国内优秀的基金经理可以在一年内,做到 2~3 左右的夏普率。大家可以调整 alpha 这个参数,看看 alpha 与夏普率的关系。

2.1. The legacy: np.random module

迄今为止,我们在网上看到的多数关于 numpy random 的教程,都是使用的 np.random module 下面的函数。除了 normal 方法之外,random 包中还有以下函数:


函数 说明
randint(a,b,shape) 生成在区间 (a,b) 之间,形状为 shape 的随机整数数组。
rand(shape) 生成 shape 形状的随机数组,使用 [0,1) 区间的均匀分布来填充。
random(shape) 生成 shape 形状的随机数组,使用均匀分布填充
randn(d1, d2, ...) 生成 shape 形状的随机数组,使用正态分布来填充。
standard_normal(shape) 生成 shape 形状的随机数组,使用标准正态分布来填充。
normal(loc,scale,shape) 生成 shape 形状的随机数组,使用正态分布来填充,loc 是均值,scale 是标准差。
choice(a,size,replace,p) 从 a 中随机抽取 size 个元素,如果 replace=True, 则允许重复抽取,否则不允许重复抽取。p 表示概率,如果 p=None, 则表示每个元素等概率抽取。
shuffle(a) 将 a 中的元素随机打乱。
seed(seed) 设置随机数种子,如果 seed=None, 则表示使用系统时间作为随机数种子。

可以看出,numpy 为使用同一个功能,往往提供了多个方法。我们记忆这些方法,首先是看生成的随机数分布。最朴素的分布往往有最朴素的名字,比如,rand, randint 和 random 都用来生成均匀分布,而 normal, standard_normal 和 randn 用来生成正态分布。

除了均匀分布之外,Numpy 还提供了许多著名的分布的生成函数,比如 f 分布、gama 分布、hypergeometric(超几何分布),beta, weibull 等等。

在同一类别中,numpy 为什么还要提供多个函数呢?有一些是为了方便那些曾经使用其它知名库(比如 matlab) 的人而提供的。randn 就是这样的例子,它是 matlab 中一个生成正态随机分布的函数,现在被 numpy 移植过来了。我们这里看到的另一个函数,rand 也是这样。而对应的 random,则是 Numpy 按自己的 API 风格定义的函数。

choice 方法在量化中有比较具体的应用。比如,我们可能想要从一个大的股票池中,随机抽取 10 只股票先进行一个小的试验,然后根据结果,再考虑抽取更多的股票。


seed 函数用来设置随机数生成器的种子。在进行单元测试,或者进行演示时(这两种情况下,我们都需要始终生成相同的随机数序列)非常有用。

2.2. New Style: default_rng

我们在上一节介绍了一些随机数生成函数,但没有介绍它的原理。Numpy 生成的随机数是伪随机数,它们是使用一个随机数生成器(RNG)来生成的。RNG 的输出是随机的,但是相同的输入总是会生成相同的输出。我们调用的每一个方法,实际上是在这个序列上的一个抽取动作(根据输入的 size/shape)。

在 numpy.random 模块中,存在一个全局的 RNG。在我们调用具体的随机函数时,实际上是通过这个全局的 RNG 来产生随机数的。而这个全局的 RNG,总会有人在它之上调用 seed 方法来初始化。这会产生一些问题,因为你不清楚何时、在何地、以哪个参数被人重置了 seed。

由于这个原因,现在已经不推荐直接使用 numpy.random 模块中的这些方法了。更好的方法是,为每一个具体地应用创建一个独立的 RNG,然后在这个对象上,调用相应的方法:

1
2
rng = np.random.default_rng(seed=123)
rng.random(size=10)

rng 是一个 Random Generator 对象,在初始化时,我们需要给它传入一个种子。如果省略,那么 Numpy 会使用系统时间作为种子。

rng 拥有大多数前一节中提到的方法,比如 normal, f, gamma 等;但从 matlab 中移植过来的方法不再出现在这个对象上。


另外,randint 被 rng.integers 替代。

除此之外,default_rng 产生的随机数生成器对象,在算法上采用了 PCG64 算法,与之前版本采用的算法相比,它不仅能返回统计上更好的随机数,而且速度上也会快 4 倍。

Warning

在 numpy 中还存在一个 RandomState 类。它使用了较慢的梅森扭曲器生成伪随机数。现在,这个类已经过时,不再推荐使用。

2.3. 数据集平衡示例

我们已经介绍了 choice 的功能,现在我们来举一个例子,如何使用 choice 来平衡数据集。

在监督学习中,我们常常遇到数据不平衡的问题,比如,我们希望训练一个分类器,但是训练集的类别分布不均衡。我们可以通过 choice 方法对数据集进行 under sampling 或者 over sampling 来解决这个问题。

为了便于理解,我们先生成一个不平衡的训练数据集。这个数据集共有 3 列,其中前两列是特征(你可以想像成因子特征),第三列则是标签。

1
2
3
4
5
6
import pandas as pd
import numpy as np

rng = np.random.default_rng(seed=42)
x = rng.random((10,3))
x[:,-1] = rng.choice([0,1], len(x), p=[0.2, 0.8])

我们通过下面的方法对这个数据集进行可视化,以验证它确实是一个不平衡的数据集。

1
2
df = pd.DataFrame(x, columns=['factor1', 'factor2', 'label'])
df.label.value_counts().plot(kind='bar')

运行结果为:

要在此基础上,得到一个新的平衡数据集,我们有两种思路,一种是 under sampling,即从多数类的数据中抽取部分数据,使得它与最小分类的数目相等;另一种是 over sampling,即从少数类的数据中复制部分数据,使得它与最大的类的数目相等。

下面的例子演示了如何进行 under sampling:


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
labels, counts = np.unique(x[:,-1], return_counts=True)

# 最小分类的标签
min_label = labels[np.argmin(counts)]

# 最小分类样本的数量,作为 UNDER SAMPLING 的样本数量
min_label_count = np.min(counts)

# 最小分类无须抽取,全部提取
results = [
    x[x[:,-1] == min_label]
]

# 对其它分类标签进行遍历,需要先剔除最小分类
for label in np.delete(labels, np.argmin(counts)):
    sampled = rng.choice(x[x[:,-1]== label], min_label_count)
    results.append(sampled)

np.concatenate(results)

这段代码先是找到最小分类及它的数量,然后遍历每个标签,再通过 rng.choice 对其它分类随机抽取最小分类的数量,最后把所有的子集拼接起来。

这段示例代码可用以多个标签的情况。如果要进行 over sampling,只要把其中的 min 换成 max 就可以了。

3. IO 操作

我们直接使用 Numpy 读写文件的场合并不多。提高 IO 读写性能一直都不是 Numpy 的重点,我们也只需要稍加了解即可。

3.1. 读写 CSV 文件

Numpy 可以从 CSV 格式的文本文件中读取数据,主要有以下方法:


api 描述
loadtxt 解析文本格式的表格数据
savetxt 将数据保存为文本文件
genfromtxt 同上,但允许数据中有缺失值,提供了更高级的用法
recfromtxt 是 genfromtxt 的快捷方式,自动推断为 record array
recfromcsv 同上,如果分隔符为逗号,无须额外指定

我们通过下面的示例简单演示一下各自的用法:

1
2
3
4
5
6
7
import io
import numpy 

buffer = io.StringIO("""1,2""")

# 默认情况下,LOADTXT 只能读取浮点数
numpy.loadtxt(buffer, delimiter=",")

这会输出数组array([1., 2.])

1
2
3
4
buffer = io.StringIO("""1,2,hello""")

# 通过指定 DTYPE 参数,可以读取其它类型
numpy.loadtxt(buffer, delimiter=",", dtype=[("age", "i4"), ("score", "f4"), ("name", "U8")])

这样我们将得到一个 Structured Array,其中第三列为字符串类型。如果我们不指定 dtype 参数,那么 loadtxt 将会解析失败。

1
2
3
4
buffer = io.StringIO("""
1,2,hello
""")
numpy.genfromtxt(buffer, delimiter=",")

这一次我们使用了 genfromtxt 来加载数据,但没有指定 dtype 参数,genfromtxt 会将非数字列解析为 nan。因此,这段代码将输出:`array([1., 2., nan])

现在,我们也给 genfromtxt 加上 dtype 参数:


1
2
3
4
5
buffer = io.StringIO("""
1,2,hello
""")

numpy.genfromtxt(buffer, delimiter=",", dtype=[("age", "i4"), ("score", "f4"), ("name", "U8")])

此时我们得到的结果是:array((1, 2., 'hello'), dtype=[('age', '<i4'), ('score', '<f4'), ('name', '<U8')])。注意它是 Structured Array。

recfromtxt 则不需要 dtype, 会自动推断数据类型。

1
2
3
4
5
buffer = io.StringIO("""
1,2,hello
""")

numpy.recfromtxt(buffer,delimiter=",")

这段代码输出为rec.array((1, 2, b'hello'), dtype=[('f0', '<i8'), ('f1', '<i8'), ('f2', 'S5')])。如果推断不准确,我们也可以自己加上 dtype 参数。

如果我们使用 recfromcsv,则连 delimiter 参数都可以省掉。

1
2
3
4
5
buffer = io.StringIO("""
age,score,name
1,2,hello
""")
numpy.recfromcsv(buffer)

输出跟上一例的结果一样。

出于速度考虑,我们还可以使用其它库来解析 CSV 文件,再转换成为 numpy 数组。比如:


1
2
3
4
5
# 利用 CSV.READER() 来解析,比 NUMPY 快 8 倍
np.asarray(list(csv.reader()))

# 利用 PANDAS 来解析,比 NUMPY 快 22 倍
pd.read_csv(buffer).to_records()

3.2. 读写二进制文件

如果我们不需要与外界交换数据,数据都是自产自销型的,也可以使用二进制文件来保存数据。

使用 numpy.save 函数来将单个数组保存数据为二进制文件,使用 numpy.load 函数来读取 numpy.save 保存的数据。这样保存的文件,文件扩展名为.npy。

如果要保存多个数组,则可以使用 savez 命令。这样保存的文件,文件扩展名为.npz。

如果有更复杂的需求,可以使用 Hdf5,pyarrow 等库来进行保存数据。

4. 日期和时间

一些第三方数据源传递给我们的行情数据,常常会用字符串形式,或者整数(从 unix epoch time 起)格式来表示行情的时间。比如,akshare 和 tushare 许多接口给出的行情数据就是字符串格式;而 QMT 很多时候,会将行情时间用整数表示。掌握这些格式与 Numpy 的日期时间格式转换、以及 Numpy 到 Python 对象的时间日期转换是非常有必要的。

但是在任何编程语言中,日期和时间的处理从来都不简单。


Info

很少有程序员/研究员了解这一点:日期和时间并不是一个数学上或者物理上的一个客观概念。时区的划分、夏令时本身就是一个政治和法律上的概念;一些地方曾经使用过夏令时,后来又取消了这种做法。其次,关于闰秒 1 的决定,也并不是有章可循的,它是由一个委员会开会来临时决定的。这种决定每年做一次。

所有这些决定了我们无法通过一个简单的数学公式来计算时间及其变化,特别是在时区之间的转换。

关于时间,首先我们要了解有所谓的 timezone aware 时间和 timezone naive 时间。当我们说到晚上 8 时开会时,这个时间实际上默认地包含了时区的概念。如果这是一个跨国会议,但你在通知时不告诉与会方时区,就会导致其它人无法准时出席 -- 他们将会在自己时区的晚上 8 时上线。

如果一个时间对象不包含时区,它就是 timezone naive 的;否则,它是 timezone aware 的。但这只是针对时间对象(比如,Python 中的 datetime.datetime)才有的意义;日期对象(比如,Python 中的 datetime.date)是没有时区的。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import pytz
import datetime

# 通过 DATETIME.NOW() 获得的时间没有时区信息
# 返回的是标准时间,即 UTC 时间,等同于调用 UTCNOW()
now = datetime.datetime.now()
print(f"now() without param: {now}, 时区信息{now.tzinfo}")

now = datetime.datetime.utcnow()
print(f"utcnow: {now}, 时区信息{now.tzinfo}")

# 构造 TIMEZONE 对象
cn_tz = pytz.timezone('Asia/Shanghai')
now = datetime.datetime.now(cn_tz)
print(f"现在时间{now}, 时区信息{now.tzinfo}")

1
2
3
4
5
6
print("现在日期:", now.date())

try:
    print(now.date().tzinfo)
except AttributeError:
    print("日期对象没有时区信息")

上述代码将依次输出:

1
2
3
4
5
now() 不带参数:2024-05-19 11:03:41.550328, 时区信息 None
utcnow: 2024-05-19 11:03:41.550595, 时区信息 None
现在时间 2024-05-19 19:03:41.550865+08:00, 时区信息 Asia/Shanghai
现在日期:2024-05-19
日期对象没有时区信息

不过,限于篇幅,我们对时间问题的介绍只能浅尝辄止。在这里,我们主要关注在 Numpy 中,日期/时间如何表示,它们彼此之间如何比较、转换,以及如何与 Python 对象进行比较和转换。

在 Numpy 中,日期/时间总是用一个 64 位整数(np.datetime64)来表示,此外,还关联了一个表示其单位(比如,纳秒、秒等)的元数据结构。np.datetime64是没有时区概念的。

1
2
3
tm = np.datetime64('1970-01-01T00:00:00')
print(tm)
print(tm.dtype)

这将显示为:

1
2
1970-01-01T00:00:00
datetime64[s]

这里的[s]就是我们前面所说的时间单位。其它常见单位还有[ms][us][ns]等等。

除了从字符串解释之外,我们还可以直接将 Python 对象转换成np.datetime64,反之亦然:


1
2
3
4
5
tm = np.datetimet64(datetime.datetime.now())
print(tm)

print(tm.item())
print(tm.astype(datetime.datetime))

下面我们来看看如何实现不同格式之间的批量转换。这在处理 akshare, tushare 或者 QMT 等第三方数据源提供的行情数据时,非常常见。

首先我们构造一个时间数组。顺便提一句,这里我们将使用np.timedelta64这个时间差分类型:

1
2
3
now = np.datetime64(datetime.datetime.now())
arr = np.array([now + np.timedelta64(i, 'm') for i in range(3)])
arr

输出结果如下:

1
2
3
4
array(['2024-05-19T12:57:47.349178', 
       '2024-05-19T12:58:47.349178',
       '2024-05-19T12:59:47.349178'], 
     dtype='datetime64[us]')

我们可以通过np.datetime64.astype()方法将时间数组转换为 Python 的时间对象:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
time_arr = arr.astype(datetime.datetime)

# 转换后的数组,每个元素都是 TIMEZONE NAIVE 的 DATETIME 对象
print(type(time_arr[0]))

# !!! 技巧
# 如何把 NP.DATETIME64 数组转换为 PYTHON DATETIME.DATE 数组?
date_arr = arr.astype('datetime64[D]').astype(datetime.date)
# 或者 -- 两者的容器不一样
date_arr = arr.astype('datetime64[D]').tolist()
print(type(date_arr[0]))

这里的关键是,我们之前生成的arr数组,其元素类型为np.datetime64[us]。它到 Python datetime.date的转换将损失精度,所以 Numpy 要求我们显式地指定转换类型。

如何将以字符串表示的时间数组转换为 Numpy datetime64 对象数组呢?答案仍然是 astype() 方法。

1
2
3
4
5
6
# 将时间数组转换为字符串数组
str_arr_time = arr_time.astype(str)
print(str_arr_time)

# 再将字符串数组转换为 DATETIME64 数组,精度指定为 D
str_arr_time.astype('datetime64[D]')

显示结果为:

1
2
3
4
5
6
7
8
9
array(['2024-05-19T12:57:47.349178', 
       '2024-05-19T12:58:47.349178',
       '2024-05-19T12:59:47.349178'], 
       dtype='datetime64[us]')

array([
    '2024-05-19', 
    '2024-05-19'],               
    dtype='datetime64[D]')

最后,我们给一个 QMT 获取交易日历后的格式转换示例。在 QMT 中,我们通过get_trading_dates来获取交易日历,该函数返回的是一个整数数组,每个元素的数值,是从 unix epoch 以来的毫秒数。

我们可以通过以下方法对其进行转换:

1
2
3
4
import numpy as np

days = get_trading_dates('SH', start_time='', end_time='', count=10)
np.array(days, dtype='datetime64[ms]').astype(datetime.date)

QMT 官方没有直接给出交易日历转换方案,但给出了如何将 unix epoch 时间戳转换为 Python 时间对象(但仍以字符串表示):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import time
def conv_time(ct):
    '''
    conv_time(1476374400000) --> '20161014000000.000'
    '''
    local_time = time.localtime(ct / 1000)
    data_head = time.strftime('%Y%m%d%H%M%S', local_time)
    data_secs = (ct - int(ct)) * 1000
    time_stamp = '%s.%03d' % (data_head, data_secs)
    return time_stamp

conv_time(1693152000000)

我们需要对每一个数组元素使用上述解析方法。官方方案的优点是不依赖任何三方库。不过,没有量化程序能离开 Numpy 库,所以,我们的方案并未增加第三方库的依赖。

字符串操作

你的数据源、或者本地存储方案很可能使用 Numpy Structured Array 或者 Rec Array 返回证券列表。很显然,证券列表中一定会包括字符串,因为它一定会存在证券代码列和证券名称列。有一些还会返回证券的地域属性和其它属性,这也往往是字符串。

对证券列表,我们常常有以下查询操作:

  1. 获取在某个板块上市的股票列表,比如,北交所、科创板和创业板与主板的个股交易规则上有一些不同,因此,我们的策略很可能需要单独为这些板块构建。这就有了按板块过滤证券列表的需要。也可能我们要排除 ST,刚上市新股。这些都可以通过字符串操作来实现。

  1. 市场上有时候会出现魔幻的名字炒作。比如龙年炒龙字头(或者含龙的个股)、炒作“东方”、炒作“中”字头。作为量化人,参与这样的炒作固然不可取,但我们要拥有分析市场、看懂市场的能力。

Numpy 中的大多数字符串操作都封装在 numpy.char 这个包下面。它主要提供了一些用于格式化的操作(比如左右填充对齐、大小写转换等)、查找和替换操作。

下面的代码展示了如何从证券列表中过滤创业板:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import numpy as np
import numpy.char as nc

# 生成 STRUCTURED ARRAY, 字段有 SYMBOL, NAME, IPO DATE
arr = np.array([('600000.SH', '中国平安', '1997-08-19'),
                ('000001.SZ', '平安银行', '1997-08-19'),
                ('301301.SZ', '川宁生物', '2012-01-01')
                ], dtype=[('symbol', 'S10'), ('name', 'S10'), ('ipo_date', 'datetime64[D]')])

def get_cyb(arr):
    mask = np.char.startswith(arr["symbol"], b"30")
    return arr[mask]

Question

我们在查找创业板股票时,使用的是 b"30"来进行匹配。为何要用 b"30"而不是"30"?

注意第 11 行,我们要通过np.char.startswith()来使用startswith函数。任何一个 numpy array 对象都没有这个方法。

".SZ"是我们的数据源给股票编制的交易所代码。不同的数据源,可能使用不同的交易所代码。比如,聚宽数据源会使用.XSHG 表示上交所,.XSHE 表示深交所。现在,如果我们要将上述代码转换为聚宽的格式,应该如何操作?


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# 生成 STRUCTURED ARRAY, 字段有 SYMBOL, NAME, IPO DATE
arr = np.array([('600000.SH', '中国平安', '1997-08-19'),
                ('000001.SZ', '平安银行', '1997-08-19'),
                ('301301.SZ', '川宁生物', '2012-01-01')
                ], dtype=[('symbol', 'U10'), ('name', 'U10'), ('ipo_date', 'datetime64[D]')])

def translate_exchange_code(arr):
    symbols = np.char.replace(arr["symbol"], ".SH", ".XSHG")
    print(symbols)
    symbols = np.char.replace(symbols, ".SZ", ".XSHE")

    arr["symbol"] = symbols
    return arr

translate_exchange_code(arr)

这一次,我们把 symbol 和 name 的定义改为 Unicode 型,以避免我们查找时,要输入像 b"30"这样的字面量。

但输出的结果可能让人意外,因为我们将得到这样的输出:

1
2
3
4
array([('600000.XSH', '中国平安', '1997-08-19'),
       ('000001.XSH', '平安银行', '1997-08-19'),
       ('301301.XSH', '川宁生物', '2012-01-01')],
      dtype=[('symbol', '<U10'), ('name', '<U10'), ('ipo_date', '<M8[D]')])

Question

发生了什么?我们得到了一堆以".XSH"结尾的 symbol,它们本应该是"600000.XSHG"这样的字符串。错在哪里,又该如何修改?

在上面的示例中,如果我们把替换字符串改为空字符串,就实现了删除操作。这里就不演示了。

char 模块还提供了字符串相等比较函数equal:

1
2
3
4
arr = array([('301301.SZ', '川宁生物', '2012-01-01')],
      dtype=[('symbol', '<U10'), ('name', '<U10'), ('ipo_date', '<M8[D]')])

arr[np.char.equal(arr["symbol"], "301301.SZ")]

在这个特殊的场景下,我们也可以直接使用以下语法:

1
arr[arr["symbol"] == "301301.SZ"]

Tip

np.char 下的函数很多,如何记忆?实际上,这些函数多数是 Python 中 str 的方法。如果你熟悉 Pandas,就会发现 Pandas 中也有同样的用法。因此,像upper, lower, strip这样的str函数,你可以直接拿过来用。

Numpy 中的字符串函数另一个比较常用的场景,就是执行格式化。你可以通过ljust, 'center', rjust在显示一个数组前,将它们的各列数据进行左右空格填充,这样,输出时就可以比较整齐。

Question

2024 年 5 月 10 日起,南京化纤走出 7 连板行情,短短 7 日,股价翻倍。市场上还有哪些名字中包含化纤的个股?它们的涨跌是否存在相关性或者跨周期相关性?

Masked Array

你可能常常在一些接近底层的库中,看到 Numpy masked array 的用法。Masked Array 是 Numpy 中很重要的概念。考虑这样的情景,你有一个数据集,其中包含了一些缺失的数据或者无效值。这些”不合格“的数据,可能以 np.nan,np.inf, None 或者其它仅仅是语法上有效的值来表示(比如,在 COVID-19 数据集中,病例数出现负数)的。


如何在保持数据集的完整性不变的前提下,仍然能对数据进行运算呢?

Note

这里有一个真实的例子。你可以在 Kaggle 上找到一个 COVID-19 的数据集,这个数据集中,就包含了累积病例数为负数的情况。该数据集由 Johns Hoopkins University 收集并提供。

很显然,我们无法直接对这些数据进行运算。请看下面的例子:

1
2
3
x = np.array([1, 2, 3, np.inf, np.nan, None])
np.mean(x)
np.nanmean(x)

只要数据中包含 np.nan, np.inf 或者 None,numpy 的函数就无法处理它们。即使数据在语法上合法,但在语义上无效,Numpy 强行进行计算,结果也是错误的。

这里有一个量化中可能遇到的真实场景,某公司有一年的年利润为零,这样使得它的 YoY 利润增长在次年变得无法计算。如果我们需要利用 YoY 数据进一步进行运算,我们就需要屏蔽这一年的无效值。否则,我们会连 YoY 利润的均值都无法计算出来。

这里有一个补救的方法,就是将原数据拷贝一份,并且将无效值替换为 np.nan。此后多数运算,都可以用np.nan*来计算。这个方法我们已经介绍过了。但是,如果你是原始数据的收集者,显然你应该照原样发布数据,任何修改都是不合适的;如果你是数据的应用者,当然应该对数据进行预处理后,才开始运算。但是,你又很可能缺少了对数据进行预处理所必须的信息 -- 你怎么能想到像-1, 0 这样看起来人畜无害的小可爱们,竟然是隐藏着的错误呢?


为了解决这个问题,Numpy 就提供了 Masked Array。但是我们不打算在 Masked Array 上过多着墨。关于 Masked Array,我们可以借用这样一句话,很多人不需要知道 Masked Array,知道 Masked Array 的人则已经精通它了。

有一点需要注意的是,仅在需要时,使用 Masked Array。因为可能与您相像的相反,Masked Array 不会提高性能,反而,它大大降低了性能:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import numpy as np

# NUMPY VERSION 1.24.4
g = np.random.random((5000,5000))
indx = np.random.randint(0,4999,(500,2))
g_nan = g.copy()
g_nan[indx] = np.nan
mask =  np.full((5000,5000),False,dtype=bool)
mask[indx] = True
g_mask = np.ma.array(g,mask=mask)

%timeit (g_mask + g_mask)**2
# 901 MS ± 52.3 MS PER LOOP ...
%timeit (g_nan + g_nan)**2
# 109 MS ± 72.2 ΜS PER LOOP ...

可以看出,Masked Array 的性能慢了接近 9 倍。

Tip

如果你不得不对含有 np.nan 的数组进行运算,那么可以尝试使用 bottleneck 库中的 nan *函数。由于并不存在 nansquare 函数,但是考虑到求方差的运算中必然包含乘方运算,因此我们可以考虑通过 nanvar 函数来评测 numpy 与 bottleneck 的性能差异。

1
2
3
4
5
6
    %timeit np.var(g_mask)
    # 587 MS ± 37.9 MS PER LOOP ...
    %timeit np.nanvar(g_nan)
    # 281 MS ± 1.46 MS PER ...
    %timeit nanvar(g_nan)
    # 61 MS ± 362 ΜS PER LOOP ...

1
bottleneck 要比 numpy 快接近 5 倍。如果你使用的 numpy 版本较旧,那么 bottleneck 还会快得更多。

ufunc

ufunc 是 Numpy 中的重要概念,它对两个输入数组同时进行逐元素的操作(比如,相加,比较大小等)。在 Numpy 中大约定义了 61 个左右的 ufunc。这些操作都是由底层的 C 语言实现的,并且支持向量化,因此,它们往往具有更快的速度。

比如,在 numpy 中,求数组中的最大值,有两个相似的函数, np.maxnp.maximum可以达成这一目标。后者是 ufunc,前者不是,两者除了用法上有所区别之外,后者的速度也要快一些。

1
2
3
4
5
6
arr = np.random.normal(size=(1_000_000,))

%timeit np.max(arr)
# 801 MS ± 54.7 MS PER LOOP ...
%timeit np.maximum.reduce(arr)
# 775 MS ± 12.1 MS PER LOOP ...

np.maximum作为 ufunc,它本来是要接收两个参数的,并不能用来求一维数组的最大值。这种情况下,我们要使用reduce操作才能得到想要的结果。

这里np.maximum是一个 ufunc,则reduce是 unfunc 对象(在 Python 中,一切都是对象,包括函数)的属性之一。ufunc的其它属性还有accumulateouterreduceat等。

accumulate是 ufunc 中的另一个常用属性,可能你之前已经有所接触。比如,在求最大回撤时,我们就会用到它:


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# 模拟一个股价序列
n = 1000
xs = np.random.randn(n).cumsum()

# 最大回撤结束期
i = np.argmax(np.maximum.accumulate(xs) - xs) 

# 最大回撤开始期
j = np.argmax(xs[:i]) 

# 最大回撤
mdd = (xs[j] - xs[i])/xs[j]

plt.plot(xs)
plt.plot([i, j], [xs[i], xs[j]], 'o', color='Red', markersize=10)

简洁就是美。在使用了accumulate之后,我们发现,计算最大回撤竟简单到只有两三行代码即可实现。

ufunc 如此好用,你可能要问,为何我用到的却不多?实际上,你很可能每天都在使用ufunc。许多二元数学操作,它们都是对 ufunc 的封装。


比如,当我们调用A + B时,实际上是调用了np.add(A, B)这个 ufunc。二者在功能和性能上都是等价的。其它的 ufunc 还有逻辑运算、比较运算等。只要某种运算接受两个数组作为参数,那么,很可能 Numpy 就已经实现了相应的 ufunc 操作。此外,一些三角函数,尽管只接受一个数组参数,但它们也是 ufunc。

因此,我们需要特别关注和学习的 ufunc 函数,可能主要就是maximumminimum等。这里再举一个在量化场景下,使用maximum的常用例子 -- 求上影线长度。

Tip

长上影线是资产向上攻击失败后留下的痕迹。它对股价后来的走势分析有一定帮助。首先,资金在这个点位发起过攻击,暴露了资金的意图。其次,攻击失败,接下来往往会有洗盘(或者溃败)。股价底部的长上影线,也被有经验的股民称为仙人指路。后面出现拉升的概率较大。上影线出现在高位时,则很可能是见顶信号。此时在较低级别的 k 线上,很可能已经出现均线拐头等比较明显的见顶信号。

现在,我们就来实现长上影线的检测。上影线的定义是:

\[ upper\_shadow = high - max(open, close) \]

下图显示了上影线:


如果 \(\(upper\_shadow > threshold\)\),则可认为出现了长上影线(当然,需要对 upper_shadow 进行归一化)。检测单日的上影线很简单,我们下面的代码将演示如何向量化地求解:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import numpy as np
import pandas as pd

rng = np.random.default_rng(seed=78)
matrix = rng.uniform(0.98, 1.02, (4, 30)).cumprod(axis=1)
opn = matrix[0]
close = matrix[-1]
high = np.max(matrix, axis=0)

upper_shadow = (high - np.maximum(opn, close))/close
np.round(upper_shadow, 2)

第 10 行的代码完全由 ufunc 组成。这里我们使用了 np.sub(减法), np.maximum, np.divide(除法)。maximum 从两个等长的数组 opn 和 close 中,逐元素比较并取出最大的那一个,组成一个新的数组,该数组也与 opn, close 等长。

如果要求下影线长度,则可以使用 minimum。


《因子投资与机器学习策略》开课啦!


目标清晰 获得感强


为什么你值得QuanTide的课程?



  1. https://zh.wikipedia.org/wiki/%E9%97%B0%E7%A7%92 

牛人太多:小市值因子之父,毕业论文被大佬狂怼

Rolf W. Banz,瑞士人,70 年代在芝加哥大学获得博士学位,并在该校任教过。此后他在伦敦经营了一家专注于小盘股投资的投资精品店,1991 年出售给 Alliance Capital。职业生涯的最后阶段,他回到了瑞士,在一家瑞士私人银行的资产管理子公司担任了高级职位。

退休后,他是一名博主。最后一篇博文,停留在 2017 年底。

小市值因子之父

Banz 因发现小市值因子一战成名。1981 年,在他的博士毕业论文中写道,基于纽交所长达 40 年的数据,小市值股票月均收益率比其他股票高 0.4%。这个发现深刻地影响了此后的市场。在后 Banz 时代,基金经理超配小盘股一度成了一种信仰。

在国内,研究表明,A 股市场在 2016 年以前,规模因子的显著性甚至超过了欧美等发达国家市场。该因子在 2017-2018 年期间,弱于大市值资产,但从 2021 年起,又跑赢了沪深 300。

这张图来自聚宽平台的回测,数据表明,两年时间里,小市值获得到 51.98%的超额收益,是相当不错的。

Rolf Banz 也因此获得了一个非正式的 title,小市值因子之父。

这个 title,不仅仅是因为小市值因子。如果说威廉. 夏普定义了第一个因子的话,那么第二个因子,小市值因子,则是 Banz 发现的。他完成了一生二,甚至是零生一的这一步。

不被大佬接受的论文

按理说,搞出了这么牛 x 的因子,作为论文答辩应该很顺利,对吧?但是,当 Banz 进行答辩时,意外发生了。

审查论文的都是大牛人,有迈伦·斯科尔斯、默顿·米勒和尤金·法玛,他们后来都获得了诺贝尔奖。但是,没有一个人喜欢 Banz 的发现。尤其是尤金·法玛(他是有效市场理论的提出者)。为什么?

有效市场的一个推论就是,只应该有一个因子,就是市场因子。因为,既然市场是有效的,就不应该存在信息差和错误定价。没有错误定价,就没有市场异象和新因子。

所以,当 Banz 在台上介绍自己的成果时,默顿·米勒就在台下问斯科尔斯,他的错误是什么?一众评委都天然地倾向否认 Banz 的结论。他们出席评审的目的,就是为了证明他是错的!

Banz 在学术成就上,是不如斯科尔斯,法玛他们这些人的,他也是一名著名的金融学家,但没能获得重要级别的大奖。但他提出的小市值因子如果用以投资,会非常有效。

Banz 自己的投资相对于其它金融学家来说比较成功。与之对比的是,斯科尔斯曾经掌管了一个名为长期资本管理的对冲基金,这个基金 1998 年倒闭,损失了数十亿美元。

简单即是美

Banz 的这篇论文并不复杂,只有 16 页纸。论文中并没有使用高深的数学,只是使用了基础的统计科学,GLS 和 OLS。如果当时有计算机可用的话,这个推导过程会显得格外简单。

Quote

if you can’t explain it simply, you don’t understand it well enough. - Einstein

Banz 发表论文时,他研究的数据来自 1926 年到 1975 年。这些数据就在芝大的 CRSP 数据库里。在漫长的 50 多年时间里,一直在静静地等待有人发现她的价值。芝大经济学人才辈出,期间不知道有多少人接触过这份数据。如此低垂和芳香的果实最终却被 Banz 采摘。

我想,这足以说明,在任何时代,在任何圈子,都会有一些低垂的果实,等待有心人去采摘

Banz 在论文中提到,他无法介绍小市值效应产生的原因。

但实际上这是容易解释的。从生物学的角度来看,幼小、成长性和风险总是相伴相随;从经济学的角度来看,一些很大的公司,往往已经做到了行业的天花板,做到了很高的市占率,如果他们也想市值 double 的话,就必须开拓一个新的赛道,并且这个赛道还得跟之前的一样宽。

但是,对大公司来说,创新和死亡都很困难。

如果在 Banz 之前的学者,能跳出计量经济学的窠臼看等市场,同时又具有扎实的统计和数学功底,发现这个因子的时间应该能早很多。

因为造原子弹的最大秘密,就是原子弹是可造的。

后 Banz 时代的 Banz

尽管小市值因子的论文使用的统计学知识非常基础,但扎实二字,其实比我们想像的要难。就这样一篇论文,尽管在形式上无懈可击,但仍然受到大量的质疑。

1999 年,Tyler Shumway and Vincent Warther 指出,Banz 的研究方法存在显著错误,他忘了减去退市股的负面效应。他们的文章发表为《 The Delisting Bias in CRSP's Nasdaq Data and Its Implications for the Size Effect》。

2011 年,Gary A. Miller and Scott A. MacKillop 认为,小市值并不优于大型股。他们使用标普 500 指数来代替大型股,使用 CRSP9-10 指数来代表小型股,两者都可以追溯到 1926 年。根据他们的研究,CRSP 9-10 和 DFA 基金在 1980 到 2000 年的回报率低于标准普尔 500 指数,虽然差距不大。

Note

实际上这个结论很容易误导人。它实际上证明了小市值因子的有效性。因为标普 500 指数的表现,远远超过了大多数基金经理。另外,根据 Miller & MacKillop,小市值因子只是在夏普上低于大型股 -- 这也是 Banz 在论文中就已经指出的事实。

Banz 对这些质疑没有进行回应。尽管他在 2012 年就开始撰写 个人博客,但就像完全没看到这些质疑一样。

[0908] QuanTide Weekly

本周要闻

  • 央行:降准有空间 利率进一步下行面临一定约束
  • 巴菲特再次减持美国银行,这次要做空自己的祖国?
  • 存量房贷下调预期落空,沪指连续三日跌破2800点

下周看点

  • 周一,8月CPI/PPI 数据发布
  • 周一,茅台业绩说明会,对白酒行业的未来发展预判是重要信号

本周精选

  • 那些必读的量化论文
  • 连载!量化人必会的 Numpy 编程 (2)

本周要闻

  • 央行货币政策司邹澜9月5日表示,是否降准降息还需观察经济走势。年初降准的政策效果还在持续显现,目前金融机构的平均法定存款准备金率大约为7%,还有一定的空间。同时,受银行存款向资管产品分流的速度、银行净息差收窄的幅度等因素影响,存贷款利率进一步下行面临一定的约束
  • 当地时间9月5日,美国证券交易委员会披露,巴菲特于9月4日前后连续三日减持美国银行套现约7.6亿美元。今年7月,他曾连续9天减持美国银行,在减持前,美国银行为其第二大重仓股,也是伯克希尔最赚钱的公司之一。
  • 沪指本周下跌2.69%,连续三日跌破2800点。目前沪指PE值处在12%分位左右,属较低位置。消息面上,上周盛传的存量房贷降息预期落空。
  • 证监会离职人员入股拟上市企业新规落地,要求更严、核查范围更广

消息来源:东方财富、财联社


那些必读的量化论文

  1. Portfolio Selection, Markovitz, 1952。在这篇文章里,马科维茨提出了现代资产组合理论(MPT),使他获得了 1990 年诺贝尔经济学奖。该论文通过将风险与回报之间的相关性纳入计算中来扩展了常见的风险回报权衡。
  2. A New Interpretation of Information Rate,Kelly, 1956。本文是当今著名的凯利准则的正式表述。该模型广泛应用于赌场游戏,特别是风险管理。作者推导出了一个公式,该公式确定了最佳分配规模,以便随着时间的推移最大化财富增长。
  3. Capital Asset Prices: A Theory of Market Equilibrium under Conditions of Risk (Sharpe, 1964)。基于Markovitz的工作,CAPM证明了只有一种有效的投资组合,即市场投资组合。这篇文章提出了著名的Beta概念。Sharpe和Markovitz是师生关系,同年获得诺贝尔经济学奖。
  4. Efficient Capital Markets: a Review of Theory and Empirical Work,Fama, 1970。这篇论文是首次提出非常流行的“有效市场假说”概念的开创性著作,尽管现在这一理论受到很大的质疑,但是它的学术价值仍然很高。
  5. The Pricing of Options and Corporate Liabilities, 1973, Black & Scholes。著名的BS公式,利用物理传热方程作为估算期权价格的起点。这也是为什么对冲基金喜欢物理系学生的原因。

  1. Does the Stock Market Overreact? 1985, Bondt & Thaler。这篇文章质疑了有效市场假说,Bondt和Thaler提出,有统计上显着的证据表明相反的情况,投资者往往会对意外的新闻事件反应过度。这也是行为金融学的经典研究之一。行为金融学在近年来是诺奖的大热门。它的底层哲学是主观价值论和与人为本的思想。Thaler也是诺奖获得者,并在电影《大空头》中扮演了自己。
  2. A closed-form GARCH option valuation model,1997, Heston & Nandi。本文提出了一种封闭式公式,用于评估现货资产,并使用广义自回归条件异方差(GARCH) 模型对其方差进行建模。由于其复杂性和实用性, GARCH 模型在 20 世纪 90 年代估计波动率方面广受欢迎,金融业也积极采用它们。
  3. Optimal Execution of Portfolio Transactions,2000,Almgren & Chriss。对于每个负责完善交易执行算法的量化开发人员来说,这篇论文绝对是必读的学术文献。文章指出,价格波动来自外生性(市场本身的波动)和内生性(自己的订单对市场造成的影响)。这是一种量子效应!作者通过最小化交易成本和波动风险的组合,形式化了一种执行和衡量交易执行绩效的方法。
  4. Incorporating Signals into Optimal Trading,2017,Lehalle。与 Almgren 和 Chris (2000) 的工作非常相似,本文讨论的是最佳交易执行。作者通过将马尔可夫信号纳入最优交易框架,进一步完善了该领域所做的工作,并针对带有漂移的随机过程(Ornstein-Uhlenbeck 过程)的资产特殊情况得出了最优交易策略。
  5. The Performance of Mutual Funds in the Period 1945-1964, Michael Jessen。Sharpe在他的文章里提出了Beta这个概念,而Alpha这个概念,就是由Jessen在这篇论文中提出的。

  1. Common risk factors in the returns on stocks and bonds, 1993, Eugene F. Fama。在这篇文章里,Famma提出了三因子。
  2. Review of Financial Studies,2017, Stambaugh, Yuan。这篇论文出现的比较晚,因此可以把之前出现的、与因子投资相关的重要论文都评述一篇,因而就成为了我们快速了解行业的一篇重要论文。这篇文章发表比较晚,但也有了952次引用。
  3. 151 Trading Strategies, 2018, Zura Kakushadze。作者来自world quant,是Alpha101的作者之一。这篇论文引用了大量的论文(2000+),也是我们进行泛读的好材料。

量化人必会的 NUMPY 编程 (2) - 核心语法

1. Structured Array

一开始,Numpy的数组只能存放同质的元素,即这些元素都必须有同样的数据类型。但对很多表格类数据,它们往往是由一条条记录组成的,而这些记录,又是由不同数据类型的数据组成的。比如,以最常见的行情数据来讲,它就必须至少包含时间、证券代码和OHLC等数据。

为了满足这种需求,Numpy扩展出一种名为Structured Array的数据格式。它是一种一维数组,每一个元素都是一个命名元组。

我们可以这样声明一个Structured Array:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
import numpy as np
import datetime
dtypes = [
        ("frame", "O"),
        ("code", "O"),
        ("open", "f4"),
        ("high", "f4"),
        ("low", "f4"),
        ("close", "f4")
    ]
secs = np.array(
    [
        (datetime.date(2024, 3, 18), "600000", 8.9, 9.1, 8.8, 9),
        (datetime.date(2024, 3, 19), "600000", 8.9, 9.1, 8.8, 9),
    ], dtype = dtypes
)

在这个数据结构中,共有6个字段,它们的名字和类型通过dtype来定义。这是一个List[Tuple]类型。在初始化数据部分,它也是一个List[Tuple]。

Warning

初学者很容易犯的一个错误,就是使用List[List]来初始化Numpy Structured Array,而不是List[Tuple]类型。这会导致Numpy在构造数组时,对应不到正确的数据类型,报出一些很奇怪的错误。
比如,下面的初始化是错误的:

1
2
3
4
secs = np.array([
    [datetime.date(2024, 3, 18), "600000", 8.9, 9.1, 8.8, 9],
    [datetime.date(2024, 3, 19), "600000", 8.9, 9.1, 8.8, 9]
], dtype=dtypes)
这段代码会报告一个难懂的"Type Error: float() argument must be a string or ..."

我们使用上一节学过的inspecting方法来查看secs数组的一些特性:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
print(f"secs的维度是{secs.ndim}")
print(f"secs的shape是{secs.shape}")
print(f"secs的size是{secs.size}")
print(f"secs的length是{len(secs)}")

print(f"secs[0]的类型是{type(secs[0])}")
print(f"secs[0]的维度是{secs[0].ndim}")
print(f"secs[0]的shape是{secs[0].shape}")
print(f"secs[0]的size是{secs[0].size}")
print(f"secs[0]的length是{len(secs[0])}")

可以看出,secs数组是一维数组,它的shape (2,)也正是一维数组的shape的表示法。前一节还介绍过这几个属性的关系,大家可以自行验证下是否仍然得到满足。


但secs的元素类型则是numpy.void,它在本质上是一个named tuple,所以,我们可以这样访问其中的任一字段:

1
2
3
4
print(secs[0]["frame"])

# 不使用列名(字段名),使用其序号也是可以的
print(secs[0][0])

我们还可以以列优先的顺序来访问其中的一个“单元格”:

1
print(secs["frame"][0])

对表格数据,遍历是很常见的操作,我们可以这样遍历:

1
2
for (frame, code, opn, high, low, close) in secs:
    print(frame, code, opn, high, low, close)

Numpy structured array在这部分的语法要比Pandas的DataFrame易用许多。我们在后面介绍Pandas时,还会提及这一点。

2. 运算类

2.1. 比较和逻辑运算

我们在上一节介绍定位、查找时,已经接触到了比较,比如:arr > 1。它的结果将数组中的每一个元素都与1进行比较,并且返回一个布尔型的数组。

现在,我们要扩充比较的指令:


函数 描述
all 如果数组中的元素全为真,返回True。可用以判断一组条件是否同时成立。
any 如果数组中至少有一个元素为真,则返回True。用以判断一组条件是否至少有一个成立
isclose 判断两个数组中的元素是否一一近似相等,返回所有的比较结果
allclose 判断两个数组中的元素是否全部近似相等
equal 判断两个数组中的元素是否一一相等,返回所有的比较结果。
not_equal 一一判断两个数组中的元素是否不相等,返回所有的比较结果
isfinite 是否为数字且不为无限大
isnan 测试是否为非数字
isnat 测试对象是否不为时间类型
isneginf 测试对象是否为负无限大
isposinf 测试对象是否为正无限大
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
# 开启多行输出模式
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

np.random.seed(78)
returns = np.random.normal(0, 0.03, size=4)
returns
# 判断是否全下跌
np.all(returns <= 0)
np.any(returns <= 0)

# 模拟一个起始价格为8元的价格序列
prices = np.cumprod(1+returns) * 8

# 对应的涨停价如下
buy_limit_prices = [8.03, 8.1, 8.1, 8.3]

# 判断是否涨停
np.isclose(prices, buy_limit_prices, atol=1e-2)

除了判断一个数组中的元素要么都为True,要么至少一个为True之外,有时候我们还希望进行模糊一点的判断,比如,如果过去20天中,超过60%的是收阳线,此时我们可以用np.count_nonzero,或者np.sum来统计数组中为真的情况:


1
2
np.count_nonzero(returns > 0)
np.sum(returns > 0)

在上一节进行比较的示例中,我们都只使用了单个条件。如果我们要按多个条件的组合进行查找,就需要依靠逻辑运算来实现。

在Numpy中,逻辑运算既可以通过函数、也可以通过运算符来完成:

函数 运算符 描述 python等价物
logical_and & 执行逻辑与操作 and
logical_or | 执行逻辑或操作 or
logical_not ~ 执行逻辑或操作 not
logical_xor '^' 执行逻辑异或操作 xor

逻辑运算有什么用呢?比如我们在选股时,有以下表格数据:

股票 pe mom
AAPL 30.5 0.1
GOOG 32.3 0.3
TSLA 900.1 0.5
MSFT 35.6 0.05

上述表格可以用Numpy的Structured Array来表示为:

1
2
3
4
5
6
tickers = np.array([
    ("APPL", 30.5, 0.1),
    ("GOOG", 32.3, 0.3),
    ("TSLA", 900.1, 0.5),
    ("MSFT", 35.6, 0.05)
], dtype=[("ticker", "O"), ("pe", "f4"), ("mom", "f4")])

现在,我们要找出求PE < 35, 动量 (mom) > 0.2的记录,那么我们可以这样构建条件表达式:


1
(tickers["pe"] < 35) & (tickers["mom"] > 0.2)

Numpy会把pe这一列的所有值跟35进行比较,然后再与mom与0.2比较的结果进行逻辑与运算,这相当于:

1
np.array((1,1,0,0)) & np.array((0, 1, 1, 0))
在Numpy中,True与1的值在做逻辑运算时是相等的;0与False也是。

如果不借助于Numpy的逻辑操作,我们就要用Python的逻辑操作。很不幸,这必须使用循环。如果计算量大,这将会比较耗时间。

在量化中使用异或操作的例子仍然最可能来自于选股。比如,如果我们要求两个选股条件,只能有一个成立时,才买入;否则不买入,就可以使用异或运算。

2.2. 集合运算

在交易中,我们常常要执行调仓操作。做法一般是,选确定新的投资组合,然后与当前的投资组合进行比较,找出需要卖出的股票,以及需要买入的股票。这个操作,就是集合运算。在Python中,我们一般是通过set语法来实现。

在Numpy中,我们可以使用通过以下方法来实现集合运算:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
import numpy as np

# 创建两个一维数组
x = np.array([1, 2, 3, 4, 5])
y = np.array([4, 5, 6, 7, 8])

# 计算交集
intersection = np.intersect1d(x, y)
print("Intersection (交集):", intersection)

# 计算并集
union = np.union1d(x, y)
print("Union (并集):", union)


1
2
diff = np.setdiff1d(x, y)
print("x - y:", diff)

此外,我们还可能使用in1d(a1, a2)方法来判断a1 中的元素是否都在 a2 中存在。比如,在调仓换股中,如果当前持仓都在买入计划中,则不需要执行调仓。

2.3. 数学运算和统计

Numpy中数学相关的运算有线性代数运算(当然还有基本代数运算)、统计运算、金融指标运算等等。

2.3.1. 线性代数

线性代数在量化中有重要用途。比如,在现代资产组合理论(MPT)中,我们要计算资产组合收益率及协方差,都要使用矩阵乘法。大家可以参考投资组合理论与实战系列文章,下面是其中的一段代码:

1
2
3
...
cov = np.cov(port_returns.T)
port_vol = np.sqrt(np.dot(np.dot(weights, cov), weights.T))

矩阵乘法是线性代数中的一个核心概念,它涉及到两个矩阵的特定元素按照规则相乘并求和,以生成一个新的矩阵。具体来说,如果有一个矩阵A 为 \(m \times n\) 维,另一个矩阵B 为 \(n \times p\) 维,那么它们的乘积 \(C = AB\)将会是一个\(m \times p\)维的矩阵。乘法的规则是A的每一行与B的每一列对应元素相乘后求和。

下面通过一个具体的例子来说明矩阵乘法的过程:

假设我们有两个矩阵A和B:


$$ A = \begin{bmatrix} 2 & 3 \ 1 & 4 \end{bmatrix} $$ 和 $$ B = \begin{bmatrix} 1 & 2 \ 3 & 1 \end{bmatrix} $$ 要计算AB,我们遵循以下步骤:

取A的第一行 \((2, 3)\) 与的第一列\((1,3)\) 相乘并求和得到\(C_{11} = [2\times1 + 3\times3 = 11]\)

同理,取A的第一行与B的第二列\((2, 1)\)相乘并求和得到\(C_{12} = [2\times2 + 3\times1 = 7]\)

取A的第二行\((1, 4)\)与B的第一列相乘并求和得到\(C_{21} = [1\times1 + 4\times3 = 13]\)

取A的第二行与B的第二列相乘并求和得到\(C_{22} = [1\times2 + 4\times1 = 5]\)

因此,矩阵C = AB为:

\[ C = \begin{bmatrix} 11 & 7 \\ 13 & 6 \ \end{bmatrix} \]

与代数运算不同,矩阵乘法不满足交换律,即一般情况下\(AB \neq BA\)

在Numpy中,我们可以使用np.dot()函数来计算矩阵乘法。


上述示例使用numpy来表示,即为:

1
2
3
4
A = np.array([[2,3],[1,4]])
B = np.array([[1,2],[3,1]])

np.dot(A, B)

最终我们将得到与矩阵C相同的结果。

除此之外,矩阵逆运算(np.linalg.inv)在计算最优投资组合权重时,用于求解方程组,特征值和特征向量(np.linalg.eig, np.linalg.svd)在分析资产回报率的主成分,进行风险分解时使用。

2.3.2. 统计运算

常用的统计运算包括:

函数 描述
np.mean 计算平均值
np.median 计算中位数
np.std 计算标准差
np.var 计算方差
np.min 计算最小值
np.max 计算最大值
np.percentile 用于计算历史数据的分位点
np.quantile 用于计算历史数据的分位数,此函数与percentile功能相同
np.corr 用于计算两个变量之间的相关性

np.percentile与np.quantile功能相同,都是用于计算分位数。两者在参数上略有区别。


当我们对同一数组,给quantile传入分位点0.25时,如果给percentile传入分位点25时,两者的结果将完全一样。也就是后者要乘以100。在量化交易中,quantile用得可能会多一些。

np.percentile(或者np.quantile)的常见应用是计算25%, 50%和75%的分位数。用来绘制箱线图(Boxplot)。

此外,我们也常用它来选择自适应参数。比如,在RSI的应用中,一般推荐是低于20(或者30)作为超卖,此时建议买入;推荐是高于80(或者70)作为超买,此时建议卖出。但稍微进行一些统计分析,你就会发现这些阈值并不是固定的。如果我们以过去一段时间的RSI作为统计,找出它的95%分位作为卖点,15%作为买点,往往能得到更好的结果。

2.3.3. 量化指标的计算

有一些常用的量化指标的计算,也可以使用Numpy来完成,比如,计算移动平均线,就可以使用Numpy提供的convolve函数。

1
2
3
import numpy as np
def moving_average(data, window_size):
    return np.convolve(data, np.ones(window_size)/window_size, 'valid')

当然,很多人习惯使用talib,或者pandas的rolling函数来进行计算。convolve(卷积)是神经网络CNN的核心,正是这个原因,我们这里提一下。

np.convolve的第二个参数,就是卷积核。这里我们是实现的是简单移动平均,所以,卷积核就是一个由相同的数值组成的数组,它们的长度就是窗口大小,它们的和为1。


如果我们把卷积核替换成其它值,还可以实现WMA等指标。从信号处理的角度看,移动平均是信号平滑的一种,使用不同的卷积核,就可以实现不同的平滑效果。

在量化中,还有一类计算,这里也提一下,就是多项式回归。比如,某两支股票近期都呈上升趋势,我们想知道哪一支涨得更好?这时候我们就可以进行多项式回归,将其拟合成一条直线,再比较它们的斜率。

下面的代码演示了如何使用Numpy进行多项式回归。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np
import matplotlib.pyplot as plt

returns = np.random.normal(0, 0.02, size=100)
alpha = 0.01
close = np.cumprod(1 + returns + alpha)

a, b = np.polyfit(np.arange(100), close, deg=1)

# 继续之前的代码

# 使用a, b构建回归线的y值
regression_line = a * np.arange(100) + b

# 绘制原始的close曲线
plt.figure(figsize=(10, 6))
plt.plot(close, label='Close Price', color='blue')

# 绘制回归线
plt.plot(regression_line, label='Regression Line', color='red', linestyle='--')

# 添加图例、标题和坐标轴标签
plt.title('Stock Close Price vs Regression Line')
plt.xlabel('Time Period')
plt.ylabel('Price')
plt.legend()

# 显示图表
plt.grid(True)
plt.show()

这将生成下图:

3. 类型转换和Typing

在不同的库之间交换数据,常常会遇到格式问题。比如,我们从第三方数据源拿到的行情数据,它们用的时间字段常常会是字符串(这是代码少写了几行吗?!)。有一些库在存储行情时,对OHLC这些字段进行了优化,使用了4个字节的浮点数,但如果要传给talib进行指标计算,就必须先转换成8个字节的浮点数,等等,这就有了类型转换的需求。

此外,我们还会遇到需要将numpy数据类型转换为python内置类型,比如,将numpy.float64转换为float的情况。


3.1. Numpy内部类型转换

Numpy内部类型转换,我们只需要使用astype函数即可。

1
2
3
4
5
6
7
8
x = np.array(['2023-04-01', '2023-04-02', '2023-04-03'])
print(x.astype(dtype='datetime64[D]'))

x = np.array(['2014', '2015'])
print(x.astype(np.int32))

x = np.array([2014, 2015])
print(x.astype(np.str_))

Tips

如何将boolean array转换成整数类型,特别是,将True转为1,False转为-1? 在涉及到阴阳线的相关计算中,我们常常需要将 open > close这样的条件转换为符号1和-1,以方便后续计算。这个转换可以用:

1
2
3
>>> x = np.array([True, False])
>>> x * 2 - 1
... array([ 1, -1])

3.2. Numpy类型与Python内置类型转换

如果我们要将Numpy数组转换成Python数组,可以使用tolist函数。

1
2
x = np.array([1, 2, 3])
print(x.tolist())

我们通过item()函数,将Numpy数组中的元素转换成Python内置类型。

1
2
3
x = np.array(['2023-04-01', '2023-04-02'])
y = x.astype('M8[s]')
y[0].item()

Warning

一个容易忽略的事实是,当我们从Numpy数组中取出一个标量时,我们都应该把它转换成为Python对象后再使用。否则,会发生一些隐藏的错误,比如下面的例子:

1
2
3
4
5
6
import json
x = np.arange(5)
print(json.dumps([0]))
print(x[0])

json.dumps([x[0]])
这里最后一行会出错。提示type int64 is not JSON serializable。把最后一行换成json.dumps([x[0].item()])则可以正常执行。

3.3. Typing

从Python 3.1起,就开始引入类型注解(type annotation),到Python 3.8,基本上形成了完整的类型注解体系。我们经常看到函数的参数类型注解,比如,下面的代码:

1
2
3
from typing import List
def add(a: List[int], b: int) -> List[int]:
    return [i + b for i in a]

从此,Python代码也就有了静态类型检查支持。

NumPy的Typing模块提供了一系列类型别名(type aliases)和协议(protocols),使得开发者能够在类型注解中更精确地表达NumPy数组的类型信息。这有助于静态分析工具、IDE以及类型检查器提供更准确的代码补全、类型检查和错误提示。

这个模块提供的主要类型是ArrayLike, NDArray和DType。


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import numpy
from numpy.typing import ArrayLike, NDArray, DTypeLike
import numpy as np

def calculate_mean(data: ArrayLike) -> float:
    """计算输入数据的平均值,数据可以是任何ArrayLike类型"""
    return np.mean(data)

def add_one_to_array(arr: NDArray[np.float64]) -> NDArray[np.float64]:
    """向一个浮点数数组的每个元素加1,要求输入和输出都是np.float64类型的数组"""
    return arr + 1

def convert_to_int(arr: NDArray, dtype: DTypeLike) -> NDArray:
    """将数组转换为指定的数据类型"""
    return arr.astype(dtype)

如果你是在像vscode这样的IDE中使用上述函数,你就可以看到函数的类型提示。如果传入的参数类型不对,还能在编辑期间,就得到错误提示。

4. 拓展阅读

4.1. Numpy的数据类型

在Numpy中,有以下常见数据类型。每一个数字类型都有一个别名。在需要传入dtype参数的地方,一般两者都可以使用。另外,别名在字符串类型、时间和日期类型上,支持得更好。比如,'S5'是Ascii码字符串别外,它除了指定数据类型之外,还指定了字符串长度。datetime64[S]除了表明数据是时间日期类型之外,还表明它的精度到秒。


类型 别名 类别 别名
np.int8 i1 np.float16 f2
np.int16 i2 np.float32 f4,还可指定结尾方式,比如'f4'表示大端字节序。其它float类型同。
np.int32 i4 np.float64 f8
np.int64 i8 np.float128 f16
np.uint8 u1 np.bool_ b1
np.uint16 u2 np.str_ U (后接长度,例如U10)
np.uint32 u4 np.bytes_ S (后接长度,例如S5)
np.uint64 u8 np.datetime64
np.timedelta64 m8和m8[D] m8[h] m8[m] m8[s]等

《因子投资与机器学习策略》开课啦!


9月8日晚,在线直播,不见不散


点击入会

金融/计量专业,硕士论文怎么确定研究课题?

面对毕业论文的压力,选择一个既具有实际应用价值又能激发研究兴趣的主题至关重要。对金融/计量专业的学生来说,在众多研究领域中,量化交易是兼顾自己的专长、又有利于未来发展的一个选择。

今天是第一期,我们介绍如何确定研究课题。

如何确定研究课题?

毕业论文最难的地方,是确定研究方向与课题。在科学研究中,提出正确的问题,就是解决了一半的问题 -- 亚里士多德。

在大师看来,有时候问题本身甚至比答案更为重要。据说有人曾建议希尔伯特去解决费马猜想问题,希尔伯特却笑了笑回答说: “我为什么要杀掉一只会下金蛋的鹅呢?” 在希尔伯特看来,一个像费马猜想那样的数学问题对数学的价值是无可估量的。

希尔伯特之墓,by By Kassandro - Own work, CC BY-SA 3.0

我们必须知道,我们必将知道

如果你有幸师从名师,确定研究方向和课题就比较容易,因为他们手中往往有一大堆问题,每一个都是一个宝库。如果没有这样的条件,这里也有一些建议。

了解行业的历史发展

确定研究课题首先要了解行业的历史发展。在量化交易领域,大概有这样几条主线。

一是资产定价研究。这是从上世纪50年代以来,由马科维茨和他的学生威廉.夏普等人开始,经过Stephen Ross,Fama等人努力,形成了因子投资与资产组合管理的理论体系。这也是经济学诺奖的一个大宝库。

二是衍生品与期权定价研究。这是由巴什利耶从1900年起开创,经过伊藤清、布莱克.费雪和迈伦.斯科尔斯等人的努力,最终以期权定价公式--BS公式提出为鼎盛时期的一个领域。

三是技术分析派。这一派以应用为主导,没有构建自己的理论体系,拿来主义比较多。比如,信息论的开山祖师香农贡献了网格交易法,他的同事,同在贝尔实验室的约翰.凯利发现了凯利公式,随后被另一名同事爱德华.索普用于赌场有股票市场中。

凯利公式-来自《孤注一掷》片头

信息论、统计、数字信号处理等学科都被用于交易,构成了量化交易很重要的一个分支。最著名的量化基金文艺复兴,早期把IBM viavoice实验室的人一锅端了,主要就是看中他们的数字信号上的研究能力。

现在则是进入了机器学习时代。通过机器学习和强化学习构建交易模型和策略的研究方兴未艾。如果说过去的金融需要的是数学和金融的复合人才,现在量化研究则更需要计算机和AI的人才。这是一个出论文研究的黄金赛道。

要完全理解这些方向,完全靠个人阅读是难以达到的。我们可以先从一些畅销书开始,再阅读发表在顶级刊物上的一些重要论文,这样可以快速掌握行业的概况。

畅销书

畅销书虽然专业性没那么强,但是对于把握研究方向、了解行业发展还是非常不错的,兼具知识性与趣味性。因为读起来没那么累,所以适合于初学者。

这些畅销书对建立自己的行业发展史观是很有帮助的:

定价未来,中信出版社

  1. 打开量化投资的黑箱,Rishi K. Narang著。本书被认为是系统性量化入门的最佳书籍之一。
  2. 宽客人生: 华尔街的数量金融大师。这本书介绍了西蒙斯、爱德华.索普等著名投资大师成功的秘诀。
  3. 高频交易员:华尔街的速度游戏。这本书介绍了颇为神秘的高频交易。看完后,你就可以放下对高频交易的执念了。因为这是一场军备竞赛,现在加入,已经没有太多机会了。
  4. 定价未来:撼动华尔街的量化金融史。这本书介绍了期权定价理论发展的历史。这是一位豆瓣网友推荐它的评语:

Quote

从第六章开始渐入佳境。怎么没有人用这样的方式讲我的大学课程呢?将人物八卦,和他们的成果,还有几代人的逐渐的智力演化、智力交锋、智力传承融合在一起。...这样写,更突出了文献检索的重要性。没有这些东西,怎么能看出智力的相互碰撞。

教科书中的众多名词、定理现在变成了活生生的人物,感觉这样的讲述方式多么生动有趣。黎曼积分、伊藤公式、马尔科夫过程、大卫希尔伯特问题、还有期权定价公式...

这也是为什么我推荐每个写论文的人,也去读这些畅销书的原因。你的大学之所以被毁掉,主要就是因为教科书写得面目可憎,读起来味同嚼蜡

现在,你应该对行业有了一个概观,大致知道了哪一个研究方向需要什么样的知识储备,未来前景如何,自己会不会感兴趣了。接下来,就要往下扎得深一些,开始接触一些专业论文了。

重要论文

这里介绍一些重要论文,你可以根据自己的方向来决定哪些是需要阅读的。

  1. Portfolio Selection, Markovitz, 1952。在这篇文章里,马科维茨提出了现代资产组合理论(MPT),使他获得了 1990 年诺贝尔经济学奖。该论文通过将风险与回报之间的相关性纳入计算中来扩展了常见的风险回报权衡。
  2. A New Interpretation of Information Rate,Kelly, 1956。本文是当今著名的凯利准则的正式表述。该模型广泛应用于赌场游戏,特别是风险管理。作者推导出了一个公式,该公式确定了最佳分配规模,以便随着时间的推移最大化财富增长。
  3. Capital Asset Prices: A Theory of Market Equilibrium under Conditions of Risk (Sharpe, 1964)。基于Markovitz的工作,CAPM证明了只有一种有效的投资组合,即市场投资组合。这篇文章提出了著名的Beta概念。Sharpe和Markovitz是师生关系,同年获得诺贝尔经济学奖。
  4. Efficient Capital Markets: a Review of Theory and Empirical Work,Fama, 1970。这篇论文是首次提出非常流行的“有效市场假说”概念的开创性著作,尽管现在这一理论受到很大的质疑,但是它的学术价值仍然很高。
  5. The Pricing of Options and Corporate Liabilities, 1973, Black & Scholes。著名的BS公式,利用物理传热方程作为估算期权价格的起点。这也是为什么对冲基金喜欢物理系学生的原因。
  6. Does the Stock Market Overreact? 1985, Bondt & Thaler。这篇文章质疑了有效市场假说,Bondt和Thaler提出,有统计上显着的证据表明相反的情况,投资者往往会对意外的新闻事件反应过度。这也是行为金融学的经典研究之一。行为金融学在近年来是诺奖的大热门。它的底层哲学是主观价值论和与人为本的思想。Thaler也是诺奖获得者,并在电影《大空头》中扮演了自己。
  7. A closed-form GARCH option valuation model,1997, Heston & Nandi。本文提出了一种封闭式公式,用于评估现货资产,并使用广义自回归条件异方差(GARCH) 模型对其方差进行建模。由于其复杂性和实用性, GARCH 模型在 20 世纪 90 年代估计波动率方面广受欢迎,金融业也积极采用它们。
  8. Optimal Execution of Portfolio Transactions,2000,Almgren & Chriss。对于每个负责完善交易执行算法的量化开发人员来说,这篇论文绝对是必读的学术文献。文章指出,价格波动来自外生性(市场本身的波动)和内生性(自己的订单对市场造成的影响)。这是一种量子效应!作者通过最小化交易成本和波动风险的组合,形式化了一种执行和衡量交易执行绩效的方法。
  9. Incorporating Signals into Optimal Trading,2017,Lehalle。与 Almgren 和 Chris (2000) 的工作非常相似,本文讨论的是最佳交易执行。作者通过将马尔可夫信号纳入最优交易框架,进一步完善了该领域所做的工作,并针对带有漂移的随机过程(Ornstein-Uhlenbeck 过程)的资产特殊情况得出了最优交易策略。

如果你的研究方向是因子投资,或者机器学习,那么还有更多、更专注于此方向的论文要阅读。这些论文在我们的课程中有介绍,数量太多,这里只介绍几个:

  1. The Performance of Mutual Funds in the Period 1945-1964, Michael Jessen。Sharpe在他的文章里提出了Beta这个概念,而Alpha这个概念,就是由Jessen在这篇论文中提出的。
  2. Common risk factors in the returns on stocks and bonds, 1993, Eugene F. Fama。在这篇文章里,Famma提出了三因子。
  3. Review of Financial Studies,2017, Stambaugh, Yuan。这篇论文出现的比较晚,因此可以把之前出现的、与因子投资相关的重要论文都评述一篇,因而就成为了我们快速了解行业的一篇重要论文。这篇文章发表比较晚,但也有了952次引用。
  4. 151 Trading Strategies, 2018, Zura Kakushadze。作者来自world quant,是Alpha101的作者之一。这篇论文引用了大量的论文(2000+),也是我们进行泛读的好材料。

每一篇论文,又会交叉引用到其它论文和资料,因此,读完这些资料,我们对自己感兴趣的领域的过去,应该有了一个比较全面的知识了解;接下来,需要读一点前沿资料,了解未来的发展。

不过,这件事比较有难度,因为每年出版的论文非常之多,根本看不过来,如何知道哪些论文是最重要的?这一部分,我们放在下一篇文章讲。

寻找课题

在我们对相关领域的历史有了比较深入的了解之后,也许你已经发现了一些尚未解决的问题;它们可能出现在你读过的论文的最后部分。但如果你仍然没有找到合适的问题,又该怎么办?

写论文实际上是一个创新的过程,毕竟,如果我们的论文全无新意,它也不可能成为一篇好的论文。

创新也是有迹可寻的。如果你读完了前面讲的那些畅销书,应该能发现一些线索。

跨界和联姻

跨界和联姻是最容易得到的创新。

BS公式部分来源于物理学中的热传导方程;爱德华.索普把凯利公式运用到赌场和股票,结果大获成功 -- 这本来是为研究长途电话线噪声而建立的理论。用在通讯方面的很多信号处理技术成为了量化算法;马科维茨是率先把统计学用在金融领域的人(这个说法可能不太正确,对很多人来说,保罗·萨缪尔森才是先驱。不过,这主要取决于如何区分金融学与经济学),仅靠均值方差理论就获得了诺贝尔奖。现在在量化中广泛使用的遗传算法,很显然来自于生物学的一些基本原理的启发。

Lcy++@cnblogs

现在,很多人都有一些机器学习和神经网络方面的知识储备。如果你一时想不到其它方向的联姻,那么可以试一试将AI与金融联系起来,并且尝试加密货币这样的新方向。

追逐新的风向和新的技术

在一些老的方向上,可能获得资料,建立知识体系比较容易,但是,这些方向上创新会比较困难。

建议尽可能找一些新的概念。加密货币尽管已经出现了快20年了,但它仍然是当红炸子鸡,是金融领域最新的方向。

在技术方向上,大模型和强化学习显然是AI上比较新、也可以嫁接到金融领域的方向之一。

寻求产业界前辈指导

通过前面两点建议,独立找到新的方向,可能对很多人来讲颇有难度。毕竟在研究生阶段,对业界的发展了解的不很太深,知识面也窄。要寻找新的问题,可以寻求产业界人士帮助。

产业界即使不能自己解决问题,他们也绝对是最先提出问题的人。实际上,在探索新课题上,学校往往都不如产业界领先。比如,计算机图形学的发展,并不是出于学术上的规划,而是一帮要打游戏的人玩出来的,后来驱动了人工智能的蓬勃发展。在量化领域也一样。

那么问题来了,作为一个普通的研究生,如何能得到前辈的指导呢? 在这个问题之前,还有一问,到哪里去寻找这些前辈呢?下面就介绍一个小技巧。

有一个六度空间理论,意思说,你与这地球上的任何一个人之间,最多只隔着5个人。要找到某个特定的人,需要找到恰当的节点。

在这个时代,知识博主就是非常好的结点。以我的小某书为例,我曾经解读过某个国外知名学者的论文,结果就在笔记评论区,遇到了他的学生。所以,如果有人要找这位学者,就有可能通过我的笔记找到他的学生,再接触到他。又比如,在我的课程中,据不完全统计(有的学员不愿意分享身份,我们也不会多问),有四五家私募(含基金)的老总级的人物(有的是为自己公司员工培训报名),所以在我的群里,就能链接上业界人士。

这只是拿我自己举例。我们总结一下方法,就是根据你要研究的方向,通过标签搜索,找到这个方向的所有的知识博主,快速过一下他们发布的内容,与你前面建立起来的知识体系对应,如果概念(关键词)重合度比较高,就关注这个博主,想办法进群。

这是第一步,现在,你很可能与你要认识的人在同一个群里了,能在公共场合说话,但无法进一步取得联系,因为你们之间的代差太大。

这里我不建议所谓的“向上社交”,我更推荐“价值互换”,平等社交。

我认识一位同学,她跟人互换的方式是帮大佬打理公众号,主要是排版校对方面的工作。在做这些工作的同时,自然也就有机会向大佬提出问题,得到指导。这就是价值互换,是健康的社交方式。

在群里不知道有没有大佬,不认识大佬怎么办?那就先从自己能帮到的人开始价值互换,认识了,再通过他认识更上一层的人。

我们还是从知识博主这个结点说起。知识博主需要的是流量,这就是我们可以帮到他的地方。我有一个券商的朋友,每天帮我点赞、转发,坚持了两个多月,我还没来得及帮到他,没想到他离职了。但我对他一直是很感谢的,如果他有任何要求,我能帮上忙的,肯定会帮。

一旦跟博主成了朋友,再提些要求,他肯定能帮忙的。

这个过程不能速成,人与人建立链接需要时间。但如果你是求知型的人,应该可以好好享受这个过程。

在今天的文章里,我介绍了量化方向上,最重要的一些论文。下一篇文章我将会介绍我平时使用的资源和工具。

带你读论文:PCA、离散小波和 XGBoost构建交易策略

这是 Nobre, Neves 发表于 2019 年的 一篇论文。在论文一起,生成了一个机器学习交易策略,取得了比 Buy-and-Hold 策略及另一个对照策略更好的回报。本文正文部分为原论文的摘要,在最后的 QuanTide 评论中,我提供了一些点评。

文章介绍了一种应用于金融领域的专家系统,该系统融合了主成分分析 (PCA)、离散小波变换 (DWT)、极限梯度提升 (XGBoost) 以及多目标优化遗传算法 (MOO-GA),旨在为投资者提供最佳的买卖点信号,以期在较低的风险水平下实现较高的投资回报。

PCA 用于缩减金融输入数据集的维度,DWT 用于对每个特征进行降噪处理,经过处理的数据集随后输入到 XGBoost 二元分类器中,其超参数由 MOO-GA 进行优化。

结果显示,PCA 显著提升了系统性能;而将 PCA 与 DWT 联合应用后,该系统能够在五个金融市场中的三个中超越传统的买入持有策略,实现了平均 49.26%的投资组合回报率,相比之下,买入持有策略的平均回报率为 32.41%。

相关工作

主成分分析

在相关工作中,主成分分析(PCA)被用于减少高维数据集的维度,同时保留数据的主要特征。PCA 通过对原始数据集进行线性组合来形成一组新的主成分,这些主成分最大限度地保留了原始数据的信息,并且具有高方差的特点。通过去除那些变化较小的维度,PCA 不仅简化了数据集,还提高了计算效率。

离散小波变换

离散小波变换(DWT)是一种强大的工具,它提供了处理现实世界问题中时间序列数据的时间变化特性的方式,与传统傅立叶变换相比,尤其适用于非平稳信号,如股票市场价格。

传统傅立叶变换假设信号是周期性的并且在整个时间区间上定义,而小波变换则允许信号表示为时间局部化的基函数的叠加。这使得小波变换能够同时捕捉信号的时间和频率信息。

小波变换的基础是将任何函数表示为一组构成小波变换基函数的小波的叠加。这些小波是有限长度振荡波形(称为母小波)的缩放和平移副本,即子小波。选择最佳的小波基依赖于待分析的原始信号特性和预期的分析目的。最终的结果是一组具有不同分辨率的时间-频率表示,这就是为什么小波变换被称为多分辨率分析的原因。

在本文提出的方法中,仅使用离散小波变换(DWT)。DWT 将信号分解为一组正交的小波,与连续小波变换相比,它更适合于实际应用,因为它可以更容易地实现数字计算机上的离散处理。

小波变换在金融时间序列数据预处理中的应用主要是为了去噪,从而帮助提高后续模型的预测精度。

遗传算法

遗传算法是一种元启发式优化方法,受到自然界生物进化过程的启发,用于解决复杂空间中的优化问题。它通过模拟自然选择和遗传机制来寻找问题的近似最优解。

遗传算法从一个由随机生成的个体组成的初始群体开始,这些个体代表了问题解的候选者。每个个体,或者说染色体,包含了一系列的参数或变量,这些变量被称为基因。

在每一代中,根据个体适应度(即个体在给定优化任务中的性能评分),选择个体进行繁殖。适应度较高的个体更有可能被选择,进而通过配对交换基因片段(交叉操作)产生新的后代。

此外,还会随机地对某些后代执行基因变异操作,以此引入新的遗传信息。

这些操作反复进行多代,直到满足预定的停止标准,例如达到最大迭代次数或群体收敛于稳定状态。最终的目标是希望经过一系列演化过程后,能够获得足够好的解。

XGBoost

XGBoost(Extreme Gradient Boosting)是一种优化的分布式梯度提升库,设计目的是为了实现高效、灵活和便携。它在传统梯度提升决策树的基础上做了改进,增加了一些增强性能的特性。具体来说,XGBoost 引入了正则化项来简化模型,帮助减少过拟合的可能性,并且支持并行处理,从而大大加快了训练的速度。

XGBoost 的一个重要特性是它允许用户自定义损失函数,并且能够自动处理缺失数据。它通过并行地构造多个树来提高计算效率,这与传统的梯度提升模型不同,传统的梯度提升模型通常是串行地建立每个树。此外,XGBoost 还可以自定义优化目标函数和评估指标,使得它非常适合于各种机器学习任务,包括在金融市场上预测股票价格的方向。

研究表明,当应用于股票市场价格方向预测时,XGBoost 作为分类器的表现超过了诸如支持向量机(SVM)和人工神经网络(ANN)等非集成方法。特别是在 60 天和 90 天的预测周期内,XGBoost 展示了更好的长期预测准确率。

提议中的方案

系统架构

目标方程

在本文中,目标是预测第 t+1 天的收盘价Closet+1相对于第 t 天的收盘价Closet是否会有一个正向或者负向的变化。因此,提出了一种监督学习解决方案,具体来说是一个二分类器。要预测的目标变量y是第 t+1 天相对于第 t 天收盘价变化的信号,它遵循一个二项式概率分布y ∈ {0, 1},其中: - 如果收盘价变化为正,则y取值为1; - 如果收盘价变化为负,则y取值为0

这个目标可以数学地定义如下:

\[ y_t = \begin{cases} 1 & \text{if } Closet_{t+1} - Closet_t \geq 0 \\ 0 & \text{if } Closet_{t+1} - Closet_t < 0 \end{cases} \]

其中: - Closet_{t+1} 是第 t+1 天的收盘价。 - Closet_t 是第 t 天的收盘价。

所有目标变量组成的数组命名为Y。金融输入数据集X是通过数据预处理模块输出的数据集,在其中应用了 PCA 和 DWT 技术于规范化后的数据集,该数据集包含了原始金融数据和技术指标。

技术分析模块

技术分析模块接收来自 Financial Data Module 的输出,并向其应用多个技术指标。使用技术指标的主要目的是,每个指标以彼此不同的方式提供有关过去原始数据的基本信息,因此将不同的技术指标组合在一起有助于检测金融数据中的模式,从而提高金融数据的性能预测系统。

该模块将创建数据集,该数据集将用作数据预处理模块的输入。该数据集由所使用的 26 个技术指标集和 5 个原始金融数据特征之间的组合组成,产生下表中所示的具有 31 个特征的数据集。技术指标的计算是使用 Python 库 TA-Lib 完成的。

数据预处理模块

数据标准化

略。

主成分分析

PCA 模块接收归一化输入数据集,包含 26 个技术指标和 5 个原始金融数据特征,总共 31 个归一化特征,为了降低数据过度拟合的风险并降低系统的计算成本,PCA 模块会将具有 31 个特征的数据集转换为较低维度的数据集,同时仍保留大部分原始数据集方差。

PCA 首先将其模型与归一化训练集进行拟合,以确定代表最大方差方向的分量。然后,主成分按照它们解释的方差量进行排序,并且仅保留那些加起来至少达到原始训练集方差 95% 的主成分。最后,将数据投影到主成分上,得到一个比原始数据集维度更低的数据集,因为它只保留了更好地解释特征之间关系的数据样本。然后将缩减后的数据集馈送到小波模块。

PCA 模块,使用了 Scikit-learn python 库。

小波变换

虽然数据集在 PCA 模块中已经进行了简化,降低了维度,只保留了能更好地解释特征之间关系的数据,但一些不相关的数据样本可能会产生负数,对系统训练和预测性能的影响可能仍然存在。

PCA 技术删除了特征子集中不相关的数据点,而 DWT 技术将对 PCA 减少的数据集中存在的每个特征在时域中执行降噪。这个过程减少了数据集中噪声的影响,同时尽可能保留每个特征的重要组成部分。

该系统中针对每个金融市场测试的小波基有:Haar 小波、Daubechies 小波和 Symlet 小波。对于 Daubechies 和 Symlet 小波,测试的阶数为 3、6、9、15 和 20。

尽管分解级别越高,可以消除的噪声就越多,从而可以更好地表示每个特征的趋势,但这也可能导致消除带有市场特征的波动。因此,在该系统中,测试的分解级别为 2、5 和 7,以找到每个金融市场的最优去噪分解级别。

DWT 首先指定小波基、阶数和所使用的分解级别。然后,对于训练集的每个特征,DWT 执行多级分解,这将产生一个近似系数和 j 个细节系数,其中 j 是所选的分解级别。

为了计算验证集和测试集的近似系数和细节系数,每次将一个数据点添加到训练集中,计算新信号的系数并保存与添加的数据点对应的系数,以便避免考虑未来信息的系数。执行此过程直到验证集和测试集中的所有点都计算出各自的系数。然后对获得的细节系数进行阈值处理并重建信号,从而产生每个原始数据集特征的去噪版本。应用 DWT 后,数据集被馈送到 XGBoost 模块。

PyWavelets 和 scikit-image 库被用于开发该系统中的 DWT 模块。

XGBoost 模块

XGBoost 二元分类器负责系统的分类过程。

二元分类器

分类器的输出 ^y 是给定当前观测值 x 的预测值,对应于实际日期 t。变量 ˆy 在 [0,1] 范围内,所有 ˆy 的集合对应于预测交易信号,该信号指示系统在下一个交易日是否应该持有多头或空头头寸。

在 XGBoost 二元分类器算法开始之前,必须选择一组参数。定义机器学习系统架构的参数称为超参数。由于每个时间序列都有其自身的特点,因此对于每个分析的时间序列都有一组不同的最优超参数,即使模型具有良好的泛化能力,从而在样本外数据中取得良好结果的超参数。因此,为了在每个分析的金融市场中获得最佳结果,必须找到最佳的超参数集。为此,使用了 MOO-GA 方法,如下一节所述。

数据预处理模块输出的预处理数据被馈送到 XGBoost 二元分类器,以及待预测的目标变量数组 Y,并且均被分为训练集、验证集和测试集。定义数据集和 XGBoost 超参数后,训练阶段就开始了。系统的泛化能力是通过系统处理未见数据的表现来衡量的。因此,在训练阶段结束后,使用样本外验证集来测试训练阶段获得的模型,以验证所获得模型的泛化能力。

该验证集在 MOO 过程中使用,有助于使用未见过的数据选择性能最佳的解决方案。训练和验证阶段结束后,将创建最终模型,即使用 MOO-GA 找到的最佳超参数集进行训练的模型,并将生成的输出与测试集进行比较,以验证预测的质量。如前所述,输出采用数据点属于 0 或 1 类之一的概率形式。

XGBoost 库用于开发该系统中的 XGBoost 二元分类器。

多目标优化遗传算法

建机器学习模型时,模型架构有很多设计选择。大多数时候,用户事先并不知道给定模型的架构应该是什么样子,因此探索一系列可能性是可取的。

XGBoost 二元分类器的超参数就是这种情况,这些参数定义了分类器的架构。这个寻找理想模型架构的过程称为超参数优化。

采用多目标优化方法而不是单目标优化方法,是因为我们的最终目标要实现一个高回报同时低风险的交易系统。因此,自然必须使用统计度量来评估系统相对于所做预测的性能 -- 在本例中是 Accuracy;但另一方面,也需要使用度量来评估系统实现良好回报的能力,同时最大限度地减少损失 -- 在这种情况下是 Sharpe Ratio。

因此,根据两个选定的目标函数来评估候选解决方案:Accuracy 和 Sharpe Ratio。每个解的集合代表了 MOO-GA 要优化的适应度函数,其目标是最大化每个目标函数。因此,给定 XGBoost 二元分类器、金融数据集 X 和目标变量数组 Y,MOO-GA 将搜索和优化 XGBoost 二元分类器超参数集,目标是最大化所获得的预测的准确性和夏普比率。

为了找到旨在最大化前面提到的两个目标函数的最佳超参数集,MOO-GA 方法基于非支配排序遗传算法-II (NSGA-II)。

由于 XGBoost 二元分类器中存在许多超参数,因此只有那些对二元分类器的架构有重大影响并因此对其整体性能影响较大的超参数才会被优化。选择要优化的超参数对偏差-方差权衡也有较大影响,它们是:学习率、最大树深度、最小子体重和子样本,这些超参数中的每一个都构成一个基因 MOO-GA 中的染色体。

在提出的 MOO-GA 中,使用了两点交叉算子,其中在父级字符串上选择两个点,并且两个选定点之间的所有内容都在父级之间交换。选择的突变率为 0.2,这意味着交叉算子生成的每个新候选解都有 20% 的概率发生突变。

我们还使用了超突变,这是一种在进化算法中将多样性重新引入群体的方法。在该系统中,在进化过程中调整变异率,以帮助算法跳出局部最优。

下图显示了重要的超参数:

DEAP python 库用于开发该系统中的 MOO-GA 模块。

交易模块

交易模块负责模拟真实的金融市场交易。其主要功能包括:

  • 接收输入:从 XGBoost 模块获取交易信号和金融数据。
  • 模拟市价订单:基于输入的交易信号,在金融市场中模拟市价订单。
  • 状态机设计:采用多头、空头和持有三种状态的状态机来执行交易。具体而言,给定交易信号后,状态机会根据信号中的市场订单执行相应的操作。XGBoost 二元分类器的预测结果 ^y 表示类别 p(^y) 的预测概率。由于选定的两个类别分别表示第 t+1 天与第 t 天收盘价的变化情况,因此可以利用这些预测结果构建交易信号。

交易信号构建方法如下:

  • 如果 p( ˆy) ≥ 0.5,且所选类别为 1,这意味着股票在 t + 1 天的收盘价(收盘价)预计将出现正变化,从而代表在 t 天有买入机会,即多头头寸被采取。该动作由训练信号中相应日期的位置 1 表示;
  • 相反,如果 p( ˆy) < 0.5,且所选类别为 0,这意味着股票在 t + 1 天的收盘价(收盘价)预计将出现负变化,从而代表在 t 天的卖出机会,即采取空头头寸。该动作由相应日期的训练信号中的位置 0 表示。

因此,交易信号的值在 [0,1] 范围内,交易模块负责解释这些值并将其转换为交易动作。

实验结果

下图显示了训练过程的示意。

案例 1 使用 PCA 的效果

在第一个案例研究中,分析了 PCA 技术对所实现系统性能的影响。

对照系统分为两种配置:

基本系统:输入数据集包含 31 个归一化的金融特征。 改进系统:输入数据集经过标准化处理。 在改进系统中,应用 PCA(主成分分析)技术来降低数据集的维度。具体来说,对包含 31 个金融特征的标准化数据集进行 PCA 处理后,每个金融市场保留了 6 个主成分。这意味着原始的 31 个金融特征被缩减为 6 个低维特征。

下图列出了每个系统的实验结果以及买入并持有策略的结果。每个金融市场获得的最佳结果用粗体突出显示。

通过检查获得的结果可以得出结论,PCA 技术的使用在提高系统性能方面发挥着重要作用,因为它不仅可以获得更高的回报,而且还可以获得更高的准确度值。

PAC 降维使得 XGBoost 二元分类器能够生成复杂度较低的模型,具备良好的泛化能力,从而避免过拟合训练数据。这样,系统能够在玉米期货合约和埃克森美孚股票上实现比 Buy & Hold 策略更高的回报。

在其他三个分市场中,使用降维技术对于系统获得正回报至关重要,因为在基本系统中,泛化能力不足以正确分类测试集。如果没有应用 PCA,系统只能在埃克森美孚公司股票上实现比买入并持有策略更高的回报。

图 4 显示了在测试期间使用基本系统和带有 PCA 的系统所得回报的对比图。可以看出,正如前面提到的,引入 PCA 后,系统能够获得更高的回报。

案例 2 使用离散小波分析的效果

在这个案例中,我们研究将 DWT 技术与 PCA 结合使用,以同时实现对输入数据的降维和降噪,看这两种技术一起应用是否能获得更好的结果。

下图列出了每个系统的结果及买入并持有策略的结果。每个金融市场获得的最佳结果用粗体突出显示。由于所有组合分析过于复杂,图中仅展示了两种最佳组合。

通过检查获得的结果可以得出结论,PCA 和 DWT 去噪技术的结合使系统比仅使用 PCA 的系统获得更好的结果。

这是因为在该系统中,PCA 降低了金融输入数据集的维度,并且 DWT 对这个降低的数据集进行了去噪,这不仅有助于避免过度拟合训练数据,也有助于去除一些可能损害系统性能的不相关样本,从而帮助系统的学习过程并提高其泛化能力。

然而,这种性能的提高只有在使用适当的小波基、阶数和分解级别时才能得到验证,因为与仅使用 PCA 相比,错误地使用 DWT 会导致更差的系统性能,甚至更低的回报。

案例 3 性能评估

在本案例研究中,将所提出的系统的性能与 Nadkarni,Neves 在 2018 年提取的一个策略进行比较,以验证使用所提出的方法与其他方法的优缺点。

下图显示了两个系统的性能对比。

结论

本文提出了一种结合 PCA 和 DWT 分别进行降维和降噪的系统,以及使用 MOO-GA 优化的 XGBoost 二元分类器,目的是实现一个尽可能的最大回报,同时最小化风险水平交易策略。实验结果支持了此方法的有效性。

QuanTide 评论

本文提出了基于机器学习构建交易策略的一个完整框架,并创造性地使用了DWT和GA等先进算法。这个交易系统使用时序特征作为训练特征,对股票、期货和加密货币都有较好的覆盖。注意它不是一个资产定价类的策略。

论文中有几处地方值得商榷,也是构建机器学习交易策略的难点

数据标注问题

论文中实现的是一个监督学习算法,数据标注不可或缺。论文使用的标注方法是:

对于t期数据,如果pnl > 0,则标记为1,否则为0。

这里的问题在于,如果一支股票的 Pnl 绝对值是 0.5%以下,它实际上并不有很强的标签意义 -- 换言之,这个结果并不一定反映主力的操作意图,它很可能只是主力缺席之后,股价随市场正常波动的一个状态。这种情况下,无论是归为1,还是归为0,都可能是错误的。为了抵消这种错误标注的负面影响,往往需要我们生产更多正确的标注才行。

遗传算法

在寻找模型的超参数时,论文使用了遗传算法。使用遗传并非必须,它只是快速找到超参数的一个算法。我们也可以使用 GridSearch 或者 RandomSearch 的方法。

我们注意到论文中并没有像PCA和DWT那样,给出使用MOO-GA与不使用MOO-GA的对比。这也验证了我们的观点。

使用DWT来去除噪声

这部分有两点可以拿出来讨论。一是论文运用DWT的阶段。它是放在TA之后。如果你熟悉技术分析(TA)的话,就知道原始的行情数据经过TA处理成技术指标之后,它已经把原始数据在时序上的特征提取出来了,而这些技术指标序列,就不见得是一种波信号了。所以,此时使用DWT进行滤波,意义何在?难以解释。

文中也提到,不是每一次使用DWT都能获得好的结果,可能原因就在这里。

除了使用DWT的阶段之外,使用DWT的用处也值得讨论。多数论文都是使用DWT来进行滤波 -- 在这件事情上,无论他们给出的数据有多好,我始终认为,都不可能比得过SMA -- SMA有明确的金融含义,其它方法只有数学上的精巧。

但是,如果你认为 DWT 能够去除噪声 -- 那么一个显而易见的推论就是,它知道哪些是信号,哪些是噪声。既然如此,为什么我们不使用小波分析来提取信号特征来进行学习?

所以,更好的方案应该是将DWT作为与TA模块并生的一个模块,用来生成特征。

关于归一化算法

论文中提到,他们使用了Min-Max来进行归一化。这一点上是个灾难的选择。我们只能对取值范围固定的数据进行Min-Max归一化,而股价的值域从理论上看,是\((0, +\infty)\)

在我们学习量化的过程中,参考已发表的论文无疑是成长的捷径。但是,有些学者并身没有太多的实际交易经验,导致论文中也可能存在各种瑕疵,这无疑也会导致论文的结果无法在实战中运用。因此,我们开设了《因子投资与机器学习策略》的课程,如果对这一领域感兴趣,但不得其门而入的,可以加入我们一起学习。

该论文可在这里下载

[0901] QuanTide Weekly

本周要闻

  • 市场传闻存量房贷利率下调,房地产 ETF 大涨,但尾盘多股炸板
  • 中国 8 月官方制造业 PMI 为 49.1% 比上月下降 0.3 个百分点
  • 国家市监总局宣布阿里整改完成
  • 半年报第一股!桐昆股份同比增长 911.35%,为已发布半年报公司中净利润增速最高。

下周看点

  • 下调存量房贷利率传闻能否兑现?
  • 周一财新制造业 PMI 发布
  • 周五,2024 低空经济发展大会将于 9 月 6 日至 8 日在芜湖举行

本周精选

  • 主力正在进场?快速傅里叶变换与股价预测研究
  • 连载!量化人必会的 Numpy 编程 (1)

本周要闻

  • 市场传出的信息,有关方面正在考虑进一步下调存量房贷利率,允许规模高达 38 万亿元人民币的存量房贷寻求转按揭,以降低居民债务负担、提振消费。截至周六,上述传闻并未获得官方证实。
  • 8 月份,PMI 为 49.1%,环比下降 0.3 个百分点。从企业规模看,大型企业 PMI 为 50.4%,仍高于临界点;中、小型企业 PMI 分别为 48.7%和 46.4%,环比下降 0.7 和 0.3 个百分点。8 月 PMI 受高温天气影响,并未显著超预期
  • 8 月 30 日下午,国家市场监督管理总局宣布阿里巴巴完成三年整改,取得良好成效。阿里巴巴股价从 2020 年峰值以来,跌去七成。
  • 美东时间周五,美股三大指数集体收涨,道指收涨 0.55%,报 41563.08 点,创历史新高。标普 500 上涨 1.01%,纳斯达克上涨 1.13%。美联储青睐的 7 月份 PCE 通胀数据基本符合预期,市场押注 9 月美联储大幅降息的可能性减少,但市场仍然预计 11 月或 12 月可能会有大幅降息。
  • 易方达纳斯达克 100ETF 发布 2024 年中期报告,常州投资集团持有 5.92%份额,成为第一大持有人。该 ETF 在去掉净值上涨 49.21%的情况下,今年仍上涨 14.91%。该 ETF2017 年发行,现净值 3.16。
  • 半年报第一股!桐昆股份发布半年报,公司实现净利润 10.65 亿元,同比增长 911.35%,为已发布半年报公司中净利润增速最高。报告期内,涤纶长丝行业下游需求较去年同期边际改善显著,产品销量与价差有所增大,行业整体处于复苏状态。整体看,电子行业成大赢家,与上年同期相比,电子行业营收增速高居首位,上半年整体营收为 1.59 万亿元,同比增长 17.3%。

信息来源:东方财富网站


下周看点

  • 周五市场传闻,下调存量房贷利率正在考虑中,随即房地产 ETF 大涨,银行股大跌以回应传闻,但尾盘回落,涨停个股纷纷炸板。下周,此传闻是被证实还是被证伪?或将对市场有重要影响。
  • 周一财新制造业 PMI 发布
  • 周五,2024 低空经济发展大会将于 9 月 6 日至 8 日在芜湖举办
  • 周五,美国 8 月失业率和非农报告公布

主力正在进场?快速傅里叶变换与股价预测研究

一个不证自明的事实:经济活动是有周期的。但是,这个事实似乎长久以来被量化界所忽略。无论是在资产定价理论,还是在趋势交易理论中我们几乎都找不到周期研究的位置 -- 在后者语境中,大家宁可使用“摆动”这样的术语,也不愿说破“周期”这个概念。

这篇文章里,我们就来探索股市中的周期。我们将运用快速傅里叶变换把时序信号分解成频域信号,通过信号的能量来识别主力资金,并且根据它们的操作周期进行一些预测。最后,我们给出三个猜想,其中一个已经得到了证明。

FFT - 时频互转

(取数据部分略)。

我们已经取得了过去一年来的沪指。显然,它是一个时间序列信号。傅里叶变换正好可以将时间序列信号转换为频域信号。换句话说,傅里叶变换能将沪指分解成若干个正弦波的组合。

1
2
3
4
5
6
7
8
# 应用傅里叶变换
fft_result = np.fft.fft(close)
freqs = np.fft.fftfreq(len(close))

# 逆傅里叶变换
filtered = fft_result.copy()
filtered[20:] = 0
inverse_fft = np.fft.ifft(filtered)

1
2
3
4
5
# 绘制原始信号和分解后的信号
plt.figure(figsize=(14, 7))
plt.plot(close, label='Original Close')
plt.plot(np.real(inverse_fft), label='Reconstructed from Sine Waves')
plt.legend()

我们得到的输出如下:

在数字信号处理的领域,时间序列被称为时域信号,经过傅里叶变换后,我们得到的是频域信号。时域信号与频域信号可以相互转换。Numpy 中的 fft 库提供了 fft 和 ifft 这两个函数帮我们实现这两种转换。

np.fft.fft 将时域信号变换为频域信号,转换的结果是一个复数数组,代表信号分解出的各个频率的振幅 -- 也就是能量。频率由低到高排列,其中第 0 号元素的频率为 0,是直流分量,它是信号的均值的某个线性函数。

np.ff.ifft 则是 fft 的逆变换,将频域信号变换为时域信号。

将时域信号变换到频域,就能揭示出信号的周期性等基本特征。我们也可以对 fft 变换出来的频域信号进行一些操作之后,再变换回去,这就是数字信号处理。


高频滤波和压缩

如果我们把高频信号的能量置为零,再将信号逆变换回去,我们就会得到一个与原始序列相似的新序列,但它更平滑 -- 这就是我们常常所说的低通滤波的含义 -- 你熟悉的各种移动平均也都是低通滤波器。

在上面的代码中,我们只保留了前 20 个低频信号的能量,就得到了与原始序列相似的一个新序列。如果把这种方法运用在图像领域,这就实现了有损压缩 -- 压缩比是 250/20。

在上世纪 90 年代,最领先的图像压缩算法都是基于这样的原理 -- 保留图像的中低频部分,把高频部分当成噪声去掉,这样既保留了图像的主要特征,又大大减少了要保存的数据量。


当时做此类压缩算法的人都认识这位漂亮的小姐姐 -- Lena,这张照片是图像算法的标准测试样本。在漫长的进化中,出于生存的压力,人类在识别他人表情方面进化出超强的能力。所以相对于其它样本,一旦压缩造成图像质量下降,肉眼更容易检测到人脸和表情上发生的变化,于是人脸图像就成了最好的测试样本。

Lena

Lena 小姐姐是花花公子杂志的模特,这张照片是她为花花公子 1972 年 11 月那一期拍摄的诱人照片的一小部分 -- 在原始的照片中,Lena 大胆展现了她诱人的臀部曲线,但那些不正经的科学家们只与我们分享了她的微笑 -- 从科研的角度来讲,这也是信息比率最高的部分。

无独有偶,在 Lena 成为数字图像处理的标准测试样本之前,科学家们一直使用的是另一位小姐姐的照片,也出自花花公子。

好,言归正传。我们刚刚分享了一种方法,去掉信号中的高频噪声,使得信号本身的意义更加突显出来。我们也希望在证券分析中使用类似的技法,使得隐藏在 K 线中的信号显露出来。

但如果我们照搬其它领域这一方法,这几乎就不算研究,也很难获得好的结果。实际上,在证券信号中,与频率相比,我们更应该关注信号的能量,毕竟,我们要与最有力量的人站在一边。


愿原力与你同在 -- 星球大战

所以,我们换一个思路,把分解后的频域信号中,能量最强的部分保留下来,看看它们长什么样。

过滤低能量信号

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# 保留能量最强的前 5 个信号
amp_threshold = np.sort(np.abs(fft_result))[-11]

# 绘制各个正弦波分量
plt.figure(figsize=(14, 7))

theforce = []
for freq in freqs:
    if freq == 0:  # 处理直流分量
        continue
    elif freq < 0:
        continue
    else:
        amp = np.abs(fft_result[np.where(freqs == freq)])
        if amp < amp_threshold:
            continue

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
        sine_wave = amp * np.sin(2 * np.pi * freq * np.arange(len(close)))
        theforce.append(sine_wave)
        plt.plot(dates, sine_wave, label=f'Frequency={freq:.2f}')

plt.legend()
plt.title('Individual Sine Wave Components')
ticks = np.arange(0, len(dates), 20)
labels_to_show = [dates[i] for i in ticks]
plt.xticks(ticks=ticks, labels=labels_to_show, rotation=45)
plt.show()

FFT 给出的频率总是一正一负,我们可以简单地认为,负的频率对我们没有意义,那是一种我们看不到、也无须关心的暗能量。所以,在代码中,我们就忽略了这一部分。

我们看到,对沪指走势影响最强的波(橙色)的周期是 7 个月左右:从峰到底要走 3 个半月,从底到峰也要走 3 个半月。由于它的能量几乎是其它波的一倍,所以它是主导整个叠加波的走势的:如果其它波与它同相,叠加的结果就会使得趋势加强;反之,则是趋势抵消。其它几个波的能量相仿,但频率不同。


这些波倒底是什么呢?它可以是一种经济周期,但是说到底,经济周期是人推动的,或者反应了人的判断。因此,我们可以把波动的周期,看成资金的操作周期

从这个分解图中,我们可以猜想,有一个长线资金(对应蓝色波段),它一年多调仓一次。有一个中线资金(对应橙色波段),它半年左右调一次仓。其它的资金则是短线资金,三个月左右就会做一次仓位变更。还有无数被我们过滤掉的高频波段,它们操作频繁,可能对应着散户,但是能量很小,一般都可以忽略;只有在极个别的时候,才能形成同方向的叠加,进而影响到走势。

现在,我们把这几路资金的操作合成起来,并与真实的走势进行对比,看看情况如何:

在大的周期上都基本吻合,也就是这些资金基本上左右了市场的走势。而且,我们还似乎可以断言,在 3 月 15 到 5 月 17 这段时间,出现了股价与主力资金的背离趋势:主力资金在撤退了,但散户还在操作,于是,尽管股价还在上涨,但最终的方向,由主力资金来决定。


Tip

黑色线是通过主力资金波段合成出来的(对未来有预测性),在市场没有发生根本性变化之前,主力的操作风格是相对固定的,因此,它可能有一定的短时预测能力。如果我们认可这个结论的话。那么就应该注意到,末端部分还存在另一个背离 -- 散户还在离场,但主力已经进场。当然,关于这一点,请千万不要太当真。

关于直流分量的解释

我过去一直以为直流分量表明资产价格的趋势,但实际上所有的波都是水平走向的 -- 但只有商品市场才是水平走向的,股票市场本质上是向上的。所以,直流分量无法表明资产价格的趋势。

直到今天我突然涌现一个想法:如果你把一个较长的时序信号分段进行 FFT 分解,这样会得到若干个直流分量。这些直流分量的回归线才是资产价格的趋势。

这里给出三个猜想:

  1. 如果分段分解后,各个频率上的能量分布没有显著变化,说明投资者的构成及操作风格也没有显著变化,我们可以用 FFT 来预测未来短期走势,直到条件不再满足为止。

  2. 沪指 30 年来直流分量应该可以比较完美地拟合成趋势线,它的斜率就等于沪指 20 年回归线的斜率。

  3. 证券价格是直流分量趋势线与一系列正弦波的组合。

下面我们来证明第二个猜想(过程略)。最终,我们将直流分量及趋势线绘制成下图:

75%

而 2005 年以来的 A 股年线及趋势线是这样的:

75%

不能说十分相似,只能说几乎完全一致。

趋势线拟合的 p 值是 0.055 左右,也基本满足 0.05 的置信度要求。


这篇文章是我们《因子投资与机器学习策略》中的内容,出现在如何探索新的因子方法论那一部分。对 FFT 变换得到的一些重要结果,将成为机器学习策略中用以训练的特征。更多内容,我们课堂上见!


量化人必会的 NUMPY 编程 (1) - 核心语法

1. 基本数据结构

NumPy 的核心数据结构是 ndarray(即 n-dimensional array,多维数组)数据结构。这是一个表示多维度、同质并且大小固定的数组对象。

ndarray 只能存放同质的数组对象,这样使得它无法表达记录类型的数据。因此,numpy 又拓展出了名为 structured array 的数据结构。它用一个 void 类型的元组来表示一条记录,从而使得 numpy 也可以用来表达记录型的数据。因此,在 Numpy 中,实际上跟数组有关的数据类型主要是两种。

前一种数组格式广为人知,我们将以它为例介绍多数 Numpy 操作。而后一种数据格式,在量化中也常常用到,比如,通过聚宽 [^聚宽] 的 jqdatasdk 获得的行情数据,就允许返回这种数据类型,与 DataFrame 相比,在存取上有不少简便之处。我们将在后面专门用一个章节来介绍。

在使用 Numpy 之前,我们要先安装和导入 Numpy 库:

1
2
# 安装 NUMPY
pip install numpy

一般地,我们通过别名np来导入和使用 numpy:

1
import numpy as np

为了在 Notebook 中运行这些示例时,能更加醒目地显示结果,我们首先定义一个 cprint 函数,它将原样输出提示信息,但对变量使用红色字体来输出,以示区别:

1
2
3
4
5
6
7
8
from termcolor import colored

def cprint(formatter: str, *args):
    colorful = [colored(f"{item}", 'red') for item in args]
    print(formatter.format(*colorful))

# 测试一下 CPRINT
cprint("这是提示信息,后接红色字体输出的变量值:{}", "hello!")

接下来,我们将介绍基本的增删改查操作。

1.1. 创建数组

1.1.1. 通过 Python List 创建

我们可以通过np.array的语法来创建一个简单的数组:

1
2
arr = np.array([1, 2, 3])
cprint("create a simple numpy array: {}", arr)

在这个语法中,我们可以提供 Python 列表,或者任何具有 Iterable 接口的对象,比如元组。

1.1.2. 预置特殊数组

很多时候,我们希望 Numpy 为我们创建一些具有特殊值的数组。Numpy 也的确提供了这样的支持,比如:


函数 描述
zeros
zeros_like
创建全 0 的数。zeros_like 接受另一个数组,并生成相同形状和数据类型的 zeros 数组。常用于初始化。以下*_like 类推。
ones
ones_like
创建全 1 的数组
full
full_like
创建一个所有元素都填充为n的数组
empty
empty_like
创建一个空数组
eye
identity
创建单位矩阵
random.random 创建一个随机数组
random.normal 创建一个符合正态分布的随机数组
random.dirichlet 创建一个符合狄利克雷分布的随机数组
arange 创建一个递增数组
linspace 创建一个线性增长数组。与 arange 的区别在于,此方法默认生成全闭区间数组。并且,它的元素之间的间隔可以为浮点数。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# 创建特殊类型的数组
cprint("全 0 数组:\n{}", np.zeros(3))
cprint("全 1 数组:\n{}", np.ones((2, 3)))
cprint("单位矩阵:\n{}", np.eye(3))
cprint("由数字 5 填充的矩阵:\n{}", np.full((3,2), 5))

cprint("空矩阵:\n{}", np.empty((2, 3)))
cprint("随机矩阵:\n{}",np.random.random(10))
cprint("正态分布的数组:\n{}",np.random.normal(10))
cprint("狄利克雷分布的数组:\n{}",np.random.dirichlet(np.ones(10)))
cprint("顺序增长的数组:\n{}", np.arange(10))
cprint("线性增长数组:\n{}", np.linspace(0, 2, 9))

Warning

尽管 empty 函数的名字暗示它应该生成一个空数组,但实际上生成的数组,每个元素都是有值的,只不过这些值既不是 np.nan,也不是 None,而是随机值。我们在使用 empty 生成的数组之前,一定要对它进行初始化,处理掉这些随机值。


生成正态分布数组很有用。我们在做一些研究时,常常需要生成满足某种条件的价格序列,再进一步研究和比较它的特性。

比如,如果我们想研究上升趋势和下降趋势下的某些指标,就需要有能力先构建出符合趋势的价格序列出来。下面的例子就演示了如何生成这样的序列,并且绘制成图形:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import numpy as np
import matplotlib.pyplot as plt

returns = np.random.normal(0, 0.02, size=100)

fig, axes = plt.subplots(1, 3, figsize=(12,4))
c0 = np.random.randint(5, 50)

for i, alpha in enumerate((-0.01, 0, 0.01)):
    r = returns + alpha
    close = np.cumprod(1 + r) * c0
    axes[i].plot(close)

绘制的图形如下:

示例中还提到了 Dirichlet(狄利克雷)分布数组。这个数组具有这样的特点,它的所有元素加起来会等于 1。比如,在现代投资组合理论中的有效前沿优化中,我们首先需要初始化各个资产的权重(随机值),并且满足资产权重之和等于 1 的约束(显然!),此时我们就可以使用 Dirichlet[^dirichlet] 分布。


Info

狄利克雷,德国数学家。他对数论、傅里叶级数理论和其他数学分析学领域有杰出贡献,并被认为是最早给出现代函数定义的数学家之一和解析数论创始人之一。

1.1.3. 通过已有数组转换

我们还可以从已有的数组中,通过复制、切片、重复等方法,创建新的数组

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
# 复制一个数组
cprint("通过 np.copy 创建:{}", np.copy(np.arange(5)))

# 复制数组的另一种方法
cprint("通过 arr.copy: {}", np.arange(5).copy())

# 使用切片,提取原数组的一部分
cprint("通过切片:{}", np.arange(5)[:2])

# 合并两个数组
arr = np.concatenate((np.arange(3), np.arange(2)))
cprint("通过 concatenate 合并:{}", arr)

# 重复一个数组
arr = np.repeat(np.arange(3), 2)
cprint("通过 repeat 重复原数组:{}", arr)

# 重复一个数组,注意与 NP.REPEAT 的差异
# NP.TILE 的语义类似于 PYTHON 的 LIST 乘法
arr = np.tile(np.arange(3), 2)
cprint("通过 tile 重复原数组:{}", arr)

Question

np.copy 与 arr.copy 有何不同?在 Numpy 中还有哪些类似函数对,有何规律?


注意在 concatenate 函数中,axis 的作用:

1
2
3
4
5
6
arr = np.arange(6).reshape((3,2))

# 在 ROW 方向上拼接,相当于增加行,默认行为
cprint("按 axis=0 拼接:\n{}", np.concatenate((arr, arr), axis=0))
# 在 COL 方向上拼接,相当于扩展列
cprint("按 axis=1 拼接:\n{}", np.concatenate((arr, arr), axis=1))

1.2. 增加/删除和修改元素

Numpy 数组是固定大小的,一般我们不推荐频繁地往数组中增加或者删除元素。但如果确实有这种需求,我们可以使用下面的方法来实现增加或者删除:

函数 使用说明
append values添加到arr的末尾。
insert obj(可以是下标、slicing)指定的位置处,插入数值value(可以是标量,也可以是数组)
delete 删除指定下标处的元素

示例如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
arr = np.arange(6).reshape((3,2))
np.append(arr, [[7,8]], axis=0)
cprint("指定在行的方向上操作、n{}", arr)

arr = np.arange(6).reshape((3,2))
arr = np.insert(arr.reshape((3,2)), 1, -10)
cprint("不指定 axis,数组被扁平化:\n{}", arr)

arr = np.arange(6).reshape((3,2))
arr = np.insert(arr, 1, (-10, -10), axis=0)
cprint("np.insert:\n{}", arr)

arr = np.delete(arr, [1], axis=1)
cprint("deleting col 1:\n{}", arr)

Tip

请一定运行一下这里的代码,特别是关于 insert 的部分,了解所谓的扁平化是怎么回事。

有时候我们需要修改个别元素的值,应该这样操作:

1
2
3
arr = np.arange(6).reshape(2,3)

arr[0,2] = 3

这里涉及到如何定位一个数组元素的问题,也正是我们下一节的内容。

1.3. 定位、读取和搜索

1.3.1. 索引和切片

Numpy 中的索引和切片语法大致类似于 Python,主要区别在于对多维数组的支持:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
arr = np.arange(6).reshape((3,2))
cprint("原始数组:\n{}", arr)

# 切片语法
cprint("按行切片:{}", arr[1, :])
cprint("按列切片:{}", arr[:, -1])
cprint("逆排数组:\n {}", arr[: : -1])

# FANCY INDEXING
cprint("fancy index: 使用下标数组:\n {}", arr[[2, 1, 0]])

上述切片语法在 Python 中也存在,但只能支持到一维,因此,对下面的 Python 数组,类似操作会出错:


1
2
3
arr = np.arange(6).reshape((3,2)).tolist()

arr[1, :]

提示 list indices must be integers or slices, not tuple。

1.3.2. 查找、筛选和替换

在上一节中,我们是通过索引来定位一个数组元素。但很多时候,我们得先通过条件运算把符合要求的索引找出来。这一节将介绍相关方法。

函数 使用说明
np.searchsorted 在有序数组中搜索指定的数值,返回索引。
np.nonzero 返回非零元素的索引,用以查找数组中满足条件的元素。
np.flatnonzero 同 nonzero,但返回输入数组的展平版本中非零的索引。
np.argwere 返回满足条件的元素的索引,相当于 nonzero 的转置版本
np.argmin 返回数组中最小元素的索引(注意不是返回满足条件的最小索引)
np.argmax 返回数组中最大元素的索引
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# 查找
arr = [0, 2, 2, 2, 3]
pos = np.searchsorted(arr, 2, 'right')
cprint("在数组 {} 中寻找等于 2 的位置,返回 {}, 数值是 {}", 
        arr, pos, arr[pos - 1])

arr = np.arange(6).reshape((2, 3))
cprint("arr[arr > 1]: {}", arr[arr > 1])

# NONZERO 的用法
mask = np.nonzero(arr > 1)
cprint("nonzero 返回结果是:{}", mask)
cprint("筛选后的数组是:{}", arr[mask])

# ARGWHERE 的用法
mask = np.argwhere(arr > 1)
cprint("argwere 返回的结果是:{}", mask)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# 多维数组不能直接使用 ARGWHERE 结果来筛选
# 下面的语句不能得到正确结果,一般会出现 INDEXERROR
arr[mask]

# 但对一维数组筛选我们可以用:
arr = np.arange(6)
mask = np.argwhere(arr > 1)
arr[mask.flatten()[0]]

# 寻找最大值的索引
arr = [1, 2, 2, 1, 0]
cprint("最大值索引是:{}", np.argmax(arr))

使用 searchsorted 要注意,数组本身一定是有序的,不然不会得出正确结果。

第 10 行到第 21 行代码,显示了如何查找一个数组中符合条件的数据,并且返回它的索引。

argwhere 返回值相当于 nonzero 的转置,在多维数组的情况下,它不能直接用作数组的索引。请自行对比 nonzero 与 argwhere 的用法。

在量化中,有很多情况需要实现筛选功能。比如,在计算上下影线时,我们是用公式\((high - max(open, close))/(high - low)\)来进行计算的。如果我们要一次性地计算过去 n 个周期的所有上影线,并且不使用循环的话,那么我们就要使用 np.where, np.select 等筛选功能。

下面的例子显示了如何使用 np.select 来计算上影线:

1
2
3
4
5
6
7
8
9
import pandas as pd
import numpy as np

bars = pd.DataFrame({
    "open": [10, 10.2, 10.1],
    "high": [11, 10.5, 9.3],
    "low": [9.8, 9.8, 9.25],
    "close": [10.1, 10.2, 10.05]
})

1
2
3
4
5
6
7
max_oc = np.select([bars.close > bars.open, 
                    bars.close <= bars.open], 
                    [bars.close, bars.open])
print(max_oc)

shadow = (bars.high - max_oc)/(bars.high - bars.low)
print(shadow)

np.where 是与 np.select 相近的一个函数,不过它只接受一个条件。

1
2
arr = np.arange(6)
cprint("np.where: {}", np.where(arr > 3, 3, arr))

这段代码实现了将 3 以上的数字截断为 3 的功能。这种功能被称为 clip,在因子预处理中是非常常用的一个技巧,用来处理异常值 (outlier)。

但它没有办法实现两端截断。此时,但 np.select 能做到,这是 np.where 与 np.select 的主要区别:

1
2
arr = np.arange(6)
cprint("np.select: {}", np.select([arr<2, arr>4], [2, 4], arr))
其结果是,生成的数组,小于 2 的被替换成 2,大于 4 的被替换成 4,其它的保持不变。

1.4. 审视 (inspecting) 数组

当我们调用其它人的库时,往往需要与它们交换数据。这时就可能出现数据格式不兼容的问题。为了有能力进行查错,我们必须掌握查看 Numpy 数组特性的一些方法。

我们先如下生成一个简单的数组,再查看它的各种特性:

1
2
3
4
arr = np.ones((3,2))
cprint("dtype is: {}", arr.dtype)
cprint("shape is: {}", arr.shape)
cprint("ndim is: {}", arr.ndim)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
cprint("size is: {}", arr.size)
cprint("'len' is also available: {}", len(arr))

# DTYPE
dt = np.dtype('>i4')
cprint("byteorder is: {}", dt.byteorder)
cprint("name of the type is: {}", dt.name)
cprint('is ">i4" a np.int32?: {}', dt.type is np.int32)

# 复杂的 DTYPE
complex = np.dtype([('name', 'U8'), ('score', 'f4')])
arr = np.array([('Aaron', 85), ('Zoe', 90)], dtype=complex)
cprint("A structured Array: {}", arr)
cprint("Dtype of structured array: {}", arr.dtype)

正如 Python 对象都有自己的数据类型一样,Numpy 数组也有自己的数据类型。我们可以通过arr.dtype来查看数组的数据类型。

从第 3 行到第 6 行,我们分别输出了数组的 shape, ndim, size 和 len 等属性。ndim 告诉我们数组的维度。shape 告诉我们每个维度的 size 是多少。shape 本身是一个 tuple, 这个 tuple 的 size,也等于 ndim。

size 在不带参数时,返回的是 shape 各元素取值的乘积。len 返回的是第一维的长度。

1.5. 数组操作

我们在前面的例子中,已经看到过一些引起数组形状改变的例子。比如,要生成一个\(3×2\)的数组,我们先用 np.arange(6) 来生成一个一维数组,再将它的形状改变为 (2, 3)。

另一个例子是使用 np.concatenate,从而改变了数组的行或者列。

1.5.1. 升维

我们可以通过 reshape, hstack, vstack 来改变数组的维度:


1
2
3
4
5
6
7
8
9
cprint("increase ndim with reshape:\n{}", 
        np.arange(6).reshape((3,2)))

# 将两个一维数组,堆叠为 2*3 的二维数组
cprint("createing from stack: {}", 
        np.vstack((np.arange(3), np.arange(4,7))))

# 将两个 (3,1)数组,堆叠为(3,2)数组
np.hstack((np.array([[1],[2],[3]]), np.array([[4], [5], [6]])))
1.5.2. 降维

通过 ravel, flatten, reshape, *split 等操作对数组进行降维。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
cprint("ravel: {}", arr.ravel())

cprint("flatten: {}", arr.flatten())

# RESHAPE 也可以用做扁平化
cprint("flatten by reshape: {}", arr.reshape(-1,))

# 使用 HSPLIT, VSPLIT 进行降维
x = np.arange(6).reshape((3, 2))
cprint("split:\n{}", np.hsplit(x, 2))

# RAVEL 与 FLATTEN 的区别:RAVEL 可以操作 PYTHON 的 LIST
np.ravel([[1,2,3],[4, 5, 6]])

这里一共介绍了 4 种方法。ravel 与 flatten 用法比较接近。ravel 的行为与 flatten 类似,只不过 ravel 是 np 的一个函数,可作用于 ArrayLike 的数组。

通过 reshape 来进行扁平化也是常用操作。此外,还介绍了 vsplit, hsplit 函数,它们的作用刚好与 vstack,hstack 相反。

1.5.3. 转置

此外,对数组进行转置也是此类例子中的一个。


比如,在前面我们提到,np.argwhere 的结果,实际上是 np.nonzero 的转置,我们来验证一下:

1
2
3
4
5
x = np.arange(6).reshape(2,3)
cprint("argwhere: {}", np.argwhere(x > 1))

# 我们再来看 NP.NONZERO 的转置
cprint("nonzero: {}", np.array(np.nonzero(x > 1)).T)

两次输出结果完全一样。在这里,我们是通过.T来实现的转置,它是一个语法糖,正式的函数是transpose

当然,由于 reshape 函数极其强大,我们也可以使用它来完成转置:

1
2
3
cprint("transposing array from \n{} to \n{}", 
    np.arange(6).reshape((2,3)),
    np.arange(6).reshape((3,2)))

快速傅里叶变换与股价预测研究

一个不证自明的事实:经济活动是有周期的。但是,这个事实似乎长久以来被量化界所忽略。无论是在资产定价理论,还是在趋势交易理论中我们几乎都找不到周期研究的位置 -- 在后者语境中,大家宁可使用“摆动”这样的术语,也不愿说破“周期”这个概念。

这篇文章里,我们就来探索股市中的周期。我们将运用快速傅里叶变换把时序信号分解成频域信号,通过信号的能量来识别主力资金,并且根据它们的操作周期进行一些预测。最后,我们给出三个猜想,其中一个已经得到了证明。

FFT - 时频互转

(取数据部分略)。

我们已经取得了过去一年来的沪指。显然,它是一个时间序列信号。傅里叶变换正好可以将时间序列信号转换为频域信号。换句话说,傅里叶变换能将沪指分解成若干个正弦波的组合。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
# 应用傅里叶变换
fft_result = np.fft.fft(close)
freqs = np.fft.fftfreq(len(close))

# 逆傅里叶变换
filtered = fft_result.copy()
filtered[20:] = 0
inverse_fft = np.fft.ifft(filtered)

# 绘制原始信号和分解后的信号
plt.figure(figsize=(14, 7))
plt.plot(close, label='Original Close')
plt.plot(np.real(inverse_fft), label='Reconstructed from Sine Waves')
plt.legend()

我们得到的输出如下:

在数字信号处理的领域,时间序列被称为时域信号,经过傅里叶变换后,我们得到的是频域信号。时域信号与频域信号可以相互转换。Numpy 中的 fft 库提供了 fft 和 ifft 这两个函数帮我们实现这两种转换。

np.fft.fft 将时域信号变换为频域信号,转换的结果是一个复数数组,代表信号分解出的各个频率的振幅 -- 也就是能量。频率由低到高排列,其中第 0 号元素的频率为 0,是直流分量,它是信号的均值的某个线性函数。

np.ff.ifft 则是 fft 的逆变换,将频域信号变换为时域信号。

将时域信号变换到频域,就能揭示出信号的周期性等基本特征。我们也可以对 fft 变换出来的频域信号进行一些操作之后,再变换回去,这就是数字信号处理。

高频滤波和压缩

如果我们把高频信号的能量置为零,再将信号逆变换回去,我们就会得到一个与原始序列相似的新序列,但它更平滑 -- 这就是我们常常所说的低通滤波的含义 -- 你熟悉的各种移动平均也都是低通滤波器。

在上面的代码中,我们只保留了前 20 个低频信号的能量,就得到了与原始序列相似的一个新序列。如果把这种方法运用在图像领域,这就实现了有损压缩 -- 压缩比是 250/20。

在上世纪 90 年代,最领先的图像压缩算法都是基于这样的原理 -- 保留图像的中低频部分,把高频部分当成噪声去掉,这样既保留了图像的主要特征,又大大减少了要保存的数据量。

当时做此类压缩算法的人都认识这位漂亮的小姐姐 -- Lena,这张照片是图像算法的标准测试样本。在漫长的进化中,出于生存的压力,人类在识别他人表情方面进化出超强的能力。所以相对于其它样本,一旦压缩造成图像质量下降,肉眼更容易检测到人脸和表情上发生的变化,于是人脸图像就成了最好的测试样本。

Lena 小姐姐是花花公子杂志的模特,这张照片是她为花花公子 1972 年 11 月那一期拍摄的诱人照片的一小部分 -- 在原始的照片中,Lena 大胆展现了她诱人的臀部曲线,但那些不正经的科学家们只与我们分享了她的微笑 -- 从科研的角度来讲,这也是信息比率最高的部分。无独有偶,在 Lena 成为数字图像处理的标准测试样本之前,科学家们一直使用的是另一位小姐姐的照片,也出自花花公子。

好,言归正传。我们刚刚分享了一种方法,去掉信号中的高频噪声,使得信号本身的意义更加突显出来。我们也希望在证券分析中使用类似的技法,使得隐藏在 K 线中的信号显露出来。

但如果我们照搬其它领域这一方法,这几乎就不算研究,也很难获得好的结果。实际上,在证券信号中,与频率相比,我们更应该关注信号的能量,毕竟,我们要与最有力量的人站在一边。

所以,我们换一个思路,把分解后的频域信号中,能量最强的部分保留下来,看看它们长什么样。

过滤低能量信号

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# 保留能量最强的前 5 个信号
amp_threshold = np.sort(np.abs(fft_result))[-11]

# 绘制各个正弦波分量
plt.figure(figsize=(14, 7))

theforce = []
for freq in freqs:
    if freq == 0:  # 处理直流分量
        continue
    elif freq < 0:
        continue
    else:
        amp = np.abs(fft_result[np.where(freqs == freq)])
        if amp < amp_threshold:
            continue
        sine_wave = amp * np.sin(2 * np.pi * freq * np.arange(len(close)))
        theforce.append(sine_wave)
        plt.plot(dates, sine_wave, label=f'Frequency={freq:.2f}')

plt.legend()
plt.title('Individual Sine Wave Components')
ticks = np.arange(0, len(dates), 20)
labels_to_show = [dates[i] for i in ticks]
plt.xticks(ticks=ticks, labels=labels_to_show, rotation=45)
plt.show()

FFT 给出的频率总是一正一负,我们可以简单地认为,负的频率对我们没有意义,那是一种我们看不到、也无须关心的暗能量。所以,在代码中,我们就忽略了这一部分。

我们看到,对沪指走势影响最强的波(橙色)的周期是 7 个月左右:从峰到底要走 3 个半月,从底到峰也要走 3 个半月。由于它的能量几乎是其它波的一倍,所以它是主导整个叠加波的走势的:如果其它波与它同相,叠加的结果就会使得趋势加强;反之,则是趋势抵消。其它几个波的能量相仿,但频率不同。

这些波倒底是什么呢?它可以是一种经济周期,但是说到底,经济周期是人推动的,或者反应了人的判断。因此,我们可以把波动的周期,看成资金的操作周期

从这个分解图中,我们可以猜想,有一个长线资金(对应蓝色波段),它一年多调仓一次。有一个中线资金(对应橙色波段),它半年左右调一次仓。其它的资金则是短线资金,三个月左右就会做一次仓位变更。还有无数被我们过滤掉的高频波段,它们操作频繁,可能对应着散户,但是能量很小,一般都可以忽略;只有在极个别的时候,才能形成同方向的叠加,进而影响到走势。

现在,我们把这几路资金的操作合成起来,并与真实的走势进行对比,看看情况如何:

在大的周期上都基本吻合,也就是这些资金基本上左右了市场的走势。而且,我们还似乎可以断言,在 3 月 15 到 5 月 17 这段时间,出现了股价与主力资金的背离趋势:主力资金在撤退了,但散户还在操作,于是,尽管股价还在上涨,但最终的方向,由主力资金来决定。

Tip

黑色线是通过主力资金波段合成出来的(对未来有预测性),在市场没有发生根本性变化之前,主力的操作风格是相对固定的,因此,它可能有一定的短时预测能力。如果我们认可这个结论的话。那么就应该注意到,末端部分还存在另一个背离 -- 散户还在离场,但主力已经进场。当然,关于这一点,请千万不要太当真。

关于直流分量的解释

我过去一直以为直流分量表明资产价格的趋势,但实际上所有的波都是水平走向的 -- 但只有商品市场才是水平走向的,股票市场本质上是向上的。所以,直流分量无法表明资产价格的趋势。

直到今天我突然涌现一个想法:如果你把一个较长的时序信号分段进行 FFT 分解,这样会得到若干个直流分量。这些直流分量的回归线才是资产价格的趋势。

这里给出三个猜想:

  1. 如果分段分解后,各个频率上的能量分布没有显著变化,说明投资者的构成及操作风格也没有显著变化,我们可以用 FFT 来预测未来短期走势,直到条件不再满足为止。
  2. 沪指 30 年来直流分量应该可以比较完美地拟合成趋势线,它的斜率就等于沪指 20 年回归线的斜率。
  3. 证券价格是直流分量趋势线与一系列正弦波的组合。

下面我们来证明第二个猜想(过程略)。最终,我们将直流分量及趋势线绘制成下图:

而 2005 年以来的 A 股年线及趋势线是这样的:

不能说十分相似,只能说几乎完全一致。

趋势线拟合的 p 值是 0.055 左右,也基本满足 0.05 的置信度要求。

本文复现代码已上传到我们的课程环境,加入我们的付费投研圈子即可获得。

这篇文章是我们《因子分析与机器学习策略》中的内容,出现在如何探索新的因子方法论那一部分。对 FFT 变换得到的一些重要结果,将成为机器学习策略中用以训练的特征。更多内容,我们课堂上见!

[0825] QuanTide Weekly

本周要闻

  • 美联储主席鲍威尔表示,美联储降息时机已经到来
  • 摩根大通港股仓位近日大量转仓,涉及市值超1.1万亿港元

下周看点

  • 广发明星基金经理刘格菘的首只“三年持有基”即将到期,亏损超58%
  • 周四A50指数交割日、周五本月收官日
  • 周六发布8月官方制造PMI

本周精选

  • 如何实现Alpha 101?
  • 高效量化编程 - Mask Array and Find Runs
  • 样本测试之外,我们还有哪些过拟合检测方法?

要闻详情

  • 美联储主席鲍威尔表示,通货膨胀率仅比美联储2%的目标高出半个百分点,失业率也在上升,“政策调整的时机已经到来”。财联社
  • 摩根大通港股仓位近日大量转仓,涉及市值超1.1万亿港元。转仓后,券端持股市值排名由第4名下跌至14名,持股市值不足2000亿港元。1个月前,摩根大通亦有超6000亿元转仓。金融界
  • 广发基金明星基金经理刘格菘的首只“三年持有基”即将到期,初期募资148.70亿元,截至今年8月22日,该基金(A/C)成立以来亏损超58%。近期,三年持有期基金集中到期。回溯来看,三年持有期主动权益基金在2021年——2022年间,公募基金行业密集推出了至少73只三年持有期主动权益基金。打开封闭之后,基民会不会巨量赎回?这是下周最重要的波动因素之一。相信有关方面已经做好了准备。新浪财经
  • 8月25日,北京商报发表《外资今天对A投爱答不理,明天就让他们高攀不起》一周年。在一年前的这篇评论中,北京商报指出,在目前股票具有极高投资价值的阶段,有一些外资流出A股,可能就是他们所谓的技术派典型代表,对指数患得患失,但他们最终一定会后悔,等想再回来的时候,势必要支付更高的价格,正所谓今天对A股爱答不理,明天就让他们高攀不起。该评论发布次日,沪指开盘于3219点。一年之后,沪指收盘于2854点。

如何实现Alpha 101?

2015 年,World Quant 发布了报告 《101 Formulaic Alphas》,它包含了 101 种不同的股票选择因子,这些因子中,有 80%是当时正在 World Quant 交易中使用的因子。该报告发表之后,在产业界引起较大反响。

目前,根据 Alpha101 生成的因子库,已几乎成为各量化平台、数据提供商和量化机构的必备。此外,一些机构受此启发,还在此基础上构建了更多因子,比如国泰君安推出的 Alpha 191 等。这两个因子库都有机构进行了实现。比如 DolphinDB聚宽 都提供了这两个因子库。

这篇文章就来介绍如何读懂 《101 Formulaic Alphas》 并且实现它。文章内容摘自我们的课程《因子分析与机器学习策略》的第8课,篇幅所限,有删节。

Alpha 101 因子中的数据和算子

在实现 Alpha101 因子之前,我们首先要理解其公式中使用的数据和基础算子。

Alpha101 因子主要是基于价格和成交量构建,只有少部分 Alpha 中使用了基本面数据,包括市值数据和行业分类数据 [^fundmental_data]。

Tip

在 A 股市场,由于财报数据的可信度问题 [^fraut],由于缺乏 T+0 和卖空等交易机制,短期内由交易行为产生的价格失效现象非常常见。因此,短期量价因子在现阶段的有效性高于基本面因子。


在价量数据中,Alpha101 依赖的最原始数据是 OHLC, volume(成交额), amount(成交量),turnover(换手率),并在此基础上,计算出来 returns(每日涨跌幅)和 vwap(加权平均成交价格)。

returns 和 vwap 的计算方法如下:

1
2
3
4
5
6
7
# THE NUMPY WAY
vwap = bars["amount"] / bars["volume"] / 100
returns = bars["close"][1:]/bars["close"][:-1] - 1

# THE PANDAS WAY
vwap = df.amount / df.volume / 100
returns = df.close.pct_change()

除此之外,要理解 Alpha101,重要的是理解它的公用算子。在 Alpha101 中,总共有约 30 个左右的算子,其中有一些像 abs, log, sign, min, max 以及数学运算符(+, -, *, /)等都是无须解释的。

下面,我们就先逐一解释需要说明的算子。

三目运算符

三目运算符是 Python 中没有,但存在于 C 编程语言的一个算子。这个运算符可以表示为:"x ? y : z",它相当于 Python 中的:

1
2
3
4
5
6
expr_result = None

if x:
    expr_result = y
else:
    expr_result = z

rank

在 Alpha101 中,存在两个 rank,一个是横截面上的,即对同一时间点上 universe 中所有的股票进行排序;另一个是时间序列上的,即对同一股票在时间序列上的排序。


在横截面上的 rank 直接调用 DataFrame 的 rank。比如,

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import pandas as pd

data = {
    'asset': ["000001", "000002", "000004", "000005", "000006"],
    'factor': [85, 92, 78, 92, 88],
    'date': [0] * 5
}
df = pd.DataFrame(data).set_index('date').pivot(index=None, columns="asset", values="factor")

def rank(df):
    return df.rank(axis=1, pct=True, method='min')

在上面这段代码中,date 为索引,列名字为各 asset,factor 为其值,此时,我们就可以通 rank(axis=1) 的方法,对各 asset 的 factor 值在截面上进行排序。当我们使用 axis=1 参数时,索引是不参与排序。pct 为 True 表示返回的是百分比排名,False 表示返回的是排名。

有时候我们也需要在时间序列上进行排序,在 Alpha101 中,这种排序被定义为 ts_rank,通过前缀 ts_来与截面上的 rank 相区分。此后,当我们看到 ts_前缀时,也应该作同样理解。

1
2
3
4
from bottleneck import move_rank

def ts_rank(df, window=10, min_count=1):
    return move_rank(df, window, axis=0, min_count=min_count)

在这里我们使用的是 bottleneck 中的 move_rank,它的速度要显著高于 pandas 和 scipy 中的同类实现。如果使用 pandas 来实现,代码如下:

1
2
3
4
5
def rolling_rank(na):
    return rankdata(na,method='min')[-1]

def ts_rank(df, window=10):
    return df.rolling(window).apply(rolling_rank)

注意第 3 行中的 [-1] 是必须的。


rank 和 ts_rank 的使用在 alpha004 因子中的应用最为典型。这个因子是:

1
2
3
# ALPHA#4    (-1 * TS_RANK(RANK(LOW), 9))
def alpha004(low):
    return -1 * ts_rank(rank(low), 9)

在这里,参数 low 是一个以 asset 为列、日期为索引,当日最低价为值的 dataframe,是一个宽表。下面,我们看一下对参数 low 依次调用 rank 和 ts_rank 的结果。通过深入几个示例之后,我们就很快能够明白 Alpha101 的因子计算过程。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from bottleneck import move_rank
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

df = pd.DataFrame(
       [[6.18, 19.36, 33.13, 14.23,  6.8 ,  6.34],
       [6.55, 20.36, 32.69, 14.52,  6.4,  6.44 ],
       [7.  , 20.79, 32.51, 14.56,  6.0 ,  6.54],
       [7.06, 21.35, 33.13, 14.47,  6.5,  6.64],
       [7.03, 21.56, 33.6 , 14.6 ,  6.5,  6.44]], 
       columns=['000001', '000002', '000063', '000066', '000069', '000100'], 
       index=['2022-01-04', '2022-01-05', '2022-01-06', '2022-01-07', '2022-01-10'])

def rank(df):
    return df.rank(axis=1, pct=True, method='min')

def ts_rank(df, window=10, min_count=1):
    return move_rank(df, window, axis=0, min_count=min_count)

df
rank(df)

-1 * ts_rank(rank(df), 3)

示例 将依次输出三个 DataFrame。我们看到,rank 是执行在行上,它将各股票按最低价进行排序;ts_rank 执行在列上,对各股票在横截面上的排序再进行排序,反应出最低位置的变化。

比如,000100 这支股票,在 2022 年 1 月 4 日,横截面的排位在 33%分位,到了 1 月 10 日,它在横截面上的排位下降到 16.7%。


通过 ts_rank 之后,最终它在 1 月 10 日的因子值为 1,反应了它在横截面上排位下降的事实。同理,000001 这支股票,在 1 月 4 日,它的横截面上的排位是 16.7%(最低),而在 1 月 5 日,它的排序上升到 50%,最终它在当日的因子值为-1,反应了它在横截面排序上升的事实。

Tip

通过 Alpha004 因子,我们不仅了解到 rank 与 ts_rank 的用法,也知道了横截面算子与时序算子的区别。此外,我们也了解到,为了方便计算 alpha101 因子,最佳的数据组织方式可能是将基础数据(比如 OHLC)都组织成一个个以日期为索引、asset 为列的宽表,以方便在两个方向上(横截面和时序)的计算。

ts_*

这一组算子中,除了之前已经介绍过的 ts_rank 之外,还有 ts_max, ts_argmax, ts_argmin, ts_min。这一些算子都有两个参数,首先时时间序列,比如 close 或者 open,然后是滑动窗口的长度。

注意这一类算子一定是在滑动窗口上进行的,只有这样,才不会引入未来数据。

除此之外,其它常用统计函数,比如 min, max, sum, product, stddev 等,尽管没有使用 ts_前缀,它们也是时序算子,而不是截面算子。考虑到我们已经通过 ts_rank 详细介绍了时序算子的用法,而这些算子的作用大家也都非常熟悉,这里就从略。

delay

在 Alpha101 中,delay 算子用来获取 n 天前的数据。比如,


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def delay(df, n):
    return df.shift(n)

data = {
    'date': pd.date_range(start='2023-01-01', periods=10),
    'close': [100, 101, 102, 103, 104, 105, 106, 107, 108, 109]
}
df = pd.DataFrame(data)

delay(df, 5)

如此一来,我们在计算第 5 天的因子时,使用的 close 数据就是 5 天前的,即原来索引为 0 处的 close。

correlation 和 covariance

correlation 就是两个时间序列在滑动窗口上的皮尔逊相关系数,这个算子可以实现为:

1
2
3
4
5
def correlation(x, y, window=10):
    return x.rolling(window).corr(y).fillna(0).replace([np.inf, -np.inf], 0)

def covariance(x, y, window=10):
    return x.rolling(window).cov(y)

注意在这里,尽管我们只对 x 调用了 rolling,但在计算相关系数时,经验证,y 也是按同样的窗口进行滑动的。

scale

按照 Alpha101 的解释,这个算子的作用,是将数组的元素进行缩放,使之满足 sum(abs(x)) = a,缺省情况下 a = 1。它可以实现为:

1
2
def scale(df, k=1):
    return df.mul(k).div(np.abs(df).sum())

decay_linear

这个算子的作用是将长度为 d 的时间序列中的元素进行线性加权衰减,使之总和为 1,且越往后的元素权重越大。

1
2
3
4
def decay_linear(df, period=10):
    weights = np.array(range(1, period+1))
    sum_weights = np.sum(weights)
    return df.rolling(period).apply(lambda x: np.sum(weights*x) / sum_weights)

delta

相当于 dataframe.diff()。

adv{d}

成交量的 d 天简单移动平均。

signedpower

signedpower(x, a) 相当于 x^a

Alpha 101 因子解读

此部分略

开源的Alpha101因子分析库

完整探索Alpha101中的定义的因子的最佳方案是,根据历史数据,计算出所有这些因子,并且通过Alphalens甚至backtrader对它们进行回测。popbo就实现了这样的功能。


运行该程序库需要安装alphalens, akshare,baostock以及jupyternotebook。在进行研究之前,需要先参照其README文件进行数据下载和因子计算。然后就可以打开research.ipynb,对每个因子的历年表现进行分析。

在我们的补充材料中,提供了该项目的全部源码并且可以在我们的课程环境中运行。


高效量化编程 - Mask Array and Find Runs

在很多量化场景下,我们都需要统计某个事件连续发生了多少次,比如,连续涨跌停、N连阳、计算Connor's RSI中的streaks等等。

比如,要判断下列收盘价中,最大的连续涨停次数是多少?最长的N连涨数是多少?应该如何计算呢?

1
2
a = [15.28, 16.81, 18.49, 20.34, 21.2, 20.5, 22.37, 24.61, 27.07, 29.78, 
     32.76, 36.04]

假设我们以10%的涨幅为限,则可以将上述数组转换为:

1
2
pct = np.diff(a) / a[:-1]
pct > 0.1

我们将得到以下数组:

1
flags = [True, False, True, False, False, False, True, False, True, True, True]

这仍然不能计算出最大连续涨停次数,但它是很多此类问题的一个基本数据结构,我们将原始的数据按条件转换成类似的数组之后,就可以使用下面的神器了:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
from numpy.typing import ArrayLike
from typing import Tuple
import numpy as np

def find_runs(x: ArrayLike) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Find runs of consecutive items in an array.
    """

    # ensure array
    x = np.asanyarray(x)
    if x.ndim != 1:
        raise ValueError("only 1D array supported")
    n = x.shape[0]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
    # handle empty array
    if n == 0:
        return np.array([]), np.array([]), np.array([])

    else:
        # find run starts
        loc_run_start = np.empty(n, dtype=bool)
        loc_run_start[0] = True
        np.not_equal(x[:-1], x[1:], out=loc_run_start[1:])
        run_starts = np.nonzero(loc_run_start)[0]

        # find run values
        run_values = x[loc_run_start]

        # find run lengths
        run_lengths = np.diff(np.append(run_starts, n))

        return run_values, run_starts, run_lengths


pct = np.diff(a) / a[:-1]
v,s,l = find_runs(pct > 0.099)
(v, s, l)

输出结果为:

1
(array([ True, False,  True]), array([0, 3, 6]), array([3, 3, 5]))

输出结果是一个由三个数组组成的元组,分别表示:

value: unique values start: start indices length: length of runs 在上面的输出中,v[0]为True,表示这是一系列涨停的开始,s[0]则是对应的起始位置,此时索引为0; l[0]则表示该连续的涨停次数为3次。同样,我们可以知道,原始数组中,最长连续涨停(v[2])次数为5(l[2]),从索引6(s[2])开始起。

所以,要找出原始序列中的最大连续涨停次数,只需要找到l中的最大值即可。但要解决这个问题依然有一点技巧,我们需要使用第4章中介绍的 mask array。

1
2
3
v_ma = np.ma.array(v, mask = ~v)
pos = np.argmax(v_ma * l)
print(f"最大连续涨停次数{l[pos]},从索引{s[pos]}:{a[s[pos]]}开始。")

在这里,mask array的作用是,既不让 v == False 的数据参与计算(后面的 v_ma * l),又保留这些元素的次序(索引)不变,以便后面我们调用 argmax 函数时,找到的索引跟v, s, l中的对应位置是一致的。

我们创建的v_ma是一个mask array,它的值为:

1
2
3
masked_array(data=[True, --, True],
            mask=[False,  True, False],
            fill_value=True)
当它与另一个整数数组相乘时,True就转化为数字1,这样相乘的结果也仍然是一个mask array:

1
2
3
masked_array(data=[3, --, 5],
             mask=[False,  True, False],
            fill_value=True)

当arg_max作用在mask array时,它会忽略掉mask为True的元素,但保留它们的位置,因此,最终pos的结果为2,对应的 v,s,l中的元素值分别为: True, 6, 5。

如果要统计最长N连涨呢?这是一个比寻找涨停更容易的任务。不过,这一次,我们将不使用mask array来实现:

1
2
3
4
v,s,l = find_runs(np.diff(a) > 0)

pos = np.argmax(v * l)
print(f"最长N连涨次数{l[pos]},从索引{s[pos]}:{a[s[pos]]}开始。")

输出结果是:最长N连涨次数6,从索引5:20.5开始。

这里的关键是,当Numpy执行乘法时,True会被当成数字1,而False会被当成数字0,于是,乘法结果自然消除了没有连续上涨的部分,从而不干扰argmax的计算。

当然,使用mask array可能在语义上更清楚一些,尽管mask array的速度会慢一点,但正确和易懂常常更重要。


计算 Connor's RSI中的streaks Connor's RSI(Connor's Relative Strength Index)是一种技术分析指标,它是由Nirvana Systems开发的一种改进版的相对强弱指数(RSI)。

Connor's RSI与传统RSI的主要区别在于它考虑了价格连续上涨或下跌的天数,也就是所谓的“连胜”(winning streaks)和“连败”(losing streaks)。这种考虑使得Connor's RSI能够更好地反映市场趋势的强度。

在前面介绍了find_runs函数之后,计算streaks就变得非常简单了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def streaks(close):
    result = []
    conds = [close[1:]>close[:-1], close[1:]<close[:-1]]

    flags = np.select(conds, [1, -1], 0)
    v, _, l = find_runs(flags)
    for i in range(len(v)):
        if v[i] == 0:
            result.extend([0] * l[i])
        else:
            result.extend([v[i] * x for x in range(1, (l[i] + 1))])

    return np.insert(result, 0, 0)

这段代码首先将股价序列划分为上涨、下跌和平盘三个子系列,然后对每个子系列计算连续上涨或下跌的天数,并将结果合并成一个新的数组。

在streaks中,连续上涨天数要用正数表示,连续下跌天数用负数表示,所以在第5行中,通过np.select将条件数组转换为[1, 0, -1]的序列,后面使用乘法就能得到正确的连续上涨(下跌)天数了。


样本测试之外,我们还有哪些过拟合检测方法?

在知乎上看到一个搞笑的贴子,说是有人为了卖策略,让回测结果好看,会在代码中植入大量的if 语句,判断当前时间是特定的日期,就不进行交易。但奥妙全在这些日期里,因为在这些日期时,交易全是亏损的。

内容的真实性值得怀疑。不过,这却是一个典型的过拟合例子。

过拟合和检测方法

过拟合是指模型与数据拟合得很好,以至于该模型不可泛化,从而不能在另一个数据集上工作。从交易角度来说,过拟合“设计”了一种策略,可以很好地交易历史数据,但在新数据上肯定会失败。

过拟合是我们在回测中的头号敌人。如何检测过拟合呢?

一个显而易见的检测方法是样本外测试。它是把整个数据集划分为互不重叠的训练集和测试集,在训练集上训练模型,在测试集上进行验证。如果模型在测试集上也表现良好,就认为该模型没有拟合。

在样本本身就不足的情况下,样本外测试就变得困难。于是,人们发明了一些拓展版本。

其中一种拓展版本是 k-fold cross-validation,这是在机器学习中常见的概念。

它是将数据集随机分成 K 个大小大致相等的子集,对于每一轮验证,选择一个子集作为验证集,其余 K-1 个子集作为训练集。模型在训练集上训练,在验证集上进行评估。这个过程重复 K 次,最终评估指标通常为 K 次验证结果的平均值。


这个过程可以简单地用下图表示:

k-fold cross validation,by sklearn

但在时间序列分析(证券分析是其中典型的一种)中,k-fold方法是不适合的,因为时间序列分析有严格的顺序性。因此,从k-fold cross-validation特化出来一个版本,称为 rolling forecasting。你可以把它看成顺序版本的k-fold cross-validation。

它可以简单地用下图表示:

rolling forecasting, by tsfresh


从k-fold cross-validation到rolling forecasting的两张图可以看出,它们的区别在于一个是无序的,另一个则强调时间顺序,训练集和验证集之间必须是连续的。

有时候,你也会看到 Walk-Forward Optimization这种说法。它与rolling forecasting没有本质区别。

不过,我最近从buildalpha网站上,了解到了一种新颖的方法,这就是噪声测试。

新尝试:噪声测试

buildalpha的噪声测试,是将一定比率的随机噪声叠加到回测数据上,然后再进行回测,并将基于噪声的回测与基于真实数据的回测进行比较。

L50

它的原理是,在我们进行回测时,历史数据只是可能发生的一种可能路径。如果时间重演,历史可能不会改变总的方向,但是偶然性会改变历史的步伐。而一个好的策略,应该是能对抗偶然性、把握历史总的方向的策略。因此,在一个时间序列加上一些巧妙的噪声,就可能会让过拟合的策略失效,而真正有效的策略仍然闪耀。

buildalpha是一个类似tradingview的平台。要进行噪声测试,可以通过图形界面进行配置。

通过这个对话框,buildalpha修改了20%左右的数据,并且对OHLC的修改幅度都控制在用ATR的20%以内。最下面的100表明我们将随机生成100组带噪声的数据。


我们对比下真实数据与叠加噪声的数据。

左图为真实数据,右图为叠加部分噪声的数据。叠加噪声后,在一些细节上,引入了随机性,但并没有改变股价走势(叠加是独立的)。如果股价走势被改变,那么这种方法就是无效的甚至有害的。

最后,在同一个策略上,对照回测的结果是:

75%

从结果上看,在历史的多条可能路径中,没有任何一条的回测结果能比真实数据好。


换句话说,真实回测的结果之所以这么好,纯粹是因为制定策略的人,是带着上帝视角,从未来穿越回去的。

参数平原与噪声测试

噪声测试是稍稍修改历史数据再进行圆滑。而参数平原则是另一种检测过拟合的方法,它是指稍微修改策略参数,看回测表现是否会发生剧烈的改变。如果没有发生剧烈的改变,那么策略参数就是鲁棒的。

Build Alpha以可视化的方式,提供了参数平原检测。

在这个3D图中,参数选择为 X= 9和Y=4,如黑色简单所示。显然,这一区域靠近敏感区域,在其周围,策略的性能下降非常厉害。按照传统的推荐,我们应该选择参数 X=8和Y=8,这一区域图形更为平坦。

在很多时候,参数平原的提示是对的 -- 因为我们选择的参数,其实价格变化的函数;但它毕竟不是价格变化。最直接的方法是,当价格发生轻微变化时,策略的性能如果仍然处在一个平坦的表面,就更能说明策略是鲁棒的。

不过,这种图很难绘制,所以,Build Alpha绘制的仍然是以参数为n维空间的坐标、策略性能为其取值的三维图,但它不再是基于单个历史数据,而是基于一组历史数据:真实历史数据和增加了噪声的数据。

高效量化编程: Mask Array应用和find_runs

在很多量化场景下,我们都需要统计某个事件连续发生了多少次,比如,连续涨跌停、N连阳、计算Connor's RSI中的streaks等等。

比如,要判断下列收盘价中,最大的连续涨停次数是多少?最长的N连涨数是多少?应该如何计算呢?

1
2
a = [15.28, 16.81, 18.49, 20.34, 21.2, 20.5, 22.37, 24.61, 27.07, 29.78, 
     32.76, 36.04]

假设我们以10%的涨幅为限,则可以将上述数组转换为:

1
2
pct = np.diff(a) / a[:-1]
pct > 0.1

我们将得到以下数组:

1
flags = [True, False, True, False, False, False, True, False, True, True, True]

这仍然不能计算出最大连续涨停次数,但它是很多此类问题的一个基本数据结构,我们将原始的数据按条件转换成类似的数组之后,就可以使用下面的神器了:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
from numpy.typing import ArrayLike
from typing import Tuple
import numpy as np

def find_runs(x: ArrayLike) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Find runs of consecutive items in an array.

    Args:
        x: the sequence to find runs in

    Returns:
        A tuple of unique values, start indices, and length of runs
    """

    # ensure array
    x = np.asanyarray(x)
    if x.ndim != 1:
        raise ValueError("only 1D array supported")
    n = x.shape[0]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
    # handle empty array
    if n == 0:
        return np.array([]), np.array([]), np.array([])

    else:
        # find run starts
        loc_run_start = np.empty(n, dtype=bool)
        loc_run_start[0] = True
        np.not_equal(x[:-1], x[1:], out=loc_run_start[1:])
        run_starts = np.nonzero(loc_run_start)[0]

        # find run values
        run_values = x[loc_run_start]

        # find run lengths
        run_lengths = np.diff(np.append(run_starts, n))

        return run_values, run_starts, run_lengths


pct = np.diff(a) / a[:-1]
v,s,l = find_runs(pct > 0.099)
(v, s, l)

输出结果为:

1
(array([ True, False,  True]), array([0, 3, 6]), array([3, 3, 5]))

输出结果是一个由三个数组组成的元组,分别表示:

value: unique values start: start indices length: length of runs 在上面的输出中,v[0]为True,表示这是一系列涨停的开始,s[0]则是对应的起始位置,此时索引为0; l[0]则表示该连续的涨停次数为3次。同样,我们可以知道,原始数组中,最长连续涨停(v[2])次数为5(l[2]),从索引6(s[2])开始起。

所以,要找出原始序列中的最大连续涨停次数,只需要找到l中的最大值即可。但要解决这个问题依然有一点技巧,我们需要使用第4章中介绍的 mask array。

1
2
3
v_ma = np.ma.array(v, mask = ~v)
pos = np.argmax(v_ma * l)
print(f"最大连续涨停次数{l[pos]},从索引{s[pos]}:{a[s[pos]]}开始。")

在这里,mask array的作用是,既不让 v == False 的数据参与计算(后面的 v_ma * l),又保留这些元素的次序(索引)不变,以便后面我们调用 argmax 函数时,找到的索引跟v, s, l中的对应位置是一致的。

我们创建的v_ma是一个mask array,它的值为:

1
2
3
masked_array(data=[True, --, True],
            mask=[False,  True, False],
            fill_value=True)
当它与另一个整数数组相乘时,True就转化为数字1,这样相乘的结果也仍然是一个mask array:

1
2
3
masked_array(data=[3, --, 5],
             mask=[False,  True, False],
            fill_value=True)

当arg_max作用在mask array时,它会忽略掉mask为True的元素,但保留它们的位置,因此,最终pos的结果为2,对应的 v,s,l中的元素值分别为: True, 6, 5。

如果要统计最长N连涨呢?这是一个比寻找涨停更容易的任务。不过,这一次,我们将不使用mask array来实现:

1
2
3
4
v,s,l = find_runs(np.diff(a) > 0)

pos = np.argmax(v * l)
print(f"最长N连涨次数{l[pos]},从索引{s[pos]}:{a[s[pos]]}开始。")

输出结果是:最长N连涨次数6,从索引5:20.5开始。

这里的关键是,当Numpy执行乘法时,True会被当成数字1,而False会被当成数字0,于是,乘法结果自然消除了没有连续上涨的部分,从而不干扰argmax的计算。

当然,使用mask array可能在语义上更清楚一些,尽管mask array的速度会慢一点,但正确和易懂常常更重要。


计算 Connor's RSI中的streaks Connor's RSI(Connor's Relative Strength Index)是一种技术分析指标,它是由Nirvana Systems开发的一种改进版的相对强弱指数(RSI)。

Connor's RSI与传统RSI的主要区别在于它考虑了价格连续上涨或下跌的天数,也就是所谓的“连胜”(winning streaks)和“连败”(losing streaks)。这种考虑使得Connor's RSI能够更好地反映市场趋势的强度。

在前面介绍了find_runs函数之后,计算streaks就变得非常简单了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def streaks(close):
    result = []
    conds = [close[1:]>close[:-1], close[1:]<close[:-1]]

    flags = np.select(conds, [1, -1], 0)
    v, _, l = find_runs(flags)
    for i in range(len(v)):
        if v[i] == 0:
            result.extend([0] * l[i])
        else:
            result.extend([v[i] * x for x in range(1, (l[i] + 1))])

    return np.insert(result, 0, 0)

这段代码首先将股价序列划分为上涨、下跌和平盘三个子系列,然后对每个子系列计算连续上涨或下跌的天数,并将结果合并成一个新的数组。

在streaks中,连续上涨天数要用正数表示,连续下跌天数用负数表示,所以在第5行中,通过np.select将条件数组转换为[1, 0, -1]的序列,后面使用乘法就能得到正确的连续上涨(下跌)天数了。