In [1]:
import numpy as np
import scipy as sp
import scipy.signal as ssig
import bokeh as bk
import bokeh.plotting as bplt
import bokeh.palettes as bpal
In [2]:
# 工具配置
TOOLS = "pan,wheel_zoom,box_zoom,box_select"
# 作图配置
plot_config = dict(plot_width=800, plot_height=500, tools=TOOLS)
bplt.output_notebook()
假设我们不断的投资年利率永远为5%的国债。
这句话的意思是说，该种国债在0时刻购买时付出1元钱，而在1年后(时刻为1)按约定取回1.05元。
我们可以反复的在第T年买入1元钱的国债，再在第T+1年获得相应的0.05元收益，那么我们的资金曲线就是线性增长的如下图所示：
In [3]:
capital_curve_1 = np.linspace(1,2,21)
print capital_curve_1
p1 = bplt.figure(title="等量购买国债的资金曲线",
background_fill="#FFFFFF",**plot_config)
p1.line(np.arange(21), capital_curve_1, line_width=3, color='steelblue', legend='等差数列')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以年为单位记)'
p1.yaxis.axis_label = '资金量'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [4]:
capital_curve_2 = np.array([1.05**i for i in range(21)])
print capital_curve_2
p1 = bplt.figure(title="再投资式购买国债的资金曲线",
background_fill="#FFFFFF",**plot_config)
p1.line(np.arange(21), capital_curve_1, line_width=3, color='steelblue', legend='等差数列')
p1.line(np.arange(21), capital_curve_2, line_width=3, color='darkred', legend='等比数列')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以年为单位记)'
p1.yaxis.axis_label = '资金量'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [5]:
capital_curve_1 = np.linspace(1,6,101)
capital_curve_2 = np.array([1.05**i for i in range(101)])
p1 = bplt.figure(title="再投资式购买国债的资金曲线（100年）",
background_fill="#FFFFFF",**plot_config)
p1.line(np.arange(101), capital_curve_1, line_width=3, color='steelblue', legend='等差数列')
p1.line(np.arange(101), capital_curve_2, line_width=3, color='darkred', legend='等比数列')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以年为单位记)'
p1.yaxis.axis_label = '资金量'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [6]:
p1 = bplt.figure(title="再投资式购买国债的资金曲线（100年，对数坐标）",
background_fill="#FFFFFF",y_axis_type="log",**plot_config)
p1.line(np.arange(101), capital_curve_1, line_width=3, color='steelblue', legend='等差数列')
p1.line(np.arange(101), capital_curve_2, line_width=3, color='darkred', legend='等比数列')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以年为单位记)'
p1.yaxis.axis_label = '资金量(对数坐标)'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [7]:
for x in range(3,24):
print '{: >3}% -- {: >6.2f}{: >6.2f}'.format(x,72.0/x,np.log(2)/np.log(1.0+x/100.0))
当你对某种属性不变资产进行投资，花了一定时间从10万元积累到了100万元，那么在同样条件下花同样的时间可以从100万元累计到1000万元（比例相等）
从10万元累计到1000万元翻了100倍，在同样外界条件下需要花的时间是2倍。这个计算等价于，以10倍为底求100的对数就是2。
In [8]:
np.log10(100)
Out[8]:
In [9]:
gauss_series = np.arange(100)+1
print gauss_series.sum(),(1+100)*100/2
等比数列求和公式：前后项之比为x时可以利用公式：(x不等于1时始终成立)
$$\sum_{k=0}^n x^k = 1+x+x^2+...x^n = \frac{1-x^{n+1}}{1-x}$$刚才我们提到的数列的性质还是比较简单，因为:
具有第一个性质的数列，我们统称为"一阶递推数列"。
而我们把第二个性质稍为改一改，如果同时乘以前一项的常数，再加上第二个常数，把乘法和加法混在一起后，事情就复杂了。
这时数列不再是等差数列或者等比数列那么简单了...但只要我们能理解递推公式，也不是很复杂。
In [10]:
A = np.float64(np.zeros_like(np.arange(21)))
A[0] = 1.0
for i in range(1,len(A)):
A[i] = A[i-1] + 0.05
print A
我们换个写法，把递推公式改成通项公式，也就是适合于每一项的公式：$A_{n} = A_{0}+0.05n$，我们可以自己验证结果的正确性。
同时我们如果去做如下的验证：（如下两点便称为数学归纳法）
验证通项公式代入递推公式：
$$A_{n+1}= A_0+0.05(n+1)\\ A_n +0.05 = A_0+0.05n+0.05 = A_0 + 0.05(n+1)$$这时我们发现，第0项，第1项，第2项...一直递推下去，所有项都满足条件。我们不使用循环一样可以得到任意长的关于A的序列。
可见，从递推公式出发，使用数学归纳法，求得通项公式是
In [11]:
print np.array([0.05*i for i in range(21)])+1.0
同理我们可以照猫画虎的对等比数列使用同样的伎俩：
等比数列$1,1.05,1.1025,...$的递推公式是$A_{n+1} = A_{n}*1.05$，如果用递推公式仍然需要循环：
In [12]:
A = np.float64(np.zeros_like(np.arange(21)))
A[0] = 1.0
for i in range(1,len(A)):
A[i] = A[i-1]*1.05
print A
但是经过简单的对通项公式的推理之后，得到$A_{n} = A_0 * 1.05^n$ ：
In [13]:
print np.array([1.05**i for i in range(21)])
In [14]:
A = np.float64(np.zeros_like(np.arange(361)))
A[0] = 0.0
for i in range(1,len(A)):
A[i] = A[i-1]*1.004+1.0
print A[:20]
In [24]:
p1 = bplt.figure(title="理财定投购买货币基金的资金曲线（30年）",
background_fill="#FFFFFF",**plot_config)
p1.line(np.arange(361), A, line_width=3, color='steelblue', legend='定投资金曲线')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以月为单位记)'
p1.yaxis.axis_label = '资金量(线性坐标)'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [16]:
p1 = bplt.figure(title="理财定投购买货币基金的资金曲线（30年，对数坐标）",
background_fill="#FFFFFF",y_axis_type="log",**plot_config)
p1.line(np.arange(1,361), A[1:], line_width=3, color='steelblue', legend='定投资金曲线')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以月为单位记)'
p1.yaxis.axis_label = '资金量(对数坐标)'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
这种奇葩的一阶线性递推数列，我们习惯性的问：怎么求通项？
处理这种问题最好的方法就是，两侧加上一个常数，使两侧的常数项之比刚好等于1.004。
$$A_0 = 0\\ A_{n+1} = A_{n}*1.004+1\\ A_{n+1}+250 = A_{n}*1.004+251 = 1.004*(A_{n}+250)$$这时这个奇葩露出了真面目：A+250是一个以250为初始项的等比数列。
$A_{n} = 250*1.004^n-250$
比较一下。
In [17]:
print (250.0*np.array([1.004**i for i in range(361)])-250.0)[-20:]
print A[-20:]
举个简单的例子，如果我们已经有购房的打算，并且准备每个月还款1万(本金+利息)，分30年还清。
如果简单假设月利率一直保持0.4%，并且要达到等额本息贷款（也就是每个月还的利息＋本金为定值）的效果，每个月该还多少钱呢？
听上去十分复杂，实际我们可以通过建立简单的金融知识体系来理解这一过程：
现金流：一组在不同时间点上发生的现金收入或者支出称为现金流。
首先我们换个思考问题的方式：假如我们没有做贷款，我们也要每个月投入等额的钱以0.4%的利息利滚利，最后在30年之后获得一笔钱，这不就是定投吗？
这说明一件事情，如下三种现金流中，所有的现金支出和现金收入的价值相等
In [18]:
t = np.arange(361)
cf1 = np.array([A[-1]*(1.004**(-360))]+[-1]*360)
cf2 = np.array([-1]*360+[A[-1]])
cf3 = np.array([-A[-1]*(1.004**(-360)),A[-1]])
In [23]:
p1 = bplt.figure(title="30年贷款并等额本息还款（可放大）",
background_fill="#FFFFFF",**plot_config)
p1.set(x_range=bk.models.Range1d(-5,20),y_range=bk.models.Range1d(-2,4))
for tp in t:
if cf1[tp]>0:
p1.line(np.array([tp,tp]), np.array([0,cf1[tp]]), line_width=3, color='steelblue', legend='获得贷款')
else:
p1.line(np.array([tp,tp]), np.array([0,cf1[tp]]), line_width=3, color='darkred', legend='等额本息还款')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以月为单位记)'
p1.yaxis.axis_label = '现金收支'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [21]:
p1 = bplt.figure(title="30年定投理财现金流（可放大）",
background_fill="#FFFFFF",**plot_config)
p1.set(x_range=bk.models.Range1d(340,365),y_range=bk.models.Range1d(-2,4))
for tp in t:
if cf2[tp]>0:
p1.line(np.array([tp,tp]), np.array([0,cf2[tp]]), line_width=3, color='steelblue', legend='银行账户结余')
else:
p1.line(np.array([tp,tp]), np.array([0,cf2[tp]]), line_width=3, color='darkred', legend='每月定投')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以月为单位记)'
p1.yaxis.axis_label = '现金收支'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
In [22]:
p1 = bplt.figure(title="30年定期按月复利储蓄",
background_fill="#FFFFFF",**plot_config)
p1.line([0,360],[0,0], line_width=1, color='black')
for tp in [0,360]:
if cf3[tp/360]>0:
p1.line(np.array([tp,tp]), np.array([0,cf3[tp/360]]), line_width=5, color='steelblue', legend='银行账户结余')
else:
p1.line(np.array([tp,tp]), np.array([0,cf3[tp/360]]), line_width=5, color='darkred', legend='储蓄金额')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = '时间(以月为单位记)'
p1.yaxis.axis_label = '现金收支'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
在这个例子中，如下几个定义分别有了具体的概念：
这三者拥有相同的价值！
我们在第三个例子中反推老爸给的钱：一笔未来值通过利率折算成为现值，就叫做折现。给定等比数列的第360项，求第0项：
In [65]:
A[-1]*(1.004**(-360))
Out[65]:
同样我们还可以在第一个例子中反推一组现金流折算的现值。360笔1元钱分别是不同等比数列的第1、第2、...第360项时，那些数列中第0项的总和：
In [66]:
np.array([1.004**(-i) for i in range(1,361)]).sum()
Out[66]:
可以看出两种折算确实是等价的。
考虑一个经典问题：Fibonacci（斐波那契）数列。
假设小张养的兔子都能连续生2年小兔，每一年生一对。假设新出生的小兔子到1岁的时候也能生了。
第0年时，小张有一对刚出生的小兔。于是第1年时，这对小兔生了一对刚出生的小兔。
第2年时，第0年的老兔和第1年的小兔同时生下了刚出生的小兔一共2对，以此类推。
第n年时，刚出生的小兔有多少对？如果想不清楚，我们仍然要先写好递推公式：
$F_{0} = 1 \\ F_{1} = 1 \\ F_{n+2} = F_{n+1} + F_{n}$
我们先用循环实现一遍：
In [23]:
N = 21
F = np.zeros_like(np.arange(N))
F[0] = 1
F[1] = 1
for i in range(2,N):
F[i] = F[i-1]+F[i-2]
print F
我们使用迭代函数实现一遍：
In [24]:
def Fibonacci(n):
N = np.int64(n)
if N<0:
return None
elif N in [0,1]:
return 1
else:
return Fibonacci(N-1)+Fibonacci(N-2)
F = [Fibonacci(i) for i in range(21)]
print F
当然，在Python语言中，为了将这个功能函数化我不会推荐第一种方法；
而如果关心性能，第二种方法又明显不合适（太多次调用自己引起栈溢出），当n太大的时候程序就会报错了。
对于这种迭代类型的函数还可以进行进一步的优化，直接不加修改地给出一个好思路，原始接链在http://code.activestate.com/recipes/474088/
In [42]:
import sys
class TailRecurseException:
def __init__(self, args, kwargs):
self.args = args
self.kwargs = kwargs
def tail_call_optimized(g):
def func(*args, **kwargs):
f = sys._getframe()
if f.f_back and f.f_back.f_back \
and f.f_back.f_back.f_code == f.f_code:
raise TailRecurseException(args, kwargs)
else:
while 1:
try:
return g(*args, **kwargs)
except TailRecurseException, e:
args = e.args
kwargs = e.kwargs
func.__doc__ = g.__doc__
return func
@tail_call_optimized
def fib(i, current = 1, next = 1):
if i == 0:
return current
else:
return fib(i - 1, next, current + next)
print [fib(i) for i in range(21)]
通过程序解决问题当然是我们强有力的一种手段，但往往结合数学知识与计算机知识会让事情变的简单。
根据问题定义：
$F_{0} = 1 \\ F_{1} = 1 \\ F_{n+2} = F_{n+1} + F_{n}$
我们称$F_{n}$为$F_{n+1}$的一阶滞后项，称$F_{n}$为$F_{n+2}$的二阶滞后项。
$$F_{n+2} = B F_{n+2} + B^2 F_{n+2}$$简化：$$F = B*F + B^2*F$$
求解初中数学：$$(1-B-B^2)F = 0\\ (1-\frac{1-\sqrt{5}}{2}B)(1-\frac{1+\sqrt{5}}{2}B)F = 0$$
然后我们就得到了两个仅含一阶滞后算子的解，每个解跟等比数列的形式都一样：
$$(1-\frac{1-\sqrt{5}}{2}B)F = 0, (1-\frac{1+\sqrt{5}}{2}B)F = 0$$通项公式：
$$F_{n} = a*(\frac{1-\sqrt{5}}{2})^n + b*(\frac{1＋\sqrt{5}}{2})^n$$这时如果我们知道
$$F_{0} = a＋b ＝ 1\\ F_{1} = a*(\frac{1-\sqrt{5}}{2})+b*(\frac{1＋\sqrt{5}}{2})＝1 $$根据这个二元一次方程组，我们先用一些作弊的方式求解（更具体的方式会在矩阵基础部分介绍）
In [30]:
A = np.mat([[1,1],[(1-np.sqrt(5))/2,(1+np.sqrt(5))/2]])
y = np.array([1,1])
res = A.I.dot(y)
print res
好吧我知道看着比较不带感是因为它毕竟不是分数而是小数，总之我们还是有方法知道：
$a = \frac{1-1/\sqrt{5}}{2}\\ b = \frac{1＋1/\sqrt{5}}{2}$
顺手验证一下：
In [43]:
print (1-1/np.sqrt(5))/2 , (1+1/np.sqrt(5))/2
最后见证奇迹的时刻到来了：现在我们有精确表达式了！
$F_{n} = a*(\frac{1-\sqrt{5}}{2})^n + b*(\frac{1＋\sqrt{5}}{2})^n$
把$a,b$代入之后究竟是什么呢？
In [44]:
def fakefib(n):
a = (1-1/np.sqrt(5))/2
b = (1+1/np.sqrt(5))/2
return a*(((1-np.sqrt(5))/2)**n)+b*(((1+np.sqrt(5))/2)**n)
print [fakefib(i) for i in range(21)]
因为计算精度问题，显得不够整洁！那么四舍五入以及取整。
In [35]:
print [np.int64(np.round(fakefib(i))) for i in range(21)]
In [46]:
1/((1+np.sqrt(5))/2)
Out[46]:
这就是我们熟悉的黄金分割数。
在时间序列的模型中，如果我们再遇到有一阶或者高阶线性递推数列，相关的简单分析我们已经能够在滞后算子B的辅助下找到数列的通项公式了。
能够良好的理解滞后算子B，对于时间序列分析是相当有帮助的。
一个序列$C = [C_1,C_2,...]$已经准备好了，而计算相应的K序列时，
$$K_0 = 50\\ K_{n+1}= \frac{2}{3}K_{n}+\frac{1}{3}C_{n+1}$$这样，KDJ指标的递推公式与初始值就被定义好了。使用滞后算子，会得到更好的形式：
$$(1-\frac{2}{3}B)K = \frac{1}{3} C$$如果你觉得最后一个公式十分眼熟的话，你可以和SciPy部分signal(信号处理)的内容进行类比，因为这完全满足ARMA模型。
我们把$K_0 = 50$的初始条件直接替换为$C_0 = 150$，我们尝试一下原有C序列用迭代法或时间序列法能否得到一致的结果：
In [53]:
C = np.array([50,70,32,58,0,23,99])
# 方法1
K1 = np.float64(np.zeros_like(C))
K1[0] = 2.0/3*50+1.0/3*C[0]
for i in range(1,len(C)):
K1[i] = K1[i-1]*2.0/3+1.0/3*C[i]
print K1
在方法2里我们把K和C序列关系的参数，以及增补过的C序列放入滤波器里，得到K[-1]=50的序列
In [57]:
# 方法2
C2 = np.array([150]+list(C))
print ssig.lfilter(np.array([1.0/3]),np.array([1.0,-2.0/3]),C2)
可以看到，这两种方式的计算，除了序列头部差了一个值，其他都一样。
在本次内容中，我们从等差数列、等比数列等最基础的内容开始复习，结合现金流等内容理解了数列问题在金融方面的应用。
同时，对于一阶或者高阶的递推数列的重要概念以及应用有了更深的认识。