前2章已经整理如下：

第1章NumPy入门

第2章NumPy基础

在本章中，我们将学习NumPy的常用函数。具体来说，我们将以分析历史股价为例，介绍怎样从文件中载入数据，以及怎样使用NumPy的基本数学和统计分析函数。这里还将学习读写文件的方法，并尝试函数式编程和NumPy 线性代数运算。

本章涵盖以下内容：

• 数组相关的函数；
• 从文件中载入数据；
• 将数组写入文件；
• 简单的数学和统计分析函数。

3.1 文件读写

首先，我们来学习使用NumPy读写文件。通常情况下，数据是以文件形式存储的。学会读写文件是深入学习NumPy的基础。

3.2 动手实践：读写文件

作为文件读写示例，我们创建一个单位矩阵并将其存储到文件中，并按照如下步骤完成。

(1) 单位矩阵，即主对角线上的元素均为1，其余元素均为0的正方形矩阵。在NumPy中可以用eye函数创建一个这样的二维数组，我们只需要给定一个参数，用于指定矩阵中1的元素个数。例如，创建2×2的数组：

回顾面试题： 创建ndarray有哪些基本的方法？

1)使用np.array()由python list创建

1
2
list1 = [1, 2, 3, 5] 
n = np.array(list1)

1
2
3
# 由ndarray变回list
list2 = n.tolist()
list2

注意：

- numpy默认ndarray的所有元素的类型是相同的(即同质)

- 如果传进来的列表中包含不同的类型，则统一为同一类型，优先级：str>float>int

例如：

2)使用np的routines函数创建—-需要记住11个函数

1
2
3
4
5
6
7
8
9
10
11
np.ones(shape, dtype=None, order='C')
np.zeros(shape, dtype=float, order='C')
np.full(shape, fill_value, dtype=None, order='C')
np.eye(N, M=None, k=0, dtype=float)      对角线为1其他的位置为0
np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)
np.arange([start, ]stop, [step, ]dtype=None)
np.random.randint(low, high=None, size=None, dtype='l')
np.random.randn(d0, d1, ..., dn)  
np.random.normal(loc=0.0, scale=1.0, size=None)
np.random.random(size=None)  
np.random.rand(d0, d1, d2, d3...,dn)

面试记忆方法：

一零满：np.ones, np.zeros, np.full（都是填充的）

角线均：np.eye, np.linspace, np.arange（对角线为1，均匀间隔，区间步长）

步长整：np.random.randint（整数）

正态浮：np.random.randn, np.random.normal, np.random.random, np.random.rand（正态分布和随机浮点数）

---

下面来看具体用法：

下面是对这些NumPy方法的详细介绍，包括它们的用法和需要注意的点：

---

1. np.ones(shape, dtype=None, order='C')

• 功能：生成一个由1填充的数组。
• 参数： > > - shape：数组的形状，可以是整数或元组。比如，(3, 4)表示一个3行4列的矩阵。
  • dtype：数据类型，默认为float。可以指定为任何NumPy支持的类型，如int，float64等。
  • order：决定数组的内存布局。'C'表示按行优先（C风格），'F'表示按列优先（Fortran风格）。
• 注意：如果不指定dtype，默认会生成浮点类型的数组。

2. np.zeros(shape, dtype=float, order='C')

• 功能：生成一个由0填充的数组。
• 参数：和np.ones()类似。
• 注意：常用于初始化一个空数组或矩阵。

3. np.full(shape, fill_value, dtype=None, order='C')

• 功能：生成一个指定形状并填充指定值的数组。
• 参数： > > - shape：数组的形状。
  • fill_value：数组中每个元素的值。
  • dtype：数据类型，默认为None，将根据fill_value自动推导。
  • order：内存布局，默认为'C'。
• 注意：可以用来生成任何值的矩阵，不仅仅是0或1。

4.np.eye(N, M=None, k=0, dtype=float)

• 功能：生成一个单位矩阵（对角线为1，其他位置为0）。
• 参数： > > - N：矩阵的行数。
  • M：矩阵的列数。如果不指定，默认为N，生成一个方阵。
  • k：对角线的位置，k=0是主对角线，k>0是上对角线，k<0是下对角线。
  • dtype：数据类型。
• 注意：如果需要生成其他对角线矩阵（如上对角线或下对角线），可以调整k的值。

5. np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)

• 功能：生成一个指定范围内的均匀间隔的数字序列。
• 参数： > > - start：序列的起始值。
  • stop：序列的结束值。
  • num：生成的数字个数，默认是50。
  • endpoint：是否包含stop值，默认是True。
  • retstep：如果为True，返回间隔的大小，默认为False。
  • dtype：数据类型。
• 注意：如果endpoint=False，序列会停止在stop之前，确保不会包含stop值。

6. np.arange([start, ]stop, [step, ]dtype=None)

• 功能：生成一个指定范围内的数字序列，类似于Python的range()函数。
• 参数： > > - start：序列的起始值，默认为0。（可选）
  • stop：序列的结束值（不包含）。（必须提供）
  • step：步长，默认为1。（可选）
  • dtype：数据类型。
• 注意：stop不包括在内。如果step为负数，start必须大于stop。

7. np.random.randint(low, high=None, size=None, dtype='l')

生成指定范围内的随机整数。

• 参数： > > - low：随机整数的下界。
  • high：随机整数的上界。
  • size：输出的数组形状。
  • dtype：数据类型，默认为整数类型。
• 注意：返回的数组包含的是指定区间内的整数，且 high 不包含在内。

8. np.random.randn(d0, d1, ..., dn)

生成标准正态分布的随机数（均值为 0，标准差为 1）。

• 参数： > > - d0, d1, ..., dn：各维度的大小。
• 注意：返回的是一个服从标准正态分布的数组，形状由参数指定。

9. np.random.normal(loc=0.0, scale=1.0, size=None)

生成符合正态分布的随机数，可以自定义均值和标准差。

• 参数： > > - loc：均值，默认为 0。
  • scale：标准差，默认为 1。
  • size：输出的数组形状。
• 注意：生成的是符合指定均值和标准差的随机数。

10. np.random.random(size=None)

生成 0 到 1 之间的随机浮点数。

• 参数： > > - size：输出的数组形状。
• 注意：返回的是一个包含 0 到 1 之间的随机浮点数的数组。

11. np.random.rand(d0, d1, ..., dn)

生成均匀分布在 [0, 1) 区间的随机数。

• 参数： > > - d0, d1, ..., dn：各维度的大小。
• 注意：类似于 np.random.random()，但接受的是维度参数而不是 size。

1
2
arr = np.eye(2, dtype = "int")
arr

注意：不指定dtype，默认创建的是浮点类型。

(2) 使用savetxt函数将数据存储到文件中，当然我们需要指定文件名以及要保存的数组。

1
np.savetxt("eye.txt", arr)

上面的代码会创建一个eye.txt文件，你可以检查文件内容和设想的是否一致。

刚才做了些什么 :读写文件是数据分析的一项基本技能。我们用savetxt函数进行了写文件的操作。在此之前，我们还用eye函数创建了一个单位矩阵。

3.3 CSV 文件

CSV（Comma-Separated Value，逗号分隔值）格式是一种常见的文件格式。通常，数据库的转存文件就是CSV格式的，文件中的各个字段对应于数据库表中的列。众所周知，电子表格软件（如Microsoft Excel）可以处理CSV文件。

3.4 动手实践：读入 CSV 文件

我们应该如何处理CSV文件呢？幸运的是，NumPy中的loadtxt函数可以方便地读取CSV文件，自动切分字段，并将数据载入NumPy数组。

下面，我们以载入苹果公司的历史股价数据为例展开叙述。

股价数据存储在CSV文件中，第一列为股票代码以标识股票（苹果公司股票代码为AAPL），第二列为dd-mm-yyyy格式的日期，第三列为空，随后各列依次是开盘价、最高价、最低价和收盘价，最后一列为当日的成交量。

下面为一行数据：

1
AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800 

从现在开始，我们只关注股票的收盘价和成交量。在上面的示例数据中，收盘价为336.1，成交量为21144800。我们将收盘价和成交量分别载入到两个数组中，如下所示：

注：下面的data.csv可以从这个网址下载 https://github.com/sundaygeek/numpy-beginner-guide/blob/master/ch3code/data.csv

1
c,v = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,7))

np.loadtxt 是 NumPy 中用于从文本文件加载数据的函数，通常用于读取包含数字的文件，并将其转换为 NumPy 数组。它非常适用于读取结构化的文本数据，如 CSV 文件、空格分隔的文件等。

1
np.loadtxt(fname, dtype=<class 'float'>, delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0, encoding='bytes', max_rows=None)

参数：

1. fname：

   • 要加载的文件名（字符串）或文件对象。如果是字符串，表示文件路径；如果是文件对象，表示已经打开的文件。

2. dtype：

   • 要加载的数据的类型，默认为 float。可以设置为其他类型（如 int，str 等）。

3. delimiter：

   • 数据分隔符，用来分割文件中的各个数据。默认为 None，表示空格或换行符分隔。如果数据是用逗号、制表符或其他分隔符分隔的，可以通过设置 delimiter 来指定。

   • 例如：delimiter=',' 用于 CSV 文件。
4. converters：

   • 一个字典，用于将列索引映射到一个函数。例如，可以用一个自定义函数来处理某一列的数据。

5. skiprows：

   • 跳过文件开始的行数（默认为 0）。常用于跳过文件中的标题行。

6. usecols：

   • 一个整数或整数元组，指定要加载的列（从 0 开始的索引）。如果是 None，表示加载所有列。

7. unpack：

   • 如果为 True，返回一个元组，其中每个元素对应一个列。默认是 False。

8. ndmin：

   • 最小维度，指定返回的数组的最小维度。默认是 0，即自动选择。

9. encoding：

   • 如果文件是字符串类型数据，可以设置编码（如 'utf-8'）。

10. max_rows：

    1
    - 最大行数，限制加载的最大行数。如果为 `None`，则加载所有行。
    

可以看到，数据存储在data.csv文件中，我们设置分隔符为,（英文标点逗号），因为我们要处理一个CSV文件。usecols的参数为一个元组，以获取第7字段至第8字段的数据，也就是股票的收盘价和成交量数据。unpack参数设置为True，意思是分拆存储不同列的数据，即分别将收盘价和成交量的数组赋值给变量c和v。

刚才做了些什么 :CSV文件是一种经常用于数据处理的文件格式。我们用loadtxt函数读取了一个包含股价数据的CSV文件，用delimiter参数指定了文件中的分隔符为英文逗号，用usecols中的参数指定了我们感兴趣的数据列，并将unpack参数设置为True使得不同列的数据分开存储，以便随后使用。

