In [2]:
%matplotlib inline
import numpy as np
from matplotlib.pyplot import plot
from matplotlib.pyplot import show
In [3]:
# 首先读入两只股票的收盘价，并计算收益率
bhp_cp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
vale_cp = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp_cp) / bhp_cp[:-1]
vale_returns = np.diff(vale_cp) / vale_cp[:-1]
协方差描述的是两个变量共同变化的趋势，其实就是归一化前的相关系数。 使用cov函数计算股票收益率的协方差矩阵
In [4]:
covariance = np.cov(bhp_returns, vale_returns)
print 'Covariance:\n', covariance
In [5]:
# 查看协方差矩阵对角线的元素
print 'Covariance diagonal:\n', covariance.diagonal()
# 计算矩阵的迹，即对角线之和
print 'Covariance trace:\n', covariance.trace()
# 计算相关系数，相关系数是协方差除以各自标准差的乘积
print 'Correlation coefficient:\n', covariance / (bhp_returns.std() * vale_returns.std())
用相关系数来度量两只股票的相关程度。相关系数的取值范围在-1到1之间，一组数据域自身的相关系数为1.使用corrcoef函数计算相关系数。
In [6]:
# 使用corrcoef计算更加精确
print 'Correlation coefficient:\n', np.corrcoef(bhp_returns, vale_returns)
相关系数矩阵关于对角线对称，BHP与VALE的相关系数等于VALE和BHP的相关系数。看起来0.68的相关系数表示他们的相关程度似乎不是很强。
In [7]:
difference = bhp_cp - vale_cp
avg = np.mean(difference)
dev = np.std(difference)
# 检查最后一次收盘价是否在同步状态
print "Out of sync : ", np.abs(difference[-1] - avg) > 2*dev
这说明，最后一次收盘价不再同步状态，我们暂时不能进行交易
In [8]:
# 绘制收益率曲线
t = np.arange(len(bhp_returns))
plot(t, bhp_returns, lw=1)
plot(t, vale_returns, lw=2)
show()
In [9]:
# 用三次多项式去拟合两只股票收盘价的差价
t = np.arange(len(bhp_cp))
poly = np.polyfit(t, bhp_cp-vale_cp, 3)
print "Polynomial fit\n", poly
In [10]:
# 用刚才得到的多项式对象，推断下一个值
print "Next value: ", np.polyval(poly, t[-1]+1)
理想情况下，BHP和VALE股票收盘价的差价越小越好。在极限情况下，差值可以在某个点为0。用roots函数找到拟合多项式函数在什么时候达到0。
In [11]:
print "Roots: ", np.roots(poly)
In [12]:
# 极值位于导数为0的点
der = np.polyder(poly)
print "Dervative:\n", der
# 得到多项式导函数的系数
In [13]:
# 求出导函数的根，即找出原多项式函数的极值点
print "Extremas: ", np.roots(der)
In [14]:
# 通过argmax和argmin函数找到最大最小值点来检查结果
vals = np.polyval(poly, t)
print "Maximum index: ", np.argmax(vals)
print "Minimum index: ", np.argmin(vals)
In [15]:
plot(t, bhp_cp-vale_cp)
plot(t, vals)
show()
我们需要在成交量前面乘上一个有收盘价变化决定的正负号。
In [16]:
cp, volume = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,7), unpack=True)
change = np.diff(cp)
print "Change:", change
使用NumPy的sign函数返回每个元素的正负号。
In [17]:
signs = np.sign(change)
print "Signs:\n", signs
使用Numpy的piecewise函数获取数组元素的正负。piecewise(分段)，可以根据给定取值，得到分段。
In [18]:
pieces = np.piecewise(change, [change<0, change>0], [-1,1])
print "Pieces:\n", pieces
In [19]:
# 检查两次输出是否一致
print "Arrays equal?", np.array_equal(signs, pieces)
In [20]:
# OBV值的计算依赖于前一日的收盘价
print "On balance volume: \n", volume[1:]*signs
使用vectorize函数可以减少你的程序中使用循环的次数，NumPy中的vectorize函数相当于Python中的map函数。我们用它来计算单个交易日的利润。
In [21]:
# 读入数据
# op is opening price,hp is the highest price
# lp is the lowest price, cp is closing price
op, hp, lp, cp = np.loadtxt('BHP.csv', delimiter=',', usecols=(3,4,5,6), unpack=True)
我们尝试以比开盘价稍低一点的价格买入股票。如果这个价格不在当日的股价范围内，则尝试买入失败，没有获利，也没有亏损，我们返回0。否则，我们将以当日收盘价卖出，所获得的利润即买入卖出的差价。
In [22]:
def calc_profit(op, high, low, close):
# 以开盘价买入,这里不考虑买入多少股
buy = op
if low < buy < high:
return (close-buy) / buy
else:
return 0
In [23]:
# 矢量化一个函数，这样可以避免使用循环
func = np.vectorize(calc_profit)
profits = func(op, hp, lp, cp)
print 'Profits:\n', profits
In [24]:
# 我们选择非零利润的交易日并计算平均值
real_trades = profits[profits != 0]
print 'Number of trades:\n', len(real_trades), round(100.0 * len(real_trades)/len(cp), 2),"%"
print "Average profit/loss % :", round(np.mean(real_trades) * 100, 2)
In [25]:
# 选择正盈利的交易日并计算平均利润
winning_trades = profits[profits > 0]
print "Number of winning trades", len(winning_trades), round(100.0*len(winning_trades)/len(cp),2),"%"
print "Average profit %", round(np.mean(winning_trades) *100, 2)
In [26]:
# 选择负盈利的交易日并计算平均损失
losing_trades = profits[profits < 0]
print "Number of winning trades", len(losing_trades), round(100.0*len(losing_trades)/len(cp),2),"%"
print "Average profit %", round(np.mean(losing_trades) *100, 2)
hanning函数式一个加权余弦的窗函数，我们使用hanning函数平滑股票收益率的数组。
离散卷积运算函数convolve的文档convolve函数文档
In [27]:
# 调用hanning函数计算权重，生成一个长度为N的窗
# 这里N为8
N = 8
weights = np.hanning(N)
print "Weights:\n", weights
In [28]:
# 首先读入两只股票的收盘价，并计算收益率
bhp_cp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
vale_cp = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp_cp) / bhp_cp[:-1]
vale_returns = np.diff(vale_cp) / vale_cp[:-1]
# convolve函数，离散线性卷积运算
smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns)[N-1 : -N+1]
smooth_vale = np.convolve(weights/weights.sum(), vale_returns)[N-1 : -N+1]
In [29]:
from matplotlib.pyplot import legend
t = np.arange(N-1, len(bhp_returns))
plot(t, bhp_returns[N-1:], lw=1.0, label='bhp returns')
plot(t, smooth_bhp, lw=2.0, label='smooth bhp')
plot(t, vale_returns[N-1:], lw=1.0, label='vale returns')
plot(t, smooth_vale, lw=2.0, label='smooth vale')
legend(loc='best')
show()
图中折线有交叉，这些交叉点可能是股价趋势的转折点，至少可以表明BHP和VALE之间的股价关系发生了变化。这些转折点可能会经常出现，我们可以利用他们预测未来的股价走势。
In [30]:
import matplotlib.pyplot as plt
# 使用多项式拟合平滑后的数据
K = 5
t = np.arange(N-1, len(bhp_returns))
poly_bhp = np.polyfit(t, smooth_bhp, K)
poly_vale = np.polyfit(t, smooth_vale, K)
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(t, smooth_bhp, label="smooth bhp")
poly_bhp_value = np.polyval(poly_bhp, t)
ax1.plot(t, poly_bhp_value, label='poly bhp')
plt.legend()
ax2 = fig.add_subplot(212)
ax2.plot(t, smooth_vale, label="smooth vale")
poly_vale_value = np.polyval(poly_vale, t)
ax2.plot(t, poly_vale_value, label='poly vale')
plt.legend()
show()
In [31]:
# 得到交叉点的x坐标
# 通过求多项式函数差，再求根
poly_sub = np.polysub(poly_bhp, poly_vale)
xpoints = np.roots(poly_sub)
print "Intersection points:", xpoints
In [32]:
# 判断是否为实数
# select选出实数
reals = np.isreal(xpoints)
print "Real number?",reals
xpoints = np.select([reals], [xpoints])
xpoints = xpoints.real
print "Real intersection points:", xpoints
In [33]:
# 去除0元素
# trim_zeros函数可以去掉一维数组中开头和末尾为0的元素
print "Sans 0s", np.trim_zeros(xpoints)
In [ ]: