简介一套基于Python语言实现的期权定价蒙特卡洛模拟源代码重点围绕雪球型期权设计定价算法与数值计算框架面向计算机类课程设计、期末项目开发的在校学生及自学者用于掌握随机模拟在金融衍生品定价中的实际应用。压缩包内共十三个文件以八个Python代码模块为核心涵盖期权定价框架、随机过程模拟、参数预处理、定价函数与对冲模拟等完整环节另有说明文档和备份文件辅助学习整体约十九KB。目前已有59人学习或下载。实现方案曾作为课程设计项目获得指导教师优秀评价成绩接近满分架构清晰、注释详尽除基础定价模型外特别针对路径依赖的雪球期权提供了专用算法并附有保本雪球、东兴等多个测试案例方便对照验证。适合希望深入理解蒙特卡洛方法与金融计算交叉应用的读者也可直接作为课程设计、期末项目的参考模板与实践起点。1. 雪球产品定价为什么绕不开蒙特卡洛条款里的敲入敲出让解析解失效做场外衍生品定价的同行应该都有体会普通欧式期权用 Black-Scholes 公式几秒钟就算完了但雪球产品一上来解析解基本失效。雪球的收益取决于标的价格是否跌破敲入线、是否在观察日涨过敲出线以及这个事件发生在哪一天。这种路径依赖结构没有封闭表达式只能靠蒙特卡洛模拟把未来可能走出的路径一条条画出来再计算每条路径上的现金流折现取平均。基于 Python 的期权定价蒙特卡洛模拟实现是个人开发者能快速上手验证雪球产品定价算法的路径本文用一个可运行的源码框架把雪球票息、敲入敲出判定和折现逻辑逐层拆开让你既能看明白原理也能自己改参数跑出结果。下面直接从蒙特卡洛的定价原理讲起再进入雪球条款的算法映射。2. 蒙特卡洛模拟定价原理风险中性测度、GBM路径与收敛速度2.1 风险中性测度为什么定价不能用历史涨跌概率蒙特卡洛模拟期权价格本质上是在算一个期望值生成大量标的价格路径把每条路径期末或中途的现金流折现到当前再对所有路径求平均。但这个期望里的概率不是指数明天的真实涨跌概率而是风险中性测度下的概率。在风险中性世界所有风险资产的期望收益率都等于无风险利率 r标的价格过程满足 dS r S dt σ S dW。如果你脑子里还是那个历史平均收益 8%的念头把漂移项设成 8%那算出来的不是公允价值而是带个人预期的估值。关于这个坑做过量化的人应该不陌生。雪球产品的敲入敲出概率对漂移项非常敏感年化漂移差 2%敲入概率可能就变化好几个百分点。所以定价引擎里用的漂移必须是 r - q其中 q 是标的分红收益率。为什么期权定价可以用这个看起来不真实的测度因为期权收益可以被标的资产动态复制复制成本在连续时间完备市场下唯一所以无论投资者喜不喜欢风险价格都等于风险中性期望折现。雪球虽然是带双障碍的场外合约复制并非完备但主流券商依旧先跑风险中性框架再叠加随机波动率或跳跃模型做修正。个人用 Python 验证雪球定价算法时第一步就是要确保漂移项写对。2.2 几何布朗运动与离散化从连续时间到每日路径风险中性下标的路径的离散格式是S_{tdt} S_t * exp( (r - q - 0.5 * sigma^2) * dt sigma * sqrt(dt) * epsilon )其中 epsilon 是标准正态随机数。注意那 0.5 * sigma^2 来自伊藤引理不是随便写的。很多人第一次模拟 GBM 时直接套 exp(r dt sigma sqrt(dt) epsilon)结果路径均值会随时间指数膨胀长期模拟完全失真。正确写法是让对数收益的漂移项等于 r - q - 0.5 * sigma^2。下面这段是我常用的路径生成函数用 NumPy 向量化一次生成所有路径避免逐路径循环import numpy as np def simulate_gbm_paths(S0, r, q, sigma, T, steps, n_paths, seed20240601): 生成风险中性测度下的GBM标的路径。 S0 : 期初标的价格 r : 无风险利率年化连续复利 q : 标的分红收益率年化连续复利 sigma : 年化波动率 T : 产品期限年 steps : 总时间步数通常按交易日数 n_paths: 模拟路径数 seed : 随机种子保证结果可复现 dt T / steps drift (r - q - 0.5 * sigma**2) * dt vol sigma * np.sqrt(dt) rng np.random.default_rng(seed) z rng.standard_normal((n_paths, steps)) # 累加对数收益得到每个时间步的对数价格偏移 log_returns np.cumsum(drift vol * z, axis1) # 第一列填期初价格后续列按路径展开 paths np.empty((n_paths, steps 1)) paths[:, 0] S0 paths[:, 1:] S0 * np.exp(log_returns) return paths这个函数的输出形状是 (n_paths, steps1)每一行就是一条从 S0 出发的完整价格路径。dt是每步代表的年份比如 252 个交易日对应一年dt就是 1/252。np.cumsum累加的是独立同分布的正态增量等价于对数收益率累加这是 GBM 离散化的标准做法。参数上最需要注意的是r、q、sigma必须统一为年化连续复利口径如果拿到的是百分比形式的离散利率先做一次换算再进模型。2.3 收敛速度与方差缩减路径数加一倍误差能少多少蒙特卡洛标准误等于样本标准差除以路径数的平方根。误差要减半路径数要翻到四倍这个平方根关系让很多新手把路径数堆到几百万然后抱怨速度太慢。正确的做法是先算标准误再决定路径数。比如 Snowball 定价里每条路径的现金流折现后标准差一般在 0.3~0.510 万条路径的标准误大约是 0.0014精度在小数点后第三位如果想看到第四位稳定就需要 50 万条左右。方差缩减是让 N 翻倍更有效的手段。最常用的是对偶变量法生成一半路径的随机数然后把随机数取负生成另一半两条路径自然负相关配对后均值方差大幅下降。在 GBM 模拟里实现非常容易只要把随机数矩阵拼成上下两半def simulate_antithetic_paths(S0, r, q, sigma, T, steps, n_paths, seed): assert n_paths % 2 0, n_paths 必须是偶数 dt T / steps drift (r - q - 0.5 * sigma**2) * dt vol sigma * np.sqrt(dt) rng np.random.default_rng(seed) z_half rng.standard_normal((n_paths // 2, steps)) z np.vstack([z_half, -z_half]) # 对偶路径 log_returns np.cumsum(drift vol * z, axis1) paths np.empty((n_paths, steps 1)) paths[:, 0] S0 paths[:, 1:] S0 * np.exp(log_returns) return paths但要注意对偶变量法对线性收益效果明显雪球这种带敲入敲出障碍的收益高度非线性对偶性可能被事件判定破坏。我在实践中会先用独立抽样跑一个标准误再跑对偶路径对比如果标准误没有显著下降就退回独立抽样不要为了高级而硬上。3. 雪球产品定价算法拆解从条款到现金流路径3.1 雪球条款拆解敲出线、敲入线与票息规则标准雪球产品可以这样理解投资者卖出一份带双障碍的奇异期权换取一个较高的票息。以挂钩中证500指数为例期初价格 S0 为 6000 点敲出线是期初价格的 103%敲入线是期初价格的 75%。产品期限 12 个月年化票息 15%。先看敲出在每个月末的观察日如果指数收盘价不低于敲出线产品立刻提前结束投资者拿回本金并按实际持有时间拿到年化 15% 的票息。再看敲入在存续期内通常每个交易日如果指数收盘价跌破敲入线就记一次敲入事件。敲入本身不结束产品但会影响最终结局。四种终局需要记清楚终局情形现金流单位本金折现年限观察日敲出无论之前是否敲入1 coupon * t_outt_out到期未敲出、未敲入1 coupon * TT到期未敲出、曾敲入S_T / S0T期初直接敲出极端测试1 coupon * 00注意最后一行的极端测试是后面验证引擎用的常用手段。另外票息是年化利率乘以实际持有年份不是乘以剩余期限。如果产品在第 6 个月敲出票息就是 15% * 0.5 7.5%现金流为 1.075。折现一律用连续复利 exp(-r * 持有年数)。3.2 事件判定状态机一条路径的五种终局从程序视角看每条路径就是一个状态机。状态分为初始、已敲入、已敲出、到期。判定顺序不能乱在每个敲出观察日先检查敲出一旦敲出立刻结束后续价格不再关心敲入的检查贯穿全路径但只影响到期未敲出且曾敲入的情况。所以敲入状态要用整个路径判断敲出则要看最早的观察日命中事件。这里有个容易出错的点敲入检查通常按收盘价但也有产品按日内最低价来判。若要支持日内价需要在每条路径里额外记录每个时间步的最低模拟价。本文的代码默认只用收盘价in_level直接和路径每个时间步的价格比较。如果你拿到的合同写的是盘中价低于敲入线即触发那就需要在 GBM 路径生成时用一个更细的时间步来捕捉最低价否则会低估敲入概率。3.3 从路径到现金流折现因子与票息累计现金流映射确定后蒙特卡洛定价就变成两步先生成路径再对每条路径做事件判定和现金流折现最后取平均。折现因子统一用 exp(-r * 持有年限)。敲出时持有年限是敲出时间乘以dt未敲出时持有年限就是 T。票息累计按年化线性比例不按月因为产品条款常用连续计息按实际天数除以 365。如果用每月固定票息那要另写一个票息支付表折现因子也要按实际支付日逐笔折现。下面的伪代码展示了判定顺序实际 Python 实现会在第 4 章给出for path in paths: knocked_in any(path[t] in_level for all t) knocked_out None for obs_day in out_obs_days: if path[obs_day] out_level: knocked_out obs_day break if knocked_out: payoff 1 coupon * knocked_out * dt elif not knocked_in: payoff 1 coupon * T else: payoff path[-1] / S0 pv payoff * exp(-r * held_years) accumulate pv / n_paths第二行any(...)的复杂度是 O(steps)每条路径都要全扫描。对 50 万条路径、252 步来说这个循环在纯 Python 里会很慢所以后面会改成 NumPy 的向量化比较。这也是源码解析里最值得看的部分如何把逐路径的for循环压缩成数组操作。4. 基于Python实现雪球蒙特卡洛定价源码分模块解析与参数说明4.1 环境准备与模块划分numpy向量化与随机数流运行这个源码不需要复杂依赖Python 3.8 以上装好 NumPy 就够了。如果你还在为 python 安装和 pycharm 配置 python 环境折腾直接把下面的代码存成一个.py文件命令行用python snowball_mc.py跑就行。模块划分建议把路径生成、现金流计算、主参数入口拆成三个函数方便后面做单元测试和参数扫描。不要把所有逻辑堆在主线里否则改一个条款要重跑全局。这个工程的模块划分snowball_mc.py ├── simulate_gbm_paths() # 路径生成支持独立采样/对偶采样 ├── price_snowball() # 事件判定 现金流折现 └── main() # 参数设置、结果输出和标准误报告我在量化项目里习惯先定好输入输出边界simulate_gbm_paths()只负责生成价格矩阵不关心任何雪球条款price_snowball()只接收路径矩阵和条款参数不关心路径是怎么来的。这样历史行情回放时也能复用price_snowball()直接喂真实价格路径。4.2 核心代码路径生成与事件判定先放路径生成函数比第 2 章多一个可选的antithetic开关便于后面做方差缩减对比。import numpy as np def simulate_gbm_paths(S0, r, q, sigma, T, steps, n_paths, seed20240601, antitheticFalse): GBM路径生成支持独立采样和对偶采样。 dt T / steps drift (r - q - 0.5 * sigma**2) * dt vol sigma * np.sqrt(dt) if antithetic: assert n_paths % 2 0, 对偶采样要求 n_paths 为偶数 rng np.random.default_rng(seed) z_half rng.standard_normal((n_paths // 2, steps)) z np.vstack([z_half, -z_half]) else: rng np.random.default_rng(seed) z rng.standard_normal((n_paths, steps)) log_returns np.cumsum(drift vol * z, axis1) paths np.empty((n_paths, steps 1)) paths[:, 0] S0 paths[:, 1:] S0 * np.exp(log_returns) return paths这个函数用np.random.default_rng(seed)创建独立的随机数生成器而不是全局的np.random.seed。两种写法的区别在于default_rng每次调用会维持自己的状态不会污染外部随机数流全局seed则会影响整个程序后续所有随机操作调试时很容易踩到奇怪的不复现问题。函数内部的dt、drift、vol参数我已经在第 2 章解释过这里要注意antithetic模式下路径数必须为偶数否则vstack会报错。4.3 核心代码雪球现金流计算与折现逻辑接下来是雪球定价函数把第 3 章的事件状态机用向量化数组操作实现避免逐路径循环。def price_snowball(paths, S0, r, T, out_level, in_level, out_obs_indices, coupon, dt): 输入 paths : shape (n_paths, steps1) 的价格路径矩阵 S0 : 期初价格 r : 无风险利率年化连续复利 T : 产品期限年 out_level : 敲出价格阈值 in_level : 敲入价格阈值 out_obs_indices : 敲出观察日索引数组例如 [21, 42, ..., 252] coupon : 年化票息率 dt : 每个时间步代表的年数 返回 (公允价格, 蒙特卡洛标准误) n_paths paths.shape[0] # 敲入路径中任意一个时间步的收盘价低于敲入线 knocked_in np.any(paths[:, 1:] in_level, axis1) # 敲出逐个观察日检查只记录最早命中 knocked_out np.zeros(n_paths, dtypebool) knocked_out_year np.full(n_paths, T) # 默认到期 for obs_idx in out_obs_indices: obs_idx int(obs_idx) if obs_idx paths.shape[1]: continue hit paths[:, obs_idx] out_level new_hit hit (~knocked_out) knocked_out[new_hit] True knocked_out_year[new_hit] obs_idx * dt # 现金流映射敲出或未敲入未敲出 本金票息敲入未敲出 期末价格比例 has_coupon knocked_out | (~knocked_in) payoff np.where( has_coupon, 1.0 coupon * knocked_out_year, paths[:, -1] / S0 ) discount np.exp(-r * knocked_out_year) present_values payoff * discount price present_values.mean() std_err present_values.std(ddof1) / np.sqrt(n_paths) return price, std_err逐段解释。敲入判断用np.any(paths[:, 1:] in_level, axis1)注意paths[:, 0]是期初价格不参与敲入判断因为期初价格不可能低于敲入线除非期初就爆产品不会成立。敲出判断循环观察日new_hit hit (~knocked_out)保证只记录第一次敲出因为已经敲出的路径后面再命中也不需要更新了。knocked_out_year默认值是 T表示未敲出路径的持有年限是到期日敲出路径会在循环里被覆盖成实际敲出时间乘以dt。现金流映射里has_coupon覆盖了两种有票息的情形敲出、和到期未敲入未敲出。np.where根据是否拿票息选择两种 payoff 公式。敲入未敲出的路径payoff 是期末价格除以期初价格这表示投资者承担标的跌幅也可能有涨幅但雪球一般不会在这种情形下让投资者赚钱不对如果敲入后指数回升但未超过敲出线期末价格可能略低于期初也可能高于期初按合同一般结算标的涨跌幅所以S_T/S0是通用写法。这段代码里np.where和np.exp都是向量化操作10 万条路径瞬间完成。标准误std(ddof1)用样本标准差这是无偏估计分母不要错。4.4 主参数入口路径数、时间步和观察日怎么设最后是主函数把所有参数串起来并验证收敛性。def main(): # ---- 标的参数 ---- S0 1.0 # 归一化期初价格 r 0.03 # 无风险利率 q 0.02 # 股息率 sigma 0.25 # 年化波动率 T 1.0 # 一年期 steps 252 # 按交易日离散 # ---- 雪球条款 ---- out_level 1.03 # 敲出价为期初的103% in_level 0.75 # 敲入价为期初的75% coupon 0.15 # 年化票息率15% # ---- 观察日 ---- # 简化每月末敲出观察按每月21个交易日 out_obs_indices np.arange(21, steps 1, 21) dt T / steps # ---- 模拟 ---- n_paths 200_000 seed 20240601 paths simulate_gbm_paths(S0, r, q, sigma, T, steps, n_paths, seed) price, stderr price_snowball( paths, S0, r, T, out_level, in_level, out_obs_indices, coupon, dt ) print(f雪球公允价格: {price:.6f}) print(f蒙特卡洛标准误: {stderr:.6f}) # 额外的收敛性检查分段计算标准误是否随路径数下降 half_price, half_stderr price_snowball( paths[: n_paths // 2], S0, r, T, out_level, in_level, out_obs_indices, coupon, dt ) print(f前一半路径价格: {half_price:.6f}, 标准误: {half_stderr:.6f}) if __name__ __main__: main()路径数n_paths我一般先设 20 万跑完看标准误。如果标准误大于 0.002再翻倍如果已经到 0.001 以下就不必再加路径了。out_obs_indices np.arange(21, steps 1, 21)是按每月 21 个交易日简化出来的实盘中必须用真实交易日历把每个月最后一个交易日对应的时间步索引放到这个数组里。用归一化S01.0的好处是所有条款阈值和 payoff 都在 1 附近避免大数运算的浮点误差结果也更容易对比。5. 雪球定价源码的五处常见踩坑随机数、时间步与方差缩减5.1 现象每次运行结果都不一样回测无法复现原因全局随机数没有固定种子或者不同 Python 进程的随机状态不同。很多人用np.random.seed()但如果你的代码里某个地方调用了第三方库偷偷消耗了全局随机数后续结果就对不上了。解决统一用np.random.default_rng(seed)并且把seed作为函数参数传入同一个 seed 在任何机器上产生完全相同的随机数流。这给后面的 Greeks 有限差分留了关键条件——微扰 S0 时必须使用完全相同的随机数序列否则差分结果会被蒙特卡洛噪声淹没。我一般固定一个种子列表例如[20240101, 20240102, ...]每个参数组合跑 10 组取均值和标准差既保证复现又削弱种子偏见。5.2 现象敲出频率比真实合约低价格偏高原因敲出观察日和时间步不匹配。如果你在代码里把每个交易日都当敲出观察日那敲出概率会显著高于真实合约真实通常每月观察一次定价就偏低反过来如果只在期末观察一次漏掉中间可能发生的敲出定价就偏高。解决把out_obs_indices定义为一个显式的列表严格对应合约里的敲出观察日历。合约上写每月第 15 个交易日就按时序把对应步索引算出来别偷懒用range(21, 252, 21)这种均匀近似。另一个相关坑是敲入观察频率写错合约写每日观察你的paths[:, 1:] in_level已经覆盖了每日但如果合约写的是收盘价而你的模拟步长是每周也会漏掉。5.3 现象路径的长期均值膨胀票息被高估原因GBM 离散化写成了S * (1 (r-q)*dt sigma*sqrt(dt)*z)或者对数收益里漏了-0.5*sigma^2。算术布朗运动的形式在dt小的时候看着没问题但线性化误差会在 252 步累积导致路径均值随模拟时间指数偏离理论值。解决使用第 4 章的log_returns cumsum(drift vol * z)写法其中drift (r - q - 0.5*sigma**2) * dt。验证方法模拟完计算所有期末价格的均值理论值应接近S0 * exp((r-q) * T)偏差如果在 0.5% 以上检查漂移项。这也是蒙特卡洛定价引擎最隐蔽的一个数学坑。5.4 现象雪球价格带着明显的市场观点不同人算出的价格天差地别原因把漂移项设成了历史平均收益率比如用指数过去五年的年化收益 10% 作为mu。这在风险中性测度下是错误的因为定价不能包含主观预测正确做法是mu r - q。如果你一定要用真实测度做场景分析比如压力测试那算出来的叫 真实概率下的预期收益不能直接与市场无套利价格混为一谈。很多刚从 python 量化交易策略代码转到期权定价的新手会在这里翻车建议代码里变量名直接写risk_free_rate不要写mu从命名上提醒自己。5.5 现象用了方差缩减后结果和独立抽样差出 3 倍标准误原因对偶变量法在非线性障碍产品上失效。雪球的现金流在敲出点附近是不连续的对偶路径之间的负相关可能被事件判定打破甚至变成正相关。解决不要盲目相信方差缩减保留独立抽样结果作为基准。我在实践中常用控制变量——以标的期末价格为控制变量因为交易员关心 Delta而雪球价格与期末价格高度相关。但控制变量的系数需要用最小二乘估计估计时要注意不能使用与定价相同的路径数据否则会出现样本内偏差。最稳妥的办法先用 10 万条独立路径跑一个基准价格和标准误再 5 万对偶路径跑一次两者在标准误 2 倍以内才认为缩减有效。6. 从源码到验证用历史行情回放与解析解校准蒙特卡洛引擎6.1 用历史行情做路径回放验证事件判定逻辑写好的price_snowball()不仅吃蒙特卡洛路径也能直接吃真实历史行情。我会从行情软件取中证500过去三年的日收盘价构造一条长度为steps1的真实路径扔进price_snowball()然后手工核对输出找到历史上实际发生的敲入敲出事件验证函数判定的敲出观察日、票息比例和折现系数是不是和按合同条款手工算的一致。这个方法虽然不能验证概率分布但能验证事件判定状态机和现金流映射里有没有低级错误。具体操作把真实价格序列归一化到 S01然后调用price_snowball(paths_real, S0, r, T, out_level, in_level, out_obs_indices, coupon, dt)。注意这里steps要和真实序列长度一致dt T / steps敲出观察日索引必须按真实交易日历重新对一遍。如果输出价格明显不合理多半是观察日索引和时间轴错位。6.2 用Delta和Greeks校验定价引擎的敏感性蒙特卡洛价格只能告诉你公允价格但交易台报价还要看风险敞口。Delta 可以用有限差分法近似但必须保证两个价格来自同一套随机数。做法是把simulate_gbm_paths()的S0参数改为S0_up和S0_down其他参数包括 seed 全部相同重新生成路径后分别定价def compute_delta(S0, r, q, sigma, T, steps, n_paths, seed, out_level, in_level, out_obs, coupon): dt T / steps eps 0.01 * S0 paths_up simulate_gbm_paths(S0 eps, r, q, sigma, T, steps, n_paths, seed) paths_dn simulate_gbm_paths(S0 - eps, r, q, sigma, T, steps, n_paths, seed) price_up, _ price_snowball(paths_up, S0 eps, r, T, out_level, in_level, out_obs, coupon, dt) price_dn, _ price_snowball(paths_dn, S0 - eps, r, T, out_level, in_level, out_obs, coupon, dt) delta (price_up - price_dn) / (2 * eps) return delta有限差分法在敲出点附近会有偏导数不连续的问题但宏观上可以作为 sanity check雪球 Delta 通常在 0.3~0.5 之间敲出概率高时 Delta 会变小如果算出来是负数或大于 1基本可以确定事件判定有 bug。更准确的路径wise导数需要对敲出事件做平滑处理这里不展开。6.3 一个自查习惯用Black-Scholes解析解给蒙特卡洛做基准测试我在改完任何雪球代码后最先跑的不是真实参数而是两个极端测试。第一个把敲入线设成 0永不敲入敲出线设成一个非常大的数永不敲出那么所有路径都落到未敲入未敲出现金流恒等于1 coupon * T折现后理论值是(1 coupon * T) * exp(-r * T)。第二个把敲出线设成 0那么期初立即可判定敲出持有年限为 0所有路径现金流为 1价格理论值为 1。这两个测试能瞬间暴露折现因子、票息公式或观察日索引的常规错误。更完整的校准是用普通欧式看涨期权删掉敲入敲出逻辑把 payoff 改成max(S_T - K, 0)用同一个 GBM 引擎跑蒙特卡洛再和 Black-Scholes 公式对比。误差应该在 2 倍标准误以内。这个测试能验证路径生成和随机数质量是否正确因为欧式期权的价值只取决于期末价格分布不依赖路径中间事件。我在做雪球定价的时候这几乎成了肌肉记忆先校准引擎再谈条款。每次把新条款丢给price_snowball()之前我都会先跑一遍这两个极端测试确认现金流计算器没被上次改参数弄坏。这是我从无数次翻车里养出来的习惯希望你也能少踩几个坑。希望帮到你。本文还有配套的精品资源点击获取