3.5 成交量加权平均价格（VWAP）

VWAP（Volume-Weighted Average Price，成交量加权平均价格）是一个非常重要的经济学量，它代表着金融资产的“平均”价格。某个价格的成交量越高，该价格所占的权重就越大。VWAP 就是以成交量为权重计算出来的加权平均值，常用于算法交易。

3.6 动手实践：计算成交量加权平均价格

我们将按如下步骤计算。 (1) 将数据读入数组。 (2) 计算VWAP。

补充：计算VWAP的公式\(VWAP = \frac{\sum(收盘价*成交量)}{\sum 成交量}\)

利用公式np.average计算。

np.average 是 NumPy 中用于计算数组元素的加权或非加权平均值的函数。它不仅能计算普通的算术平均值，还可以根据指定的权重计算加权平均值。

1
np.average(a, axis=None, weights=None, returned=False)

参数：

1. a：

   • 输入数组，表示要计算平均值的数组。

2. axis：

   • 指定沿哪个轴计算平均值。如果是 None（默认值），则对整个数组计算平均值。如果是整数或元组，表示沿着指定轴或多个轴计算。

3. weights：

   • 指定权重的数组。默认是 None，表示所有元素的权重相等。如果提供了 weights 数组，则按该数组提供的权重进行加权平均。

4. returned：

   • 如果为 True，则返回加权平均值以及权重和的总和（仅在 weights 非 None 时有效）。默认是 False，只返回平均值。

1
2
3
c,v = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,7))
vwap = np.average(c, weights = v)
vwap

刚才做了些什么 :这很容易，不是吗？我们仅仅调用了average函数，并将v作为权重参数使用，就完成了 VWAP的计算。此外，NumPy中也有计算算术平均值的函数。

3.6.1 算术平均值函数

NumPy中的mean函数很友好，一点儿也不mean（该词有“尖酸刻薄”的意思）。这个函数可以计算数组元素的算术平均值。具体用法如下：

1
np.mean(a, axis=None, dtype=None, out=None, keepdims=False)

np.mean 是 NumPy 中用于计算数组的算术平均值（即均值）的函数。它可以计算整个数组的平均值，也可以按指定的轴计算平均值。它与 np.average 的功能类似，但不支持加权平均。

参数：

1. a：
   • 输入数组，表示要计算平均值的数组。
2. axis：
   • 指定沿哪个轴计算平均值。如果是 None（默认值），则对整个数组计算平均值。如果是整数或元组，表示沿着指定轴或多个轴计算。
3. dtype：
   • 结果的类型。如果提供了 dtype，则计算结果将转换为该类型。默认为 None，即根据输入数组的类型自动选择。
4. out：
   • 用于存储结果的数组。默认是 None，即返回一个新的数组。
5. keepdims：
   • 如果为 True，则在计算平均值后，保留原数组的维度。默认是 False。

例如：

3.6.2 时间加权平均价格

在经济学中，TWAP（Time-Weighted Average Price，时间加权平均价格）是另一种“平均” 价格的指标。既然我们已经计算了VWAP，那也来计算一下TWAP吧。其实TWAP只是一个变种而已，基本的思想就是最近的价格重要性大一些，所以我们应该对近期的价格给以较高的权重。最简单的方法就是用arange函数创建一个从0开始依次增长的自然数序列，自然数的个数即为收盘价的个数。当然，这并不一定是正确的计算TWAP的方式。事实上，本书中关于股价分析的大部分示例都仅仅是为了说明问题。计算TWAP的代码如下。

1
2
3
t = np.arange(len(c))
TWAP = np.average(c, weights = t)
TWAP

在这个例子中，TWAP的值甚至比算术平均值还要高。

突击测验：计算加权平均值 问题1 以下哪个函数可以返回数组元素的加权平均值？

(1) weighted average (2) waverage (3) average (4) avg

答案：（3），np.average() 函数可以计算加权平均值。通过 weights 参数，你可以为每个元素指定权重，以计算加权平均值。

其他选项 weighted average、waverage 和 avg 不是 NumPy 中的有效函数名称。

勇敢出发：计算其他数据的平均值 请尝试对开盘价做相同的计算，并计算成交量和其他各种价格的算术平均值。

此答案，略，上面几乎都有了，无非换个对象。

3.7 取值范围

通常，我们不仅仅想知道一组数据的平均值，还希望知道数据的极值以及完整的取值范 围——最大值和最小值。我们的股价示例数据中已经包含了每天的股价范围——最高价和最低价。但是，我们还需要知道最高价的最大值以及最低价的最小值。不然，我们怎样才能知道自己的股票是赚了还是赔了呢？

3.8 动手实践：找到最大值和最小值

min函数和max函数能够满足需求。我们按如下步骤来找最大值和最小值。

回顾数据：

股价数据存储在CSV文件中，第一列为股票代码以标识股票（苹果公司股票代码为AAPL），第二列为dd-mm-yyyy格式的日期，第三列为空，随后各列依次是开盘价、最高价、最低价和收盘价，最后一列为当日的成交量。

1
AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800 

(1) 首先，需要再次读入数据，将每日最高价和最低价的数据载入数组：

1
max_price, min_price = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols = (4,5))

唯一需要修改的就是usecols中的参数，因为最高价和最低价与之前的数据在不同的列中。

(2) 下方的代码即可获取价格区间：

1
2
highest = np.max(max_price)
lowest = np.min(min_price)

现在，计算区间中点就很容易了，留给读者自己尝试练习。

计算区间中点：

区间的中点可以通过最高值和最低值的平均值来计算。\(midpoint = \frac{highest+lowest}{2}\)

1
midpoint = np.mean([highest, lowest])

(3) NumPy中有一个ptp函数可以计算数组的取值范围。该函数返回的是数组元素的最大值和最小值之间的差值。也就是说，返回值等于max(array) - min(array)。调用ptp函数：

刚才做了些什么 :我们定义了股价的最大值和最小值的取值范围。最大值是在每日最高价数据上使用max函数得到的，而最低价是在每日最低价数据上使用min函数得到的。我们还用ptp函数计算了极差，即最大值和最小值之间的差值。

3.9 统计分析

股票交易者对于收盘价的预测很感兴趣。常识告诉我们，这个价格应该接近于某种均值。算数平均值和加权平均值都是在数值分布中寻找中心点的方法。然而，它们对于异常值（outlier）既不鲁棒也不敏感。举例来说，如果我们有一个高达100万美元的收盘价，这将影响到我们的计算结果。

算术平均值和加权平均值受异常数据（数据偏差特别大的数据）的影响较大。也就是说，算数平均值和加权平均值对极端值（如非常大的或非常小的数值）过于敏感，可能会导致计算结果不准确。

3.10 动手实践：简单统计分析

我们可以用一些阈值来除去异常值，但其实有更好的方法，那就是中位数。将各个变量值按大小顺序排列起来，形成一个数列，居于数列中间位置的那个数即为中位数。例如，我们有1、2、3、4、5这5个数值，那么中位数就是中间的数字3。下面是计算中位数的步骤。

(1) 计算收盘价的中位数。你已经知道如何从CSV文件中读取数据到数组中了，因此只需要复制一行代码并确保只获取收盘价数据即可，如下所示.

(2) 一个叫做median的函数将帮助我们找到中位数。我们调用它并立即打印出结果。

1
2
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
np.median(c)

(3) 既然这是我们首次使用median函数，我们来检查一下结果是否正确。这可不是因为我们多疑！当然，我们可以将整个数据文件浏览一遍并人工找到正确的答案，但那样太无趣了。我们将对价格数组进行排序，并输出排序后居于中间位置的值，这也就是模拟了寻找中位数的算法。msort函数可以帮我们完成第一步。我们将调用这个函数，获得排序后的数组，并输出结果。

使用msort时，有提示：

1
msort is deprecated, use np.sort(a, axis=0) instead

所以还是用sort。

太好了，代码生效了！现在，我们来获取位于中间的那个数字：

(4) 咦，这个值和median函数给出的值不一样，这是怎么回事？经过仔细观察我们发现，median函数返回的结果甚至根本没有在我们的数据文件里出现过。这就更奇怪了！在给NumPy 团队提交bug报告之前，我们先来看下文档。原来这个谜团很容易解开。原因就在于我们的简单算法模拟只对长度为奇数的数组奏效。对于长度为偶数的数组，中位数的值应该等于中间那两个数的平均值。因此，输入如下代码：

