《Python数据分析基础教程：NumPy学习指南（第2版）》笔记10：第四章 便捷函数2——

本章通过计算股票收益率相关性的案例演示NumPy数据分析。

第四章 便捷函数

4.5 净额成交量

成交量（volume）是投资中一个非常重要的变量，它可以表示价格波动的大小。 OBV（On-Balance Volume，净额成交量或叫能量潮指标）是最简单的股价指标之一，它可以由当日收盘价、前一天的收盘价以及当日成交量计算得出。这里我们以前一日为基期计算当日的OBV值（可以认为基期的OBV值为0）。若当日收盘价高于前一日收盘价，则本日OBV等于基期OBV加上当日成交量。若当日收盘价低于前一日收盘价，则本日OBV等于基期OBV减去当日成交量。若当日收盘价相比前一日没有变化，则当日成交量以0计算。

4.6 动手实践：计算OBV

换言之，我们需要在成交量前面乘上一个由收盘价变化决定的正负号。在本节教程中，我们将学习该问题的两种解决方法，一种是使用NumPy中的sign函数，另一种是使用NumPy的piecewise函数。

• (1) 把BHP数据分别加载到收盘价和成交量的数组中：

import numpy as np
c, v=np.loadtxt('BHP.csv', delimiter=',', usecols=(6, 7), unpack=True)

为了判断计算中成交量前的正负号，我们先使用diff函数计算收盘价的变化量。 diff函数可以计算数组中两个连续元素的差值，并返回一个由这些差值组成的数组：

change = np.diff(c)
print("Change", change)

收盘价差值的计算结果如下：

Change [ 1.92 -1.08 -1.26 0.63 -1.54 -0.28 0.25 -0.6 2.15 0.69 -1.33 1.16
1.59 -0.26 -1.29 -0.13 -2.12 -3.91 1.28 -0.57 -2.07 -2.07 2.5 1.18
-0.88 1.31 1.24 -0.59]

• (2) NumPy中的**sign函数可以返回数组中每个元素的正负符号**，数组元素为负时返回-1，为正时返回1，否则返回0。对change数组使用sign函数：

signs = np.sign(change)
print("Signs", signs)

change数组中各元素的正负符号如下所示：

Signs [ 1. -1. -1. 1. -1. -1. 1. -1. 1. 1. -1. 1. 1. -1. -1. -1. -1. -1. -1. -1. -1.
1. 1. 1. -1. 1. 1. -1.]

另外，我们也可以使用piecewise函数来获取数组元素的正负。顾名思义， piecewise函数可以分段给定取值。使用合适的返回值和对应的条件调用该函数：

pieces = np.piecewise(change, [change < 0, change > 0], [-1, 1])
print("Pieces", pieces)

再次输出数组元素的正负，结果如下：

Pieces [ 1. -1. -1. 1. -1. -1. 1. -1. 1. 1. -1. 1. 1. -1. -1. -1. -1. -1. -1. -1. -1.
2. 1. 1. -1. 1. 1. -1.]

检查两次的输出是否一致：

print("Arrays equal?", np.array_equal(signs, pieces))

结果如下：

Arrays equal? True

• (3) OBV值的计算依赖于前一日的收盘价，所以在我们的例子中无法计算首日的OBV值：

print("On balance volume", v[1:] * signs)

计算结果如下：

[2620800. -2461300. -3270900. 2650200. -4667300. -5359800. 7768400.
-4799100. 3448300. 4719800. -3898900. 3727700. 3379400. -2463900.
-3590900. -3805000. -3271700. -5507800. 2996800. -3434800. -5008300.
-7809799. 3947100. 3809700. 3098200. -3500200. 4285600. 3918800.
-3632200.]

小结

我们刚刚计算了OBV值，它依赖于收盘价的变化量。我们分别使用了NumPy中的sign函数和piecewise函数这两种不同的方法来判断收盘价变化量的正负。

示例完整代码如下：

import numpy as np

c, v=np.loadtxt('BHP.csv', delimiter=',', usecols=(6, 7), unpack=True)

change = np.diff(c)
print("Change", change)

signs = np.sign(change)
print("Signs", signs)

pieces = np.piecewise(change, [change < 0, change > 0], [-1, 1])
print("Pieces", pieces)

print("Arrays equal?", np.array_equal(signs, pieces))

print("On balance volume", v[1:] * signs)

4.7 交易过程模拟

你可能经常想尝试干一些事情，做一些实验，但又不希望造成任何不良后果。而NumPy就是用于实验的完美工具。我们将使用NumPy来模拟一个交易日，当然，这不会造成真正的资金损失。许多人喜欢抄底，也就是等股价下跌后才买入。类似的还有当股价比当日开盘价下跌一小部分（比如0.1%）时买入。

4.8 动手实践：避免使用循环

使用vectorize函数可以减少你的程序中使用循环的次数。我们将用它来计算单个交易日的利润。

• (1) 首先，读入数据：

import numpy as np
o, h, l, c = np.loadtxt('BHP.csv', delimiter=',', usecols=(3, 4, 5, 6), unpack=True)

• (2) NumPy中的vectorize函数相当于Python中的map函数。调用vectorize函数并给定calc_profit函数作为参数，尽管我们还没有编写这个函数：

func = np.vectorize(calc_profit)

• (3) 我们现在可以先把func当做函数来使用。对股价数组使用我们得到的func函数:

profits = func(o, h, l, c)

• (4) calc_profit函数非常简单。首先，我们尝试以比开盘价稍低一点的价格买入股票（取值为1.001）。如果这个价格不在当日的股价范围内，则尝试买入失败，没有获利，也没有亏损，我们均返回0。否则，我们将以当日收盘价卖出，所获得的利润即买入和卖出的差价。事实上，计算相对利润更为直观：

def calc_profit(open, high, low, close):
    #buy just below the open
    buy = open * float(1.001)

    # daily range
    if low <  buy < high:
        return (close - buy)/buy
    else:
        return 0

• (5) 在所有交易日中有两个零利润日，即没有利润也没有损失。我们选择非零利润的交易日
  并计算平均值：

real_trades = profits[profits != 0]
print("Number of trades", len(real_trades), round(100.0 * len(real_trades)/len)(c), 2), "%"
print("Average profit/loss %", round(np.mean(real_trades) * 100, 2))

交易结果如下：

Number of trades 28 93.33 %
Average profit/loss % -0.01

(6) 乐观的人们对于正盈利的交易更感兴趣。选择正盈利的交易日并计算平均利润：

winning_trades = profits[profits > 0]
print("Number of winning trades", len(winning_trades), round(100.0 * len(winning_trades)/len(c), 2), "%")
print("Average profit %", round(np.mean(winning_trades) * 100, 2))

正盈利交易的分析结果如下：

Number of winning trades 16 53.33 %
Average profit % 0.7

• (7) 悲观的人们对于负盈利的交易更感兴趣，选择负盈利的交易日并计算平均损失：

losing_trades = profits[profits < 0]
print("Number of losing trades", len(losing_trades), round(100.0 * len(losing_trades)/len(c), 2), "%")
print("Average loss %", round(np.mean(losing_trades) * 100, 2))

负盈利交易的分析结果如下：

Number of losing trades 12 40.0 %
Average loss % -0.96

小结

我们矢量化了一个函数，这是一种可以避免使用循环的技巧。我们使用一个能返回当日相对利润的函数来模拟一个交易日，并分别打印出正盈利和负盈利交易的概况。

示例完整代码如下：

import numpy as np

o, h, l, c = np.loadtxt('BHP.csv', delimiter=',', usecols=(3, 4, 5, 6), unpack=True)

def calc_profit(open, high, low, close):
    #buy just below the open
    buy = open * float(1.001)

    # daily range
    if low <  buy < high:
        return (close - buy)/buy
    else:
        return 0

func = np.vectorize(calc_profit)
profits = func(o, h, l, c)
print("Profits", profits)

real_trades = profits[profits != 0]
print("Number of trades", len(real_trades), round(100.0 * len(real_trades)/len(c), 2), "%")
print("Average profit/loss %", round(np.mean(real_trades) * 100, 2))

winning_trades = profits[profits > 0]
print("Number of winning trades", len(winning_trades), round(100.0 * len(winning_trades)/len(c), 2), "%")
print("Average profit %", round(np.mean(winning_trades) * 100, 2))

losing_trades = profits[profits < 0]
print("Number of losing trades", len(losing_trades), round(100.0 * len(losing_trades)/len(c), 2), "%")
print("Average loss %", round(np.mean(losing_trades) * 100, 2))

4.9 数据平滑

噪声数据往往很难处理，因此我们通常需要对其进行平滑处理。除了用计算移动平均线的方法，我们还可以使用NumPy中的一个函数来平滑数据。
**hanning函数是一个加权余弦的窗函数。**在后面的章节中，我们还将更为详细地介绍其他窗函数。

4.10 动手实践：使用 hanning 函数平滑数据

我们将使用hanning函数平滑股票收益率的数组，步骤如下。

• (1) 调用hanning函数计算权重，生成一个长度为N的窗口（在这个示例中N取8）：
  import numpy as np

N = 8
weights = np.hanning(N)
print(“Weights”, weights)
得到的权重如下：

Weights [ 0. 0.1882551 0.61126047 0.95048443 0.95048443 0.61126047 0.1882551 0. ]

