美式期权baw定价的python实现

这一年多，写小程序的语言，逐步从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)