Numpy 确定稳态向量|极客教程

确定稳态向量，马尔科夫链（Markov chain）被用来表示至少有两个状态的系统。有关马尔科夫链的详细信息，请参阅 http://en.wikipedia.org/wiki/Markov_chain。此类系统t时刻的状态仅取决于t-1时刻的状态。系统的当前状态在这些状态之间随机地切换。如果把一支股票的股价变动情况定义为一个马尔科夫链，并定义持平F、上涨U和下跌D这三个状态，则我们可以根据当日收盘价确定其稳态。

在未来某个时刻之后，或者从理论上讲经过无限长时间之后，马尔科夫链系统的状态将不再改变。这个状态被称为稳态。随机矩阵包含了状态之间转移的概率。把随机矩阵A应用于稳态，状态将保持不变，用数学符号表示如下：
确定稳态向量
也可以把稳态看作特征值为1的特征向量。

具体步骤

现在我们需要先获取数据。

获取一年的数据。

获取数据的一种办法是使用Matplotlib（如有必要，请参阅1.5节）。我们将获取从现在算起一年内的数据，相关代码如下。

today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo('AAPL', start, today)

选取收盘价。

我们已经从雅虎财经频道获得了需要的历史数据。这些数据被表示为一个由元组构成的列表，但我们只关心元组中包含的收盘价。历史数据列表的例子如下。

[(734744.0, 675.25, 673.47000000000003, 677.66999999999996, 
672.60000000000002, 7228500.0), …,(734745.0,670.63999999999999, 
663.87, 671.54999999999995, 662.85000000000002, 10799600.0)]

元组中的第一个元素代表日期，后面依次是开盘价、最高价、最低价和收盘价，最后一个元素是成交量。可以用如下方式选择收盘价。

close =  [q[4] for q in quotes]

收盘价是每个元组中的第五个数。现在我们获得了一个由大约253个收盘价构成的列表。

确定状态数组。

使用NumPy的diff函数获得两个相邻交易日的收盘价差，进而确定状态数组。状态数组中的每个状态由价差的符号决定。当NumPy的sign函数的输入参数为负值、正值和零时，其返回结果分别为-1、1和0。

states = numpy.sign(numpy.diff(close))

初始化随机矩阵。

对于每一次状态转换，其起始状态和结束状态都有三种可能。例如，如果起始状态为U，则可能切换到：

状态U
状态F
状态D

使用NumPy的zeros函数，对随机矩阵进行初始化。

NDIM = 3
SM = numpy.zeros((NDIM, NDIM))

选取每一种符号对应的起始状态索引。

代码现在变得有点复杂了，我们不得不使用循环语句。对每一种可能的符号，使用NumPy的where函数选取其对应的起始状态索引。这里的k是一个平滑常数，后面将讨论其作用。

signs = [-1, 0, 1]
k = int(sys.argv[2])

for i in xrange(len(signs)):
    #从特定符号对应的状态开始转换
    start_indices = numpy.where 
        (states[:-1] == signs[i])[0]

平滑处理和随机矩阵。

给定起始状态，可以对每种类型转换发生的次数进行计数。把计数结果除以该起始状态对应的总的转换次数，就能获得随机矩阵中相应的转换概率。顺便说一下，这不是最佳的做法，因为可能会造成过拟合。
在现实生活中，相邻的两个交易日的收盘价有可能相同。但对于易变的股票市场而言，这种可能性很小。对于转换次数为零的情况，一种应对方法是使用加法平滑。基本思路是在转换次数上加一个常数，从而避免出现概率为0的情况。计算随机矩阵中各个元素的代码如下。

N = len(start_indices) + k * NDIM

# 跳过总转换次数为零的情况
if N == 0:
    continue

#获取结束状态的状态值
end_values = states[start_indices + 1]

for j in xrange(len(signs)):
    # 特定类型转换发生的次数 
    occurrences = len 
        (end_values[end_values == signs[j]])
    SM[i][j] = (occurrences + k)/float(N)
print SM

上述代码的功能是，对每一种可能的转换类型，基于转换发生的次数和加法平滑技术，计算其转换概率。
使用AAPL（苹果公司）的数据和平滑常数k=1，得到如下随机矩阵。

[[ 0.50925926  0.00925926  0.48148148]
 [ 0.33333333  0.33333333  0.33333333]
 [ 0.35135135  0.00675676  0.64189189]]

特征值和特征向量。

为了能得到特征值和特征向量，我们需要使用NumPy的linalg模块和eig函数。

eig_out = numpy.linalg.eig(SM)
print eig_out

eig函数返回的是一个特征值数组和一个特征向量数组。

(array([ 1.        ,  0.15817566,  0.32630882]),  
 array([[ 0.57735027,  0.74473695,  0.00297158],
        [ 0.57735027, -0.39841481, -0.99983179],
        [ 0.57735027, -0.53538072,  0.01809841]]))

选择特征值1对应的特征向量。
现在我们只对特征值1对应的特征向量感兴趣。实际上，特征值可能不会精确地等于1，因此应该设立一个误差边界。可以寻找取值介于0.9和1.1之间的特征值对应的特征向量的索引，代码如下。

idx_vec = numpy.where 
    (numpy.abs(eig_out[0] - 1)

本章的完整代码如下。

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy
import sys

if (len(sys.argv) != 3):
    print "Usage python %s SYMBOL k" % (sys.argv[0])
    print "For instance python %s AAPL 1" % (sys.argv[0])
    sys.exit()

today = date.today()
start = (today.year - 1, today.month, today.day)

quotes = quotes_historical_yahoo(sys.argv[1], start, today)
close =  [q[4] for q in quotes]

states = numpy.sign(numpy.diff(close))

NDIM = 3
SM = numpy.zeros((NDIM, NDIM))

signs = [-1, 0, 1]
k = int(sys.argv[2])

for i in xrange(len(signs)):
    #从特定符号对应的状态开始转换
    start_indices = numpy.where(states[:-1] == signs[i])[0]

    N = len(start_indices) + k * NDIM

    # 跳过总转换次数为零的情况
    if N == 0:
        continue

    #获取结束状态的状态值
    end_values = states[start_indices + 1]

    for j in xrange(len(signs)):
        # 特定类型转换发生的次数 
        occurrences = len(end_values[end_values == signs[j]])
        SM[i][j] = (occurrences + k)/float(N)

print SM
eig_out = numpy.linalg.eig(SM)
print eig_out

idx_vec = numpy.where(numpy.abs(eig_out[0] - 1)

这段代码的输出结果如下所示。

[[ 0.4952381   0.00952381  0.4952381 ]
 [ 0.33333333  0.33333333  0.33333333]
 [ 0.34210526  0.00657895  0.65131579]]
(array([ 1.        ,  0.15328174,  0.32660547]),  
 array([[ 0.57735027,  0.7424435 ,  0.00144451],
        [ 0.57735027, -0.44112882, -0.99982343],
        [ 0.57735027, -0.50416566,  0.01873551]]))
Index eigenvalue 1 (array([0]),)
Steady state vector [ 0.57735027  0.57735027  0.57735027]
Check [ 0.57735027  0.57735027  0.57735027]