• (2) 使用convolve函数计算BHP和VALE的股票收益率，以归一化处理后的weights作为参数：

bhp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp) / bhp[ : -1]
smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns)[N-1:-N+1]

vale = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
vale_returns = np.diff(vale) / vale[ : -1]
smooth_vale = np.convolve(weights/weights.sum(), vale_returns)[N-1:-N+1]

• (3) 用Matplotlib绘图：

t = np.arange(N - 1, len(bhp_returns))

plt.plot(t, bhp_returns[N-1:], lw=1.0)
plt.plot(t, smooth_bhp, lw=2.0)

plt.plot(t, vale_returns[N-1:], lw=1.0)
plt.plot(t, smooth_vale, lw=2.0)
plt.show()

绘制的折线图如下。

图中的细线为股票收益率，粗线为平滑处理后的结果。如你所见，图中的折线有交叉。这些交叉点很重要，因为它们可能就是股价趋势的转折点，至少可以表明BHP和VALE之间的股价关系发生了变化。这些转折点可能会经常出现，我们可以利用它们预测未来的股价走势。

• (4) 使用多项式拟合平滑后的数据（设K=N）：

K = 8
t = np.arange(N - 1, len(bhp_returns))
poly_bhp = np.polyfit(t, smooth_bhp, K)
poly_vale = np.polyfit(t, smooth_vale, K)

(5) 现在，我们需要解出上面的两个多项式何时取值相等，即在哪些地方存在交叉点。这等价于先对两个多项式函数作差，然后对所得的多项式函数求根。使用polysub函数对多项式作差：

poly_sub = np.polysub(poly_bhp, poly_vale)
xpoints = np.roots(poly_sub)
print("Intersection points", xpoints)

解出的交叉点如下：

Intersection points [ 27.73321597+0.j 27.51284094+0.j 24.32064343+0.j
18.86423973+0.j 12.43797190+1.73218179j 12.43797190-1.73218179j
6.34613053+0.62519463j 6.34613053-0.62519463j]

• (6) 得到的结果为复数，这不利于我们后续处理，除非时间也有实部和虚部。因此，这里需要用isreal函数来判断数组元素是否为实数：

reals = np.isreal(xpoints)
print("Real number?", reals)

结果如下：

Real number? [ True True True True False False False False]

可以看到有一部分数据为实数，因此我们用select函数选出它们。 select函数可以根据一组给定的条件，从一组元素中挑选出符合条件的元素并返回数组：

xpoints = np.select([reals], [xpoints])
xpoints = xpoints.real
print("Real intersection points", xpoints)

得到的实数交叉点如下所示：

Real intersection points [ 27.73321597 27.51284094 24.32064343 18.86423973 0. 0. 0. 0.]

• (7) 我们需要去掉其中为0的元素。 trim_zeros函数可以去掉一维数组中开头和末尾为0的元素：

print("Sans 0s", np.trim_zeros(xpoints))

去掉0元素后，输出结果如下所示：

Sans 0s [ 27.73321597 27.51284094 24.32064343 18.86423973]

小结

我们使用hanning函数对股票收益率数组进行了平滑处理，使用polysub函数对两个多项式作差运算，以及使用isreal函数判断数组元素是否为实数，并用select函数选出了实数元素。最后，我们用trim_zeros函数去掉数组首尾的0元素。

除了hanning函数之外，还使用其他的平滑函数，如hamming、 blackman、 bartlett以及kaiser。它们的使用方法和hanning函数类似。

示例完整代码如下：

import numpy as np
import matplotlib.pyplot as plt

N = 8

weights = np.hanning(N)
print("Weights", weights)

bhp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp) / bhp[ : -1]
smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns)[N-1:-N+1]

vale = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
vale_returns = np.diff(vale) / vale[ : -1]
smooth_vale = np.convolve(weights/weights.sum(), vale_returns)[N-1:-N+1]

K = 8
t = np.arange(N - 1, len(bhp_returns))
poly_bhp = np.polyfit(t, smooth_bhp, K)
poly_vale = np.polyfit(t, smooth_vale, K)

poly_sub = np.polysub(poly_bhp, poly_vale)
xpoints = np.roots(poly_sub)
print("Intersection points", xpoints)

reals = np.isreal(xpoints)
print("Real number?", reals)

xpoints = np.select([reals], [xpoints])
xpoints = xpoints.real
print("Real intersection points", xpoints)

print("Sans 0s", np.trim_zeros(xpoints))

plt.plot(t, bhp_returns[N-1:], lw=1.0)
plt.plot(t, smooth_bhp, lw=2.0)

plt.plot(t, vale_returns[N-1:], lw=1.0)
plt.plot(t, smooth_vale, lw=2.0)
plt.show()