1
(sorted_arr[N//2] + sorted_arr[(N-1)//2])/2

成功了！ (5) 另外一个我们关心的统计量就是方差。方差能够体现变量变化的程度。在我们的例子中，方差还可以告诉我们投资风险的大小。那些股价变动过于剧烈的股票一定会给持有者制造麻烦。在NumPy中，计算方差只需要一行代码，看下面：

1
np.var(c)

(6) 既然我们不相信NumPy的函数，那就再次根据文档中方差的定义来复核一下结果。注意，这里方差的定义可能与你在统计学的书中看到的不一致，但这个定义在统计学上更为通用。

方差是指各个数据与所有数据算术平均数的离差平方和除以数据个数所得到的值。

这段话的解释是关于方差（Variance）的定义。方差是描述数据集分布的一个统计量，它衡量的是数据点相对于其均值（算术平均数）的离散程度。具体来说，方差度量的是每个数据点与均值之间的差异有多大，并且该差异是以平方的形式来衡量的。

对于一个数据集\(x_1,x_2,…,x_n\)，其方差\(\sigma^2\)（对于总体数据）可以通过以下公式计算：\(\sigma^2 = \frac{1}{n}\sum_{i = 1}^n (x_i- \mu)^2\)

• \(x_i\)是数据集中的第i个数据点。
• \(\mu\) 是数据集的均值（算术平均数），即\(\mu = \frac{1}{n}\sum_{i = 1}^n {x_i}\)

一些书里面告诉我们，应该用数据个数减1去除离差平方和:注意样本方差和总体方差在计算上的区别。总体方差是用数据个数去除离差平方和，而样本方差则是用样本数据个数减1去除离差平方和，其中样本数据个数减1（即n-1）称为自由度。之所以有这样的差别，是为了保证样本方差是一个无偏估计量。

1
np.mean((c - c.mean())**2)

这正是我们希望得到的结果！

刚才做了些什么 ：你可能已经注意到了一些新东西。我们直接在数组c上调用了mean方法c.mean()，是的，没有写错。ndarray对象有mean方法，这将给你带来便利。从现在开始，记住这种用法是正确无误的。

3.11 股票收益率

在学术文献中，收盘价的分析常常是基于股票收益率和对数收益率的。

简单收益率是指相邻两个价格之间的变化率，而对数收益率是指所有价格取对数后两两之间的差值。我们在高中学习过对数的知识，“a”的对数减去“b”的对数就等于“a除以b”的对数。

因此，对数收益率也可以用来衡量价格的变化率。

注意，由于收益率是一个比值，例如我们用美元除以美元（也可以是其他货币单位），因此它是无量纲的。

总之，投资者最感兴趣的是收益率的方差或标准差，因为这代表着投资风险的大小。

3.12 动手实践：分析股票收益率

按照如下步骤分析股票收益率。 (1) 首先，我们来计算简单收益率。NumPy中的diff函数可以返回一个由相邻数组元素的差值构成的数组。这有点类似于微积分中的微分。为了计算收益率，我们还需要用差值除以前一天的价格。不过这里要注意，diff返回的数组比收盘价数组少一个元素。经过仔细考虑，我们使用如下代码：

1
returns = np.diff(arr)/arr[:-1]

• np.diff(arr) 用于计算一个数组 arr 中相邻元素之间的差值。也就是说，它会返回一个新数组，其中每个元素是原数组中相邻两个元素的差。
• arr[:-1] 是对数组 arr 的切片操作，表示从数组的第一个元素到倒数第二个元素（不包括最后一个元素）。

注意，我们没有用收盘价数组中的最后一个值做除数。接下来，用std函数计算标准差：

1
np.std(returns)

(2) 对数收益率计算起来甚至更简单一些。我们先用log函数得到每一个收盘价的对数，再对结果使用diff函数即可。

1
logreturns = np.diff(np.log(arr))

一般情况下，我们应检查输入数组以确保其不含有零和负数。否则，将得到一个错误提示。 不过在我们的例子中，股价总为正值，所以可以将检查省略掉。

(3) 我们很可能对哪些交易日的收益率为正值非常感兴趣。在完成了前面的步骤之后，我们只需要用where函数就可以做到这一点。where函数可以根据指定的条件返回所有满足条件的数组元素的索引值。输入如下代码：

1
posretindices = np.where(returns > 0)

即可输出该数组中所有正值元素的索引。

(4) 在投资学中，波动率（volatility）是对价格变动的一种度量。历史波动率可以根据历史价格数据计算得出。计算历史波动率（如年波动率或月波动率）时，需要用到对数收益率。年波动率等于对数收益率的标准差除以其均值，再除以交易日倒数的平方根，通常交易日取252天。我们用std和mean函数来计算，代码如下所示：

1
2
annual_volatility = np.std(logreturns)/np.mean(logreturns)
annual_volatility = annual_volatility/np.sqrt(1/252)

(5) 请注意sqrt函数中的除法运算。在Python中，整数的除法和浮点数的除法运算机制不同， 我们必须使用浮点数才能得到正确的结果。与计算年波动率的方法类似，计算月波动率如下：

1
monthly_volatility = annual_volatility * np.sqrt(1/12.)

刚才做了些什么 : 我们用计算数组相邻元素差值的diff函数计算了简单收益率，用计算数组元素自然对数的log函数计算了对数收益率。最后，我们计算了年波动率和月波动率。

3.13 日期分析

你是否有时候会有星期一焦虑症和星期五狂躁症？想知道股票市场是否受上述现象的影响？我认为这值得深入研究。

3.14 动手实践：分析日期数据

首先，我们要读入收盘价数据。

数据是这样的：

股价数据存储在CSV文件中，第一列为股票代码以标识股票（苹果公司股票代码为AAPL），第二列为dd-mm-yyyy格式的日期，第三列为空，随后各列依次是开盘价、最高价、最低价和收盘价，最后一列为当日的成交量。

下面为一行数据：

1
AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800 

随后，根据星期几来切分收盘价数据，并分别计算平均价格。

最后，我们将找出一周中哪一天的平均收盘价最高，哪一天的最低。在我们动手之前，有一个善意的提醒：你可能希望利用分析结果在某一天买股票或卖股票，然而我们这里的数据量不足以做出可靠的决策，请先咨询专业的统计分析师再做决定！

程序员不喜欢日期，因为处理日期总是很烦琐。NumPy是面向浮点数运算的，因此需要对日期做一些专门的处理。请自行尝试如下代码，单独编写脚本文件或使用本书附带的代码文件：

1
dates,closing_price = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(1,6))

执行以上代码后会得到一个错误提示：

1
ValueError: could not convert string '28-01-2011' to float64 at row 0, column 2.

按如下步骤处理日期。

(1) 显然，NumPy尝试把日期转换成浮点数。我们需要做的是显式地告诉NumPy怎样来转换日期，而这需要用到loadtxt函数中的一个特定的参数。这个参数就是converters，它是一本数据列和转换函数之间进行映射的字典。

为此，我们必须写出转换函数：

1
2
3
# 将一个日期字符串转换成对应星期几的数字，例如   星期一  0       星期二  1
def dateset2num(s):
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()

datetime.datetime.strptime(s, "%d-%m-%Y"):

• 这一部分是将字符串 s 解析为一个 datetime 对象。strptime 是一个用于将日期字符串转换为 datetime 对象的函数。这里的格式 "%d-%m-%Y" 表示：

  • %d：表示日期（01到31），例如 05。

  • %m：表示月份（01到12），例如 09。

  • %Y：表示年份（四位数），例如 2025。

• .date():

  • 从 datetime 对象中提取出日期部分，返回一个 date 对象，丢弃时间信息。

• .weekday():

  方法会返回该日期是星期几。返回值是一个整数，表示星期几：

  • 0：星期一
  • 1：星期二
  • 2：星期三
  • 3：星期四
  • 4：星期五
  • 5：星期六
  • 6：星期天

测试：

我们将日期作为字符串传给datestr2num函数，如“28-01-2011”。这个字符串首先会按照指定的形式”%d-%m-%Y“转换成一个datetime对象。补充一点，这是由Python标准库提供的功能，与NumPy无关。随后，datetime对象被转换为date对象。最后，调用weekday方法返回一个数字。如同你在注释中看到的，这个数字可以是0到6的整数，0代表星期一，6代表星期天。当然，具体的数字并不重要，只是用作标识。

(2) 接下来，我们将日期转换函数挂接上去，这样就可以读入数据了。

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

# 将一个日期字符串转换成对应星期几的数字
def dateset2num(s):
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()

# 读取数据时，确保使用正确的编码
dates, closing_price = np.loadtxt(
    "data.csv",  # CSV 文件名
    delimiter=",",  # 分隔符
    unpack=True,  # 将结果拆分为不同的数组
    usecols=(1, 6),  # 读取第2列和第7列
    converters={1: dateset2num},  # 对第2列应用日期转换函数
    dtype="float64",  # 确保其他列的类型为 float64
)

注意：上面的代码会报错：ValueError: could not convert string '28-01-2011' to float64 at row 0, column 2.

明明写了dateset2num函数，处理还是有问题。书中也没有说明为什么。

根据错误信息，问题出在 strptime() 的 argument 1 must be str, not bytes，即尝试将一个字节（bytes）类型的数据传递给 strptime，而 strptime 期望的是一个字符串（str）类型的输入。

这是因为 np.loadtxt 在读取文件时，默认情况下将所有的数据按 bytes 格式读取。这导致你在 dateset2num 函数中传入的参数是字节类型（bytes），而 datetime.strptime() 只能处理字符串类型（str），因此报错。

解决办法：在dateset2num函数中，传入strptime之前，将s进行编码，即s = s.decode('utf-8')

1
2
3
4
5
6
7
8
9
10
11
12
13
# 读取数据时使用 converters 来转换日期列
def dateset2num(s):
    # s 是字节类型，转换成字符串后再进行处理
    s = s.decode('utf-8')  # 解码为字符串
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()
dates, closing_price = np.loadtxt(
    "data.csv",  # CSV 文件名
    delimiter=",",  # 分隔符
    unpack=True,  # 将结果拆分为不同的数组
    usecols=(1, 6),  # 读取第2列和第7列
    converters={1: dateset2num},  # 对第2列应用日期转换函数
    dtype="float64"  # 确保其他列的类型为 float64
)

然后就可以成功了：

如你所见，没有出现星期六和星期天。股票交易在周末是休市的。

(3) 我们来创建一个包含5个元素的数组，分别代表一周的5个工作日。数组元素将初始化为0。

1
averages = np.zeros(5)

这个数组将用于保存各工作日的平均收盘价。

(4) 我们已经知道，where函数会根据指定的条件返回所有满足条件的数组元素的索引值。take函数可以按照这些索引值从数组中取出相应的元素。我们将用take函数来获取各个工作日的收盘价。在下面的循环体中，我们将遍历0到4的日期标识，或者说是遍历星期一到星期五，然后用where函数得到各工作日的索引值并存储在indices数组中。在用take函数获取这些索引值相应的元素值。最后，我们对每个工作日计算出平均值存放在averages数组中。代码如下：

1
2
3
4
5
6
for i in range(5):
    indices = np.where(dates == i)
    prices = np.take(closing_price, indices)
    avg = np.mean(prices)
    print(avg)
    averages[i] = avg

(5) 如果你愿意，还可以找出哪个工作日的平均收盘价是最高的，哪个是最低的。这很容易做到，用max和min函数即可，代码如下：

1
2
3
4
np.max(averages)
np.argmax(averages)
np.min(averages)
np.argmin(averages)

刚才做了些什么 : argmin函数返回的是averages数组中最小元素的索引值，这里是4，也就是星期五。而argmax函数返回的是averages数组中最大元素的索引值，这里是2，也就是星期三。

勇敢出发：计算VWAP和TWAP

这一节的内容很有趣！在示例数据中，你的苹果股票在星期五是最便宜的一天而星期三是最值钱的一天。先不管我们的数据量有多小，思考一下是否有更合理的计算平均值的方式？我们是不是应该考虑成交量呢？如果计算时间加权的平均值是不是更有意义呢？快尝试一下吧！请计算VWAP和TWAP，你可以在本章的开头部分找到一些相关的提示。

【计算VWAP】 首先注意：成交量是最后一列

股价数据存储在CSV文件中，第一列为股票代码以标识股票（苹果公司股票代码为AAPL），第二列为dd-mm-yyyy格式的日期，第三列为空，随后各列依次是开盘价、最高价、最低价和收盘价，最后一列为当日的成交量。

下面为一行数据：

1
AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800 

计算VWAP的公式\(VWAP = \frac{\sum(收盘价*成交量)}{\sum 成交量}\)

这个已经学习过了：

首先读取需要的数据：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 读取日期,收盘价和成交量，其中收盘价和成交量用来计算VWAP
# 读取数据时使用 converters 来转换日期列
def dateset2num(s):
    # s 是字节类型，转换成字符串后再进行处理
    s = s.decode('utf-8')  # 解码为字符串
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()
dates, closing_price, turnover = np.loadtxt(
    "data.csv",  # CSV 文件名
    delimiter=",",  # 分隔符
    unpack=True,  # 将结果拆分为不同的数组
    usecols=(1, -2,-1),  # 读取第2列和第7列
    converters={1: dateset2num},  # 对第2列应用日期转换函数
    dtype="float64"  # 确保其他列的类型为 float64
)
dates

1
2
# 计算VWAP
np.average(closing_price, weights = turnover)

或者：

1
vwap = np.sum(closing_price * turnover) / np.sum(turnover)

和上面的结果相同。

计算TWAP：

1
twap = np.mean(closing_price)

所以不读取dates也是可以的。

3.15 周汇总

在之前的“动手实践”教程中，我们用的是盘后数据。也就是说，这些数据是将一整天的交易数据汇总得到的。如果你对棉花市场感兴趣，并且有数十年的数据，你可能希望对数据做进一步的汇总和压缩。开始动手吧。我们来把苹果股票数据按周进行汇总。

3.16 动手实践：汇总数据

我们将要汇总整个交易周中从周一到周五的所有数据。数据覆盖的时间段内有一个节假日：2月21日是总统纪念日。这天是星期一，美国股市休市，因此在我们的示例数据中没有这一天的数据记录。数据中的第一天为星期五，处理起来不太方便。按照如下步骤来汇总数据。

(1) 为了简单起见，我们只考虑前三周的数据，这样就避免了节假日造成的数据缺失。你可以稍后尝试对其进行拓展。

1
2
dates = dates[:16]
closing_price = closing_price[:16]

代码基于3.14节的教程。

(2) 首先我们来找到示例数据中的第一个星期一。回忆一下，在Python中星期一对应的编码是0，这可以作为where函数的条件。接着，我们要取出数组中的首个元素，其索引值为0。但where函数返回的结果是一个多维数组，因此要用ravel函数将其展平。

1
np.ravel(np.where(dates == 0))[0]

(3) 下面要做的是找到示例数据的最后一个星期五，方法和找第一个星期一类似。星期五相对应的编码是4。此外，我们用-1作为索引值来定位数组的最后一个元素。

1
last_friday = np.ravel(np.where(dates == 4))[-1]

接下来创建一个数组，用于存储三周内每一天的索引值。

1
weeks_indices = np.arange(first_monday, last_friday + 1)

(4) 按照每个子数组5个元素，用split函数切分数组：

1
np.split(weeks_indices, 3)

注意：原书此处写的参数是5，应该是3.

(5) 在NumPy中，数组的维度也被称作轴。现在我们来熟悉一下apply_along_axis函数。 这个函数会调用另外一个由我们给出的函数，作用于每一个数组元素上。目前我们的数组中有3 个元素，分别对应于示例数据中的3个星期，元素中的索引值对应于示例数据中的1天。在调用apply_along_axis时提供我们自定义的函数名summarize，并指定要作用的轴或维度的编号（如取1）、目标数组以及可变数量的summarize函数的参数。

np.apply_along_axis() 是 NumPy 中的一个强大函数，它可以对 数组的某一维度 应用一个指定的函数。

np.apply_along_axis(func1d, axis, arr, *args, **kwargs)

参数解析

• func1d：应用到 1D 数组 的函数（必须是作用于 1D 数组的函数）。
• axis：沿着哪个轴（axis=0 代表列（每一列作为 1D 处理），axis=1 代表行）。
• arr：输入的 numpy 数组。
• *args, *kwargs：额外传递给 func1d 的参数。

举例：

1
2
n = np.arange(1, 13).reshape(3, 4)
n1 = np.apply_along_axis(np.sum, axis = 1, arr = n)

1
weeksummary = np.apply_along_axis(summarize, 1, weeks_indices,opening_price, high, low, closing_price) 

(6) 编写summarize函数。该函数将为每一周的数据返回一个元组，包含这一周的开盘价、最高价、最低价和收盘价，类似于每天的盘后数据。

1
2
3
4
5
6
def summarize(a, o, h, l, c):
    monday_open = o[a[0]]
    week_high = np.max( np.take(h, a) )  
    week_low = np.min( np.take(l, a) )
    friday_close = c[a[-1]] 
    return("APPL", monday_open, week_high, week_low, friday_close) 

注意，我们用take函数来根据索引值获取数组元素的值，并用max和min函数轻松计算出一周的最高股价和最低股价。一周的开盘价即为周一的开盘价，而一周的收盘价即为周五的收盘价。

到此为止，完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
dates = dates[:16]
first_monday = np.ravel(np.where(dates == 0))[0]
last_friday = np.ravel(np.where(dates == 4))[-1]
weeks_indices = np.arange(first_monday, last_friday + 1)
weeks_indices = np.split(weeks_indices, 3)
opening_price, high, low, closing_price = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(3,4,5,6))

