【NumPy学习指南】笔记--第三章 常用函数

概述

这一章学习了NumPy的常用函数，主要涵盖以下内容：

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

下表列举这一章用到的NumPy常用函数以及ndarray对象的函数：

接下来过一遍本章的程序代码：

首先导入需要用到的包

import numpy as np

import matplotlib.pyplot as plt

% matplotlib inline

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签

plt.rcParams['axes.unicode_minus']=False #用来正常显示负号

读写文件

创建一个单位矩阵并将其存储到文件中

i2 = np.eye(2) # 创建一个二维单位矩阵数组

print(i2)

np.savetxt('eye.txt',i2) # 用savetxt函数进行写文件的操作

读入CSV文件

载入苹果公司的历史股价数据为例。股价数据存储在CSV中，第一列为股票代码，第二列为dd-mm-yyyy格式日期，第三列为空，随后各列依次为开盘价、最高价、最低价和收盘价，最后一列为当日的成交量。

# 用loadtxt函数读取csv文件，用delimiter参数指定文件中的分隔符

# 用usecols参数指定要导入的数据列，将unpack参数设置为True使不同列的数据分开存储

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

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

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

vmap = np.average(c, weights=v) 

print("VWAP = " + str(vmap))

算术平均值

mean = np.mean(c)

print ("算术平均值mean = " + str(mean))

时间加权平均价格（TWAP）

在经济学中，TWAP是另一种“平均”价格的指标。基本思想就是最近的价格重要性大一些，所以应对近期的价格给以较高的权重。最简单的方法是用arange函数创建一个从0开始依次增长的自然数序列，自然数的个数即为收盘价的个数。这并不一定是正确的计算TWAP的方式，示例仅仅是为了说明问题。

t = np.arange(len(c))

twap = np.average(c, weights=t)

print("TWAP = " + str(twap))

最大值和最小值

# 载入每日最高价和最低价的数据

h, l = np.loadtxt('data.csv', delimiter=',', usecols=(4,5), unpack=True)

# 获取价格区间

highest = np.max(h)

lowest = np.min(l)

print("最高价格highest =  " + str(highest))

print("最低价格lowest =  " + str(lowest))

# ptp函数计算数组的取值范围，返回的是数组元素的最大值和最小值之间的差值

spread_high_price = np.ptp(h)

spread_low_price = np.ptp(l)

print("最高价之间的差值spread_high_price =  " + str(spread_high_price))

print("最低价之间的差值spread_low_price =  " + str(spread_low_price))

简单统计分析

# 调用median函数计算收盘价的中位数

median = np.median(c)

print("收盘价的中位数median = " + str(median))

数学上常规方法计算中位数：

      首先对数组进行从小到大排序；

      如果数组长度为奇数，中位数就是中间位的数；

      如果数组长度为偶数，中位数的值等于中间那两个数的平均值。

# 对价格数组进行排序

sorted_c = np.msort(c)

print("排序后的价格数组sorted = " + str(sorted_c))

# 根据中位数定义求中位数

N = len(c)

if N % 2 == 1:   # 如果长度为奇数

    middle = sorted_c[(N - 1)/2]

else:   # 长度为偶数

    middle = (sorted_c[N / 2] + sorted_c[(N-1) / 2]) / 2

print("中位数middle = " + str(middle))

方差：能体现变量变化的程度。在示例中，方差体现为投资风险的大小。

numpy文档中方差的定义：各个数据与所有数据算术平方数的离差平方和除以数据个数所得到的值。

一般统计学书上样本方差的定义：离差平方和/(n-1) ,其中n为样本个数。

之所以有这样的差别，是为了保证样本方差是一个无偏估计量，但numpy文档中方差的定义在机器学习领域更为通用。

# 根据方差函数计算方差

variance = np.var(c)

print("根据函数计算收盘价的方差variance = " + str(variance))

# 根据定义计算方差

variance_from_definition = np.mean((c - np.mean(c)) ** 2)

print("根据定义计算收盘价的方差variance_from_definition = " + str(variance_from_definition))

股票收益率

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

简单收益率是指相邻两个价格之间的变化率。

对数收益率是指所有价格取对数后两两之间的差值。对数收益率也可以用来衡量价格的变化率。

