这一年多，写小程序的语言，逐步从matalb转到了python。实际上，我觉得MATLAB对于我来说，要友好很多，帮助文件论述清晰，还有详尽的示例，而且相关工具函数齐备，临时做些计算非常方便。不过，知乎上提及金融分析，基本上都是python了，而且社区论坛python也要活跃很多。虽然这些论坛主要是搞程序的在做，但对新进入的人来说，肯定优选python作为入门工具了。Matlab在金融工程里面的份额肯定要减少很多了。

# -*- coding: utf-8 -*-

"""

__author__ = 'laofish'

__title__=pyderivatives_fun.py

__mtime__ = 2017-01-07

"""

import numpy as np

import time

import pandas as pd

from optbkfun import *

from scipy.stats import norm

# from statsmodels import *

import scipy.optimize as opt

from math import * # cmath支持负数

class opttypeerr(Exception):

"""自定义期权类型错误"""

# 上面这种语法实际上是一个类的继承语法

def __init__(self, value):

Exception.__init__(self)

self.value = "期权类型错误"

def bsprice(S,K,sigma,T,r):

"""普通期权的bs定价公式"""

d1=(log(S/K)+(r+sigma**2/2)*T)/(sigma*sqrt(T))

d2=d1-sigma*sqrt(T)

c=S*norm.cdf(d1)-K*exp(-r*T)*norm.cdf(d2)

p=K*exp(-r*T)*norm.cdf(-d2)-S*norm.cdf(-d1)

return [c,p]

def bsmprice(S,K,sigma,T,r,q):

"""普通期权的bsm定价公式，带股息率"""

d1=(log(S/K)+(r-q+sigma**2/2)*T)/(sigma*sqrt(T))

d2=d1-sigma*sqrt(T)

c=S*exp(-q*T)*norm.cdf(d1)-K*exp(-r*T)*norm.cdf(d2)

p=K*exp(-r*T)*norm.cdf(-d2)-S*exp(-q*T)*norm.cdf(-d1)

return [c,p]

def futOptPrice(F,K,sigma,T,r):

"""black76 期货期权定价公式"""

# sigma为期货波动率

# 期货和期权的到期日不需要一致

d1=(log(F/K)+sigma**2*T/2)/(sigma*sqrt(T))

d2=d1-sigma*sqrt(T)

c=exp(-r*T)*(F*norm.cdf(d1)-K*norm.cdf(d2))

p=exp(-r*T)*(K*norm.cdf(-d2)-F*norm.cdf(-d1))

return [c,p]

def imp_br(c,p,S,K,r,T):

"""计算隐含股息率"""

q=-1/T*log((c-p+K*exp(-r*T))/S)

return q

def findSx(Sx,K,r,q,v,T,PutCall):

"""找到美式期权提前行权的Sx"""

n = 2 * (r - q) / v ** 2

k = 2 * r / v ** 2 / (1 - exp(-r * T))

sigma=v

if Sx<0:

y=1e1000

elif PutCall=='C':

q2 = (1-n+sqrt((n-1)**2+4*k))/2

y = (bsmprice(Sx,K,sigma,T,r,q)[0]

+ (1 - exp(-q * T) * norm.cdf((log(Sx / K) + (r - q + v **2 / 2)) / v / sqrt(T))) * Sx / q2

- Sx + K) ** 2

elif PutCall=='P':

q1 = (1 - n - sqrt((n - 1) ** 2 + 4 * k)) / 2

y = (bsmprice(Sx,K,sigma,T,r,q)[1]

- (1 - exp(-q * T) * norm.cdf(-(log(Sx / K) + (r - q + v ** 2 / 2)) / v / sqrt(T))) * Sx / q1

+ Sx - K) ** 2

else:

raise opttypeerr

return y

# def findSx_baw(x)

def bawPrice(S,K,sigma,T,r,q):

"""美式期权baw定价近似解法"""

v=sigma

[c,p]=bsmprice(S,K,sigma,T,r,q)

start = S

func=lambda s:findSx(s, K, r, q, v, T, 'C')

data = opt.fmin(func, start)

Sx=data[0]

# print(data)

d1 = (log(Sx / K) + (r - q + v ** 2 / 2)) / v / sqrt(T)

n = 2 * (r - q) / v ** 2

k = 2 * r / v ** 2 / (1 - exp(-r * T))

q2 = (1 - n + sqrt((n - 1) ** 2 + 4 * k)) / 2

A2 = Sx * (1 - exp(-q * T) * norm.cdf(d1)) / q2

if S < Sx:

ac = c + A2 * (S / Sx) ** q2

else:

ac = S - K

print(Sx)

func = lambda s: findSx(s, K, r, q, v, T, 'P')

data= opt.fmin(func, start)

Sx=data[0]

d1 = (log(Sx/K) + (r - q + v ** 2 / 2)) / v / sqrt(T)

n = 2 * (r - q) / v ** 2

k = 2 * r / v ** 2 / (1 - exp(-r * T))

q1 = (1 - n - sqrt((n - 1) ** 2 + 4 * k)) / 2

A1 = -Sx * (1 - exp(-q * T) * norm.cdf(-d1)) / q1

if S > Sx:

ap = p + A1 * (S / Sx) ** q1

else:

ap = K - S

print(Sx)

return [ac,ap]

if __name__=='__main__':

try:

bawPrice(100, 150, 0.2, 0.3, 0.04, 0.01)

print(bsprice(100,100,0.2,1,0.05))

print(futOptPrice(20, 18, 0.25, 4/12, 0.09))

except opttypeerr as e:

print(e.value)