def summarize(a, o, h, l, c):
    monday_open = o[a[0]]
    week_high = np.max( np.take(h, a) )  
    week_low = np.min( np.take(l, a) )
    friday_close = c[a[-1]] 
    return("APPL", monday_open, week_high, week_low, friday_close) 
  
weeksummary = np.apply_along_axis(summarize, 1, weeks_indices,opening_price, high, low, closing_price) 

(7) 使用NumPy中的savetxt函数，将数据保存至文件。

1
np.savetxt("weeksummary.csv", weeksummary, delimiter=",", fmt = "%s")

如代码中所示，我们指定了文件名、需要保存的数组名、分隔符（在这个例子中为英文标点逗号）以及存储浮点数的格式。

格式字符串以一个百分号开始。接下来是一个可选的标志字符：-表示结果左对齐，0表示左端补0，+表示输出符号（正号+或负号-）。第三部分为可选的输出宽度参数，表示输出的最小位数。第四部分是精度格式符，以”.”开头，后面跟一个表示精度的整数。最后是一个类型指定字符，在我们的例子中指定为字符串类型。

刚才做了些什么 : 我们刚刚完成的事情在其他编程语言中几乎是无法完成的。我们定义了一个函数summarize，把它作为参数传给了apply_along_axis函数，而summarize的参数也通过apply_along_axis的参数列表便捷地传递了过去。

勇敢出发：优化代码

改进代码，使其能够处理节假日休市的情况。对代码进行性能分析，考察apply_along_ axis函数的使用能带来多少速度上的提升。

首先整理完整代码：

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
import numpy as np
import datetime

def dateset2num(s):
    # s 是字节类型，转换成字符串后再进行处理
    s = s.decode('utf-8')  # 解码为字符串
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()

def summarize(a, o, h, l, c):
    monday_open = o[a[0]]
    week_high = np.max( np.take(h, a) )  
    week_low = np.min( np.take(l, a) )
    friday_close = c[a[-1]] 
    return("APPL", monday_open, week_high, week_low, friday_close) 

dates, opening_price, high, low, closing_price = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(1,3,4,5,6), converters={1:dateset2num})

dates = dates[:16]
first_monday = np.ravel(np.where(dates == 0))[0]
last_friday = np.ravel(np.where(dates == 4))[-1]
weeks_indices = np.arange(first_monday, last_friday + 1)
weeks_indices = np.split(weeks_indices, 3)

weeksummary = np.apply_along_axis(summarize, 1, weeks_indices,opening_price, high, low, closing_price) 

np.savetxt("weeksummary.csv", weeksummary, delimiter=",", fmt = "%s")

【处理节假日休市】

假设有一个包含 假期日期 的列表或文件（例如，holidays）。在处理时，我们可以过滤掉这些日期的数据，以确保分析时不包含节假日的数据。

1
2
3
# 假设你有一个包含节假日日期的列表（例如 [datetime.date(2023, 1, 1), datetime.date(2023, 12, 25)]）
# holidays = {datetime.date(2011, 2, 21), datetime.date(2023, 12, 25)}  # 示例假期
holidays = {datetime.date(2011, 2, 21)} 

【优化性能】np.apply_along_axis() 在执行时其实是对每一维度的数组调用一个函数，可能导致性能较差，尤其在大数据集时。我们可以通过 列表推导式 和 预先计算索引 来优化。

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
import datetime

# 将日期字符串转换为星期几的数字
def dateset2num(s):
    s = s.decode('utf-8')  # 解码为字符串
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()

# 假设你有一个包含节假日日期的列表（例如 [datetime.date(2023, 1, 1), datetime.date(2023, 12, 25)]）
holidays = {datetime.date(2011, 2, 21)}   # 示例假期

# 检查日期是否是节假日
def is_holiday(date):
    return date in holidays

# 过滤掉节假日的数据
def filter_holidays(dates, open_prices, high_prices, low_prices, close_prices):
    valid_indices = []
    for i, date in enumerate(dates):
        if not is_holiday(date):
            valid_indices.append(i)
    return dates[valid_indices], open_prices[valid_indices], high_prices[valid_indices], low_prices[valid_indices], close_prices[valid_indices]

# 计算每周的数据汇总
def summarize(a, o, h, l, c):
    monday_open = o[a[0]]
    week_high = np.max(np.take(h, a))  
    week_low = np.min(np.take(l, a))
    friday_close = c[a[-1]] 
    return ("APPL", monday_open, week_high, week_low, friday_close)

# 读取数据
dates, opening_price, high, low, closing_price = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(1, 3, 4, 5, 6), converters={1: dateset2num})

# 过滤节假日数据
dates, opening_price, high, low, closing_price = filter_holidays(dates, opening_price, high, low, closing_price)

# 获取每周的索引
dates = dates[:16]  # 假设你只关心前16天的数据
first_monday = np.ravel(np.where(dates == 0))[0]  # 找到第一个星期一
last_friday = np.ravel(np.where(dates == 4))[-1]  # 找到最后一个星期五

# 获取每周的索引范围
weeks_indices = np.arange(first_monday, last_friday + 1)
weeks_indices = np.split(weeks_indices, 3)  # 假设我们分成三周

# 用列表推导代替 np.apply_along_axis 来加速处理
weeksummary = [summarize(week_indices, opening_price, high, low, closing_price) for week_indices in weeks_indices]

# 保存结果到 CSV
np.savetxt("weeksummary_new.csv", weeksummary, delimiter=",", fmt="%s")

性能比较：

1
2
3
4
5
6
7
8
9
10
11
import timeit

# 计时 np.apply_along_axis
apply_time = timeit.timeit(lambda: np.apply_along_axis(summarize, 1, weeks_indices, opening_price, high, low, closing_price), number=10)

# 计时 列表推导式
list_comprehension_time = timeit.timeit(lambda: [summarize(week_indices, opening_price, high, low, closing_price) for week_indices in weeks_indices], number=10)

print(f"np.apply_along_axis 执行时间: {apply_time:.4f} 秒")
print(f"列表推导式执行时间: {list_comprehension_time:.4f} 秒")

3.17 真实波动幅度均值（ATR）

ATR（Average True Range，真实波动幅度均值）是一个用来衡量股价波动性的技术指标。 ATR的计算并不是重点，只是作为演示几个NumPy函数的例子，包括maximum函数。

ATR 计算的核心在于对市场波动性进行平滑处理，具体是通过考虑市场的最高价、最低价以及前一天的收盘价来计算波动幅度。

3.18 动手实践：计算真实波动幅度均值

提前准备数据：读取高价、最低价和收盘价，用于后续计算

1
h, l, c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(4, 5, 6))

按照如下步骤计算真实波动幅度均值。 (1) ATR是基于N个交易日的最高价和最低价进行计算的，通常取最近20个交易日。

1
2
3
4
import sys
N = int(sys.argv[1])
h = h[-N:]
l = l[-N:]

注：在jupyter notebook中运行以上几行代码会报错ValueError: invalid literal for int() with base 10: '-f'，因为试图从命令行参数中读取一个整数值，而 sys.argv[1] 提供的值是 -f，它不能转换为整数。

在 Jupyter Notebook 或某些 IDE 中运行代码时，可能会自动传递一些参数给脚本，通常 -f 是指文件路径或者其他参数。

在某些环境（如 Jupyter Notebook）中，sys.argv 可能会包含额外的参数，导致无法按预期提取正确的参数。

解决方案：修改代码以处理默认参数： 如果你并不一定要从命令行获取 N 的值，你可以将其改为通过代码直接赋值，或者在没有提供命令行参数时使用一个默认值。

一会儿改成N = 20即可。

(2) 我们还需要知道前一个交易日的收盘价。

1
previousclose = c[-N -1:-1]

注：c[-N - 1: -1] 表示从数组 c 中取出从 -N-1 到 -1 索引范围的元素。负索引是从数组末尾开始倒数的索引

对于每一个交易日，计算以下各项。

• h – l 当日股价范围，即当日最高价和最低价之差。
• h – previousclose 当日最高价和前一个交易日收盘价之差。
• previousclose – l 前一个交易日收盘价和当日最低价之差。

(3) max函数返回数组中的最大值。基于上面计算的3个数值，我们来计算所谓的真实波动幅度，也就是这三者的最大值。现在我们想在一组数组之间按照元素挑选最大值——也就是在所有的数组中第一个元素的最大值、第二个元素的最大值等。为此，需要用NumPy中的maximum函数，而不是max函数。

1
truerange = np.maximum(h - l, h - previousclose, previousclose - l)

(4) 创建一个长度为N的数组atr，并初始化数组元素为0。

1
atr = np.zeros(N)

(5) 这个数组的首个元素就是truerange数组元素的平均值。

1
atr[0] = np.mean(truerange)

用如下公式计算其他元素的值： \(\frac{((N-1)PATR+TR)}{N}\)

这里，PATR表示前一个交易日的ATR值，TR即当日的真实波动幅度。

1
2
3
for i in range(1, N):
    atr[i] = (N - 1) * atr[i-1] + truerange[i]
    atr[i] /= N

刚才做了些什么 : 我们生成了3个数组，分别表示3种范围——当日股价范围（h-l），当日最高价和前一个交易日收盘价之差(h - previousclose)，以及前一个交易日收盘价和当日最低价之差(previousclose - l)。这告诉我们股价波动的范围，也就是波动性的大小。ATR的算法要求我们找出三者的最大值。而之前使用的max函数只能给出一个数组内的最大元素值，并非这里所需要的。我们要在一组数组之间挑选每一个元素位置上的最大值，也就是在所有的数组中第一个元素的最大值、第二个元素的最大值等。在这一节的“动手实践”教程中，我们了解到maximum函数可以做到这一点。最终，我们根据每一天的真实波动幅度值计算出一个移动平均值。

在随后的教程中，我们将学习更好的计算移动平均值的方法。

勇敢出发：尝试minimum函数

除了maximum函数，NumPy中还有一个minimum函数。你可能已经猜到了这个函数的功能。请创建一小段脚本或在IPython中创建一个会话，以证明你的猜想。

1
2
3
4
5
# randint(low, high=None, size=None, dtype=int)
n1 = np.random.randint(0, 10, size=(3, 4))
n2 = np.random.randint(0, 10, size=(3, 4))
n3 = np.random.randint(0, 10, size=(3, 4))
np.minimum(n1, n2, n3)

3.19 简单移动平均线

简单移动平均线（simple moving average，SMA）通常用于分析时间序列上的数据。为了计算它， 我们需要定义一个N个周期的移动窗口，在我们的例子中即N个交易日。我们按照时间序列滑动这个窗口，并计算窗口内数据的均值。

3.20 动手实践：计算简单移动平均线

移动平均线只需要少量的循环和均值函数即可计算得出，但使用NumPy还有更优的选 择——convolve函数。简单移动平均线只不过是计算与等权重的指示函数的卷积，当然，也可以是不等权重的。

卷积是分析数学中一种重要的运算，定义为一个函数与经过翻转和平移的另一个函数的乘积的积分。

卷积的数学定义：

卷积通常定义为两个函数\(f(t)\)和\(g(t)\)的乘积，记作\((f*g)(t)\)，它的定义如下：\((f*g)(t)=\int_{- \infty}^{+ \infty}f(\tau)g(t- \tau)d \tau\)。

其中：

• \(f(t)\)和\(g(t)\)是两个函数。
• \(g(t- \tau)\)表示将\(g(t)\)函数沿时间轴翻转，并对其进行平移.
• \(\tau\)是一个积分变量，表示\(f(t)\)和翻转平移后的\(g(t)\)函数在不同位置上的重叠程度。

举例：使用 NumPy convolve 计算 SMA

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

# 示例数据
prices = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

# 计算 3 天的 SMA
N = 3
weights = np.ones(N) / N  # 生成等权重窗口 [1/N, 1/N, ..., 1/N]
print(weights)

# 使用 NumPy 的 `convolve` 计算移动平均
sma = np.convolve(prices, weights, mode='valid')

print(sma)  # 输出: [2. 3. 4. 5. 6. 7. 8.]

按照如下步骤计算简单移动平均线。 (1) 使用ones函数创建一个长度为N的元素均初始化为1的数组，然后对整个数组除以N，即可得到权重。如下所示：

1
2
3
4
5
import sys

# N = int(sys.argv[1])
N = 5
weights = np.ones(N)/ N

在N = 5时，输出结果如上。

(2) 使用这些权重值，调用convolve函数：

1
2
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
np.convolve(weights, c)[N-1 : -N+1]

注意：默认情况下，np.convolve(weights, c) 返回比 c 更长的数组：len(output)=len(c)+len(weights)−1

所以，必须裁剪结果。这就是N-1:-N+1的作用。

N-1：去掉 前 N-1 个无效值（初始点缺少足够的数据）。

-N+1：去掉 最后 N-1 个无效值（卷积计算超出 c 的范围）。

例如：

1
2
3
4
c = np.array([103, 104, 105, 106, 107])
N = 3
weights = np.ones(N) / N  # [1/3, 1/3, 1/3]
np.convolve(weights, c, mode='full')

(3) 我们从convolve函数返回的数组中，取出中间的长度为N的部分(即两者做卷积运算时完全重叠的区域)。下面的代码将创建 一个存储时间值的数组，并使用Matplotlib进行绘图。我们会在后续章节学习这个绘图库。

1
2
3
4
5
6
7
8
import matplotlib 
from matplotlib.pyplot import plot  
from matplotlib.pyplot import show 

t = np.arange(N-1, len(c))
plot(t, c[N-1:], lw=1.0) 
plot(t, sma, lw=2.0) 
show()

在下图中，相对较平滑的粗线描绘的是5日移动平均线，而锯齿状的细线描绘的是每天的收盘价。

刚才做了些什么 : 我们计算出了收盘价数据的简单移动平均线。是的，你掌握了很重要的知识，那就是简单移动平均线可以用信号处理技术求解——与1/N的权重进行卷积运算，N为移动平均窗口的大小。我们还学习了ones函数的用法，即可以创建元素均为1的数组，以及convolve函数，计算一组数据与指定权重的卷积。

完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import sys
import matplotlib 
from matplotlib.pyplot import plot  
from matplotlib.pyplot import show 

# N = int(sys.argv[1])
N = 5
weights = np.ones(N)/N
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
sma = np.convolve(weights, c)[N-1 : -N+1]
t = np.arange(N-1, len(c))
plot(t, c[N-1:], lw=1.0) 
plot(t, sma, lw=2.0) 
show()

3.21 指数移动平均线

除了简单移动平均线，指数移动平均线（exponential moving average）也是一种流行的技术指标。指数移动平均线使用的权重是指数衰减的。对历史上的数据点赋予的权重以指数速度减小，但永远不会到达0。我们将在计算权重的过程中学习exp和linspace函数。

3.22 动手实践：计算指数移动平均线

给定一个数组，exp函数可以计算出每个数组元素的指数。例如，看下面的代码：

1
2
n = np.arange(5)
np.exp(n)

linspace函数需要一个起始值和一个终止值参数，以及可选的元素个数的参数，它将返回一个元素值在指定的范围内均匀分布的数组。如下所示：(下面的0是包含的)

1
arr = np.linspace(-1, 0, 5)

下面我们来对示例数据计算指数移动平均线。 (1) 还是回到权重的计算——这次使用exp和linspace函数。

1
2
3
4
5
import sys
# N = int(sys.argv[1]) # ValueError: invalid literal for int() with base 10: '-f'
N = 5
weights = np.exp(np.linspace(-1, 0, N))
weights

(2) 对权重值做归一化处理。我们将用到ndarray对象的sum方法。

1
2
3
# 归一化
weights /= weights.sum()
weights

(3) 接下来就很容易了，我们只需要使用在简单移动平均线一节中学习到的convolve函数即可。同样，我们还是将结果绘制出来。