收益率是无量纲的。投资者感兴趣的是收益率的方差或标准差，因为这代表这投资风险的大小。

# 计算简单收益率。

# diff函数返回一个由相邻数组元素的差值构成的数组

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

# 用std函数计算标准差

standard_deviation = np.std(returns)

print("简单收益率标准差standard_deviation = " + str(standard_deviation))

# 计算对数收益率

# 一般情况下需要检查输入数组以确保其不含零和负数，在该例中股价总为正值，可省略检查。

logreturns = np.diff(np.log(c)) 

# 我们对哪些交易日的收益率为正值感兴趣

# where函数可以根据指定的条件返回所有满足条件的数组元素的索引值。

posretindices = np.where(returns > 0) 

print("收益率为正的索引： " + str(posretindices))

波动率是对价格变动的一种度量。年波动率等于对数收益率的标准差除以其均值，再除以交易日倒数的平方根，通常交易日取252天。

# 计算波动率

annual_volatility = np.std(logreturns) / np.mean(logreturns)

annual_volatility = annual_volatility / np.sqrt(1. / 252.)

print("年波动率annual_volatility = " + str(annual_volatility))

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

print("月波动率monthly_volatility = " + str(monthly_volatility))

日期分析

Numpy是面向浮点数运算的，因此需要对日期做一些专门的处理。我们需要显式地告诉Numpy怎样来转换日期，需要用到loadtxt函数中converters参数，该参数是数据列和转换函数之间进行映射的字典。

写出转换函数：

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

import datetime

# 将日期作为字符串传给datestr2num函数，按照指定的形式"%d-%m-%Y"转换成一个datetime对象，

# 随后datetime对象被转换为date对象，最后调用weekday方法返回一个数字

def datestr2num(s):

    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()

# 将日期转换函数挂接上去，读入数据

dates, o = np.loadtxt('data.csv', delimiter=',', usecols=(1,3), converters={1:datestr2num}, unpack=True)

print("日期Dates = " + str(dates))

# 创建一个包含5个元素的数组，分别代表一周的5个工作日，数组初始化为0，这个数组将用于保存各工作日的平均收盘价。

averages = np.zeros(5)

# 遍历星期一到星期五

for i in range(5):

    indices = np.where(dates == i)  # 用where函数得到各工作日的索引值并存储在indices数组中

    prices = np.take(c,indices)     # 再用take函数获取这些索引值相应的元素值

    avg = np.mean(prices)

    print("Day "+str(i)+" prices "+str(prices)+" Average "+str(avg))

    averages[i] = avg  # 对每个工作日计算出的平均值存放在averages数组中

# 找出哪个工作日的平均收盘价最高，哪个最低

top = np.max(averages)

print("最高工作日平均收盘价top = " + str(top))

# argmax函数返回的是averages数组中最大元素的索引值

print("平均收盘价最高的工作日是： " + str(np.argmax(averages)))   

bottom = np.min(averages)

print("最低工作日平均收盘价bottom = " + str(bottom))

# argmin函数返回的是averages数组中最小元素的索引值

print("平均收盘价最低的工作日是： " + str(np.argmin(averages)))

汇总数据

汇总整个交易周中从周一到周五的所有数据。

# 只考虑前三周的数据，避免节假日造成的数据缺失。

c = c[:16]

dates = dates[:16]

# 找到第一个星期一

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

print('第一个星期一的索引：' + str(first_monday))

# 找到最后一个星期五

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

print('最后一个星期五的索引：' + str(last_friday))

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

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

print("周索引初始化： " + str(weeks_indices))

# 按照每个子数组5个元素，用split函数切分数组

weeks_indices = np.split(weeks_indices, 3)

print("分组之后的周索引： " + str(weeks_indices))

# 编写summarize函数。该函数将为每一周的数据返回一个元组，包含这一周的开盘价，最高价，最低价和收盘价

def summarize(a, o, h, l, c):

    monday_open = o[a[0]]  # 一周的开盘价即为周一的开盘价

    # 用take函数来根据索引值获取数组元素的值，并用max函数计算出一周的最高股价

    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)

