能源价格风险管理尤其是电力价格的风险建模是我近期用MATLAB从零复现的一个项目。外行看到“电力价格”会觉得它和普通大宗商品差不多但实际上电力价格的波动性远远大于传统商品刚拿到数据时那种日内冲高几倍、隔几天又跌到负值的情况会让你怀疑是不是统计口径出了问题。而所谓“按高水平文章复现源代码”核心就是把论文里那几行优雅的随机微分方程落地成可运行、可校验、可解释的MATLAB代码这个过程中要补的细节远比想象中多。文章会按完整的项目推进顺序来写先说电力价格为什么难建模再讲模型怎么选、为什么这么选然后给出整套源代码的实现思路和关键代码片段最后专门讲复现过程中踩过的坑以及怎么证明你的代码是对的。适合正在做电价预测、电力市场风控、或者论文和毕设需要复现能源价格模型的同学参考。1. 电力价格风险管理难在哪里搞清楚建模基础再上手1.1 电力价格波动的三个特殊性和两个必须处理的特征我在项目一开始就吃了“惯性思维”的亏总想着套用股票价格里的几何布朗运动结果被真实电价数据教做人。先说三个特殊性。第一电力几乎不能大规模存储发出来就要用掉所以供需必须实时平衡任何一侧的突然波动都会直接反映到价格上。第二电力需求的短期价格弹性极低用户不会因为一小时涨五倍就不开空调这意味着价格调整几乎全靠供给侧响应。第三电力市场存在物理网络约束输电阻塞会让不同节点的价格出现明显差异这也就是人们常说的“节点电价”或者“分区电价”做风险管理时必须想清楚手里数据属于哪个区域。有了这三个特殊性电价序列自然呈现出两个必须建模的特征一是均值回归电价在大部分时间会围绕发电成本、供需关系形成的某个长期均值上下波动偏离后会被市场拉回来二是跳跃和尖峰极端天气、机组故障、需求突增都会让价格在短时间内飙升这种跳跃不是正态分布能描述的传统的连续扩散模型根本接不住。做模型设计之前建议先对你的数据集做个简单统计算出日收益率的峰度、偏度画一张直方图和QQ图。如果峰度明显大于3偏度明显右偏那基本可以判定纯正态扩散模型是不适用的。这一步花十分钟能帮你省下后面几天和参数纠缠的时间。1.2 复现高水平论文之前先梳理代码模块清单“按高水平文章复现”听起来像是个技术活实际上更像是把一篇论文拆成一个一个能跑的模块。我在动手写代码前先画了一张纸上的流程图把整个项目分成五个模块后来项目全程都是按这个清单走的非常有帮助。第一个模块是数据预处理包括缺失值处理、异常值识别、对数价格转换。第二个模块是参数估计这一步最核心因为论文里只会给你一个连续时间模型你需要自己推导离散化形式再选择OLS、极大似然或者其他估计方法。第三个模块是情景生成也就是蒙特卡洛模拟用估计出的参数生成大量未来价格路径。第四个模块是风险度量计算把模拟路径转化为VaR、CVaR、最坏损失这类风险指标。第五个模块是回测验证把模型预测结果和真实数据对比验证覆盖率和模型稳定性。这五个模块之间的接口一定要清晰。我的习惯是每个模块写成独立的函数文件输入输出都用结构体打包这样后面改参数估计方法或者换一种模型结构不需要动其他模块的代码。尤其是参数估计和模拟这两个模块数据接口一定要统一比如模拟函数接收一个params结构体里面存kappa、mu、sigma、lambda等字段这样新模型只要往结构体里加字段就行主脚本基本不用改。2. 选对模型骨架均值回归、跳扩散与机制转换怎么权衡2.1 均值回归OU过程刻画电价向成本回归的基本盘电力市场的长期经济逻辑决定了电价会向发电边际成本回归这是均值回归模型的经济学基础。最常用的是Ornstein-Uhlenbeck过程通常建模在对数价格上表达式为d lnP κ(μ - lnP)dt σdW其中κ是回归速度表示价格偏离均值后以多快的速度被拉回μ是长期均值σ是波动率。κ越大价格回归越快κ趋近于0时模型退化成几何布朗运动。在实际写代码时我们需要把连续时间方程离散化。对OU过程做精确离散化后可以得到一阶自回归形式lnP_t a b * lnP_{t-1} ε_t其中b e^{-κΔt}a μ(1 - e^{-κΔt})ε_t的方差是σ²(1 - e^{-2κΔt})/(2κ)。这个形式的好处是直接用OLS回归就能估计a和b再反解出κ、μ、σ实现成本极低。不过要注意这个模型只能刻画电价的正常波动区间对尖峰几乎没有解释能力。在电价数据里跳跃往往只占百分之几的时间却贡献了很大的风险值。只用OU过程算出来的VaR很可能是被严重低估的。2.2 跳扩散模型把尖峰做进随机过程为了处理尖峰我选择了跳扩散模型一种在高水平能源价格模型中比较常见的结构。它是在均值回归基础上叠加一个泊松跳跃项dlnP κ(μ - lnP)dt σdW JdN其中dN是泊松过程λ表示跳跃强度即单位时间平均发生多少次跳跃J是跳跃幅度可以设定为服从正态分布或双指数分布。在电力市场中双指数分布往往能更好地刻画正负跳跃的不对称性但正态分布在实现上更简单对初学者更友好我建议先从正态假设入手。参数估计是跳扩散模型的难点。论文里可以在连续时间框架下推导出漂亮的条件但落到代码里你需要先判断哪些观测点是“跳跃点”再对跳跃部分和扩散部分分别估计。比较常用的做法是阈值法计算对数收益率的滚动标准差凡是对数收益率偏离均值超过一定倍数比如三倍滚动标准差的点就视为跳跃点剩下的数据用来估计均值回归参数。估计跳跃强度λ时直接统计跳跃点的数量除以总时间长度跳跃幅度的均值和标准差则用跳跃点收益率的均值和标准差来近似。这个方法虽然朴素但我在项目里实测下来对模拟数据的还原度令人满意且代码简单不容易出错。2.3 三种常见模型的优缺点对比与选型建议很多新手一上来就想用最复杂的模型总觉得模型越复杂越能发好文章或体现水平但实际项目里往往相反。我把常用模型的适用场景整理成了表格方便你对照选择。模型刻画能力优点缺点适用场景均值回归OU回归特性参数少、稳定、易实现忽略尖峰VaR低估价格稳定期初步建模跳扩散回归特性尖峰能刻画跳跃精度明显提升参数估计相对复杂含极端事件的常规市场马尔可夫机制转换多种状态切换区分正常态和尖峰态参数多容易出现识别问题结构性变化明显的市场我的建议是先老老实实把OU模型跑通拿到一组可解释的参数再在这个基础上加跳跃项。这样每一步都有对照即使后面模型出问题你也知道问题大概率出在哪儿。3. MATLAB源代码完整实现从数据清洗到风险指标输出3.1 数据准备对数价格与异常电价处理这一步看似基础实际上是决定参数估计质量的关键因为真实电价数据往往比教科书脏得多。我拿到的数据集里存在三个典型问题一是某些小时的数据缺失可能是市场停机或者数据上传失败二是出现明显的录入错误比如价格写成0.01又立刻变成1200三是极端尖峰到底是真实市场事件还是数据错误需要结合交易日志判断。数据处理的第一原则是任何合理模型都建立在数据是真实、连续、干净的假设上如果原始数据有问题后面的复现再漂亮也没有意义。我的处理顺序是先用线性插值补齐短期缺失值对于连续缺失超过三天的区间直接标记并剔除该阶段避免插值引入伪信息。对于负价格问题这里必须展开说。部分电力市场允许出现负电价这时对数价格变换完全失效。如果你的数据里有负值可选方案有三种一是把绝对值整体平移确保序列全部为正再取对数但解释变量时要注意偏移量二是改用算术OU过程直接对原始价格建模三是如果负值只占极少数可以考虑剔除后插值平滑。我比较推荐第二种直接扩展成算术模型毕竟在真实市场中负电价本身就是一个重要风险点直接剔除会低估尾部风险。预处理代码示例如下% price: a column vector containing spot prices % timestamps: corresponding datetime vector % Step 1: fill short gaps by linear interpolation nanIdx isnan(price); priceFill fillmissing(price, linear); % Step 2: detect negative prices, print for manual check negIdx priceFill 0; if any(negIdx) fprintf(Found %d negative prices, need modelling decision.\n, sum(negIdx)); end % Step 3: if no negative price, take log transform logPrice log(priceFill);这里用到了MATLAB内置的fillmissing函数这个函数比手动循环快得多而且能自动处理数据首尾的NaN非常省心。如果你用的是老版本MATLAB没有这个函数可以自己写一个基于线性插值的简单版本逻辑并不复杂。3.2 参数估计OLS和极大似然在MATLAB中的落地写法参数估计我分两条路来讲。第一条路是OLS适合OU模型代码短、速度快、结果稳定第二条路是极大似然适合跳扩散模型或者你想在论文里用更“标准”的方法时使用。OU过程的OLS估计代码非常简洁dt 1/365; % daily data, assuming 365 days per year Y logPrice(2:end); X [ones(length(Y), 1), logPrice(1:end-1)]; beta (X * X) \ (X * Y); a beta(1); b beta(2); residuals Y - X * beta; sigma2 var(residuals); kappa -log(b) / dt; mu a / (1 - b); sigma sqrt(2 * kappa * sigma2 / (1 - b^2)); fprintf(kappa %.4f\n, kappa); fprintf(mu %.4f\n, mu); fprintf(sigma %.4f\n, sigma);需要注意的是这里使用的是对数价格的“价格水平值”做回归而不是收益率。很多人第一次写会直接对logPrice(2:end)和logPrice(1:end-1)做差做成收益率回归那样得到的是扩散项参数不是均值回归参数两者含义完全不同。跳扩散模型用极大似然估计会更严谨。由于跳跃项的加入似然函数不再有解析形式需要用数值优化求解% define negative log-likelihood for jump-diffusion model function nll jumpDiffusionNLL(params, data, dt) kappa params(1); mu params(2); sigma exp(params(3)); % ensure positivity lambda exp(params(4)); muJ params(5); sigmaJ exp(params(6)); n length(data); nll 0; for t 2:n r data(t) - data(t-1); meanPart kappa * (mu - data(t-1)) * dt; diffVar sigma^2 * dt; % mixture density: no jump one jump p0 exp(-lambda*dt) * normpdf(r, meanPart, sqrt(diffVar)); p1 (1 - exp(-lambda*dt)) * ... normpdf(r, meanPart muJ, sqrt(diffVar sigmaJ^2)); nll nll - log(p0 p1 1e-12); end end然后调用fmincon或fminsearch做优化注意初始值要合理。我在项目里是从OLS估计结果出发kappa初值取0.5左右mu取历史均值sigma取收益率的滚动标准差lambda取0.05muJ取0附近sigmaJ取和sigma差不多的值。优化时要加上参数边界约束特别是sigma、lambda、sigmaJ必须为正。使用极大似然时有一个容易踩的坑对参数做log变换可以避免在优化过程中产生负值也就是我上面写的exp(params(3))这样的形式。这种技巧可以在不引入复杂约束的情况下保证参数始终为正优化速度也更快。3.3 蒙特卡洛模拟生成价格路径随机数种子的正确打开方式参数估计完之后下一步是生成海量未来价格路径。这里我用的是蒙特卡洛模拟原因在于非线性风险度量、跳跃项和路径依赖问题用解析公式很难处理而模拟的办法能把所有复杂度都吸收进去。模拟的核心代码思路如下rng(20240601); % set seed for reproducibility nScen 10000; % number of simulated paths nDays 30; % risk horizon in days dt 1/365; S0 logPrice(end); % start from latest observed log price simLogPrice zeros(nScen, nDays 1); simLogPrice(:, 1) S0; for t 1:nDays dW sqrt(dt) * randn(nScen, 1); jumpCount poissrnd(lambda * dt, nScen, 1); jumpSize normrnd(muJ, sigmaJ, nScen, 1); jumpTerm jumpCount .* jumpSize; simLogPrice(:, t1) simLogPrice(:, t) ... kappa * (mu - simLogPrice(:, t)) * dt ... sigma * dW jumpTerm; end simPrice exp(simLogPrice); % convert back to price level第一行rng固定随机数种子这是几乎所有复现类项目里最关键的细节之一。没有固定种子你每次跑出来的VaR都不一样就没办法排查代码改动带来的真实影响。我在项目里会把种子写成日期加版本号的形式比如rng(20240601)这样别人拿到代码也能完全复现你的结果。还有一点需要提醒的是poissrnd生成的是跳跃次数意味着一个时间步里可能发生多次跳跃。这和真实电力市场的情况比较符合极端天气下几个小时内的连续冲击可能形成一串尖峰。但如果你的步长比较大比如按天模拟λ*dt通常都小于0.1一次以上跳跃的概率很低影响并不大。模拟完成后一定要做基本合理性检查比如画出前100条路径看看尖峰分布是否符合你对电力市场的直觉统计一下模拟结果中出现负价格的占比如果占比过高而历史数据几乎没出现过说明模型设定有问题需要回头检查。3.4 VaR和CVaR计算把模拟路径变成风控指标有了未来价格路径风险指标的计算就顺理成章了。这里的核心概念有三层一是价格收益二是价格相对当前水平的损失三是组合或头寸的损失敞口。如果你需要管理的是一个实物购电合同持仓量是固定MW那么损失就是价格差乘以电量。假设持有100MW的购电仓位持有期为30天计算方式和代码如下exposureMW 100; loss (simPrice(:, end) - S0) * exposureMW; VaR95 -prctile(loss, 5); CVaR95 -mean(loss(loss -VaR95)); fprintf(95%% VaR %.2f\n, VaR95); fprintf(95%% CVaR %.2f\n, CVaR95);VaR的含义是在95%置信水平下未来30天内的最大可能损失不会超过这个数值CVaR是损失超过VaR的那些尾部场景下的平均损失。相比VaRCVaR对极端尾部更敏感这也是巴塞尔协议逐渐从VaR转向ESExpected Shortfall的原因如果你要给风控委员会汇报建议把两个指标都列出来。需要注意这里的loss用的是期末价格和当前价格的价差没有考虑中途每天的盯市实际应用中还要根据头寸的结算规则调整。如果是期货逐日盯市就要把每一天的损失路径都算出来然后再汇总逻辑会复杂一些。4. 复现正确性如何保证踩坑记录和验证套路4.1 新手最容易翻车的三个隐性错误第一类错误是时间单位不统一。论文里通常写年化参数比如kappa5表示一年回归5次sigma0.3表示年化波动率30%但数据是日度频率。如果在代码里dt用了1就会把年化参数当成日参数导致模拟结果完全失真。我的建议是一开始就统一所有参数的时间基准全部用年化单位代码里只出现一次dt的定义后面所有地方都引用这个变量。第二类错误是跳跃点和异常值的处置不当。用阈值法识别跳跃时阈值取太大会漏掉真实跳跃取太小又会把正常波动当跳跃导致均值回归参数估计失真。我的经验是先画出对数收益率的分布图观察尾部形态再取一个合理的分位数作为阈值比如99.5%分位数。不用一来就上3倍标准差因为真实数据的尖峰尾部比正态分布厚得多3倍标准差可能把很多跳跃都漏掉了。第三类错误是忽视负价格导致的对数变换失效。如果数据里存在负电价直接log会得到NaN程序可能不会报错但你的参数估计结果会全是NaN或无穷大。我建议在数据预处理模块里就增加一个显式检查一旦发现负价格立刻打印提示而不是等到后面模拟的时候才发现结果完全不对。下面这个表是我在项目中整理的高频错误速查表适合贴在代码旁边。现象常见根因排查方向参数估计结果异常大/小时间单位混乱dt用错检查dt定义和年化换算模拟路径全部为NaN数据里有负价格或零值检查对数变换前是否有非正值VaR结果忽大忽小没有固定随机数种子在模拟前固定rng种子模拟路径波动剧烈失真跳跃参数初始值不合理调整lambda初值限制sigmaJ范围极大似然优化不收敛参数空间过宽/目标过复杂用OLS结果做初值加参数边界约束4.2 用合成数据反验代码让模型自己证明自己我在项目中最坚持的一个习惯是用已知参数生成合成数据再用自己的估计程序把参数估回来对不上就找问题对上了才用真实数据。这种方法叫反验能在代码开发阶段就把大部分bug暴露出来而不是等风控报告出了问题才回头debug。反验的思路很简单rng(42); trueKappa 0.8; trueMu 3.2; trueSigma 0.4; trueParams [trueKappa, trueMu, trueSigma]; % generate synthetic price data using the same simulation function % then call the estimation function or script [kappaEst, muEst, sigmaEst] estimateOU(syntheticLogPrice, dt); fprintf(kappa true%.2f est%.4f\n, trueKappa, kappaEst); fprintf(mu true%.2f est%.4f\n, trueMu, muEst); fprintf(sigma true%.2f est%.4f\n, trueSigma, sigmaEst);如果估计出来的参数和真实参数差别很大不要急着改参数先检查数据生成函数和估计函数是否使用了相同的离散化公式。我在第一次做反验时就发现模拟时用的是精确离散化估计时用的是欧拉离散化两者在小步长下差距不大但在大步长下会导致估计结果系统性地偏离。后来我把两边统一成同一种离散化方式参数就基本上能回归到真实值了。反验通过以后还要做一次模拟均值与理论均值的对照。OU过程的条件期望是E[lnP_t | lnP_0] lnP_0 * e^{-κΔt} μ(1 - e^{-κΔt})你可以对比模拟10000条路径的平均值和这个理论值的差异通常应该在千分位级别如果偏差过大大概率是模拟代码的累积误差或者随机数生成方式有问题。4.3 性能优化当模拟路径从1万变成100万条初期验证阶段用1万条路径足够但到了正式风险报告阶段尤其是尾部风险指标1万条路径的抽样误差太大。想要把CVaR估计得比较稳通常需要50万到100万条模拟路径这时MATLAB代码的性能就变成了一个很现实的问题。第一招是向量化。上面的模拟代码里循环内部几乎所有的操作都是矩阵运算已经天然向量化了。但如果你还在用下面的循环方式就要特别注意% slow version to avoid for s 1:nScen for t 1:nDays ... end end这种双重循环在小规模测试时无所谓路径一多直接卡死。MATLAB的向量化能力很强应该把场景维度作为矩阵的行时间维度作为列用一次randn(nScen,1)生成所有场景在同一时刻的随机扰动然后利用矩阵运算符一次性推进所有路径。第二招是预分配数组。模拟前先zeros好整个矩阵避免在循环里不断扩展数组的长度。这个习惯在MATLAB里极其重要数组在循环里每次扩展都会触发内存重新分配性能会急剧下降。我在上面的代码里已经提前分配好了simLogPrice这是标准做法。第三招是当路径数量特别大、内存吃紧时可以考虑把模拟拆成多块比如每次模拟20万条循环5次最后合并结果。这样每个块所需的内存可控也能利用parfor做简单并行化。并行化的时候要注意随机数流问题每个worker要设置独立的随机子流避免并行计算产生重复的随机序列MATLAB的RandStream支持这个功能值得研究。5. 模型结果怎么用风险报告、压力测试和扩展方向5.1 把模拟结果转成决策可读的风险报告模型算出来的是一堆概率分布但管理层想知道的是“最坏情况下我们亏多少”“要不要买套保合约”。所以我习惯把输出整理成一张简洁的风险报告表这也是风控团队、交易团队都能看懂的格式。我的风险报告通常包含三块内容一是风险指标汇总列出不同置信水平下的VaR和CVaR二是极端情景描述比如历史上最严重的尖峰事件在未来模拟中出现的概率三是敏感性分析展示kappa、sigma、lambda分别变化百分之十时VaR的变化幅度。最后这一块特别重要因为管理层通常会问“你的模型如果参数估计有一点误差结果还可靠吗”用敏感性分析可以回答这个问题。在实际操作中我还会额外做一次压力测试。比如人为地设定未来一周内出现三次极端跳跃事件然后用模型模拟这种情况下的组合损失。压力测试和蒙特卡洛模拟的区别在于模拟是在自然概率下进行压力测试则是人为改参数看最坏情况会不会击穿公司的风险限额。电力市场里历史尖峰对VaR的影响往往滞后压力测试是对VaR这类统计指标的重要补充。5.2 后续扩展机制转换、机器学习回归和高维情景项目跑通之后我个人的建议是不要停止在基础模型上因为真实电力市场的结构性变化比任何单一模型能刻画的都要多。比如有些市场有明显的峰谷时段差异价格在不同时段的均值回归速度完全不同这时候可以考虑在模型中引入时间虚拟变量或者改用马尔可夫机制转换模型让模型自动识别“正常状态”和“尖峰状态”。此外如果你手头有负荷数据、新能源出力数据、燃料价格数据可以考虑把它们作为外生变量引入模型。比如在均值回归项中加入负荷偏离度因子模型会更有业务解释力。我见过比较实用的做法是用MATLAB优化工具箱训练BP神经网络建立“负荷—电价”映射然后把这个映射作为跳跃强度的驱动程序这样模型既能捕捉价格尖峰又能解释尖峰的来源。最后还想说一点关于代码工程化的体会。能源价格风险管理不是一次性分析而是需要反复跑、反复更新的流程。所以代码从一开始就要按工程标准来写模块化、加注释、固定随机种子、输出可复现的结果文件这些基础工作会给你后面每一次更新节约大量时间。我自己最大的进步不是学会了更多高级模型而是学会了让每一段代码都经得起回头检验。最后分享一个项目里的习惯所有源代码文件名带版本号每次跑完就保存一份estimate_results.mat里面存参数估计值、随机数种子、数据区间和代码版本号。风险管理模型不怕复杂怕的是算完不知道对不对。数据和代码的可追溯性永远比模型本身的花哨程度更重要。