1
2
3
4
5
6
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
ema = np.convolve(weights, c)[N-1:-N+1]
t = np.arange(N-1, len(c))
plot(t, c[N-1:], lw=1.0) 
plot(t, ema, lw=2.0) 
show()

刚才做了些什么 : 我们对收盘价数据计算了指数移动平均线。首先，我们使用exp和linspace函数计算出指数衰减的权重值。linspace函数返回的是一个元素值均匀分布的数组，随后我们计算出它们的指数。为了将这些权重值归一化，我们调用了ndarray对象的sum方法。最后，我们再次应用了在前面简单移动平均线一节中学习到的convolve函数，最终计算出指数移动平均线。

完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import sys
import matplotlib 
from matplotlib.pyplot import plot  
from matplotlib.pyplot import show 

N = 5
weights = np.exp(np.linspace(-1, 0, N))
# 归一化
weights /= weights.sum()
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
ema = np.convolve(weights, c)[N-1:-N+1]
t = np.arange(N-1, len(c))
plot(t, c[N-1:], lw=1.0) 
plot(t, ema, lw=2.0) 
show()

3.23 布林带

布林带（Bollinger band）又是一种技术指标。是的，股票市场的确有成千上万种技术指标。 布林带是以发明者约翰·布林格（John Bollinger）的名字命名的，用以刻画价格波动的区间。布林带的基本型态是由三条轨道线组成的带状通道（中轨和上、下轨各一条）。

• 中轨 简单移动平均线。
• 上轨 比简单移动平均线高两倍标准差的距离。这里的标准差是指计算简单移动平均线所用数据的标准差。
• 下轨 比简单移动平均线低两倍标准差的距离。

3.24 动手实践：绘制布林带

我们已经掌握了计算简单移动平均线的方法。如有需要，请复习3.20节“动手实践：计算简单移动平均线”的内容。接下来的例子将介绍NumPy中的fill函数。fill函数可以将数组元素的值全部设置为一个指定的标量值，它的执行速度比使用array.flat = scalar或者用循环遍历数组赋值的方法更快。按照如下步骤绘制布林带。

(1) 我们已经有一个名为sma的数组，包含了简单移动平均线的数据。因此，我们首先要遍历和这些值有关的数据子集。数据子集构建完成后，计算其标准差。注意，从某种意义上来说，我们必须去计算每一个数据点与相应平均值之间的差值。如果不使用NumPy，我们只能遍历所有的数据点并逐一减去相应的平均值。幸运的是，NumPy中的fill函数可以构建元素值完全相同的数组。这可以让我们省去一层循环，当然也就省去了这个循环内作差的步骤。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import sys

# N = int(sys.argv[1]) # ValueError: invalid literal for int() with base 10: '-f'
N = 5
weights = np.ones(N)/N

for i in range(N - 1, C):
    if i + N < C:
        dev = c[i: i + N]
    else:
        dev = c[-N:]
    sma = np.convolve(weights, c)[N-1 : -N+1]
    averages = np.zeros(N)
    averages.fill(sma[i - N - 1])
    dev = dev - averages     
    dev = dev ** 2     
    dev = np.sqrt(np.mean(dev))     
    deviation.append(dev) 
    
deviation = 2 * np.array(deviation)  
upperBB = sma + deviation  
lowerBB = sma - deviation

(2) 使用如下代码绘制布林带（不必担心，我们将在第9章中学习绘图方面的知识）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import matplotlib 
from matplotlib.pyplot import plot  
from matplotlib.pyplot import show 

c_slice = c[N-1:]
between_bands = np.where((c_slice < upperBB) & (c_slice > lowerBB))

between_bands = len(np.ravel(between_bands))

t = np.arange(N - 1, C)  
plot(t, c_slice, lw=1.0)  
plot(t, sma, lw=2.0)  
plot(t, upperBB, lw=3.0)  
plot(t, lowerBB, lw=4.0)  
show() 

上图是用我们的示例数据绘制出来的布林带。中间锯齿状的细线描绘的是每天的收盘价，而稍微粗一点也平滑一点的穿过它的曲线即为简单移动平均线。

刚才做了些什么 ： 我们在示例数据上计算得到了布林带。更重要的是，我们还了解了NumPy中fill函数的用法。该函数可以用一个指定的标量值填充数组，而这个标量值也是fill函数唯一的参数。

完整代码如下：(补充了注释，可以详细看看)

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import sys
import matplotlib 
from matplotlib.pyplot import plot  
from matplotlib.pyplot import show 

deviation = []
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
# 计算数据点的总数
C = len(c)

# N = int(sys.argv[1]) # ValueError: invalid literal for int() with base 10: '-f'
N = 5  # 移动平均的窗口大小设定为 5
weights = np.ones(N)/N  # 创建一个长度为 N 的均等权重数组 [1/N, 1/N, ..., 1/N]，用于计算移动平均

for i in range(N - 1, C):  # 这个循环从 N-1（即第一个完整窗口）到 C（数据长度）遍历,dev 选取窗口内的 N 个数据
    if i + N < C:
        dev = c[i: i + N]  # 取 N 个窗口内的数据
    else:
        dev = c[-N:] # 取最后 N 个数据
    sma = np.convolve(weights, c)[N-1 : -N+1]  # 计算移动平均 （SMA）,进行一维卷积，相当于对 c 进行滑动窗口均值计算,[N-1 : -N+1] 提取有效部分，去掉边界影响
    averages = np.zeros(N)
    averages.fill(sma[i - N - 1])  # averages.fill(sma[i - N - 1]) 创建一个 N 个相同值的数组，每个值都是该窗口对应的 sma 值
    dev = dev - averages  # 计算窗口内每个点与 SMA 之间的差距    
    dev = dev ** 2      # 计算平方
    dev = np.sqrt(np.mean(dev))     # 计算标准差（均方根） 
    deviation.append(dev) # 将标准差存入 deviation 列表
    
deviation = 2 * np.array(deviation)  # 计算 2 倍标准差（标准布林带的默认系数）
upperBB = sma + deviation  # 计算上轨
lowerBB = sma - deviation  # 计算下轨

# 计算布林带范围内的数据点

c_slice = c[N-1:]   # 对齐数据，去掉前面 N-1 个点
between_bands = np.where((c_slice < upperBB) & (c_slice > lowerBB))  # 找出落在布林带之间的点的索引

between_bands = len(np.ravel(between_bands))  # 将索引数组展平,计算落在布林带之间的数据点数量

# 绘制布林带
t = np.arange(N - 1, C)  # 生成横坐标数据
plot(t, c_slice, lw=1.0)   # 绘制原始数据
plot(t, sma, lw=2.0)  # 绘制移动平均线
# 绘制布林带上轨和下轨
plot(t, upperBB, lw=3.0)  
plot(t, lowerBB, lw=4.0)  
show() 

勇敢出发：转换为指数移动平均线

人们通常将简单移动平均线作为布林带的中轨线。而以指数移动平均线作为中轨线也是一种流行的做法，因此我们将它留作练习。如果需要提示，你可以在本章中找到合适的示例。 验证一下fill函数的执行速度是否真的比使用array.flat = scalar或者用循环遍历数组赋值的方法更快。

是的，通常布林带（Bollinger Bands, BB）采用 简单移动平均线（SMA） 作为中轨线。但在一些场景下，使用 指数移动平均线（EMA, Exponential Moving Average） 作为中轨线也是一种常见的变体。这是因为 EMA 对最近的数据点赋予更高的权重，使得它对价格变化更加敏感。

将 SMA 替换为 EMA，然后计算标准差并绘制布林带。

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
import sys
import matplotlib 
from matplotlib.pyplot import plot  
from matplotlib.pyplot import show 

deviation = []
c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6,))
# 计算数据点的总数
C = len(c)

# N = int(sys.argv[1]) # ValueError: invalid literal for int() with base 10: '-f'
N = 5  # 移动平均的窗口大小设定为 5
# weights = np.ones(N)/N  # 创建一个长度为 N 的均等权重数组 [1/N, 1/N, ..., 1/N]，用于计算移动平均
weights = np.exp(np.linspace(-1, 0, N))# 指数移动平均线的权重
# 权重归一化
weights /= weights.sum()

for i in range(N - 1, C):  # 这个循环从 N-1（即第一个完整窗口）到 C（数据长度）遍历,dev 选取窗口内的 N 个数据
    if i + N < C:
        dev = c[i: i + N]  # 取 N 个窗口内的数据
    else:
        dev = c[-N:] # 取最后 N 个数据
    # sma = np.convolve(weights, c)[N-1 : -N+1]  # 计算移动平均 （SMA）,进行一维卷积，相当于对 c 进行滑动窗口均值计算,[N-1 : -N+1] 提取有效部分，去掉边界影响
    ema = np.convolve(weights, c)[N-1:-N+1]   # 计算指数移动（EMA）
    averages = np.zeros(N)
    averages.fill(ema[i - N - 1])  # averages.fill(sma[i - N - 1]) 创建一个 N 个相同值的数组，每个值都是该窗口对应的 sma 值
    dev = dev - averages  # 计算窗口内每个点与 SMA 之间的差距    
    dev = dev ** 2      # 计算平方
    dev = np.sqrt(np.mean(dev))     # 计算标准差（均方根） 
    deviation.append(dev) # 将标准差存入 deviation 列表
    
deviation = 2 * np.array(deviation)  # 计算 2 倍标准差（标准布林带的默认系数）
upperBB = ema + deviation  # 计算上轨
lowerBB = ema - deviation  # 计算下轨

# 计算布林带范围内的数据点

c_slice = c[N-1:]   # 对齐数据，去掉前面 N-1 个点
between_bands = np.where((c_slice < upperBB) & (c_slice > lowerBB))  # 找出落在布林带之间的点的索引

between_bands = len(np.ravel(between_bands))  # 将索引数组展平,计算落在布林带之间的数据点数量

# 绘制布林带
t = np.arange(N - 1, C)  # 生成横坐标数据
plot(t, c_slice, lw=1.0)   # 绘制原始数据
plot(t, ema, lw=2.0)  # 绘制移动平均线
# 绘制布林带上轨和下轨
plot(t, upperBB, lw=3.0)  
plot(t, lowerBB, lw=4.0)  
show() 

注：上面蓝色的线是原始数据线，在两幅图中是一样的。

还有一个问题：验证一下fill函数的执行速度是否真的比使用array.flat = scalar或者用循环遍历数组赋值的方法更快。

下面是一个简单的例子：

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
31
32
import numpy as np
import time

# 设置数组大小
size = 10**6

# 创建一个全零数组
arr1 = np.zeros(size)
arr2 = np.zeros(size)
arr3 = np.zeros(size)

# 定义要填充的值
value = 5.0

