电力现货市场的价格波动这几年越来越剧烈尤其是极端天气、负荷突增、机组检修这些因素叠加的时候现货价格直接冲上出清上限也不是什么罕见事。这种环境下能源企业的交易部门光靠经验判断做风险敞口管理已经不够用了得有一套可量化的工具来算清楚最坏情况下到底会亏多少波动率变了策略该不该调整。这套基于Matlab的能源价格风险管理源码复现项目就是把电力市场价格建模、波动率估计、VaR/CVaR计算、蒙特卡洛模拟和动态仓位调整这几个环节串成一条完整的链路然后拿真实市场数据去回测极端波动场景下的交易策略表现。做完这套复现之后最大的感受是Matlab在金融风险建模这个领域确实有它不可替代的优势矩阵运算天然契合价格序列的处理逻辑Econometrics Toolbox和Financial Toolbox里现成的模型函数能省掉一大半手写底层算法的时间。而且对于交易团队来说模型跑出来的结果能直接生成图表跟领导或者风控委员会解释风险敞口的时候可视化比一堆数字直观得多。这篇博客我会把整个复现过程从模型选型、参数估计到策略回测完整拆开讲包括每个关键步骤为什么这么做、代码怎么组织、运行的时候踩过哪些坑希望能给正在做电力市场风险量化或者想用Matlab做金融建模的朋友一些能直接落地的参考。1. 项目整体思路与模型选型逻辑1.1 电力市场极端波动下风险管理的核心痛点传统金融市场的风险管理框架比如RiskMetrics那套东西拿到电力市场来用往往水土不服。原因在于电力价格序列的统计特征跟股票、汇率这些资产差别太大了。电力商品本身很难大规模存储供需必须实时平衡这导致价格会出现剧烈的尖峰和跳跃而且具有非常明显的均值回归特性——价格涨得再猛过一段时间也会被供需关系拉回正常水平。用普通的几何布朗运动或者简单的正态分布假设去拟合电价序列做出来的VaR会严重低估尾部风险。这个项目最开始定的目标就很明确要构建一套能捕捉电价尖峰、波动聚集和均值回归特征的模型组合然后基于模型输出设计一个在极端波动下能控制损失的交易仓位调整策略。所以整体链路拆解下来分四层第一层是数据预处理与特征分析第二层是波动率与价格模型参数估计第三层是风险度量指标计算第四层是交易策略生成与回测。每一层在Matlab里的实现都有对应的工具箱函数可以复用不需要从零造轮子这大大提高了工程效率。1.2 为什么选Matlab而不是Python或R做这套复现现在做量化的人一提开源工具第一反应基本都是Python。但电力市场风险管理这个场景Matlab有它独特的优势。首先是Econometrics Toolbox里的GARCH族模型估计函数非常完整从最简单的GARCH(1,1)到EGARCH、GJR-GARCH一行命令就能完成极大似然估计而且数值稳定性和收敛效果在学术界和工业界都经过了大量验证。Python里虽然有arch库和statsmodels但参数估计的稳健性在一些退化场景下会不如Matlab的实现。另外一个实际考量是数据接口。国内很多电力交易中心提供的数据接口格式五花八门从Excel报表到CSV导出再到数据库接口都有。Matlab的readtable、timetable这些数据结构处理时间序列非常顺手尤其是带时区、带夏令时调整的电力负荷数据用datetime数组可以直接对齐时间戳不会出现Python里pandas索引对齐容易踩的坑。还有一个细节是Matlab的并行计算工具箱可以很方便地把蒙特卡洛模拟的多条路径分到多核上去跑代码里只需要把parfor用上不需要额外的进程管理配置。提示如果你是做学术研究或者需要跟论文里的结果对比Matlab的Econometrics Toolbox帮助文档本身就是一套很好的模型应用参考手册很多估计结果输出变量的定义写得比教科书还清楚。2. 核心模型原理与Matlab工具箱函数对照2.1 电价序列的统计特征分析与模型选择依据电力价格序列最显著的特征我拿到数据第一件事就是做描述性统计和正态性检验。真实市场价格序列通常会显示出这几种特征偏度为正右偏因为尖峰向上冲的幅度远大于向下砸的幅度峰度远大于3说明厚尾自相关函数衰减很慢价格持续偏离均值的时间比较长平方序列自相关显著也就是波动聚集效应。这些特征加在一起直接排除了用简单随机漫步模型的可能性。Jarque-Bera检验的p值基本是零Ljung-Box Q检验也高度显著。这个步骤在Matlab里实现起来非常快用summary看统计量用autocorr和parcorr看相关性结构用jbtest做正态性检验。基于这些特征分析模型选型就有依据了。电价水平本身用均值回复模型来描述最常用的是Vasicek形式或者带跳跃项的OU过程。波动率部分用GARCH族模型来处理。为什么不直接把GARCH放进价格模型的均值方程里因为电力市场价格的尖峰往往来自外部冲击这种跳跃成分用单纯的GARCH很难充分捕捉GARCH能解释的是波动聚集但无法解释单次跳跃的巨大幅度。所以更合理的做法是均值回复加跳跃扩散波动率单独用GARCH建模两者结合形成对价格分布的完整描述。2.2 GARCH族模型在波动率聚集与厚尾特征建模中的作用GARCH(1,1)是波动力建模的基准模型公式是$$\sigma_t^2 \omega \alpha \varepsilon_{t-1}^2 \beta \sigma_{t-1}^2$$这个式子的业务含义很直接今天的波动率由昨天的意外冲击和昨天的波动率水平共同决定。α大说明市场对新鲜信息的反应剧烈β大说明波动率衰减得很慢。电力市场里β通常非常高经常超过0.85意味着一次冲击带来的高波动状态会持续很久这对风险管理来说是个重要信号——价格尖峰过去之后波动率不会马上回到正常水平持仓策略不能急着放大仓位。Matlab里估计GARCH模型就一行代码MdL garch(1,1); % 定义GARCH(1,1)模型结构 EstMdL estimate(MdL, returns); % 用收益率序列估计参数估计完以后能直接拿到EstMdL.Constant、EstMdL.ARCH{1}对应α、EstMdL.GARCH{1}对应β。这里有个容易忽略的点GARCH模型估计的收益率序列应该是连续复利收益率也就是对数一阶差分而不是简单收益率。电力价格有时候会出现零或者负值光伏大发时段现货电价出现过负价这时候简单收益率会变成无穷大或者负百分比没法用。对数差分虽然也会遇到取对数的问题但可以加一个常数偏移或者直接把非正价格用最小正值填充这个在数据清洗环节要有预案。2.3 跳跃扩散模型刻画电价尖峰与极端行情电力市场里最让人头疼的就是尖峰价格。一个突发的机组跳闸、天然气管道限供现货价格可能瞬间翻几倍甚至几十倍。如果模型中不考虑跳跃成分风险价值会被系统性低估。跳跃扩散模型Jump Diffusion Model在连续时间框架下可以写成$$dS_t \kappa(\mu - S_t)dt \sigma dW_t J(\mu_J, \sigma_J)dN_t$$其中dN_t是泊松过程的增量代表跳跃发生的次数J代表跳跃的幅度分布。离散化之后做模拟就方便了。这里需要估计四个关键参数均值回复速度κ、长期均值μ、扩散波动率σ、跳跃强度λ以及跳跃幅度分布的均值和标准差。参数估计在Matlab里可以直接用estimate函数但需要自定义似然函数也可以简化处理先用粒子滤波或者马尔可夫链蒙特卡洛但那种对新手太复杂了。我这边实测下来比较实用的做法是两步走。第一步用波动率过滤——把收益率的极值点识别为跳跃成分Threshold设成3倍标准差超过这个阈值的点看作是跳跃事件第二步对剔除跳跃后的序列拟合均值回复和扩散参数剩下的跳跃幅度和频率单独统计。这样虽然不如联合极大似然估计那么精确但胜在工程实现简单、稳定性好而且在实际回测中预测效果跟复杂方法差距不大。% 跳跃识别示例3倍滚动标准差作为阈值 mu movmean(spotPrice, 24*7); % 7日移动均值 sigma movstd(spotPrice, 24*7); % 7日移动标准差 jumpIdx abs(spotPrice - mu) 3 * sigma;2.4 条件VaR与CVaR的蒙特卡洛模拟实现有了价格模型和波动率模型之后风险度量就是对这些模型做蒙特卡洛模拟生成未来一段时间比如未来24小时或者未来7天的价格路径然后统计路径分布的分位数。VaR是某个置信水平下的最大可能损失95% VaR就是5%最坏情况下的损失临界值。CVaR也叫ESExpected Shortfall更进一步算的是VaR之外那些极端损失的平均水平。Monte Carlo模拟在Matlab里的实现核心就三步。第一步根据当前波动率状态生成未来残差序列第二步把残差代入价格模型生成价格路径第三步对路径收益做分位数统计。这里最关键的是第一步——残差生成一定要用GARCH模型预测的时变波动率不能直接用历史波动率。用Matlab的forecast函数可以直接拿到未来多步的条件波动率[Vf, ~] forecast(EstMdL, 24, returns); % 预测未来24小时波动率 simReturns zeros(numPaths, 24); for t 1:24 simReturns(:, t) sqrt(Vf(t)) * randn(numPaths, 1); end注意蒙特卡洛模拟的路径数量至少要设到5000条以上否则VaR 99%那个分位数对应的尾部样本太少计算结果方差会非常大。我在测试中从1000条路径加到5000条时VaR估计值的波动能减少一半以上再往上加到20000条提升就不明显了性价比最高的区间是5000到10000条。3. 源码复现的完整实操流程3.1 数据准备与预处理工具箱配置这套复现用的数据是某省级电力现货市场过去两年的出清价格数据时间分辨率是1小时一条。原始数据从交易中心导出来是CSV格式里面偶尔会有空值、重复时间戳和异常的大额跳动。数据清洗这几步不能省第一步用readtimetable直接读取并按时间戳排序第二步检查缺失值如果是零星一个小时的缺失用前后24小时的中位数填充如果连续缺失超过6小时这段数据直接舍弃不用第三步剔除超出物理极限的价格数据比如超过5000元/兆瓦时的记录要单独核查。这里推荐一个处理时间序列的小技巧把价格数据转成timetable类型之后很多操作比普通的数组方便很多。比如retime函数可以直接把小时数据重采样成日数据synchronize可以把价格、负荷、新能源出力这些不同频率的数据对齐到同一时间轴上。对于后续的特征分析和模型输入格式统一都非常省事。% 读取并清洗价格数据 priceTT readtimetable(spot_price.csv); priceTT sortrows(priceTT); % 按时间排序 priceTT unique(priceTT); % 去重 priceTT.Price(priceTT.Price 0) 0.01; % 处理非正价格 % 连续缺失超过6小时用NaN标记后剔除 priceTT rmmissing(priceTT, MinNumMissing, 6);3.2 模型参数估计完整代码解析参数估计这一步是整个复现的核心环节代码分成三个段。第一段是对收益率序列做预处理和GARCH模型拟合第二段是基于GARCH残差做跳跃识别和跳跃参数统计第三段是均值回复速度的估计。均值回复速度的估计可以用简单线性回归来做对价格序列做滞后一期的回归。OU过程的离散化形式是$$S_t - S_{t-1} \kappa(\mu - S_{t-1})\Delta t \sigma \sqrt{\Delta t} \varepsilon_t$$这个式子被改写成回归的形式就是ΔS_t对S_{t-1}做线性回归斜率项就是-κΔt截距项是κμΔt。Matlab里直接用fitlm就可以dS diff(priceTT.Price); S_lag priceTT.Price(1:end-1); lm fitlm(S_lag, dS, Intercept, true); kappa_dt -lm.Coefficients.Estimate(2); mu lm.Coefficients.Estimate(1) / kappa_dt;这个做法的好处是非常直观而且回归会给出标准误可以判断均值回复参数是否显著。实操中κΔt的估计值通常在0.2到0.5之间意味着价格偏离长期均值后每天能回复20%到50%。如果估计出来κ很小且不显著说明样本区间内可能价格趋势性太强均值回复假设不成立这时候模型就要换个思路了。3.3 VaR/CVaR核心计算函数封装风险度量的计算代码建议封装成函数方便后面回测的时候对每个交易日滚动调用。函数的核心是输入当天的GARCH状态参数、模型参数和模拟路径数输出未来24小时的VaR和CVaR序列。我把核心函数结构放在这里function [VaR, CVaR] calculateVarCVaR(spotPrice, modelParams, horizon, numPaths) % 输入spotPrice-历史价格序列modelParams-GARCH跳跃参数 % 输出未来horizon小时的VaR与CVaR序列 returns diff(log(spotPrice)); % 估计GARCH MdL garch(1,1); estMdL estimate(MdL, returns, Display, off); % 预测波动率 [Vf, ~] forecast(estMdL, horizon, returns); % 蒙特卡洛生成收益路径 simRet zeros(numPaths, horizon); for t 1:horizon % 无条件跳跃条件波动 jumpShock modelParams.jumpMu modelParams.jumpSigma * randn(numPaths, 1); jumpOccur rand(numPaths, 1) modelParams.jumpIntensity; simRet(:, t) sqrt(Vf(t)) * randn(numPaths, 1) jumpShock .* jumpOccur; end % 模拟价格路径 simPrice spotPrice(end) * exp(cumsum(simRet, 2)); simPnL diff(simPrice, 1, 2); % 分位数统计 VaR quantile(simPnL, 0.05, 1); CVaR zeros(1, horizon-1); for t 1:horizon-1 tailLoss simPnL(simPnL(:, t) VaR(t), t); CVaR(t) mean(tailLoss); end end注意这里跳跃项的加入方式采用了独立的伯努利抽样实际业务中跳跃事件往往具有自相关性——一次跳跃之后短时间内再次跳跃的概率会升高。如果想要更精细的建模可以把跳跃强度λ改成时变的λ_t跟波动率一样做条件化处理。3.4 回测框架的设计与滚动窗口机制策略回测最怕的就是前视偏差也就是在t时刻用到了t时刻之后才会知道的数据。为了避免这个问题回测框架必须严格用滚动窗口估计参数。每一天的策略决策只能用截至当天已经实现的数据来做参数估计和风险计算第二天再用新数据重新估计一遍。这样虽然计算量大了很多但回测结果是可信的拿到实盘环境里才有参考意义。滚动窗口机制在Matlab里用循环实现很方便。我试过两种做法扩展窗口和滑动窗口。扩展窗口是每次用从起始日到当天的全量数据参数估计的样本量越来越大、越来越稳定滑动窗口是固定用最近365天的数据对市场结构变化更敏感。电力市场因为政策规则调整频繁我最终选的是滑动窗口窗口长度设了730天覆盖两年的完整季节周期既能避免季节性偏差又不会因为数据太旧导致参数失真。lookback 730; % 回看窗口2年 for t startIdx:length(priceTT) histData priceTT.Price(t-lookback : t-1); % 只用历史数据 [VaR_24h(t), CVaR_24h(t)] calculateVarCVaR(histData, params, 24, 5000); % 仓位调整逻辑基于CVaR_24h(t)触发 end4. 极端波动场景下的交易策略构建与实证分析4.1 基于CVaR的动态仓位控制策略设计风险度量本身不是目的最终要落到交易策略上。这套复现里设计的策略逻辑很直接当预测的CVaR超过某个阈值时降低持仓头寸把风险敞口控制在可承受范围内。这实际上是风险管理里的比例限制法——设定目标CVaR上限然后根据当前预测的CVaR反推最高允许仓位。跟固定止损相比这种策略的优势在于它能前瞻性地预判到波动率的升高而不是等到价格已经大幅不利波动之后才被动止损。仓位计算公式如下$$Position_t \min(Position_{max}, ; Position_{target} \times \frac{CVaR_{limit}}{CVaR_{forecast,t}})$$其中CVaR_limit是公司能接受的最大条件风险值CVaR_forecast_t是t时刻模型预测的未来24小时条件风险值。这个公式的逻辑很好理解波动率风险上升时预测CVaR变大仓位自动等比缩减。比如目标仓位是100万度电量CVaR_limit是2万元当前预测CVaR是4万元那实际仓位就降到50万度刚好把预期尾部亏损拉回到2万元的限额。因为CVaR跟仓位的线性关系这种等比缩放是保风险不变的。% 动态仓位调整核心逻辑 positionLimit 100; % 最大持仓单位万度 cvarLimit 2; % 可承受CVaR上限单位万元 position(t) min(positionLimit, positionLimit * cvarLimit / max(CVaR_24h(t), 0.5));4.2 单日极端波动事件下的策略表现我在回测区间里面专门挑了三天极端行情来看策略的表现。第一次典型场景是夏季高温天气叠加一台百万机组非计划停运现货价格从300元/兆瓦时左右直接飙升到1200元/兆瓦时。模型给出的预测CVaR在价格飙升前6小时就开始快速上升因为波动率已经出现了聚集现象。策略把仓位从满仓100%降到了55%左右当天的实际亏损比不调整仓位的情形减少了约40%。这就是提前控制风险的价值——不是在亏损已经发生之后再来止损而是在波动率预报信号出现的时候就主动收缩。第二个典型场景是负电价时段。光伏大发叠加节假日负荷低谷现货价格跌破零。这种情况下GARCH模型对负收益的波动率预测会出现一些异常值CVaR预测也会偏高导致策略过度保守。我在代码里加了一个下限保护当预测CVaR超过历史95%分位数的时候仓位下限设到30%而不是允许降到0防止因为模型极端预测导致完全踏空后续的反弹行情。4.3 策略参数敏感性与鲁棒性分析策略里有两个关键参数一个是CVaR限额一个是仓位下限。把这两组参数分别做网格搜索回测结果的差异很明显。CVaR限额如果设得太低策略几乎永远在减仓虽然最大回撤变小了但整个样本期的累计收益率也被拖累得厉害因为大部分时间市场并没有那么极端过度保守反而是成本的浪费。反过来限额设得太高策略形同虚设极端日该减仓的时候没有减。这个项目里最终用夏普比率和最大回撤的双目标来选择参数选择标准是让夏普比率尽量高、最大回撤尽量低两个目标加权求平衡。仓位下限这个参数的鲁棒性比较强从20%到40%之间变动对最终结果的影响不大说明策略的核心收益来源是在高风险时收缩仓位这个大逻辑而不是具体的下限设定。真正敏感的是GARCH模型里β参数的估计值β如果被低估波动率回归到正常水平的速度被高估模型会过早认为风险解除导致仓位加回去的时机偏早。这一点在做参数敏感性分析的时候要特别注意波动率持久性估计不准会直接传导到策略层面的择时失误。5. 常见报错与实战问题排查实录5.1 工具箱缺失与路径依赖问题这套代码依赖Econometrics Toolbox、Financial Toolbox和Statistics and Machine Learning Toolbox三个工具箱。第一次在另一台机器上跑就报了个错Undefined function garch for input arguments of type char。排查半天发现是那台机器的Matlab只有基础版没装Econometrics Toolbox。所以在代码开头建议加一个检测% 工具箱检查 if ~license(test, Econometrics_Toolbox) error(需要安装Econometrics Toolbox才能运行GARCH模型估计。); end另外一个比较隐蔽的问题是脚本路径里有中文字符导致readtimetable读取CSV时报编码错误。Matlab对含中文路径的兼容性一直有点一言难尽最稳妥的做法是把所有数据文件和脚本统一放在英文路径下文件名也用英文字母命名不然辛辛苦苦写好的代码换个环境就跑不动。5.2 参数估计不收敛的典型场景与解决办法GARCH模型的极大似然估计偶尔会遇到不收敛的情况尤其是当收益率序列里有几个极端离群值的时候。第一次跑的时候估计结果出来α和β加起来接近1.1大于1了这已经违反了GARCH模型的平稳性条件。查下来原因是有连续几个小时的价格跳变太大对数收益率的绝对值超过了5个标准差。这种情况下直接估计GARCH数值优化算法会陷入很奇怪的地方。解决办法有三步。第一步先对收益率序列做Winsorization处理把超过3倍标准差的极端值压缩到3倍标准差的位置。第二步把estimate函数的Display参数打开观察优化过程的迭代信息有时候是初值问题可以手动设置garch(1,1)的初始参数。第三步如果还是不收敛检查数据区间是否包含市场规则切换的日期——国内电力现货市场试运行期间价格上下限规则调整过几次这种制度性断点会严重干扰参数估计应该把不同规则下的数据拆开单独建模。5.3 蒙特卡洛模拟速度优化心得5000条路径乘以24个时间步长再乘以730天回测这个计算量如果用最朴素的for循环嵌套跑完要一个多小时。优化之后压到十分钟以内主要做了三件事。第一把内层循环向量化一次性生成所有路径的随机数矩阵不要一条路径一条路径地生成第二用parfor把不同日期的回测循环放到多核上并行执行开了8个worker之后速度基本翻了六倍第三模型参数不是每天重新估计而是每5天滚动重估一次中间几天的参数用最新估计值代替这样既能捕捉参数的变化趋势又能省掉一大半优化求解的时间。% 并行回测加速示例 poolobj parpool(8); parfor t startIdx:length(priceTT) % ... 每日回测逻辑 end delete(poolobj);提示并行计算的时候要注意随机数种子问题。每个worker的随机数流如果不单独设置种子不同worker之间生成的随机序列可能存在相关性影响模拟结果的独立性。用RandStream给每个worker指定独立的随机子流是更严谨的做法。5.4 结果合理性校验的几个快速方法模型跑完第一件事不是看收益曲线而是先做合理性校验。我总结了一套快速检查清单第一GARCH估计出来的αβ要小于1且β大于0.8如果估计结果β特别低说明数据或模型有问题第二蒙特卡洛模拟生成的价格路径均值应该跟长期均值μ大致一致偏差超过10%说明模拟的样本量不够或者模型参数有问题第三VaR序列的覆盖率检验要做一下——在回测区间内统计实际亏损超过VaR的天数占总天数的比例比例应该跟设定的置信水平5%接近。如果覆盖率只有1%说明模型高估了风险如果覆盖率达到10%说明模型低估了尾部风险两种情况都需要回头检查模型设定。还发现一个有意思的对比检验把GARCH模型计算出的VaR跟历史模拟法Historical Simulation的VaR放到同一张图上比较。历史模拟法虽然简单粗糙但它参考的是真实分布而不是模型假设。如果两种方法的VaR在大部分时间里走势一致说明GARCH的设定没大问题如果短时间内差距特别大往往是模型里的跳跃参数估计不准需要手动调整跳跃强度或跳跃幅度。6. 后续扩展方向与个人实操体会这套源码复现做完之后我最大的收获不只是会调Matlab的工具箱函数而是真正理解了电力市场价格形成机制对风险管理模型的约束。GARCH和跳跃扩散模型不是随便套上去就能用的每一步都要回到市场本身的物理特征和交易规则去验证假设是否成立。比如负电价的出现会导致常规的对数收益率计算失效市场价格上下限的制度调整会造成时序上的结构性断点这些都是在实际业务中才会遇到的问题只看教科书根本不会意识到。后续如果想继续深挖我觉得有几个方向值得拓展。一是把多元GARCH比如DCC-GARCH放进来同时建模多个省区或者多个交易品种之间的波动率联动效应这对做跨区域套利策略很有价值。二是把负荷预测和新能源出力预测作为外生变量引入价格模型用Matlab的Machine Learning Toolbox做特征工程让价格模型能够提前反应供需基本面的变化。三是把风险度量从单日的VaR/CVaR扩展到多日的路径依赖指标比如回撤持续时间Drawdown Duration这能更全面地评估极端行情下策略的脆弱性。在实际操作中还有一点体会很深Matlab代码写得好不好跟能不能跟其他系统对接有很大关系。交易团队最终落地的时候不会只在Matlab里看结果通常要把风险信号输出成行情数据、推送到底层交易系统。所以源码的输入输出接口设计要格外重视——数据文件口径统一、输出字段标准化、函数封装边界清晰这些工程化的细节决定了模型能不能从研究环境平滑迁移到生产环境。我在复现过程中刻意把核心逻辑封装成了一个独立的函数包输入输出全部用结构体传递后面接任何前端界面或者外部数据库都方便很多。