# apply_along_axis函数会调用summarize函数，分别作用于3个星期数组，数组元素中的索引值对应于示例中的一天

weeksummary = np.apply_along_axis(summarize, 1, weeks_indices, o, h, l, c)

print("周总结：" + str(weeksummary))

# 用savetxt函数，指定文件名，需要保存的数组名，分隔符以及存储浮点数的格式

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

格式字符串的说明

格式字符串的格式为：%[标志字符][输出宽度参数][精度格式符]字符类型

其中[]部分为可选部分，以下对各个部分进行说明。

标志字符：

输出宽度参数：表示输出的最小位数。

精度格式符：以“.”开头，后面跟一个表示精度的整数。

字符类型：

真实波动幅度均值（ATR）

ATR是一个用来衡量股价波动性的技术指标。

ATR是基于N个交易日的最高价和最低价进行计算的，通常取最近20个交易日。

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

# 获取近N个交易日的最高价和最低价

N = int(raw_input(u"请输入一个正整数：\n"))

h = h[-N:]

l = l[-N:]

print("最高价h的长度：" + str(len(h)) + " 最低价l的长度：" + str(len(l)))

print("收盘价c：" + str(c))

# 获取前一个交易日的收盘价

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

print("前一个交易日的收盘价previousclose = " + str(previousclose))

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



h-l 当日股价范围，即当日最高价和最低价之差；



h-previousclose 当日最高价和前一个交易日收盘价之差；