# 使用 fill() 方法
start_time = time.time()
arr1.fill(value)
fill_time = time.time() - start_time

# 使用 flat 赋值方法
start_time = time.time()
arr2.flat = value
flat_time = time.time() - start_time

# 使用循环遍历赋值
start_time = time.time()
for i in range(size):
    arr3[i] = value
loop_time = time.time() - start_time

# 结果输出
fill_time, flat_time, loop_time

结论

1. fill() 方法最快，其优化可能源于 NumPy 内部使用高效的 C 代码直接操作内存。
2. flat = scalar 方法次之，但仍然比 fill() 慢，可能是由于 flat 属性涉及迭代操作。
3. 循环遍历赋值最慢，因为它涉及 Python 级别的循环，每次赋值都有额外的开销。

如果要填充 NumPy 数组，推荐使用 fill() 方法！

3.25 线性模型

许多科学研究中都会用到线性关系的模型。NumPy的linalg（linear algebra的缩写词）包是专门用于线性代数计算的。 下面的工作基于一个假设，就是一个价格可以根据N个之前的价格利用线性模型计算得出。

3.26 动手实践：用线性模型预测价格

我们姑且假设，一个股价可以用之前股价的线性组合表示出来，也就是说，这个股价等于之前的股价与各自的系数相乘后再做加和的结果，这些系数是需要我们来确定的。用线性代数的术语来讲，这就是解一个最小二乘法的问题。步骤如下。

(1) 首先，获取一个包含N个股价的向量b。

数据读取（读取收盘价这一列）：

1
2
3
import sys

b = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6, ))

1
2
3
4
# N = int(sys.argv[1])
N = 5
b = b[-N:]    # 取最新的N个数据
b = b[::-1]   # 倒着排（第一个是最新的）

(2) 第二步，初始化一个N×N的二维数组A，元素全部为0。

1
A = np.zeros((N, N))

(3) 第三步，用b向量中的N个股价值填充数组A。

1
2
for i in range(N):
    A[i, ] = c[-N-1-i:-i-1]

由于 A[i, :] 每次都被填充一个包含 c 中某个子集的数据，在 for 循环结束时，A 将会是一个 N x N 的矩阵，包含了从 c 数组中按上述方式选择的每一行数据。

(4) 我们的目标是确定线性模型中的那些系数，以解决最小平方和的问题。我们使用linalg包中的lstsq函数来完成这个任务。

1
(x, residuals, rank, s) = np.linalg.lstsq(A, b)

返回的元组中包含稍后要用到的系数向量x、一个残差数组、A的秩以及A的奇异值。

补充：np.linalg.lstsq(A, b)：这是 NumPy 中的 最小二乘法 解法，用于求解线性方程组 A * x = b，即找到向量 x，使得：\(min_x\vert \vert Ax-b \vert \vert ^2\)。

lstsq 是 Least Squares（最小二乘）的缩写。

该函数返回一个元组 (x, residuals, rank, s)，其中：

• x 是最小二乘解（我们需要的目标向量）。
• residuals 是残差，即 A * x - b 的平方和，表示解的误差。
• rank 是矩阵 A 的秩，表示矩阵的线性独立性。
• s 是矩阵 A 的奇异值，提供了矩阵的稳定性和条件信息。

(5) 一旦得到了线性模型中的系数，我们就可以预测下一次的股价了。使用NumPy中的dot函数计算系数向量与最近N个价格构成的向量的点积（dot product）。

1
np.dot(b, x)

np.dot(b, x) 计算 b 和 x 的点积，即： \(\text{dot product} = b_1 \cdot x_1 + b_2 \cdot x_2 + \cdots + b_N \cdot x_N\)，这个值通常用于表示向量 b 和 x 的相似性。

查了一下记录，下一个交易日实际的收盘价为353.56。因此，我们用N = 5做出的预测结果并没有差得很远。

刚才做了些什么 : 我们预测了明天的股价。如果这真的有效，我们就可以提早退休了！你看，买这本书是多么正确的投资！我们为股价预测建立了一个线性模型，于是这个金融问题就变成了一个线性代数问题。NumPy中的linalg包里有一个lstsq函数，帮助我们求出了问题的解——即估计线性模型中的系数。在得到解之后，我们将系数应用于NumPy中的dot函数，通过线性回归的方法预测了下一次的股价。

完整代码如下：

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

b = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(6, ))
# N = int(sys.argv[1])
N = 5
b = b[-N:]    # 取最新的N个数据
b = b[::-1]   # 倒着排（第一个是最新的）

A = np.zeros((N, N))
for i in range(N):
    A[i, ] = c[-N-1-i:-i-1]
    
(x, residuals, rank, s) = np.linalg.lstsq(A, b)
np.dot(b, x)

补充：为什么矩阵 A 是按这种方式构造的，特别是通过 从数组 c 中选取特定子集 来填充矩阵。其实，这种构造矩阵的方法通常出现在 线性回归、时间序列建模、最小二乘法拟合 或 平滑 等问题中。让我们逐步分析：

1. 为什么构造矩阵 A？

矩阵 A 是为了解线性方程 A * x = b，它用于最小二乘法（Least Squares）问题。在这个问题中，我们需要通过 A 和 b 来拟合一个模型，通常这是一个回归问题，目标是找到一个最优解向量 x，使得 A * x 尽可能接近 b。

假设：我们要用 A 来拟合一个基于过去 N 个数据点的模型，并预测当前或未来的值。因此，A 的每一行应当代表过去 N 个数据点的 一个时间窗口。

2. 为什么每行包含 c[-N-1-i:-i-1] 这种形式的切片？

时间序列建模和滞后特性

在许多时间序列分析中，特别是在 回归分析 或 平滑 中，常常需要用到 滞后数据（Lagged Data）来预测当前或未来的值。通过使用滞后数据，我们试图找到过去的趋势、模式或关联，以便对当前的输出做出预测。

假设 c 是某种时间序列数据（例如股价、温度、销售量等）。如果我们希望基于 最近的 N 个数据点 来预测当前数据，那么每个时间点的预测通常会考虑其 前几个时间点的值。

为何选取 c[-N-1-i:-i-1]？

• -N-1-i 和 -i-1 代表的就是在时间序列中 从倒数第 N+1 项到倒数第 i+1 项 的一段数据。这种选择确保了对于每一行（即时间窗口），我们都能提取过去的 N 个数据点。我们从数据的末尾（最近的）开始逆向选择，以便构建符合时序的模型。
• 逐步递减：对于矩阵 A 的每一行，i 的变化（从 0 到 N-1）表示我们通过向后滑动时间窗口来填充数据。也就是说，第一行使用最近的 N 个数据点，第二行使用倒数第二个时间点开始的 N 个数据点，依此类推。

3. 为什么矩阵 A 的每行是过去的 N 个数据点？

• 回归模型的要求：通常情况下，我们需要通过过去的数据来预测当前或未来的值。在 线性回归 中，目标是找到一组 权重，使得某些 特征（自变量） 的线性组合能尽可能准确地拟合目标值（因变量）。在这里，我们使用过去的 N 个数据点（即历史数据）作为特征，用它们来预测未来的值。
• 时间窗口（Sliding Window）：每一行代表一个时刻的 时间窗口，即过去 N 个数据点。滑动窗口的方法将允许我们利用每一时刻的历史数据进行预测。这是非常常见的一种时间序列预测方法。

4. 矩阵 A 的作用

• 矩阵 A 的每一行包含 过去 N 个数据点，这将允许我们通过最小二乘法（Least Squares）找到一个合适的线性模型来描述数据的趋势。

矩阵 A 的构建方式实际上反映了一种常见的 时间序列建模 技术，通过使用过去 N 个数据点来预测未来的值。在这个过程中：

• 每一行代表一个时间窗口，即过去 N 天的数据。
• 切片 c[-N-1-i:-i-1] 是为了从时间序列中提取合适的滞后数据，用作模型的输入特征。
• 通过这种方式，我们希望拟合一个模型，该模型可以根据过去的 N 个数据点来预测当前或未来的值。

3.27 趋势线

趋势线，是根据股价走势图上很多所谓的枢轴点绘成的曲线。顾名思义，趋势线描绘的是价格变化的趋势。过去的股民们在纸上用手绘制趋势线，而现在我们可以让计算机来帮助我们作图。在这一节的教程中，我们将用非常简易的方法来绘制趋势线，可能在实际生活中不是很奏效，但这应该能将趋势线的原理阐述清楚。

3.28 动手实践：绘制趋势线

按照如下步骤绘制趋势线。 (1) 首先，我们需要确定枢轴点(pivots)的位置。这里，我们假设它们等于最高价、最低价和收盘价的算术平均值。

data.csv有以下几列：第一列为股票代码以标识股票（苹果公司股票代码为AAPL），第二列为dd-mm-yyyy格式的日期，第三列为空，随后各列依次是开盘价、最高价、最低价和收盘价，最后一列为当日的成交量。

下面为一行数据：

1
AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800 

1
2
3
# 读数据：最高价、最低价和收盘价
h, l, c = np.loadtxt("data.csv", delimiter=",", unpack=True, usecols=(4,5,6))
pivots = (h + l + c)/3

从这些枢轴点出发，我们可以推导出所谓的阻力位和支撑位。

• 阻力位是指股价上升时遇到阻力，在转跌前的最高价格；
• 支撑位是指股价下跌时遇到支撑，在反弹前的最低价格。

需要提醒的是，阻力位和支撑位并非客观存在，它们只是一个估计量。基于这些估计量，我们就可以绘制出阻力位和支撑位的趋势线。我们定义当日股价区间为最高价与最低价之差。

(2) 定义一个函数用直线\(y= at + b\)来拟合数据，该函数应返回系数a和b。这里需要再次用到linalg包中的lstsq函数。将直线方程重写为\(y = Ax\)的形式，其中A = [t 1]，x = [a b]。使用ones_like和vstack函数来构造数组A。

lstsq 是 Least Squares（最小二乘）的缩写。

1
2
3
def fit_line(t, y):
    A = np.vstack(t, np.ones_like(t)).T
    return np.lstsq(A, y)[0]

详细解释这个步骤的数学原理。

1. 目标

本步骤的目标是 用一条直线 \(y= at + b\)来 拟合数据，并求得 系数\(a\)和\(b\)，以便用于趋势预测。这本质上是一个 线性回归 问题，即： \(y= at + b\)，需要找到最佳的\(a\)和\(b\)使得误差最小。

2. 线性方程的矩阵形式

方程\(y= at + b\)可以转换为矩阵形式：

