四变量函数最优误差求解、可视化及Excel导出问题咨询

一、修正初始误差计算代码的核心问题 #

原代码在计算call时sigma未定义会直接报错，同时数组初始化可简化为更高效的方式：

import numpy as np
import scipy.stats as si
from scipy import integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
import openpyxl

# 基础参数
r = 0.03
S0 = 1
T = 1
u1 = 1/4  # 短到期时间
strike = 1

# 初始化参数数组（替代原循环，更简洁高效）
N = np.arange(3, 31, 1)  # 积分点数：3到30
K = np.arange(0.05, 1.0, 0.05)  # 区间宽度：0.05到0.95（步长0.05）
sigma_vals = np.arange(0.21, 1.51, 0.01)  # 波动率：0.21到1.5（对应原k=20到149）

# 初始化误差数组：维度[sigma数量, N数量, K数量]
E = np.zeros((len(sigma_vals), len(N), len(K)))

# 遍历所有参数组合计算误差
for sigma_idx, sigma in enumerate(sigma_vals):
    # 计算真实期权价格call（移到sigma循环内，因依赖sigma）
    d1 = (np.log(S0 / strike) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S0 / strike) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    call = S0 * si.norm.cdf(d1) - strike * np.exp(-r * T) * si.norm.cdf(d2)
    
    for N_idx, N1 in enumerate(N):
        for K_idx, K1 in enumerate(K):
            def f(x):
                d = (np.log(x/strike) + (r + 0.5 * sigma ** 2) * (T-u1)) / (sigma * np.sqrt(T-u1))
                W1 = si.norm.pdf(d) / (x * sigma * np.sqrt(T-u1))
                d11 = (np.log(S0/x) + (r + 0.5 * sigma ** 2) * u1) / (sigma * np.sqrt(u1))
                d22 = (np.log(S0/x) + (r - 0.5 * sigma ** 2) * u1) / (sigma * np.sqrt(u1))
                call1 = S0 * si.norm.cdf(d11) - x * np.exp(-r * u1) * si.norm.cdf(d22)
                return W1 * call1
            
            hedge, _ = integrate.fixed_quad(f, strike-K1, strike+K1, n=int(N1))
            error = np.log(np.abs(hedge - call))
            E[sigma_idx, N_idx, K_idx] = error
            print(f"波动率{sigma:.2f}, 积分点数{N1}, 区间宽度{K1:.2f}的对数误差: {error:.4f}")

修正点说明： #

• 将call计算移到sigma循环内部，解决原代码sigma未定义的报错问题
• 用np.arange替代手动循环初始化参数数组，代码更简洁高效
• 用enumerate遍历数组获取索引和值，避免手动计算索引
• 解构integrate.fixed_quad返回值，忽略无意义的误差估计输出

二、误差可视化解决方案 #

E是三维数组（sigma × N × K），需分场景可视化：

1. 单sigma下的二维热力图（展示N和K对误差的影响） #

# 选择目标sigma（比如第0个，对应sigma=0.21）
target_sigma_idx = 0
sigma_val = sigma_vals[target_sigma_idx]

# 转为DataFrame
df_error = pd.DataFrame(E[target_sigma_idx], index=N, columns=K)

# 绘制热力图
plt.figure(figsize=(12,8))
heatmap = plt.imshow(df_error, cmap='viridis', aspect='auto', extent=[K.min(), K.max(), N.max(), N.min()])
plt.colorbar(heatmap, label='对数误差')
plt.xlabel('区间宽度K')
plt.ylabel('积分点数N')
plt.title(f'波动率{sigma_val:.2f}下的对数误差热力图')
plt.show()

2. 四参数下的3D可视化（展示N、K、sigma与误差的关系） #

# 生成网格数据
sigma_mesh, N_mesh, K_mesh = np.meshgrid(sigma_vals, N, K, indexing='ij')

# 绘制3D曲面图
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(sigma_mesh, N_mesh, K_mesh, E, cmap='viridis')
ax.set_xlabel('波动率sigma')
ax.set_ylabel('积分点数N')
ax.set_zlabel('区间宽度K')
ax.set_title('四参数组合下的对数误差分布')
fig.colorbar(surf, label='对数误差')
plt.show()

三、最优误差求解（找到最小误差对应的参数组合） #

# 找到最小误差的索引
min_error_idx = np.unravel_index(np.argmin(E), E.shape)
sigma_opt = sigma_vals[min_error_idx[0]]
N_opt = N[min_error_idx[1]]
K_opt = K[min_error_idx[2]]
min_error = E[min_error_idx]

print(f"最优参数组合：")
print(f"波动率sigma: {sigma_opt:.2f}")
print(f"积分点数N: {int(N_opt)}")
print(f"区间宽度K: {K_opt:.2f}")
print(f"最小对数误差: {min_error:.4f}")

四、Excel导出修正方案 #

原代码存在df未初始化、工作表名重复的问题，修正后代码如下：

with pd.ExcelWriter('error_output.xlsx') as writer:
    for sigma_idx, sigma in enumerate(sigma_vals):
        # 生成当前sigma对应的误差DataFrame
        df_current = pd.DataFrame(E[sigma_idx], index=N, columns=K)
        # 动态生成工作表名（避免重复）
        sheet_name = f"Sigma_{sigma:.2f}"
        # 导出到Excel
        df_current.to_excel(writer, sheet_name=sheet_name)
        print(f"已导出波动率{sigma:.2f}的误差矩阵到工作表{sheet_name}")

修正点说明： #

• 循环内直接创建独立的df_current，避免未初始化变量报错
• 用Sigma_xx.xx格式动态生成工作表名，避免重复和无意义名称
• 确保每个sigma对应独立的DataFrame和工作表

内容的提问来源于stack exchange，提问作者Purba Banerjee