简介PyVaR 是一套面向金融工程学习者与量化分析人员的风险价值计算代码包围绕蒙特卡罗模拟与几何布朗运动GBM模型估算股票 VaR并借助共享内存并行机制提升多核环境下的模拟效率。压缩包共 182 个文件、约 14.62MB以 98 个 m 文件、29 个 txt 说明、11 个 f90 源码及若干 mex 二进制、mat 数据、py 脚本和 readme 文档为主兼顾 MATLAB 调用、Fortran 数值计算与 Python 示例目录层次便于按模块查阅。已有 858 人学习适合需要理解 VaR 建模流程、并行加速思路与跨语言实现细节的读者。包内可看到蒙特卡罗路径生成、GBM 参数设置、示例脚本、使用说明与测试验证材料能帮助读者把置信水平、持有期和价格分布等要素落到代码中并观察并行任务拆分对计算速度的影响为风险管理与投资决策场景提供可复用的实验基础。1. 收盘前十分钟风控要的那个数下午两点五十风控同事发来一句话手上这只票95% 置信度下明天最多亏多少。要的不是“大概跌三个点”而是一个带分布假设、能写进日报的数字——这就是 VaR风险价值。pyVaR 做的是这件事用几何布朗运动GBM描述股价用蒙特卡罗抽出几十万到上百万条未来路径再从损益分布的尾部取分位数顺带给出期望尾部损失 ES。整个计算走共享内存并行进程池分块跑路径结果落在同一块共享内存里省掉子进程把大数组 pickle 回主进程的来回拷贝。它适合手里有一篮子标的、每天都要批量出一版风险测算又不想为这点计算量去碰集群的人。模型口径和原版 MATLAB 脚本一致我拿它对拍过数值中间有几处对不上的坑后面细讲。2. 工程结构与运行入口从一份 toml 配置到命令行跑通2.1 技术选型numpy 向量化打底multiprocessing 做并行这类计算天生适合向量化生成标准正态随机数、乘波动率、累加、取终值四步全是 C 层循环numpy 不用写一行显式 for。真正卡住性能的从来不是算法而是“一天几十次调用、每次几百万条路径”堆出来的墙钟时间。选型上有三个明显的分岔口。第一不要用 pandas 逐行模拟DataFrame 的索引对齐和 NaN 语义在这个场景里纯属负担价格矩阵一进 numpy 就该丢掉列名。第二不要指望多线程np.random.Generator.standard_normal和np.cumsum这两个热点都不释放 GIL线程开 16 个也只是轮流排队还额外多出上下文切换。第三numba 或 Cython 确实能把单线程压到极限但环境里多一套编译链部署到别人的机器上就容易翻车收益也不如直接多进程来得实在。所以主干就是numpy出数值、multiprocessing.Pool出并行、shared_memory出结果回收。依赖上只需要numpy、pandas读行情、pyarrowparquet、scipy可选回测里的卡方分位数没有重框架。2.2 目录结构与模块职责包名用全小写pyvar避免和项目名大小写混用导致import踩坑。目录按“数据 → 模型 → 统计 → 调度”切每个文件只干一件事pyvar/ ├── pyvar/ │ ├── __init__.py │ ├── data.py # 行情读取、清洗、GBM 参数估计 │ ├── gbm.py # GBM 路径/终值生成单进程纯函数 │ ├── var.py # 分位数、VaR、ES、回测统计 │ ├── parallel.py # 共享内存 进程池调度 │ ├── config.py # toml 配置解析与参数校验 │ └── cli.py # 命令行入口 ├── configs/ │ └── default.toml ├── tests/ │ ├── test_gbm.py # 均值/方差收敛性测试 │ └── test_var.py # 与解析解、分位数方法对照 ├── requirements.txt └── README.md文件职责关键点data.py宽表读入、停牌填充、对数收益、年化 mu/sigma停牌日必须 ffill否则收益序列缺一段gbm.py纯函数无全局状态输入参数输出价格单进程可测方便和并行版对拍var.pynp.percentile取左尾ES 取尾部均值插值方法要和 MATLAB 对齐parallel.py切块、建共享内存、发任务、收结果这里是唯一有并发语义的文件config.py类型转换 范围校验confidence 必须在 (0,1)n_paths 要能整除 chunk这种切法的好处是可测性gbm.py和var.py都是纯函数单元测试和串行对照都方便parallel.py只负责“怎么切”出问题九成是并发或内存问题范围收得很窄。2.3 环境安装与最小可运行命令Python 3.10 以上主要是shared_memory和np.random.Generator的行为在旧版本上有差异。虚拟环境建好之后装依赖、装本包python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate pip install numpy pandas pyarrow scipy pip install -e . # 以可编辑模式装进当前环境配置文件用 toml比 argparse 堆一屏参数清爽也方便版本管理# configs/default.toml [market] prices data/prices.parquet # 宽表index日期columns代码 lookback_days 500 # 估计 mu/sigma 的回看窗口 confidence 0.95 # VaR 置信度 horizon_days 1 # 持有期 [model] steps_per_day 1 # 每天拆几步欧式一步即可 n_paths 2000000 # 总路径数 [runtime] workers 8 # 进程数 物理核数 chunk_paths 250000 # 每块路径数 seed 20240101 # 固定种子保证可复现跑起来只要一行python -m pyvar.cli --config configs/default.toml --out reports/var.csv输出是一张按标的排列的表每行包含ticker、mu、sigma、var、es、n_paths、elapsed_s。跑通之后先别急着调大n_paths先拿一两个标的和串行版本对一下数值确认分位数口径没问题再去压性能。2.4 输入数据长什么样行情必须是宽表索引是交易日列是标的代码值是收盘价复权后。用长表也行但每次 pivot 都是一笔额外开销批量测算时没必要。读入后的第一件事是ffill停牌日不填收益率序列就会少一段年化波动率直接偏小VaR 自然偏乐观这是最隐蔽的一类错误。# pyvar/data.py import numpy as np import pandas as pd def load_prices(path: str) - pd.DataFrame: 宽表读入index日期, columns代码, values复权收盘价 df pd.read_parquet(path) df df.sort_index().ffill() # 停牌日沿用前值 df df.dropna(howany) # 上市晚的标的直接剔除别补零 if (df 0).any().any(): raise ValueError(价格序列存在非正值检查复权方式) return df def estimate_gbm_params(prices: pd.DataFrame, lookback: int 500, ann: int 252): 返回 {ticker: (mu, sigma)}均为年化值 px prices.tail(lookback) log_ret np.log(px / px.shift(1)).dropna() sigma log_ret.std(ddof1) * np.sqrt(ann) # 年化波动率 mu log_ret.mean() * ann # 年化漂移 return {c: (float(mu[c]), float(sigma[c])) for c in prices.columns}lookback决定参数稳定性太短sigma 跟着最近一波行情剧烈跳动太长又跟不上市场结构变化。500 个交易日是常见折中约两年。ddof1用样本标准差和 MATLAB 的std默认行为一致这点对后续对拍很重要。mu的估计噪声远大于sigma实际做短期风险测算时把mu直接设成无风险利率甚至 0、只让sigma从历史数据估是更稳的做法——漂移项对一天的 VaR 影响几乎可以忽略但估计误差会实打实传进结果里。3. 核心算法落地GBM 从参数到损益分布3.1 两套写法逐步模拟和一步到位别用错场景GBM 的解析解是S_T S_0 * exp((mu - sigma^2/2)T sigma*sqrt(T)*Z)Z 是标准正态。很多入门实现会老老实实按天迭代生成一个(n_paths, steps)的完整价格矩阵。这样做不是错但内存占用是n_paths * steps * 8字节100 万条路径、252 个交易日就是 2GB机器立马吃紧。对纯欧式的 VaR 来说只需要终值分布中间路径一点用没有。一步到位只需要(n_paths, 1)内存直接降两个数量级。所以工程里我一般写两条路# pyvar/gbm.py import numpy as np def gbm_terminal_fast(S0, mu, sigma, T, n_paths, rng): 只求终值欧式 VaR 用这个内存 O(n_paths) z rng.standard_normal(n_paths) return S0 * np.exp((mu - 0.5 * sigma ** 2) * T sigma * np.sqrt(T) * z) def gbm_paths(S0, mu, sigma, T, steps, n_paths, rng): 完整路径矩阵给障碍/路径依赖产品用内存 O(n_paths * steps) dt T / steps z rng.standard_normal((n_paths, steps)) incr (mu - 0.5 * sigma ** 2) * dt sigma * np.sqrt(dt) * z log_paths np.cumsum(incr, axis1) # 沿时间维累加 return S0 * np.exp(log_paths)两个函数在数学上等价区别只在内存和用途。incr里(mu - 0.5*sigma**2)*dt是漂移修正项少了这个减号模拟出来的价格期望会系统性偏高一年期上能差出几个百分点属于那种“结果看着合理、但就是不对”的典型错误。np.cumsum(..., axis1)的axis必须是 1写成 0 就变成沿路径方向累加出来是一片乱码数字。T用年为单位horizon_days1时T 1/252steps只有需要中间路径时才大于 1欧式场景固定 1。3.2 参数估计mu 和 sigma 从哪来决定了结果可信度估值口径必须统一。历史数据是日频收盘价年化因子在 A 股常用 252商品或加密资产会是 365这个数字一改sigma 差 20%VaR 跟着差 20%。团队里最好在配置里写死别让每个人自己拍。mu的处理更微妙。用历史对数收益均值年化500 天样本估出来的值动辄 ±30%噪声比信号大得多。我一般的做法是短期风险测算里把mu设成 0 或资金成本只让sigma从数据来只有在做长期压力测试、需要模拟价格中枢时才用历史漂移。这一条在 2.4 的代码里是显式暴露的改配置就能切不用动模型层。还有一点容易忽略GBM 假设对数收益服从正态真实市场是尖峰厚尾95% 分位数上还好99% 上 GBM 会系统性低估尾部风险。所以项目里保留了分布可替换的位置第 6 章会讲怎么换成 t 分布。3.3 向量化路径生成为什么不用循环rng.standard_normal((n_paths, steps))一次调用就拿到整个随机数矩阵内部是 C 的梅森旋转或 PCG64比 Python 层循环快三到四个数量级。代价是内存峰值高所以并行版本里要配合分块一块一块生成见第 4 章。incr.sum(axis1)和log_paths[:, -1]在数值上等价前者少一次完整矩阵拷贝只求终值时用sum。如果steps1两者一样。一个细节S0 * np.exp(...)里的乘法在数组外面numpy 会广播不用写成np.exp(...) * S0没有性能差异但前者语义更清楚——先算收益率再折回价格。3.4 从损益分布到 VaR 与 ES拿到终值数组后损益pnl ST - S0VaR 就是左尾分位数取正# pyvar/var.py import numpy as np def var_es(pnl: np.ndarray, confidence: float 0.95): pnl: 损益数组正赚负亏 返回 (VaR, ES)都取正数表示损失幅度 q np.percentile(pnl, (1 - confidence) * 100, methodlinear) var -q tail pnl[pnl q] es -tail.mean() if tail.size else var return float(var), float(es)methodlinear是 numpy 1.22 之后的写法之前叫interpolation用旧版本会直接抛TypeError。这个参数是 VaR 对拍失败的头号嫌疑MATLAB 的prctile默认插值方式、pandas 的quantile默认linear、Excel 的PERCENTILE又是另一种同样一组数能差出千分之几路径数少的时候这个差异肉眼可见。tail为空表示分位点正好落在最小值上属于极小概率兜底返回var避免nan。ES 是尾部均值比 VaR 更能反映极端损失两个一起输出别只报 VaR。注意pnl的符号约定必须全项目统一。有的实现用收益率、有的用绝对金额混着传进去会让 VaR 符号翻转而且这种错在正收益占多数时看不出来只有亏钱那天才暴露。4. 共享内存并行切分、回传与调优4.1 为什么必须绕开 GIL以及怎么切numpy的部分 ufunc 会释放 GIL但随机数生成和cumsum这类归约操作不在其列。实测下来8 线程跑这段代码的加速比不到 1.3还不如单线程。所以并行单元必须是进程。切分策略就一个原则每块要够大摊平进程启动和内存映射的固定成本。我一般让chunk_paths 100000workers取物理核数超线程对纯浮点计算帮助有限反而增加调度开销。块数不用等于进程数Pool.map会自动排队块数多一点比如 2 到 4 倍 workers反而能减少长尾防止某个进程分到的块特别慢把整体拖住。但块数也不是越多越好。每多一块就多一次共享内存 attach、一次随机数生成器初始化单块 5 万条以下时这部分开销能占到总时间的 40% 以上。4.2 共享内存写回与可复现的种子管理朴素的写法是让 worker 返回自己的terminal数组主进程np.concatenate。问题是每次返回都要 pickle 一遍100 万条 float64 是 8MB8 个进程来回一趟就是几十兆的序列化开销延迟还高。共享内存的思路是主进程先划一块地worker 直接把结果写进去返回值只是一个状态码。# pyvar/parallel.py import numpy as np from multiprocessing import Pool, shared_memory def _worker(job): name, n_total, start, stop, S0, mu, sigma, T, steps, seed job n stop - start rng np.random.default_rng(seed) # 每个 worker 独立种子 dt T / steps z rng.standard_normal((n, steps)) incr (mu - 0.5 * sigma ** 2) * dt sigma * np.sqrt(dt) * z terminal S0 * np.exp(incr.sum(axis1)) # 只留终值省内存 shm shared_memory.SharedMemory(namename) # attach 已有块 out np.ndarray((n_total,), dtypenp.float64, buffershm.buf) out[start:stop] terminal # 直接写进共享区 shm.close() # worker 只 close不 unlink return 0 def simulate_terminal_parallel(S0, mu, sigma, T, steps, n_paths, workers8, chunk250_000, seed0): shm shared_memory.SharedMemory(createTrue, sizen_paths * 8) out np.ndarray((n_paths,), dtypenp.float64, buffershm.buf) starts list(range(0, n_paths, chunk)) # 用 SeedSequence 派生子种子保证可复现且互不重叠 seeds [int(s.generate_state(1)[0]) for s in np.random.SeedSequence(seed).spawn(len(starts))] jobs [(shm.name, n_paths, st, min(st chunk, n_paths), S0, mu, sigma, T, steps, sd) for st, sd in zip(starts, seeds)] try: with Pool(workers) as pool: pool.map(_worker, jobs, chunksize1) result out.copy() # 拷出来再释放共享块 finally: del out # 先删 ndarray 引用 shm.close() shm.unlink() # 主进程负责销毁 return result几个参数和顺序都不能乱。sizen_paths * 8是按 float64 算的字节数换 float32 要同步改成 4两边不一致会直接抛ValueError: buffer is smaller than requested size。SeedSequence(seed).spawn(n)是关键如果用seed i当每个 worker 的种子相邻种子在 PCG64 下的前几个输出相关性强等于人为制造了路径相关性VaR 会偏小。spawn派生的是统计独立的子序列可复现性也保住了——同一seed跑两遍结果逐位相同。写回时out[start:stop] terminal用的是切片赋值不产生中间数组如果写成out terminal只是重新绑定名字共享区里那片内存还是空的主进程拿到的是一堆 0。这是我自己早期犯过的错调了半天以为是并行没跑。finally里的顺序是del out→close→unlink。少了delndarray 的 buffer 还引用着共享区unlink在某些平台上会抛BufferError。worker 里只能close绝不能unlink——那会把还在被其他进程使用的块直接删掉。unlink只在主进程、且只做一次。提示Windows 上SharedMemory(name...)的name不要带路径分隔符用默认生成的随机名最稳跨机器或跨容器的场景不要共享内存改用np.memmap落到磁盘。4.3 和 MATLAB parfor 版对拍顺带把分块调优做掉原版是 MATLAB 脚本核心循环用parfor并行口径和这里一致。对拍不能逐条路径比因为rng(..., twister)和 numpy 的 PCG64 生成的随机序列完全不同步要比的是统计量。% 原版 MATLAB 核心片段仅用于对拍口径 rng(20240101, twister); Z randn(nPaths, steps); incr (mu - 0.5*sigma^2)*dt sigma*sqrt(dt)*Z; ST S0 * exp(sum(incr, 2)); % 注意沿第 2 维求和 varValue -prctile(ST - S0, (1-confidence)*100);对拍时固定S0、mu、sigma、T、nPaths两边各跑 20 次取均值。我能接受的一致性标准是var的相对差在 1% 以内es在 2% 以内且随nPaths增大差异单调收窄。如果差异不收敛一定是口径问题不是随机误差。对不上的地方主要有三处列出来备查差异点MATLABPython影响分位数插值prctile默认方法np.percentile(methodlinear)千分之几路径少时明显标准差自由度std默认 n-1std(ddof1)忽略会系统性低估 sigma矩阵方向列优先sum(...,2)行优先sum(axis1)写反直接结果作废调优部分我在 8 核测试机上跑了一组固定 200 万条路径、steps1workerschunk_paths耗时(s)加速比1200000012.41.00210000006.61.8845000003.53.5482500002.15.90161250002.45.178 进程之后开始回落因为测试机只有 8 个物理核超出的进程在抢资源。加速比到不了线性的 8 倍剩下的损耗来自共享内存映射和进程启动路径数越大这部分占比越小。要是把chunk_paths压到 5 万8 进程的总耗时反而会涨到 3 秒以上。5. 避坑与常见问题排查五个我真实踩过的坑5.1 进程开太多反而变慢现象把workers从 8 调到 32总耗时从 2.1 秒涨到 3.8 秒机器风扇狂转但 CPU 利用率上不去。原因物理核只有 8 个多出来的 24 个进程在做时间片轮转每个进程的共享内存 attach 和随机数生成器初始化都要完整走一遍纯开销。另外每块路径数被摊薄到 6 万多条固定成本占比飙升。解决workers上限锁死在multiprocessing.cpu_count()的物理核数配置里加断言chunk_paths保持 10 万以上。真要压榨宁可开 8 个进程每块 50 万也别开 16 个进程每块 12 万。5.2 Windows 下启动就抛 RuntimeError现象Linux 上跑得好好的换到 Windows 直接报An attempt has been made to start a new process before the current process has finished its bootstrapping phase。原因Linux 默认fork子进程继承父进程内存Windows 和 macOS 新版默认spawn子进程要重新 import 主模块。如果Pool(...)写在模块顶层import 时会再次执行到创建进程池那一行无限递归Python 用这个异常把它拦下来。解决所有进程池创建代码必须放在if __name__ __main__:下面_worker函数定义在模块顶层而不是闭包里——闭包里的函数 spawn 时 pick 不到。实在绕不开用multiprocessing.get_context(spawn)显式指定行为在三个平台上就一致了。5.3 共享内存泄漏与 BufferError现象跑几十次之后/dev/shm被占满或者程序退出时报BufferError: cannot close exported pointers exist严重时下一次运行直接FileExistsError。原因SharedMemory是操作系统级资源进程正常退出不会自动清理必须显式unlink。另一种情况是 ndarray 还引用着这块 buffer 就调了close引用计数没归零资源挂在那里。解决try/finally包住整个计算过程finally里严格按del ndarray→shm.close()→shm.unlink()三步走。不要图省事只写close也不要在 worker 里调unlink。调试阶段可以在启动时扫一遍/dev/shm把残留的psm_*清掉。5.4 VaR 数值和 MATLAB 版对不上现象同一组参数MATLAB 出 0.0231Python 出 0.0227差 1.7%而且加大路径数差异不收敛。原因查下来是三处叠加。一是np.percentile的method没指定默认linear和 MATLABprctile的插值方式不同二是年化因子一边用 252 一边用 250三是最隐蔽的——MATLAB 那边用的是简单收益Python 这边估参数用的是对数收益两条链路的前提假设不一样。解决把插值方法、年化因子、收益类型三件事写进配置并打印到日志里对拍前先核这三项。年化因子对 sigma 的影响是线性的收益类型的影响在波动率大的标的上能到 5% 以上不是小数。5.5 路径矩阵把内存吃满现象n_paths调到 500 万机器卡死dmesg里能看到 OOM killer。原因gbm_paths一次性生成(5000000, 252)的 float64 矩阵将近 10GB再加上incr和log_paths两个中间数组峰值能到 30GB。解决欧式场景直接换gbm_terminal_fast内存降到 40MB。确实要中间路径的按chunk_paths分块生成、算完即弃或者用float32把内存砍半——代价是尾部精度损失99% 分位数上要谨慎。别指望 swap 能救这类大数组一换出就是灾难。6. 二次开发加对偶变量、换厚尾分布、做回测验证拿到项目后最值得动的一刀是方差缩减。蒙特卡罗的误差按1/sqrt(n)收敛路径数翻 4 倍才能把标准误减半代价是 4 倍时间。对偶变量法不增加路径数只让一半路径用Z、另一半用-Z同样的计算量下标准误能降 30% 到 50%。# 在 gbm_terminal_fast 基础上改造 def gbm_terminal_antithetic(S0, mu, sigma, T, n_paths, rng): half n_paths // 2 z_half rng.standard_normal(half) z np.concatenate([z_half, -z_half]) # 对偶配对 return S0 * np.exp((mu - 0.5 * sigma ** 2) * T sigma * np.sqrt(T) * z)n_paths必须是偶数奇数时最后一条会被丢掉half用整除在配置校验里加一条if n_paths % 2: raise。换厚尾分布也简单把标准正态换成自由度 5 到 7 的 t 分布再标准化尾部就会明显变厚99% VaR 会往上走一截from scipy import stats z stats.t.rvs(df5, sizen_paths, random_staterng) z z / np.sqrt(5 / 3) # t(5) 方差是 df/(df-2)5/3标准化到方差 1不标准化的话波动率和 t 分布的方差叠在一起sigma的物理含义就变了VaR 会莫名偏大又是一个隐蔽的坑。验证方法别只看单次结果要跑重复实验固定参数跑 200 次统计var的均值和标准误比较改造前后的收敛速度再拿真实历史数据做 Kupiec POF 回测看例外率是否落在置信区间里def kupiec_pof(exceptions: int, n: int, confidence: float): 返回 LR 统计量与卡方分布 3.84(df1, 95%) 比较 p 1 - confidence ph exceptions / n return -2 * (np.log((1 - p) ** (n - exceptions) * p ** exceptions) - np.log((1 - ph) ** (n - exceptions) * ph ** exceptions))exceptions是实际损失超过 VaR 的天数。LR 超过 3.84 说明模型在统计上不可接受例外率偏高就是低估风险偏低则是模型过于保守——后者不一定是好事资本占用会虚高。这套实验我一般跑 250 个交易日样本太短统计量没意义。从那以后我每次改模型都会先跑一遍对拍脚本、再跑一遍回测两关都过才敢把新版本挂到日报上光看单次数字顺眼就上线迟早翻车。希望帮到你。本文还有配套的精品资源点击获取