用矩阵表示为：\(Ax=y\)，其中：\(A\)是 设计矩阵（Design Matrix），定义为：\(A = \begin{bmatrix} t_1 & 1 \\ t_2 & 1 \\ t_3 & 1 \\ . & . \\ . & . \\ . & . \\ t_n & 1 \end{bmatrix}\)，\(x\)是 需要求解的参数向量：\(x = \begin{bmatrix} a \\ b \end{bmatrix}\)，\(y\)是 目标变量：\(y = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ . \\ . \\ . \\ y_n \end{bmatrix}\)。

任务：求解\(x = \begin{bmatrix} a, b \end{bmatrix}^T\)，即最优的斜率和截距。

np.vstack([t, np.ones_like(t)])

• t 是时间序列 t = np.arange(len(c))，即 [0, 1, 2, ..., n-1]。

• np.ones_like(t) 创建一个全 1 的数组，与 t 形状相同，用于表示偏置项（截距）。

• vstack（垂直堆叠）后：

  • ```shell A = [[0 1] [1 1] [2 1] [3 1] … [n-1 1]]

  ```

• A.T 之后变成：

  • ```shell [[0 1 2 3 … n-1] [1 1 1 1 … 1 ]]

  ```

• 这样 A 就成为了 时间变量 t 与常数项 1 组成的矩阵。

np.linalg.lstsq(A, y, rcond=None)[0]使用 最小二乘法 计算线性拟合的参数 \(x\)，x[0] = a 是斜率，x[1] = b 是截距。

(3) 假设支撑位在枢轴点下方一个当日股价区间的位置，而阻力位在枢轴点上方一个当日股价区间的位置，据此拟合支撑位和阻力位的趋势线。

1
2
3
4
5
t = np.arange(len(c))  
sa, sb = fit_line(t, pivots - (h - l))  # 支撑位
ra, rb = fit_line(t, pivots + (h - l))  # 阻力位
support = sa * t + sb  
resistance = ra * t + rb 

pivots - (h - l) 计算支撑位。

pivots + (h - l) 计算阻力位。

fit_line() 用最小二乘法分别计算支撑位和阻力位的 趋势线参数 a 和 b。

support：支撑位趋势线 y = sa * t + sb。

resistance：阻力位趋势线 y = ra * t + rb。

在 股市技术分析 中：

• 支撑位 代表价格下跌时可能的反弹点。
• 阻力位 代表价格上涨时可能的回调点。
• 拟合趋势线 旨在找出股价变动的 长期趋势，而非短期波动。

通过最小二乘法：

• 计算 y = at + b 形式的直线，使其最接近计算出的支撑位和阻力位。
• 这样，我们可以用趋势线预测 未来的支撑位和阻力位。

(4) 到这里我们已经获得了绘制趋势线所需要的全部数据。但是，我们最好检查一下有多少个数据点落在支撑位和阻力位之间。显然，如果只有一小部分数据在这两条趋势线之间，这样的设定就没有意义。设置一个判断数据点是否位于趋势线之间的条件，作为where函数的参数。

1
2
condition = (c > support) & (c < resistance) 
between_bands = np.where(condition)

以下是根据条件判断的布尔值：

复查一下具体取值：

注意，where函数返回的是一个秩为2的数组，因此在使用len函数之前需要调用ravel函数。

1
between_bands = len(np.ravel(between_bands))

我们还得到了一个额外的奖励：一个新的预测模型。我们可以用这个模型来预测下一个交易日的阻力位和支撑位。

1
2
sa * (t[-1] + 1) + sb
ra * (t[-1] + 1) + rb 

此外，还有另外一种计算支撑位和阻力位之间数据点个数的方法：使用[]和intersect1d函数。在[]操作符里面定义选取条件，然后用intersect1d函数计算两者相交的结果。

1
2
a1 = c[c > support]  
a2 = c[c < resistance]

(5) 我们再次将结果绘制出来，如下所示：

1
2
3
4
plot(t, c)  
plot(t, support)  
plot(t, resistance)  
show() 

绘制结果如下图所示，其中包含了股价数据以及对应的支撑位和阻力位。

刚才做了些什么 : 我们用NumPy画出了趋势线，省去了用尺、铅笔和绘图纸的麻烦。我们定义了一个用直线拟合数据的函数，其中用到NumPy中的vstack、ones_like和lstsq函数。拟合数据是为了得到支撑位和阻力位两条趋势线的方程。随后，我们用两种不同的方法分别计算了有多少个数据点落在支撑位和阻力位之间的范围内，并得到了一致的结果。

第一种方法使用where函数和一个条件表达式。第二种方法使用[]操作符和intersect1d函数。intersect1d函数返回一个由两个数组的所有公共元素构成的数组。

3.29 ndarray 对象的方法

NumPy中的ndarray类定义了许多方法，可以在数组对象上直接调用。通常情况下，这些方法会返回一个数组。你可能已经注意到了，很多NumPy函数都有对应的相同的名字和功能的ndarray对象。这主要是由NumPy发展过程中的历史原因造成的。

NumPy 是 Python 中用于科学计算的核心库之一，它最初的实现是在 Numeric 和 Numarray 这两个库的基础上。后来，这两个库合并并演变成了今天的 NumPy。NumPy 的设计目标是提供一个高效、灵活的多维数组对象及其操作，同时保持良好的性能和兼容性。

早期的 Numeric 库定义了基础的多维数组对象，提供了基本的数学操作。随着科学计算的需求增加，Numarray 引入了更多的高级功能，特别是支持较大的数组以及更多的数据类型。在这两个库的基础上，NumPy 于 2006 年诞生，结合了两者的优点，并在很多方面做了优化和改进。

ndarray 是 NumPy 中最核心的类，它表示一个多维的数组对象。ndarray 类有大量的 实例方法，这些方法可以直接在 ndarray 对象上调用，例如 sum()、mean()、reshape()、transpose() 等。

我们注意到，许多 NumPy 的函数和 ndarray 对象的实例方法是 功能相同 的。这是因为，在 NumPy 的设计中，有两个并行的接口：

模块级函数接口：像 np.sum()、np.mean() 等是作为模块级的函数提供的。

对象方法接口：像 ndarray 对象的 sum()、mean() 等是作为 ndarray 类的方法提供的。

这两种接口实现了相同的功能，但是它们的使用方式不同。

这种设计的历史原因与 NumPy 的发展历程 有关。早期的 Numeric 和 Numarray 库分别为数组提供了不同的接口，有时会使用函数，有时会使用对象方法。随着 NumPy 的发展，NumPy 逐渐将这两者合并，但为了兼容历史和提供更灵活的接口，它保留了 函数接口 和 对象方法接口 两种调用方式。

提供两种接口（函数和方法）是为了让用户能够根据自己的需要选择使用哪一种，更加灵活。

这种重复设计的目的就是为了兼容不同的使用方式和编程风格，同时保证 代码的简洁性 和 性能的高效性。

ndarray对象的方法相当多，我们无法在这里逐一介绍。前面遇到的var、sum、std、argmax、argmin以及mean函数也均为ndarray方法。 数组的修剪和压缩请参见下一节中的内容。

3.30 动手实践：数组的修剪和压缩

这里给出少量使用ndarray方法的例子。按如下步骤对数组进行修剪和压缩操作。 (1) clip方法返回一个修剪过的数组，也就是将所有比给定最大值还大的元素全部设为给定的最大值，而所有比给定最小值还小的元素全部设为给定的最小值。例如，设定范围1到2对0到4 的整数数组进行修剪：

1
2
a = np.arange(5)
a.clip(1, 2)

(2) compress方法返回一个根据给定条件筛选后的数组。例如：

1
2
3
a = np.arange(4)
a.compress(a > 2)
a.compress(a > 0)

注：看下面的例子

compress(condition, axis=None) 方法要求 condition 必须是一维数组，但 n > 5 生成的是一个 与 n 形状相同的布尔数组，它是 二维的，导致 compress() 无法正确执行。

---

正确的例子：

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

n = np.random.randint(0, 10, (3, 4))  # 生成 3x4 的随机整数数组
print("n:\n", n)

# 生成一个长度为 3 的布尔数组，指定要保留的行
condition = np.array([True, False, True])

# 过滤行（axis=0 表示按行过滤）
result = np.compress(condition, n, axis=0)
print("Filtered Rows:\n", result)

3.31 阶乘

许多程序设计类的书籍都会给出计算阶乘的例子，我们应该保持这个传统。

3.32 动手实践：计算阶乘

ndarray类有一个prod方法，可以计算数组中所有元素的乘积。按如下步骤计算阶乘。

(1) 计算8的阶乘。为此，先生成一个1~8的整数数组，并调用prod方法。

1
2
n = np.arange(1, 9)
n.prod()

这很不错，但如果我们想知道1~8的所有阶乘值呢？ (2) 没问题！调用cumprod方法，计算数组元素的累积乘积。

1
n.cumprod()

刚才做了些什么 : 我们使用prod和cumprod方法计算了阶乘。

3.33 本章小结

本章我们学习了很多常用的NumPy函数。我们用loadtxt读文件，用savetxt写文件，用eye函数创建单位矩阵，用loadtxt函数从一个CSV文件中读取股价数据。NumPy中的average 和mean函数可以用来计算数据的加权平均数和算术平均数。

本章还提到了一些常用的统计函数。首先，我们使用min和max函数来确定股价的范围；然后，用median函数获取数据的中位数；最后，用std和var函数计算数据的标准差和方差。

diff函数可以返回数组中相邻元素的差值，因此我们用它来计算股票的简单收益率。log函数可以计算数组元素的自然对数。

loadtxt函数默认将所有数据转换为浮点数类型，它有一个特定的参数可以完成转换。这个参数就是converters，它是一个可以将数据列和所谓的转换函数连接起来的参数。

我们自定义了一个函数并将其作为参数传给了apply_along_axis函数。我们实现了一个可以在多个数组间找出每个位置上最大元素的算法。

我们了解到ones函数可以创建一个全为1的数组，而且convolve函数可以根据指定的权重计算卷积。

我们用exp和linspace函数得到了一组指数衰减的权重值。linspace可以给出一个均匀分布的数组，然后我们计算出该数组元素的指数。

我们还调用ndarray类的sum方法对权重值做归一化处理。

我们还接触到了fill函数，这个函数可以用一个指定的标量值填充数组，而这个标量值也是其唯一的参数。

结束了本章的NumPy常用函数之旅，我们将来到NumPy便捷函数的世界。下一章将学习polyfit、sign和piecewise等NumPy函数。

• Previous
  numpy(2)_numpy基础
• Next
  numpy(4)便捷函数