previousclose-l 前一个交易日收盘价和当日最低价之差'''

# 返回一组数组之间挑选每一个元素位置上的最大值

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

print("真实波动truerange：" + str(truerange))

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

atr = np.zeros(N)

# 数组首元素就是truerange数组元素的平均值

atr[0] = np.mean(truerange)

# 其他元素的计算公式为((N-1)PATR+TR)/N，PATR表示前一个交易日的ATR值，TR即当日的真实波动幅度

for i in range(1, N):

    atr[i] = (N - 1) * atr[i-1] +truerange[i]

    atr[i] /= N

print("ATR = " + str(atr))

简单移动平均线

简单移动平均线通常用于分析时间序列上的数据。

简单移动平均线可以用信号处理技术求解，即与1/N的权重进行卷积运算，N为移动平均窗口的大小，也可以不等权重。

卷积定义为一个函数与经过翻转和平移的另一个函数的乘积的积分。

# 创建一个长度为N的数组，并初始化元素的值为1，除以N，得到权重

N = int(raw_input(u"请输入一个正整数：\n"))

weights = 1.0 * np.ones(N) / N

print("权重weights = " +str(weights))

# 使用权重值，调用convolve函数，获取两者做卷积运算时完全重叠的区域

sma = np.convolve(weights, c)[N-1: -N+1]

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

plt.plot(t, c[N-1:], lw=1.0, label="close price")  # 每天的收盘价

plt.plot(t, sma, lw=2.0, label=str(N) + 'days sma')  # N日简单移动平均线

plt.legend(loc="upper left")

plt.show()

指数移动平均线

指数移动平均线使用的权重是指数衰减的。对历史上的数据点赋予的权重以指数速度减小，但永远不会到达0.

# linspace函数需要一个起始值和一个终止值参数，以及可选的元素个数的参数，返回一个元素值在指定的范围内均匀分别的数组

weights = np.exp(np.linspace(-1., 0., N))

# 对权重值做归一化处理

weights /= weights.sum()

print("权重weights = " + str(weights))

ema = np.convolve(weights, c)[N-1: -N+1]

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

plt.plot(t, c[N-1:], lw=1.0, label="close price")

plt.plot(t, ema, lw=2.0, label=str(N) + "days ema")

plt.legend(loc="upper left")

plt.show()

布林带

布林带的基本形态是由三条轨道线组成的带状通道(中轨和上、下轨各一条)。

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

deviation = []

C = len(c)

# 计算标准差

for i in range(N-1, C):

    if i+N < C:

        dev = c[i: i+N]

    else:

        dev = c[-N:]

    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

c_slice = c[N-1:]

between_bands = np.where((c_slice < upperBB) & (c_slice > lowerBB))

# print(lowerBB[between_bands])

# print(c[between_bands])

# print(upperBB[between_bands])

between_bands = len(np.ravel(between_bands))

print("布林带：" + str(float(between_bands)/len(c_slice)))

t = np.arange(N-1, C)

plt.figure(figsize=(12,6))

plt.plot(t, c_slice, lw=1.0, label="close price")

plt.plot(t, sma, lw=2.0, label="sma")

plt.plot(t, upperBB, lw=3.0, label="upperBB")

plt.plot(t, lowerBB, lw=4.0, label="lowerBB")

plt.legend(loc="upper left")

plt.show()

用线性模型预测价格

# 获取一个包含N个股价的向量b

b = c[-N:]

b = b[::-1]

print("b: " + str(b))

# 初始化一个N*N的二维数组A，元素全部为0

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

# 用b向量中的N个股价值填充数组A

for i in range(N):

    A[i,] = c[-N-1-i: -1-i]

print("A: " + str(A))

# 使用linalg包中的lstsq函数来确定线性模型中的系数,原理是最小二乘

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

print("系数向量: " + str(x))

print("残差数组: " + str(residuals))

print("A的秩: " + str(rank))

print("A的奇异值: " + str(s))

# 使用dot函数计算系数向量与最近N个价格构成的向量的点积

pre_price = np.dot(b,x)

print("预测下一个交易日的收盘价为：" + str(pre_price))

绘制趋势线

阻力位：股价上升时遇到阻力，在转跌前的最高价格；

支撑位：股价下跌时遇到支撑，在价格反弹前的最低价格。

阻力位和支撑位只是一个估计量。

# 确定枢轴点的位置，假设为最高价、最低价和收盘价的算术平均值

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

pivots = (h + l + c)/3.0

print("枢轴点：" + str(pivots))

# 定义一个函数用直线y=ax+b来拟合数据，该函数应返回系数a和b。

# 将直线方程重写y=Ax的形式，其中A =[t 1]，x=[a b]

def fit_line(t, y):

    A = np.vstack([t, np.ones_like(t)]).T

    return np.linalg.lstsq(A, y)[0]

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

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

第一种方法计算有多少数据点落在支撑位和阻力位之间的范围内：使用where函数和一个条件表达式。

# 设置一个判断数据点是否位于趋势线之间的条件

condition = (c > support) & (c < resistance)

print("条件：" + str(condition))

between_bands = np.where(condition)

# 复查具体取值

# print support[between_bands]

# print c[between_bands]

# print resistance[between_bands]

between_bands = len(np.ravel(between_bands))

print ("位于两条趋势线之间的数据点的个数：" + str(between_bands))

print ("位于两条趋势线之间的数据比例：" + str(float(between_bands)/len(c)))

# 用模型来预测下一个交易日的阻力位和支撑位

next_support = sa * (t[-1] + 1) +sb

next_resistance = ra * (t[-1] +1) + rb

print("下一个交易日的支撑位：" + str(next_support))

print("下一个交易日的阻力位：" + str(next_resistance))

另一种计算支撑位和阻力位之间数据点个数的方法：使用[]操作符和intersect1d函数。

intersect1d函数返回一个由两个数组的所有公共元素构成的数组。

a1 = c[c > support]

a2 = c[c < resistance]

between_bands2 = len(np.intersect1d(a1, a2))

print("位于两条趋势线之间的数据点的个数: " + str(between_bands2))

plt.plot(t, c)

plt.plot(t, support)

plt.plot(t, resistance)

plt.show()

数组的修剪和压缩

clip方法返回一个修剪过的数组，也就是将所有比给定最大值还大的元素全部设为给定的最大值，而所有比给定最小值还小的元素全部设为给定的最小值。

a = np.arange(5)

print("a = " + str(a))

print("修剪后: " + str(a.clip(1, 3)))

compress方法返回一个根据给定条件筛选后的数组

a = np.arange(5)

print("a = " + str(a))

print("压缩后：" + str(a.compress(a>2)))

连乘

b = np.arange(1, 9)

print("b = " + str(b))

print("8的阶乘：" + str(b.prod()))

print("数组元素的累积乘积：" + str(b.cumprod()))

以上都是在jupyter notebook中运行，chap03.常用函数.ipynb和所需数据文件可见这里