(Python代码）编程实现局部波动率模型

#导入库
import numpy as np
import math

#函数定义
def tree_simulate(S0,T,r,Vol,M): ##SO为标的价格，T是总时间，r为无风险利率，Vol为初始波动率，M是时间间隔个数
   #参数初始处理与生成
    dt = T / M  ##总的时间/时间间隔个数 = 时间间隔个数
    S = np.zeros((M+1, M+1)) ##定义期权价格矩阵
    S[0, 0] = S0 ##第一个元素为初始价格
    Sigma = np.zeros((M+1, M+1)) ##定义波动率矩阵
    Sigma[0, 0] = Vol ##第一个元素为初始波动率
    P= np.zeros((M+1, M+1)) ##定义概率矩阵
    CuP=np.zeros((M+1, M+1))##定义累计概率矩阵
    CuP[0,0]=1 ##第一个元素定义为100%，即1
    #完成局部波动率二叉树修改
    col=np.linspace(2,2*math.floor(M/2),math.floor(M/2))#依据需要重新设定的中心节点的列序数规律，生成对应列的等差数组
    n=int(len(col))
    row=col/2 #依据需要重设的行列序数之间的对应关系，生成对应行的序数数组
    for i in range(n):
        S[int(row[i]),int(col[i])]=S0
        Sigma[int(row[i]),int(col[i])]=max(Vol-((S[int(row[i])*math.e**(r*dt),int(col[i])]-S0)/S0),0.01)
    #局部波动率二叉树生成
    for t in range(1, M+1):
        if t==1:#解决第一步正常二叉树
            S[0,t]=S[0,t-1]*math.e**(Sigma[0,t-1]*math.sqrt(dt))
            Sigma[0,t]=max(Vol-(S[0,t]-S0)/S0,0.01)
            S[1,t]=S[0,t-1]*math.e**(-Sigma[0,t-1]*math.sqrt(dt))
            Sigma[1,t]=max(Vol-(S[1,t]-S0)/S0,0.01)
            P[0,0]=(S[0,0]*math.e**(r*dt)-S[1,t])/(S[0,t]-S[1,t])
            
        if (col==t).any()==True:#奇数个节点二叉树重构
            #从重构的中心节点开始划分为上行节点和下行节点，通过不同方式计算
            
            for i in range(int(row[int(t/2-1)]-1),-1,-1):#上行节点，从最接近中心节点的位置向上计算
                S[i,t]=S[i,t-1]*math.e**(r*dt)+(math.pow(S[i,t-1]*Sigma[i,t-1],2)*dt)/(S[i,t-1]*math.e**(r*dt)-S[i+1,t])
                Sigma[i,t]=max(Vol-(S[i,t]-S0)/S0,0.01)
                P[i,t-1]=(S[i,t-1]*math.e**(r*dt)-S[i+1,t])/(S[i,t]-S[i+1,t])
            
            for j in range(int(row[int(t/2-1)]+1),int(t+1)):#下行节点，从最接近中心节点的位置向下计算
                S[j,t]=S[j-1,t-1]*math.e**(r*dt)-(math.pow(S[j-1,t-1]*Sigma[j-1,t-1],2)*dt)/(S[j-1,t]-S[j-1,t-1]*math.e**(r*dt))
                Sigma[j,t]=max(Vol-(S[j,t]-S0)/S0,0.01)
                P[j-1,t-1]=(S[j-1,t-1]*math.e**(r*dt)-S[j,t])/(S[j-1,t]-S[j,t])
            
            #补足中心节点对应的上行概率
            P[int(row[int(t/2-1)]),t-1]=(S[int(row[int(t/2-1)]),t-1]*math.e**(r*dt)-S[int(row[int(t/2-1)]+1),t])/(S[int(row[int(t/2-1)]),t]-S[int(row[int(t/2-1)])+1,t])
        
        if (col==t).any()==False and t>1:#偶数个节点二叉树调整
            
            #先调整由重构后的中心节点出发到达的两个节点具体数值
            S[int(row[int(t/2-1)]),t]=S[int(row[int(t/2-1)]),t-1]*math.e**(Sigma[int(row[int(t/2-1)]),t-1]*math.sqrt(dt))
            Sigma[int(row[int(t/2-1)]),t]=max(Vol-(S[int(row[int(t/2-1)]),t]-S0)/S0,0.01)
            S[int(row[int(t/2-1)]+1),t]=S[int(row[int(t/2-1)]),t-1]*math.e**(-Sigma[int(row[int(t/2-1)]),t-1]*math.sqrt(dt))
            Sigma[int(row[int(t/2-1)]+1),t]=max(Vol-(S[int(row[int(t/2-1)]+1),t]-S0)/S0,0.01)
            
            #上行节点计算
            for i in range(int(row[int(t/2-1)]-1),-1,-1):
                S[i,t]=S[i,t-1]*math.e**(r*dt)+(math.pow(S[i,t-1]*Sigma[i,t-1],2)*dt)/(S[i,t-1]*math.e**(r*dt)-S[i+1,t])
                Sigma[i,t]=max(Vol-(S[i,t]-S0)/S0,0.01)
                P[i,t-1]=(S[i,t-1]*math.e**(r*dt)-S[i+1,t])/(S[i,t]-S[i+1,t])
            
            #下行节点计算
            for j in range(int(row[int(t/2-1)]+1),int(t+1)):
                S[j,t]=S[j-1,t-1]*math.e**(r*dt)-(math.pow(S[j-1,t-1]*Sigma[j-1,t-1],2)*dt)/(S[j-1,t]-S[j-1,t-1]*math.e**(r*dt))
                Sigma[j,t]=max(Vol-(S[j,t]-S0)/S0,0.01)
                P[j-1,t-1]=(S[j-1,t-1]*math.e**(r*dt)-S[j,t])/(S[j-1,t]-S[j,t])
        
        #累积概率树图生成：累计概率之间存在明显规律直接计算
        CuP[0,t]=CuP[0,t-1]*P[0,t-1]
        CuP[t,t]=CuP[t-1,t-1]*(1-P[t-1,t-1])
        for i in range(1,t):
            CuP[i,t]=CuP[i-1,t-1]*(1-P[i-1,t-1])+CuP[i,t-1]*P[i,t-1]   
    
    #期权定价
    Price=sum((np.maximum(S-102,0)*CuP)[0:M+1,M])
    return S,Sigma,P,CuP,Price

#输出显示4位小数
np.set_printoptions(formatter={'float':
                               lambda x: '%8.4f' % x})