蒙特卡洛模拟定价

蒙特卡洛模拟定价（精选7篇）

蒙特卡洛模拟定价 第1篇

关键词：住房抵押贷款支持证券,定价,蒙特卡洛模拟,提前偿付,违约

1 引言

MBS偿付资金主要来源于借款人的还款资金, 借款人提前偿付或违约都必然影响MBS的现金流, 从而影响MBS的价格。在过去的研究中, 不少学者对提前偿付和违约行为进行了研究, 根据提前偿付 (违约) 模型的不同可分为两类, 结构化模型和简约化模型。结构化模型, 又称理性期权定价模型, 利用期权理论来描述借款人的提前偿付和违约行为。Dunn等[1,2]是利用结构化模型对提前偿付 (违约) 进行研究的早期代表, 认为利率变化是提前偿付 (违约) 的主要原因。Stanton[3]在模型中引入交易成本概念, 解释了在期权理论下不是所有借款人都会同时提前偿付 (违约) 的问题。Leung等[4]将房价因素引入模型, 建立了双因素模型。近年来Longstaff等[5,6,7]更进一步对有序提前偿付进行了研究, 构建了多元决策模型。

结构化模型估计结果与历史数据存在矛盾, 为此部分学者提出简约式模型, 试图利用统计方法从历史数据中发现规律。Schwartz[8]等将比例危险函数用于MBS定价, 对后续相关研究产生了深远影响。由于比例危险模型已成为一种研究提前偿付和违约行为的常用方法, 本文假设提前偿付和违约行为符合比例危险模型。

MBS定价的常用方法主要包括蒙特卡洛模拟和有限差分法, 特别在多维随机过程情况下, 几乎所有定价模型都采用蒙特卡洛方法求解[9]。Chen[10]认为由于MBS现金流存在路径依赖的同时还隐含着一个美式看涨期权和一个美式看跌期权, 因此蒙特卡洛模拟成为唯一适合于MBS定价的方法。

本文目的在于构建包含提前偿付和违约因素的定价模型, 并运用蒙特卡洛方法进行模拟定价, 希望能对我国MBS的开展提供参考。本文结构安排如下:第二部分对MBS定价模式进行分析;第三部分是定价模型的构建, 分别对MBS的现金流、提前偿付和违约模型以及利率期限结构进行分析和描述;第四部分利用蒙特卡洛方法对MBS的价格进行模拟;第五部分是全文的研究结论。

2 MBS定价模式分析

一般来说, MBS的定价模式主要包括现金流现值定价、利率模拟定价和期权调整利差定价。

2.1 现金流现值定价

(1) 静态现金流量报酬法。

静态现金流量报酬率指MBS未来现金流现值等于当前市场价格时的贴现率。静态现金流量报酬率可由公式P=∑${}_{\text{t}=1}^{\text{Τ}}$CFt/ (1+r) t解出。其中P为MBS市场价格, CFt为MBS第t期现金流量, T为MBS持续期, r为静态现金流量报酬率。由上式可以看出, 定价成功与否很大程度上取决于CFt预测的准确性。静态现金流量报酬率法简单实用, 但无法反映利率期限结构, 同时也没有考虑不同利率路径下现金流量的波动性。

(2) 静态利差定价。

静态利差法对不同期限贴现率进行了区分, 以“可比的”国债到期收益曲线上不同期限折现率加一个固定利差作为贴现率。其公式为:P=∑${}_{\text{t}=1}^{\text{Τ}}$CFt/ (1+rt+s) t。其中, rt是可比零息国债即期利率, s为可调利差。

2.2 利率模拟定价

利率是影响MBS现金流的主要因素之一, 同时又对MBS贴现率产生影响。因此可通过给出利率期限结构, 来预测MBS未来现金流和贴现率, 从而求出MBS的理论价值。

2.3 期权调整利差定价

静态利差定价假定再融资利率在整个期间内不变, 设定了一定提前偿付的可能性, 由它预测出的各期现金流按利率变化路径上相应收益率进行贴现, 各期现金流的现值和即是MBS的价格。静态利差定价没有考虑通行利率可能存在的多种走势, 因此有学者对其进行改进提出期权调整利差法 (Option-Adjusted Spread, OAS法) 。OAS法考虑了通行利率的多种可能走势, 同时也考虑了由此引起的再融资利率变化和相应提前偿付的各种预测。

OAS法通常利用蒙特卡洛模拟产生足够数量的利率路径, 在各利率路径下, 确定MBS可能的现金流, 然后假定一个利差, 在每条路径上都使用该利差和该路径利率之和作为现金流贴现率, 将各路径下的现金流进行贴现可以得到各路径下MBS价格的一种可能情形。将上述计算出来的MBS价格加以平均得到的价格称为理论价格, 而使得理论价格与市场价格相等的利差即为期权调整利差。

3 定价模型构建

3.1 MBS现金流分析

本文仅考虑以固定利率每月等额偿付贷款本息的转付型MBS, 为研究方便, 不考虑服务费用。假设资产池初始加权平均贷款 (年) 利率为R0, 债券持续期为T。令Bt表示无提前偿付和违约前提下, 第t期期末未付本金余额, 其中B0为资产池初始本金, 此时MBS月现金流Flowt为:

t期末, 剩余本金为:

令It和Pt分别表示t期利息和本金的支付额, 则有:

进一步定义违约率λ${}_t^D$为t期违约本金占期初剩余本金百分比, 提前偿付率λ${}_t^{{\rm P}}$为t期提前偿付本金占期初剩余本金百分比, 则第t期期末实际剩余本金$\stackrel{}{{\displaystyle B}}{}_t$为:

t期实际现金流$\stackrel{}{{\displaystyle \text{F}\text{l}\text{o}\text{w}{}_{\text{t}}}}$为:

3.2 定价模型

假设1 (提前偿付与违约函数) :提前偿付率λ${}_t^{{\rm P}}$和违约率λ${}_t^D$符合比例危险模型, 即:

λ${}_t^l$=λ${}_0^1$exp[ (βl) TX${}_t^1$] l={P, D} (7)

其中λ${}_0^l$为基本危险函数, (X${}_t^l$) t≥0是「适应的多维随机变量, βl是变量系数。为研究方便进一步假设λl的影响因素为贷款利率 (R) , 贷款利率平方 (R2) , 贷款利率与初始贷款利率的利差 (S) , 贷款利率与初始贷款利率利差的平方 (S2) 等四个变量[11]。

假设2 (利率过程) :即期利率 (rt) 和贷款利率与即期利率的利差 (st) 满足如下的平方根均值回复模型 (CIR model) :

εt、ε't～iid N (0, 1)

其中△rt=rt-rt-1, △st=st-st-1, θ0、θ'0表示利率和利差调整速度, θ1、θ'1表示长期平均利率和利差, θ22rt-1、θ'22st-1表示利率和利差的波动率。由定义有Rt=rt+st, St=Rt-R0, (t>0) 。根据CIR模型的性质, 上述方法产生的即期利率和贷款利率在绝大多数情况下满足Rt>rt。

定理[12]:在假设1成立的前提下, 持续期为T的MBS在t时刻的均衡价格为

其中$\stackrel{}{{\displaystyle Flow}}{}_{t‚m}=E{}_m^{*}[△(t‚m)\stackrel{}{{\displaystyle \text{F}\text{l}\text{o}\text{w}}}{}_t]‚△(t‚m)=\exp (-{\displaystyle {\sum }_{j=t}^{m}r}{}_j/12)$, 这里rj为第j期的即期 (年) 利率, △ (t, m) 为从m期到t期的贴现利率, E*m () 是一种对 (rj) 的风险中性测度。

4 定价的蒙特卡洛模拟

4.1 定价思路

根据利率期限结构和利率波动性假说, 模拟足够数量的即期利率和中长期贷款利率路径, 据此确定MBS的提前偿付和违约情况, 从而进一步的确定未来现金流, 将现金流贴现后即可得到MBS一种可能的价格, 最后将各种可能价格求取期望即是所求的MBS价格。具体定价思路如下:

(1) 产生I个随机序列, 每个序列包含N个随机变量, 由假设2可得I个即期利率序列 (r${}_1^{(\text{i})}$, r${}_2^{(\text{i})}$, , r${}_{\text{Ν}}^{(\text{i})}$) 、贷款利率序列 (R${}_1^{(\text{i})}$, R${}_2^{(\text{i})}$, , R${}_{\text{Ν}}^{(\text{i})}$以及贷款利率与初始贷款利率的利差序列 (S${}_1^{(\text{i})}$, S${}_2^{(\text{i})}$, , S${}_{\text{Ν}}^{(\text{i})}$, 其中i={1, 2, , I}。对每个即期利率序列都存在一个贴现因子${}_{\text{m}‚\text{n}}^{(\text{i})}$=exp (-∑${}_{j=m+1}^n$r (i) jh) (h=1/12, n>m) (11)

这样, 基于I个序列我们估计的从第n期到第m期的贴现函数△ (m, n) 即为

(2) 由假设2估算提前偿付率λP (i) 和违约率λD (i) , 进而通过 (6) 式得到各期实际现金流$\stackrel{}{{\displaystyle \text{F}\text{l}\text{o}\text{w}}}{}_t^{(i)}$, 将未来现金流贴现汇总可以得到第i个路径对应的MBS价格。最后这I个MBS价格的期望就是所求MBS产品的最终价格。

4.2 参数设置

考虑一个30年期MBS, 资产池由按月等额偿付固定利率抵押贷款构成, 资产池初始贷款利率为5.5%, 规模为100万元, 即T=360, R0=0.055。式 (8) 中θ0=0.4, θ1=0.015, θ2=0.06, r0=0.015;式 (9) 中θ'0=0.5, θ'1=0.04, θ'2=0.05, s0=0.04;进一步的, 在提前偿付模型和违约模型中λ${}_0^{\text{Ρ}}$=0.0512 (1-exp (-0.3607 (t-1) /12) ) , λ${}_0^{\text{D}}$=0.0514 ( (t-1) /70) 0.8exp (- ( (t-1) /70) 1.8) , R, R2, S, S2的系数分别为1.8978, -0.11102, 0.04907, 0.00008和0.97596, -0.01569, 0.00794, 0.00003, 上述参数的选择主要参考了[11,12]。

4.3 结果分析

一般来说, 随着模拟路径增加, 模拟误差会越来越小。Davidson[13]的研究表明:模拟路径达到1000次左右时, 模拟误差已经趋于很小, 之后再增加模拟次数, 模拟误差不再有明显变化。因此本文选择的模拟次数为1000次。

图1反映了即期利率和贷款利率的一条模拟路径, 图2反映了与之对应的提前偿付率和违约率情况。在图3中分别显示了无提前偿付风险无违约风险 (图3a) 、有提前偿付风险无违约风险 (图3b) 、无提前偿付风险有违约风险 (图3c) 以及有提前偿付风险和违约风险 (图3d) 四种情况下实际现金流的模拟结果。结果表明:提前偿付改变了现金流结构, 资产池前期现金流增加, 后期现金流减少;在最初的几年违约对现金流的影响较明显, 后期由于资产池规模缩小以及违约率下降, 违约的影响也相应下降。图4反映了初始时刻的MBS价值分布 (含提前偿付和违约风险) , MBS的平均价格为139万, 最低价为138.2万, 最高价为139.5万, 模拟MBS价值的标准差为1855.03元。图5则反映了模拟MBS价格的收敛情况。

5 结论

MBS产品本身的性质决定了其现金流不仅具有随机性而且具有路径依赖性。传统的债券定价方法在处理路径依赖性方面存在不足, 因此本文在建立了包含提前偿付和违约风险的MBS定价模型后, 应用蒙特卡洛模拟方法对某一具体案例进行了定价分析, 反映了提前偿付和违约行为对MBS现金流的影响, 同时也说明该方法具有一定的实用性。限于篇幅, 本文并没有对提前偿付率和违约率模型及参数进行深入的研究, 仅是参考了部分国外的研究成果, 因此具有一定的局限性。在MBS的定价过程中, 构建符合我国国情的提前偿付模型和违约模型是精确定价的关键, 同时各种风险因素之间相互影响, 也是有待未来需要深入研究的一个问题。

蒙特卡洛模拟定价 第2篇

随着国际金融市场的不断发展和完善, 大量金融衍生产品应运而生, 其中就包括奇异期权, 即由标准期权经过组合、变化衍生出来的各种新产品。障碍期权是其中一种, 因其比普通标准欧式期权价格便宜, 倍受市场青睐, 被广泛用于资产风险管理。金融产品作为资本市场投资工具, 具有较高的风险收益特性, 如何正确评估产品的内在价值就显得尤为重要。而选择合理的定价模型则是正确评估产品价值的重中之重。本文主要研究了上升敲入欧式看涨障碍期权 (对于看跌期权, 可以通过期权评价公式求得) 的二叉树和蒙特卡洛模拟定价模型, 以比较二者之间的差异。

二、障碍期权的涵义及定价模型

2.1向上敲入障碍期权的涵义

障碍期权是一种路径依赖期权, 其收益依赖于标的资产的路径变化, 当期权在有效期内的资产价格触及到某一水平, 即通常所说的障碍时, 期权合约生效或失效。向上敲入障碍期权是指障碍水平高于标的资产初始价格, 且资产价格在规定时间内达到障碍水平, 期权才得以生效的金融衍生品。

2.2二叉树定价模型

假设市场是完全市场, 不存在套利行为和交易成本, 所有资产可以无限细分, 没有卖空限制, 且股票不支付红利。设股票初始时刻价格为S0, 期权初始价值为f, 无风险收益率为常数r, 期权执行价格为K, 障碍值为SB (SB>S0) , 期权的有效区间为[0, T], 将其细分为N个等长的区间[tn, tn+1] (0≤n≤N-1) , 在每个区间上股票价格或者以概率p上升为原来的u倍, 或者以概率1-p下降为原来的1/u (d=1/u) , 相应的, 期权价值变为fu或fd。

根据风险中性原理, 可以求得模型的参数p, u, d的表达式:

其中, σ是股票价格的实际波动率, △t=T/N。

由向上敲入欧式看涨障碍期权的定义知, 时刻t的期权价值为

其中, Iω (S) 是集合ω的特征函数, 可简写为Iω, 即

在风险中性世界里, 知在T-△t时刻的期权价值可由t时刻期权价值的期望值在△t时间内用无风险收益率r贴现得到, 即

同理, T-2△t时刻的期权价值可由T-△t时刻期权价值的期望值在△t时间内用无风险收益率r贴现得到, 以此类推, 可以得到在初始时刻的期权价值。

但是存在一个问题, 即设定的障碍价格水平不一定正好落在二叉树结点上, 即障碍价格 (真实障碍) 水平处于两个结点水平之间。虽然这种情况可以通过将时间步长缩短, 细化时间间隔的方法减少误差, 但是会大大增加计算复杂度。本文采用线性插值的方法进行处理。首先选定最为临近障碍价格水平的两个结点水平, 分别作为内障碍水平和外障碍水平, 然后利用线性插值公式求得 (图1) 。

2.3蒙特卡洛模拟定价模型

蒙特卡洛模拟方法的思路很简单, 首先是随机生成风险中性条件下资产价格的多条不同路径, 然后计算所有路径期权价值的期望值, 最后用无风险利率贴现即可得到所求期权价格。

障碍期权是路径依赖期权, 其收益是到期日前所有股票价格的函数, 即

, 假设股票价格遵循几何布朗运动, 由伊藤引理可得股票价格的随机路径:

其中, ε是由标准正态分布产生的随机数, 可由Box-Muller法或大数法则产生。因此, 只要随机生成一个ε, 就可以得到股票价格的一条路径。

设第i条路径的障碍期权收益为fi, 则

计算k条路径期权收益的均值作为期权价值的期望值, 即

在风险中性条件下, 用无风险利率r贴现得到障碍期权的估计价格, 即

三、实例与结论

3.1实例

选取的数据来源于姜礼尚所编的《期权定价的数学模型和方法 (第二版) 》, 设一张6个月的看涨期权股票, 现股价为145元, 无风险利率r=6%, 波动率σ=0.295。期权的执行价格为145元, 障碍价格分别为160元和190元。

通过matlab编程实现了上述实例中的向上敲入障碍期权的二叉树和蒙特卡罗模拟定价过程。表1给出了两种模型计算的结果和实际值。

3.2结论

表1中障碍价格为0时的期权价格, 表示标准期权的价格为12.87元。从中可以看出, 无论采用二叉树定价法还是蒙特卡洛模拟定价法, 都验证了障碍期权的价值要比普通期权的价值要低这一事实。

随着模拟路径的增加, 蒙特卡洛模拟方法得到的障碍期权价格与实际越接近, 这是因为样本数据增加后, 会减少估计值与实际值的误差水平, 使估计值更加接近实际值。

蒙特卡罗模拟比二叉树方法得出的结果更接近于实际价格, 更具有良好的实用性。其主要原因: (1) 二叉树定价模型中, 由于时间步长不够长, 由内外障碍价格计算出来的期权价格不能很好的精确到真实障碍水平下的障碍价格; (2) 随着障碍的增加, 障碍期权的价格偏离实际价格越大, 其主要原因是随着障碍价格的增加, 内外障碍价格差增大, 导致偏离真实障碍越大; (3) 二叉树定价时, 仅仅将到期时股票价格低于障碍水平时的期权价格赋予零值, 是存在缺陷的。因为二叉树定价法的原理是从后往前计算, 这与股票价格从前向后变动的方向是相反, 股票价格可能在到期之前已经触及到障碍价格, 而在到期时又下降到障碍价格以下, 在此情况下, 该条路径仍然是有效的, 其值并不为零。

摘要：基于对障碍期权的涵义理解以及对二叉树和蒙特卡洛模拟定价模型步骤的详细阐述, 采用两种方法分别对上升敲入看涨障碍期权进行了价格评估, 并通过实例分析两种方法的差异性。

关键词：障碍期权,二叉树定价,蒙特卡洛模拟定价

参考文献

[1]郑振龙, 陈蓉.金融工程 (第二版) [M].北京:高等教育出版社, 2008.141-233.

[2]姜礼尚.期权定价的数学模型和方法 (第二版) [M].北京:高等教育出版社, 2003.255-284.

[3]徐腾飞, 曹小龙, 胡云姣.离散障碍期权定价的蒙特卡罗模拟[J].北京化工大学学报 (自然科学版) , 2013, (3) :123-127.

[4]邢晓芳, 周圣武, 石广平等.基于跳扩散模型的欧式上升敲入期权定价研究[J].云南大学学报 (自然科学版) , 2009, (s2) :354-357.

设计投标中的蒙特卡洛模拟报价研究 第3篇

自我国改革开放以来, 勘察设计单位实行了技术经济责任制, 与建设单位 (业主) 之间逐步形成供需关系和价值的交换, 以一种行业出现在经济社会中。从1996年的建筑法到2000年国务院颁布的《建设工程勘察设计管理条例》、2003 年国家发改委等八部委联合发布的《工程建设项目勘察设计招标投标办法》, 勘察设计市场逐步走向成熟。根据《招标投标法》的规定, 对工程建设项目的勘察、设计、施工、监理以及与工程建设有关的重要设备、材料等的采购, 必须进行招标。2000年5月1日国家发展计划委员会发布了第3号令《工程建设项目招标范围和规模标准规定》, 规定了工程项目必须进行招标的范围, 其中勘察设计单项合同估算价在50万元人民币以上和项目总投资额在3000万元人民币以上的项目, 都必须进行勘察设计招标。自此, 设计投标已经成为设计咨询企业获得订单的主要途径。另外, 随着我国加入WTO满5年的到来, 国外设计咨询企业越来越多的进入中国市场, 设计咨询行业的市场竞争也越来越激烈。因此, 工程勘察设计企业只有参加投标竞争, 依靠其经济实力、技术优势、社会声誉等条件, 通过竞争占领市场、立足市场, 并不断地在市场中完善自己, 才能在激烈的市场竞争中长久地生存与发展。为此, 研究工程项目勘察设计的投标策略, 用以指导设计投标活动, 最大限度地提高设计投标中标率, 是设计咨询企业的一件大事。

勘察设计是服务性工作, 其招投标与施工招标的条件和要求都有所不同, 它不能象施工招投标那样以投标报价作为主要的定标因素, 这是由勘察设计行业特点决定的。勘察设计招标一般采用综合评估法评标, 评标委员会按照招标文件规定的评标标准和方法, 结合前一阶段设计批复文件, 对投标人的业绩、信誉和勘察设计人员的能力以及勘察设计方案的优劣进行综合评定。一般招标文件都要求投标文件分为技术标和商务标两部分, 商务标又包括企业资质、工程业绩、投入人力、获奖情况等资信条件和投标报价。归纳起来, 参与投标的设计咨询企业能否中标, 不考虑人际关系、攻关能力等方面, 决定因素主要包括技术文件、资信条件和设计费报价。技术和资信对于实力较强的设计咨询企业来说都不成为问题, 通过正确的决策、合理地组织与控制、精心设计和准备, 可以拿到较高得分, 而且对于实力相当的竞争对手来说, 在这两方面也不会有太大差距。但在报价方面, 由于一般招标文件都规定投标报价的满分与各家报价的平均值有关, 使得投标报价得分情况产生不确定性, 报价也就成为能否中标的重要影响因素。因此, 研究设计费报价问题, 是在现有技术水平和资信条件下提高中标率的关键。解决不确定性问题, 采用蒙特卡洛模拟技术是一种有效办法, 本文主要研究了设计投标中的模拟报价方法。

1 模拟报价方法研究

1.1 蒙特卡洛模拟法

模拟法即蒙特卡洛模拟 (Monte Carlo Simulation) [1], 其名称源于美国在第二次世界大战中研制原子弹的“曼哈顿计划”, 是从不同变量的分布中随机抽样进行模拟试验, 以求解各类技术问题近似解的数学方法。蒙特卡洛模拟法借助人对未来事件的主观概率估计及计算机随机模拟, 在解决难以用数学分析方法求解的动态系统复杂问题上具有极大的优越性。

用模拟法进行投标报价决策, 首先应构建报价得分的数学计算模型, 然后分析模型中每一随机变量可能变化的范围并确定这些变量的概率分布;通过模拟试验, 为各随机变量抽取随机数, 并将随机数按照概率分布模型转化为变量的抽样值;将抽样值组成一组基础数据, 计算出报价得分值;最后进行重复试验, 统计试验结果即可求出特定报价得分达到一定水平的概率, 或在给定得分下限值及报价得分高于下限值所要求的概率时通过试算求得报价区间。

1.2 构建报价得分模型

当招标文件规定以给定的设计费额度为满分报价或规定下浮率报价时, 报价决策不属于不确定性问题, 报价工作相对比较简单;但勘察设计招标一般采用综合评分法评标, 招标文件往往规定投标报价满分与各家报价的平均值有关, 所报价格比满分价格或高或低均有扣分, 因此报价得分模型是由评标办法决定的, 一般性地表示如下:

$F=\{\begin{array}{l}F{}_0-\frac{|x{}_1-g(x)|}{g(x)}100dF{}_n>F{}_{min}\\ F{}_{min}F{}_{n}F{}_{min}\end{array}(1)$

式中:

F己方报价得分

x1己方报价

F0招标文件规定的报价满分分值

Fmin招标文件规定的投标报价保底得分

d报价每高于 (或低于) 满分基准价1%的扣分值 (高或低时可以不同)

g (x) 满分对应的报价基准价。表达式由评分办法确定, 一般与各家报价有关, 例如按“各家报价的平均值与报价最低值取平均”时, $g(x)=\frac{1}{2}\left[\frac{1}{{\rm N}}{\displaystyle \sum _{\text{i}=1}^{\text{Ν}}x}{}_i+min(x{}_i)\right]$;又如按“各家报价扣除最高值和最低值后的平均值”, $g(x)=\frac{{\displaystyle \sum _{\text{i}=1}^{\text{Ν}}}x{}_i-max(x{}_i)-min(x{}_i)}{{\rm N}-2}$。式中N为投标人数量, xi (i=1, 2N) 为第i家投标人报价, 其中第一家投标人为己方。

编制投标文件进行报价决策时, 竞争对手如何报价是不确定的, 只能通过估计预测其可能的报价区间, 因此在上述得分模型中, 以g (x) 中各家报价xi为随机变量。

1.3 各家报价概率分布的确定

随机变量概率分布的确定是进行模拟分析的基础, 也是决定模拟准确性的关键。取值范围为一个区间的连续变量, 工程上常用的概率分布有正态分布、三角形分布、β分布、阶梯分布、梯形分布、直线分布等, 各有其特点和适用范围。投标报价作为一种特定的决策过程, 与评分办法、投标人的偏好、对招标人的倾向性预计、对竞争对手的报价预测等有关, 多数情况下趋向于一种不对称的偏态分布。β分布不仅具有较好的适应性, 它由参数决定分布的形状, 只要选择适当的形状参数, 就可以描述不同的分布;而且可以通过三点估计法确定模型参数, 便于实际应用, 是一种常用的偏态分布模型, 适合描述投标报价的概率分布。β分布由变量的最小值a、最大值b和形状参数r、s确定, 概率密度函数为:

$f(x)=\frac{\Gamma (r+s)}{(b-a){}^{r+s+1}\Gamma (r)\Gamma (s)}(x-a){}^{r-1}(b-x){}^{s-1}(axb,r＞0,s＞0)(2)$

其均值和方差分别为:

$v=a+(b-a)\frac{r}{r+s}$ (3)

$\delta {}^2=(b-a){}^2\frac{r\cdot s}{(r+s){}^2(r+s+1)}$ (4)

变量的分布参数通常采用历史数据推定法或专家调查法 (常用德尔菲法) 确定。一般情况下, 对某一竞争对手如何报价是难于估计的, 但根据项目特点和评分办法, 估计全部竞标人报价的整体分布趋势相对较易, 从而确定本次投标活动中各家报价服从的概率分布。采取专家调查法, 由多名熟悉投标活动和市场情况的经营人员, 估计本次报价得分最高的最可能报价值 (m) 、最小报价值 (a) 和最大报价值 (b) , 按照三点估计法确定报价服从的概率分布模型。根据经验公式[2], 三点估计法的β分布均值和方差的估算值分别为:$v=\frac{a+4m+b}{6}$;$\delta {}^2=\frac{(b-a){}^2}{36}$, 分别代入式 (3) 和式 (4) 连立求解, 即可得出形状参数r、s的值。由于Γ函数较为复杂, 手工计算很困难, 但通过计算机程序能够快速准确的得出β分布函数模型。另外, 对于经常同台竞争、熟悉对方报价策略的竞争对手, 还可以通过收集对方以往投标项目报价资料, 采用历史数据推定法或专家调查方法来估计其报价的最可能值 (m) 、最小值 (a) 和最大值 (b) , 预测每一竞争对手的报价概率分布。

1.4 随机数的抽取

常用的随机变量抽样方法有逆变换法、舍选法、近似法等。β分布函数模型较为复杂, 一般应用中常采用舍选法, 通过两个伪随机数计算特定的函数值并经反复比较验算, 最终得到满足要求的随机变量值。本课题采用计算机编程计算, 随机变量抽样参照逆变换方式通过计算程序进行。

随机抽样中, 首先要产生 (0, 1) 区间的均匀分布随机数作为伪随机数。产生伪随机数的方法很多, 本课题采用编程语言提供的随机数函数获取。然后进行随机变量抽样值的转换, 即以获取的伪随机数R作为变量的概率值, 则有p{xX}=∫${}_{\text{a}}^{\text{X}}$ (x) dx=R, 根据其概率密度函数通过计算机求解变量的抽样值X, 使伪随机数转化为随机变量抽样值。如图1所示。

1.5 模拟过程

为了保证蒙特卡洛模拟结果的准确性, 需要对相同的随机试验过程重复进行成千上万次, 然后整理全部结果进行统计计算。

(1) 以全部投标人的报价为随机变量, 模拟计算满分报价基准值的平均值, 以此作为己方初始报价的参考。即每次按照报价概率分布模型, 随机抽取每一投标方的报价值, 按照评分办法确定的报价满分基准价g (x) 函数计算满分对应的报价基准价, 重复进行若干次模拟后, 计算满分基准价的平均值。因为报价越接近满分基准价得分越高, 所以用模拟出的满分基准价平均值作为己方初始报价的参考值, 报价是容易得到高分的。

(2) 给定己方报价和期望下限得分, 以全部竞争对手报价为变量, 模拟计算己方得分大于期望下限得分的概率。设定己方报价方案后, 随机抽取每一竞标对手的报价值, 按照报价得分模型计算己方报价得分;重复进行若干次模拟后, 对得分值排序整理, 可以计算出得分大于期望下限得分的概率值。概率表达的是得高分可能性的量化值, 以百分数表示的概率值越大则得高分的可能性越高, 以此来判断己方报价方案的可行性。

(3) 设定己方得分下限值和报价得分高于下限的概率 (即置信水平) , 反向模拟试算己方报价区间。前面对特定报价模拟计算得高分的概率, 用来判断己方具体报价是否可行;但如果需要提出可行的报价建议, 则要通过反向模拟进行区间估计。这种反向模拟是在报价所有可能取值的范围内按照一定步距逐一试算, 每一次试算都是经过若干次模拟得出得分高于下限的概率并与设定概率值比较, 从而得出满足条件的报价临界值, 形成推荐的报价区间。

2 模拟报价应用实例

2.1 项目概况及报价模型

某设计招标项目共有6家设计单位参与投标, 我方投资测算结果建安工程费约50000万元, 按照国家设计费收费标准计算出的设计费为1477万元。投标报价满分为20分, 以“各家报价平均值与次低值取平均”为满分基准价, 报价低于基准价1%扣0.30分, 高于基准价1%扣0.50分, 直至扣至保底分5分。以xmin2表示6家报价中的次低价, 满分对应的报价基准价为:

$g(x)=\frac{1}{2}(\frac{1}{6}$${\displaystyle \sum _{\text{i}=1}^6}$xi+xmin2)

报价得分模型可以表示为:

$F=\{\begin{array}{l}20-\frac{|x{}_1-g(x)|}{g(x)}100d(F{}_n＞5)\\ 5(F{}_{n}5)\end{array}$

其中x1>g (x) 时d=0.5, 否则d=0.3。

2.2 确定报价概率分布

经4名主管领导和经营人员进行专家估计, 得到各家投标报价服从最可能值1060万元、最小值980万元、最大值1270万元的β分布。经计算得出形状参数r=2.5228, s=4.6734, 概率密度函数为:

$\begin{array}{l}f(x)=\frac{\Gamma (7.1962)}{290{}^{8.1962}\Gamma (2.5228)\Gamma (4.6734)}\\ \leftarrow (x-980){}^{1.5228}(1270-x){}^{3.6734}\end{array}$

2.3 模拟结果

以所有投标方报价 (含我方) 为随机变量, 进行5000次模拟计算满分报价基准值的平均值, 结果为1065.19万元, 此平均值可以作为我方报价的参考值。

以1100万元为我方报价、其他各家报价为随机变量, 进行5000次模拟计算, 得出我方得分均值为18.52分, 得分大于18分 (设定的分数下限值) 的概率为70.70%。以1080万元为我方报价时, 我方得分均值为19.18分, 得分大于18分的概率为94.44%, 相应累计概率图见图2。

设定我方得分下限值为18分、得分高于18分的概率不小于90%, 反向模拟我方报价区间, 结果为1005万元～1085万元。报价对应得分均值散点图见图3, 报价对应得分大于18分的概率散点图见图4。也就是说, 在竞争对手报价分布模型专家估计的基础上, 经过报价模拟分析, 建议我方报价在1005万元～1085万元之间取值, 可使报价得分高于18分的概率达到90%以上。

经领导决策, 我方在该项目上报价1080万元。开标后我方报价得分为最高, 但因技术标失分及其他因素, 最终未能中标。

3 结语

采用蒙特卡洛模拟技术, 以各家报价为随机变量, 通过估计竞标人报价的概率分布, 可以模拟计算出特定报价得高分的概率值, 或在一定置信水平下模拟出得高分的报价区间, 从而为投标报价决策提供依据, 以提高设计投标的中标率。当然, 投标报价模拟过程必须借助计算机程序才能进行, 而且模拟结果的可靠性依赖于报价分布模型专家估计的有效性。对于熟悉的竞争对手, 还可以分别估计每一对手报价的概率分布, 从而使模拟结果的可靠性更高。

参考文献

[1].徐钟济.蒙特卡罗方法[M].上海:上海科学技术出版社, 1985

蒙特卡洛模拟定价 第4篇

在计量经济学中, 古典线性回归模型的基本假定是对参数和模型进行显著性检验的基础, 如果违背了其中任意一条或多条假定, 假设检验的结果将不再有效, 甚至可能会得到完全错误的结论。

许多计量经济学教科书对此给出了严谨的数学证明, 这里不再赘述。笔者试图从另一个角度, 即运用数据生成过程, 以参数的显著性检验为研究对象, 探讨存在异方差和多重共线性时, 检验结果是否有效。

第二章, 简单地介绍一下古典线性回归模型的基本假定中的“同方差”假定和“无多重共线性”假定, 为后续的论述打下基础。第三章, 介绍数据生成过程, 利用一个简单的二元回归模型, 生成所有需要的变量, 用来进行回归和检验。第四章分别讨论不同条件下的检验结果和结论。

二、同方差假定和无多重共线性假定

在进行数据生成之前, 先简单地介绍一下古典线性回归模型基本假定中的“同方差”假定和“无多重共线性”假定。以一个二元回归模型为例,

其中, Y为被解释变量, X1和X2为解释变量, β是参数, 随机误差项εi通常被假设为服从均值为0, 方差为常数σ2的正态分布。

可以看出, 由于假定方差σ2为常数, 随机误差项是“同方差的”, 每个Y是以相同的方差分布在其均值周围, 这就是古典线性模型中的“同方差”假定, 如果该假定被违背, 则称“存在异方差”。

另外, 古典线性回归模型还假定, 解释变量X2与解释变量X1之间没有较强的线性关系, 即“无多重共线性”假定。如果该假定被违背, 则称“存在多重共线性”。

三、数据生成过程

第一, 生成解释变量:

X1i=1, 2, 3200, 其中i=1, 2, 3200。

可以看出, 变量X1实际上是一组简单的等差数列。

第二, 生成解释变量:

由于变量X2是由变量X1生成的, 因此, 如果随机误差项ηi的方差θ越小, X2与X1的线性关系就越强, 即存在多重共线性。

第三, 生成被解释变量Y:

任意设定以下一个二元回归模型,

将前两步中生成的解释变量X1和X2放入模型中, 即可生成被解释变量Y。

这里, 随机误差项εi的方差被设置为xα1i, 只有当α=0时, 该模型才不存在异方差, 否则就会存在异方差。

第四, 将生成的解释变量Y对解释变量X1和X2进行二元回归, 采用普通最小二乘法估计每一个参数, 并依次对各个参数进行显著性检验, 检验其在5%显著水平下是否显著。

第五, 重复以上步骤1~步骤4, 为了增加研究结论的准确性和可靠性, 重复次数设置为3000。

第六, 统计在这3000次回归分析中, 参数的显著性检验结果显著的比例。

在整个数据生成过程中, 存在异方差与否, 取决于选择的参数α是否为0, 例如当α=0时, 不存在异方差, 而当α=1时, 则存在异方差。

存在多重共线性与否, 取决于选择的参数θ是否足够大, 例如当θ=5, 不存在多重共线性, 而当θ=0.01时, 则存在多重共线性。

四、结论

由于被解释变量Y是由两个解释变量X1和X2共同生成, 也就是说, X1和X2对Y均有显著性影响, 参数的显著性检验结果应该为显著, 这是“真实情况”。

下面依次来看, 在不同条件下, 检验结果是否与“真实情况”一致, 据此推断:在不同条件下, 参数的显著性检验是否依然有效。

(一) 不存在异方差和多重共线性

为了方便比较, 首先考察古典线性回归模型基本假定成立的情形, 即, 不存在异方差也不存在多重共线性, 令α=0, θ=5。

结果见下表:

由此可以看出, 当不存在异方差和多重共线性时, 参数的显著性检验结果全部为显著, 这与“真实情况”完全吻合, 说明此时参数的显著性检验完全有效。

(二) 存在异方差

结果见下表:

由此可以看出, 当存在异方差时, 参数的显著性检验结果显著的比例非常低, 这与“真实情况”不符, 说明此时参数的显著性检验失效。

(三) 存在多重共线性

结果见下表:

由此可以看出, 当存在多重共线性时, 参数的显著性检验结果显著的比例非常低, 这与“真实情况”不符, 说明此时参数的显著性检验失效。值得注意的是, 对常数项进行显著性检验的结果似乎依然有效, 这是由于X2的生成函数中缺少截距项造成的。

综上所述, 当存在异方差或多重共线性时, 参数的显著性检验结果显著的比例大幅降低, 说明在此条件下, 参数的显著性检验失效。

摘要：通过蒙特卡洛实验, 在计量经济学软件R中实现数据生成过程, 当古典线性回归模型基本假定中的同方差假定和无多重共线性假定被违背时, 参数的显著性检验是否依然有效。

关键词：线性回归,异方差,多重共线性,参数的显著性检验

参考文献

[1]Gujarati and Porter, Basic Econometrics, 5th ed., McGraw-Hill, New York, 2009.

蒙特卡洛模拟定价 第5篇

作为全球最大的金融市场之一, 外汇市场吸引了无数投资者。现代金融市场的发展以及信息技术的兴起, 越来越多的金融分析人员依靠分析金融资产的价格趋势来建立投资决策模型。交易策略在本世纪普遍盛行。根据Aite Group的统计预测, 外汇市场上未来将超过30%的交易量是由算法完成

技术分析, 主要从市场历史数据挖掘信息, 来找到盈利的模式和趋势, 所以也被称为“图表法”或者称为技术交易策略。技术分析最主要是使用技术指标同时结合价格变动图表来预测未来价格变动趋势。此预测方法最早是Charles Dow在19世纪末期发明使用的, 曾经被投资者广泛使用, 作为投资交易决策的指导。近几年, 许多经济学家通过对外市场上使用技术分析的盈利能力文献进行梳理, 得出技术分析手段在近些年的外汇市场上的收益逐渐消失。其中大多数的解释都认为外汇市场变得更加有效。最新的金融理论适应性市场假说 (AMH) (罗闻全, 2004) 从另一个角度也解释了这一现象, AMH理论认为个人投资者依据以往经验来进行投资决策, 他们通过学习以往决策结果来强化之后的投资。因此, 以往能获得超额收益的技术策略由于广泛地被运用到市场投资中, 而不再有效。同时AMH在实际市场操作也有指导作用:首先, 适应性假说认为风险和收益之间的关系, 在时间序列上是不稳定的, 即在一定时期内, 投资的风险和收益可能并不对称。其次, 适应性市场假说认为市场的套利机会总是存在的。第三, 适应性进化假说认为套利机会存在, 并至少有一种策略在一定时期是有效的, 只不过连续使用一种策略会使收益下降, 单一的投资策略难以应付市场环境的变化。最后, 适应性市场假说认为创新是在市场中生存的核心。随着进化算法的发展, 这种方法克服了以往单个技术指标的弊端, 通过优化各种技术指标组合, 被大量运用于外汇市场。

本文从该理论结果出发, 利用蒙特卡洛模拟方法选择最优的外汇交易策略, 从而不断适应市场环境变化调整交易策略, 以获得超额的收益。

2. 建模方法

2.1 神经网络交易模型

神经网络模型是一种半参数方法。Cheng等 (1994) 从统计的角度解释了神经网络提供的信息。图1为一个包含隐含层简单的前馈网络结构。输入层有2个节点, 隐含层包含3个节点。输入层每个节点与隐含层的节点相连接, 同时这些隐含层的节点与输出层的单个节点相连。我们称这个网络为一个2-3-1前馈神经网络。

神经网络通过“活动函数”把信息从一个层传递给下一层。考虑包含一个隐含层的前馈神经网络, 隐含层的第j个节点定义为:

其中xi为第i个输入节点的值, fj (·) 为激活函数, 通常取logistic函数:

αoj为偏差, 从i→j求和是指所有输入节点反馈给j, wij表示权重。

应用神经网络包括两个步骤。第一步是训练网络 (即:建立一个网络, 其中包括确定节点的数量以及估计它们的偏差和权重) 。第二步是推断, 特别是预测。在训练阶段, 数据通常分为两个不重叠的子样本。第一子样本是用来估计给定前馈神经网络的参数。建造完成的网络再用于第二子样本进行预测, 并计算其预测的准确性。通过比较预测性能, 选择优于其他的网络, 作为“最好”的网络来进行推理。这是广泛应用于统计模型的选择一种交叉验证 (cross validation) 的想法。当然, 其他的模型选择方法也可使用。

在时间序列的应用中, 假设{ (yt, xt) |t=1, ···, T}为网络训练样本数据, 其中xt为输入向量, yt为因变量 (如资产收益率等) 。对于一个给定的网络结构, 令ot表示输入为xt时网络的输出结果。通过最小化一些拟合标准选择参数偏差和权重, 来训练神经网络, 常用的标准如最小平方法:

这是一个非线性的估计问题, 可以通过一些迭代方法解决。为了确保拟合函数的平滑性, 现有的最小化问题可以添加一些额外的约束。在神经网络文献中, 反向传播 (BP) 学习算法也是网络的常用训练方法之一。

2.2 支持向量机模型

支持向量机 (SVMs) 和神经网络一样, 可以用于回归和分类任务的建模。支持向量机模型以强大的理论基础和不同领域的成功应用获得了越来越多的关注。支持向量机是一般线性分类器的一种, 最主要的特点就是能够同时最小化经验误差和最大化几何边缘区。Vapnik (1995, 1998) 和Shawe-Taylor, Cristianini (2000) 对SVMs模型基础做了详细的论述。Smola和Scholkopf (2004, 1998) 对SVMs用于回归应用提出了一些基本的想法。

SVMs模型背后最基本的想法就是将原始数据映射到一个新的高维空间上, 从而可以使用线性模型得到一个分割平面, 对分类问题进行处理。而将原始数据映射到新的空间上主要是通过核函数 (kernel function) 来实现。所以从这个角度上, SVMs属于线性分类器的一种。

超平面的选择主要是通过最大化不同类别样本间分离边际来得到, 如下图所示 (2) 。这个最优化问题主要通过二次规划来解决。软边际方法允许小比例样本分到错误边际一侧, 这些就被称为“花费” (cost) 。

支持向量回归过程也类似, 最主要的区别在于误差的形式和花费 (cost) 的计算。这通常采用ε-insensitive损失函数

2.3 多元适应性回归样条

多元适应性回归样条 (Friedman, 1991) 是可加线性模型 (Hastie和Tibshirani, 1990) 的特例。MARS的表达形式如下:

其中, cis为常数, Bis为基础函数。

基础函数可以取多种不同的形式, 从最简单的常数到多变量之间交互作用函数。最常用的基础函数是hinge函数, 其形式为:

其中, xi为预测变量, t为变量的阈值。下图 (3) 为一个hinge函数例子:

3. 构建模型及交易策略

3.1 定义预测模型

细化目标是通过构建的模型预测外汇价格的未来值, 同时设计交易系统好用, 在该下单时下单, 该平仓时平仓, 以获得最大利润。所以, 最重要的就是要能通过模型准确预测出外汇价格未来值。其中需要面对的问题有: (1) 是否需要预测每日报价; (2) 预测未来多长时间。这些通常取决于使用的模型如何预测, 并如何用于交易决策。

假设, 如果价格相差超过p%, 则交易是有价值的 (考虑交易成本后) 。我们通过上面讨论的三种建模方法来预测未来k天的价格是否达到这个边界值。同时, 在未来的k天的总体价格动态预测, 这不是在某个时间由一个特定的价格获得。如:在时刻t+k收盘价表示一个变化比p%低得多的值, 但它可能在窗宽t...t+k内已被前面有一个时期代表价格的变化比p%要高得多。所以, 重点是要预测未来k天价格总体趋势。引入一个新的变量, 它是在报价数据的基础上计算的, 可被视为在未来k天价格趋势的指标, 用字母T表示。同时, 这个指标的值与设定未来k天的边界p%相关。p%的变化, 是指价格高于或低于目前的价格。基本思想是:正向的变化将导致买入, 而负向的变化将导致卖出行为。指标T值假设趋势作为一个单一值, 正值代表向上的趋势, 负值代表价格下调的趋势。定义每天的平均价格近似为:

其中, Ci, Hi和Li, 分别表示外汇的第i天的收盘价, 最高价, 最低价。

令Vi表示未来k天平均价格对当天收盘价的k百分比变动 (也成为算数收益率) :

定义的指标变量T为所有变动超过目标边际p%的变化率之和:

变量T的含义为:在k天的周期内所有日平均价格高于目标变化率的总和。高的正T值意味着有多个未来日均价格高于今天的收盘价p%。这种情况是发出买入指令的潜在机会, 因为有良好的预期外汇价格将上涨。另一方面, 高的负T值则建议卖出操作, 预计未来外汇的价格很可能会下降。如果T值在零附近, 可能的原因为期间内价格比较“平坦”或者是由于正和负的变化相互抵消所引起, 即价格处于上下波动状态。

通常, 技术指标方法有两种, 每一个技术指标都是从一个特定的方面对外汇市场进行观察。第一种是按照明确的数学公式产生新的数字, 这是技术指标中较为广泛的一类, 著名的KDJ指标、RSI指标、MA指标都属于这类。第二种是没有明确的数学公式, 只有处理数据的文字叙述方法, 这一类指标相对较少。技术指标的原始数据是开盘价、最高价、最低价、收盘价、成交量、成交金额和成交笔数等。由于外汇数据的对交易量的衡量不准确, 所以与交易量相关的交易指标本文没有涉及。本文采用的技术指标如下:ATR (Average True Range) 指标表示序列的波动性;SMI (Stochastic Momentum Index) 指标表示动量效应;SMI (Stochastic Momentum Index) 指标表示随机动量效应;ADX (Directional Movement Index) 表示方法变动指标;Aroon指标识别开始趋势;Bollinger带来比较一段时期内的波动性;CV波动率指标;CLV (Close Location Value) 表示收盘价在交易区间相对值;MACD指标;SAR (Parabolic Stop-and-Reverse) 表示反转过程;Volatility指标表示波动率, 共11个价格指标。

本文使用三种方法构建指标T值与上述11个指标变量关系模型, 从而得到对T值得预测模型。通过指标T值来描述未来价格的变化趋势。相应地, 构建的交易规则为:在t时刻上, 如果T值超过某个阈值, 则正确的交易信号是“买入”;如果T值低于另一个阈值, 则交易信号为“卖出”;其他情况下, 交易信号是“持有”。T值处理有两种方法, 一种为直接通过技术指标来预测T值, 即T值为连续变量, 所以是一个回归问题, 当然最后还是要根据预测的T值和设定的阈值来转化成外汇交易信号;另一种方法则是通过计算历史数据相应的T值, 依据阈值方法, 记录每一天正确的交易信号, 这种方法将问题转化为分类问题, 即通过技术指标来直接预测买卖信号。具体说:

第一个选择是使用T值作为目标变量, 并尝试使用预测信息构建模型以预测T值。这种方法与多元回归相似。即将模型预测结果“转化”为交易信号, 通过决定阈值的大小, 根据预测的T值, 将交易信号分为3种可能。如下面的公式所示:

阈值0.01和-0.01的选择是纯粹启发式的, 也可以使用其他的阈值。但0.01值的含义也表示在k取10天期间内计算出来的T值中有至少四个日均价格是超过目前的价格0.25%, (4×0.0025=0.01) 。如果使用其他阈值, 过高的阈值将减少信号产生, 而太小的阈值则会导致频繁的交易, 市场微小的变化都有可能一起交易, 从而会导致风险加剧。

第二种方法是直接预测信号, 将问题转化为处理一个分类问题。使用指标T值和上述公式中相同阈值判断交易信号, 即使用历史数据计算k去10天周期的T值, 再通过上述公式和阈值将T值转化成每一天的交易信号。对交易信号 (买, 卖, 持有) 进行建模, 这就类似于一个分类问题。回归和分类的主要区别在于目标变量的类型上, 如果目标变量时数值型的 (如T值) , 就可以采用回归方法来研究;如果目标变量时分类变量 (如买卖信号) , 就可以使用分类方法来处理。不同的方法和技术分别应用与这两类问题。

3.2 模型评价指标

信号预测可以通过错误率的度量来评估, 其定义为:

其中ŷi是模型在测试样本i的预测值, 同时, 真实值用yi, L0 1被称为0/1损失函数的预测:

通常使用精度来度量精确性, 其中精度= (1-错误率) 得到。这两个统计量用来比较模型之间预测未来k天的准确性。

精度 (或错误率) 对于这类问题的度量并不是一个很好的选择。事实上, 市场上三种可能的结果 (买, 持有, 卖) 之间具有很强的不平衡性, 持有信号出现的可能性远远大于其他的两类信号, 因为在金融市场上大变动的价格是十分罕见的现象。这意味着, 模型表现的精确性主要由最常见的结果持有信号所控制, 而这不是需要的结果。我们希望得到的模型对罕见信号 (买入和卖出) 能有准确的预测。因为这些信号会引发市场交易行为, 从而能达到获得潜在获利的最终目标。

金融市场预测是典型的罕见的事件驱动的例子。通常, 基于事件的预测任务通常是由精确度和召回率混淆矩阵来评价的, 其重点是对事件进行评估, 而忽视的常见情况下的表现 (如, 持有信号) 。精确度可以认为是模型预测信号正确的比例。召回率为模型预测出的信号数和原始数据中所有的相应信号的比率。这些指标可以很容易地通过在测试数据之间比较预测值和真实值的混乱矩阵来得到模型的结果。表1为混淆矩阵的一个例子。

从表1中, 对于这个问题我们定义的准确率和召回率, 如下所示:

同时, 还可以通过获得买卖信号的准确率和召回率来独立计算对某个特定信号的统计数据, 例如:

准确率和率也可以“合并”成单个统计量, 被称为在F-测量 (Rijsbergen, 1979) , 由下式给出:

其中0≤β≤1, 用来控制rec相对prec的重要程度。

4. 实证结果

4.1 实证数据

本文选取数据主要是为外汇市场上8个主要的外汇汇率, 数据主要来自Oanda, 时间段为2008年9月3日到2013年12月31日。其他的一些网站也有类似的数据, 如雅虎财经网站、新浪财经数据等。

所有模型都采用了三种不同的数据处理方式: (1) 单个模型应用所有的测试期数据方式; (2) 以固定更新w天的增长窗口方式; (3) 每次更新w天的移动窗口方式。图4显示的是三种不同的处理方式。

4.2 变量选择

上述各种指标变量用来预测未来的指标T值。我们构建随机森林模型 (Breiman Leo和Adele Cutler, 2001) , 来对模型中重要的指标变量进行选择。定义一个初始的指标集, 然后使用随机森林模型来评估每个指标变量的重要性。基于这些估计, 选择最相关的指标变量。随机森林也可以用来估计参与预测的任务的变量的重要性。通俗地说, 如果我们依次取出每个变量, 这种重要性可以通过计算百分比增加随机森林的错误进行估计。因为它包括在变量选择的过程的建模工具, 所以有类似过滤器的概念。然而, 这不是一个迭代搜索过程, 而且, 我们将使用其他的预测模型来预测T, 这意味着通过该方法选择的变量不为这些其他模型进行优化, 并在这个意义上, 此方法就像一个过滤器的方法。

(随机森林) 算法是通过训练多个决策树, 生成模型, 然后综合利用多个决策树进行分类。随机森林算法只需要两个参数:构建的决策树的个数t, 在决策树的每个节点进行分裂时需要考虑的输入特征的个数m。

4.2.1. 单棵决策树的构建:

(1) 令N为训练样例的个数, 则单棵决策树的输入样例的个数为N个从训练集中有放回的随机抽取N个训练样例。

(2) 令训练样例的输入特征的个数为M, 切m远远小于M, 则我们在每颗决策树的每个节点上进行分裂时, 从M个输入特征里随机选择m个输入特征, 然后从这m个输入特征里选择一个最好的进行分裂。m在构建决策树的过程中不会改变。

(3) 每棵树都一直这样分裂下去, 直到该节点的所有训练样例都属于同一类。不需要剪枝。

4.2.2 随机森林的分类结果

按照1生成t个决策树之后, 对于每个新的测试样例, 综合多个决策树的分类结果来作为随机森林的分类结果。

(1) 目标特征为数字类型:取t个决策树的平均值作为分类结果。

(2) 目标特征为类别类型:少数服从多数, 取单棵树分类结果最多的那个类别作为整个随机森林的分类结果。

后续对每个货币对汇率建模时, 就采用上表对应的技术指标与T值。如, 对EUR_USD建模时, 本文选取变量重要性大于10的指标相应的技术指标就选择:ATR, ADX, SAR, Mean和SD指标, 如图5所示。

三个模型准确度度量结果, 如下表 (3~5) 所示。从各个模型的准确性比较中可以看出, 对于EUR_USD和NAD_USD数据集支持向量机模型都有很好的预测结果, 而对于USD_CAD、GBP_USD和USD_CHF数据选用神经网络模型有更高的预测准确性, 同时AUD_USD、USD_JPY和EUR_JPY数据则是使用多元适应回归样条模型更好。

4.3 交易策略

我们的交易系统有两个交易状态:多头和空头。多头是指在t时刻以价格p购入资产, 同时在t+x卖出, 建立多头预期未来价格会上涨, 从而卖出资产获利。空头是指在t时刻以价格p卖出资产, 同时在t+x买入资产返还, 这是一种卖空的模式, 这种模式以期未来资产价格会下跌, 从而套利。一般来讲, 建立多头预示着我们期望未来的价格会上升;建立空头预示着我们期望未来的价格会下降。

使用蒙特卡洛模拟来选择具有稳定收益率的模型。其中通过改变各个模型的参数, 得到神经网络模型组成60个模型, 支持向量机模型也组成60个模型, 多元适应性回归样条模型组成120个模型。如图5所示, 我们在6年数据中随机地选择r个时点 (满足初始大于2年) , 同时选用该时点前两年数据来训练模型, 使用后1年数据来测试。对每个模型进行r次测试, 计算平均收益率等指标, 从而在选择较好的模型进行交易, 图6为某一模型策略的累积收益率变化情况。结果表明, 在EUR_USD和NAD_USD数据集上使用支持向量机模型可以得到很好的表现, 而在USD_CAD、GBP_USD和USD_CHF数据选用神经网络模型有更高的收益率, 同时AUD_USD、USD_JPY和EUR_JPY数据则是使用多元适应回归样条模型更好, 这个结果与模型预测准确度相吻合。表7、表9分别是使用两种不同策略后最优模型在测试数据上的收益结果, 相应的策略参数设置在表6和表8中给出。从中, 可以看出在每个数据集选用了最优模型之后, 都可以获得超额的收益, 比原来单一使用某个模型的收益得到大幅提高。

策略1:

其中, 各个评价指标的含义为:NTrades表示交易次数, NProf表示盈利次数, Perc Prof表示盈利百分比, PL为最终金额, Ret为收益率, Ret Over BH为超过买入持有策略收益率, Max DD为最大回撤额, Sharp Ratio为夏普比率, Avg Prof为平均收益, Avg Los为平均损失, Avg PL为平均资金额, Max Prof为最大利润, Max Loss为最大损失额。

策略2:

结论

本文使用蒙特卡洛模拟来确定不同市场上使用的模型, 从而使得在测试数据上能够获得更高的收益率。实证结果表明, 与单一的模型使用在不同市场相比较, 基于蒙特卡洛模拟的模型选择在测试样本上可以取得更好的收益水平。这个结果在一定程度上验证了AMH假设, 即单一的投资策略难以应付市场环境的变化, 只有不断地调整交易策略, 才能在市场上获得超额收益。未来的研究可以在结合更多的模型进来, 同时增加样本量来增加模型的稳定性。

摘要：为了选择不同外汇市场交易采用的模型, 本文使用蒙特卡洛模拟方法比较构建了三种模型 (神经网络模型, 支持向量机模型, 多元适应性回归样条模型) 的市场适应性。考虑外汇市场上8种主要的外汇货币, 通过方法选择出来的模型进行交易, 发现不同的外汇价格时间序列选择模型有所不同, 同时得出结论:不同的外汇市场上使用不同的模型会比单一模型获得更高的收益率, 这在一定程度上也解释了适应性市场假设 (AMH) , 即单一的投资策略难以应付市场环境的变化。

关键词：蒙特卡洛模拟,外汇交易策略,适应性市场假设 (AMH)

参考文献

[1]Park C H, Irwin S H.What do we know about the profitability of technical analysis?[J].Journal of Economic Surveys, 2007, 21 (4) :786-826.

[2]Friedman, J. (1991) .Multivariate adaptive regression splines.The Annals of Statistics, 19 (1) :1-144.

[3]Shawe-Taylor, J.and Cristianini, N. (2000) .An introduction to support vectormachines.Cambridge University Press.

[4]周东生, 唐焕文.有效市场假说与市场适应性有效[J].管理科学, 2005, 18 (3) :76-80.

蒙特卡洛模拟定价 第6篇

采用逐一试验的方法同时考虑这些因素对分离过程的影响,将是一项庞大而费时的任务。蒙特卡洛方法是一种通过随机变化的统计获得近似结果的方法,可用来解决具有不确定性的复杂问题。

本文以某型导弹级间分离为研究对象,将飞行力学和分离动力学相结合,建立考虑各种干扰因素和偏差的数学模型,并采用MATLAB/Simulink搭建了分离仿真模型。为了计算各种随机干扰因素的任意组合情况,采用蒙特卡洛方法对分离过程进行仿真计算,得到了各种干扰因素对分离的影响程度。

1 影响级间分离的因素分析

在分离过程中,受到各种干扰因素的影响,使得导弹的弹道参数偏离理想弹道而可能造成两级发生碰撞。本文考虑的干扰因素主要有:

1)初始弹道参数(质量、速度、攻角)的偏差以及分离机构工作不同步等。初始弹道参数偏差在仿真中通过改变模型仿真初值来实现。

2)风场的影响。按照空气团的运动特性,风可以分为常值风和变化风。常值风是指在一定的空间和时间范围内,风速矢量为常值的风;变化风是指在一定的空间和时间范围内,风速的大小和方向都发生变化的风,变化风即为大气紊流,其一般采用Dryden模型或Von Karman模型[1,2],具体建模过程这里不再详述。

风场的干扰使得仿真模型中气动力和气动力矩需要采用相对速度Va,相对攻角αa和相对侧滑角βa进行计算。图1给出了速度矢量关系,VW表示风速,Va表示导弹的相对速度,V表示导弹原来的速度。

由速度矢量关系可得Va在地面坐标系中的投影为:

式(1)中,(Va.gx,Va.gy,Va.gz)T是Va在地面坐标系的分量;(Vx,Vy,Vz)T为V在地面坐标系的分量;(VWx,VWy,VWz)T为VW在地面坐标系的分量。

导弹相对速度在弹体坐标系中的分量为:

式(2)中Cbg(ψ,ν,γ)为从地面坐标系到弹体坐标系的坐标转换矩阵,导弹相对速度的大小在弹体坐标系可表示为:

由攻角αa和侧滑角βa的定义可得:

从上面的公式可以看出由于大气风的存在,使得导弹出现附加的攻角和侧滑角,同时实际相对速度也可能变大或变小。

3)助推器发动机推力偏差的干扰。推力偏差[3]包括推力偏斜和推力偏心,是导致分离干扰的重要因素之一,对导弹助推级和主级飞行及级间分离过程都产生重要的影响。

设发动机推力P与弹体纵轴不重合,推力向量和弹体纵轴的关系可用倾斜角θp和偏角ψp两个角度来表示,θp是推力向量P与Obxbzb(弹体坐标系)平面的夹角,向上为正;ψp是P在Obxbzb平面的投影与Obxb轴的夹角,ψp在Obxbzb平面的右侧为正,如图2所示。

推力在弹体坐标系中的分量可表示为:

通过坐标转换,可推导出推力在弹道坐标系中的分量:

式(7)中,(Pkx,Pky,Pkz)T为推力在弹道坐标系中的分量,Ctg(ψC,θ)为从地面坐标系到弹道坐标系的坐标转换矩阵。

设推力P的作用点为Op;Op在弹体坐标系下的坐标为(xp,yp,zp)T,如图3所示。

从导弹质心Ob到Op的矢径为rp=[xp,yp,zp]T,则推力对导弹的力矩为:

推力偏斜和偏心与发动机的工艺控制水平有关。由于推力偏差的存在使得导弹在俯仰、偏航和滚转方向都产生了干扰力矩,当其偏差过大有可能导致导弹飞行失稳。

2 分离过程数学模型

2.1 基本假设

1)导弹的助推级、主级分离发生在稠密大气层内,尚未脱离大风区。

2)忽略地球自转运动和导弹弹性运动,将导弹本体及助推级、主级分离视作刚体。

3)分离解锁瞬间完成,不计分离不同步的干扰。

4)导弹轴对称,各种干扰和偏差在各方向上出现概率相同、大小相等。

2.2 分离动力学方程

为了研究导弹级间分离的运动特性,需要分别建立整体、主级导弹、助推器各自的动力学和运动学数学模型。考虑到模型的通用性,下面以导弹整体为例给出其动力学和运动学数学模型[4]。

为了便于分析导弹运动特性,将导弹质心运动动力学方程建立在弹道坐标系下。弹道坐标系下的导弹质心运动动力学微分方程为:

式(9)中,G为重力,A为气动力,P为发动机推力,可表示为F=P+R+G。式(9)在弹道坐标系下的分量形式可以表示为:

式(10)中Fx、Fy、Fz为合力F在弹道坐标系的分量,V、θ、ψC分别表示为导弹速度,弹道倾角、弹道偏角。

在弹体坐标系上建立的导弹绕质心转动动力学方程:

式(11)中:J为转动惯量,ω为转动角速度,M为导弹受到的合力矩。

因为导弹是轴对称外形,弹体的所有惯性积为零,将式(11)分解到弹体坐标系中可得导弹绕质心转动的动力学方程为:

2.3 分离运动学方程

在地面坐标系下建立的导弹质心运动学方程:

在弹体坐标系下的导弹绕质心转动的运动学方程的运动方程为:

式(14)中,ν、ψ、γ分别表示俯仰角、偏航角和滚转角。ωx、ωy、ωz为导弹绕自身三个坐标轴的旋转角速度。

2.4 碰撞分析

在级间分离研究中最关心的是主级和助推器是否发生碰撞。本文采用危险截面特征点方法来判断两级是否发生碰撞。分别在主级弹体坐标系和助推器弹体坐标系下选取n个可能发生碰撞的点,记为P=[P1…Pj…Pn],Q=[Q1…Qj…Qn],如图4所示。为了计算Pj与Qj之间的最小距离,将其转换到地面坐标系下,记为P'、Q',可求出主级与助推器之间的最小距离为:Lmin=min(‖P'i-Qj'‖i,j=1,…,n)。在分离过程中通过计算最小距离Lmin来确定主级尾段与助推器前段是否发生碰撞,当Lmin<0认为两级发生碰撞。此外,在级间分离仿真中,还需要考虑助推级和主级运动状态是否平稳、分离后主级控制系统是否达到了能力极限等。

3 蒙特卡洛方法

蒙特卡洛方法是一种通过随机变量的统计试验或随机模拟,求解数学、物理和工程技术问题近似解的数值方法[5],其通过对系统输入偏差在一定范围内多次随机抽样组合,获得输出样本的统计特性,并发现可能导致系统输出超出预期临界的偏差组合。在应用蒙特卡洛方法进行级间分离模拟仿真时,其基本步骤是:

1)确定级间分离过程中受到的干扰因素及分布规律。

2)根据级间分离动力学,建立考虑各种干扰因素的分离数学模型。

3)根据各种随机因素的分布律,构造相应的数学概率模型。

4)根据模拟精度要求确定仿真次数。

5)随机抽样进行多次模拟打靶,获得干扰作用下的弹道参数。

6)对结果处理,得到级间分离的统计特征。

4 仿真结果与分析

本文基于受扰分离动力学数学模型,在MAT-LAB/Simulink环境下编写了蒙特卡洛仿真程序。仿真初始条件为:发射高度为0 km,发射角度45°,当导弹的过载为-3时,给出级间分离指令,主级与助推器最小距离大于4~6倍弹径时认为分离结束,这里取为3 m。各种干扰因素及服从分布类型如表1所示,以表1中的干扰作为随机量进行蒙特卡洛仿真,试验次数取1 000次,所需要的机时不到30min,通过仿真发现仅有1次分离失败,分离成功率为99.8%,失败的情况为推力偏心、偏斜取最大值、逆风大气紊流模型、风速平均值达到10 m/s组合情况。表2给出了分离前后的弹道参数统计特性,由表数据对比可知,分离前后主级导弹的速度损失约为9.2 m/s,助推级的速度损失为57.02 m/s,这主要是由于助推器是一个钝头体,且处于主级喷管冷气流区域,受到的气动干扰较大的缘故。

两级之间的相对距离随时间变化如图5所示,从图中可以看出,由于分离机构工作不同步(0.01s)的影响,两级之间在开始时并没有分离,随着分离机构的完全分离,两级之间的最小距离才开始逐渐单调增大,在0.129 s的时候相对距离的最小值已经达到3 m,分离迅速,两级没有发生碰撞。

图6为蒙特卡洛仿真结果曲线。图6(a)和图6(b)分别给出了主级和助推器攻角随时间的变化。由图可知,助推级姿态变化比较剧烈,攻角由0.015°到2.524°的变化,而主级在分离后攻角变化趋势平缓。主要是由于分离机构的不同步引起了导弹姿态的变化,而主级质量是助推器质量的9倍左右,分离机构不同步引起的干扰对主级的姿态影响比较小。由图6(c)可知,导弹在助推段攻角变化较大,主要是因为助推开始段速度较低,受风干扰和推力偏心偏斜的影响较大,从而影响了分离时刻的导弹姿态。

5 结论

本文以某型导弹级间分离为研究模型,以正态分布的随机量来模拟各种干扰的偏差,提出了受扰的导弹弹道仿真模型的构建方法,利用蒙特卡洛方法对级间分离过程进行仿真模拟,得到了有益的结论:

1)蒙特卡洛方法可方便、快速地对各种偏差组合情况进行仿真,通过大量的仿真计算,可以全面地分析各种可能出现的情况和导致分离失败的最严酷组合,从而减小飞行试验的风险。

2)风场干扰、推力偏心偏斜对助推初始段弹道参数影响很大,严重时可能导致导弹在助推段失控。由于风干扰是无法控制的,因此必须对推力偏心偏斜、结构安装偏差等在生产制造过程中是可以加以控制的,以提高分离的安全性。

3)分离机构不同步时间越长,其带来的干扰力和力矩对级间分离影响越大,严重时会导致分离失败。

4)两级导弹分离过程因各种偏差和扰动的存在需采取安全余量设计。

参考文献

[1] Anjan C,Langelaan J W.Flight path planning for UAV atmospheric energy harvesting using heuristic search.AIAA Guidance,Navigation,and Control Conference,August Toronto,Ontario Canada,2010:1—18

[2] Gedeon J.Dynamic analysis of dolphin style thermal cross country flight.Technical Soaring,1973;3(1):9—19

[3] 杨涛,王中伟,张为华,等.固体发动机推力偏差对导弹级间分离的影响.固体火箭技术,2006;29(1):1—4Yang Tao,Wang Zhongwei,Zhang Weihua.Effect of solid motor thrust biason on stage separation of missile.Journal of Solid Rocket Technology,2006;29(1):1—4

[4] 李新国,方群.有翼导弹飞行动力学.西安:西北工业大学出版社,2005:56—80Li Xinguo,Fang Qun.Flight dynamics of missile.Xi'an:Northwestern Polytechnical University Press,2005:56—80

蒙特卡洛模拟定价 第7篇

关键词：连续油藏模拟,马尔可夫链,蒙特卡洛,概率预测,不确定性量化

1 背景

1.1 不确定性量化技术

对产量预测中的不确定性进行量化完全依赖于历史拟合。进行历史拟合包含三个主要步骤:第一, 根据描述油藏几何与流动特性的整套参数来定义油藏;第二, 根据预估的概率分布对油藏未知参数进行赋值;第三, 从预先的分配油藏模型可作为实际模拟例子, 并在目标函数中量化模拟结果和观察数据的差异。许多方法用于寻求合理的模型, 某些历史拟合的技术梯度法技术, 它们的目标就是最优化, 即模拟结果最符合观测结果的模型就是所需要的模型。用梯度法进行的不确定性量化也许还不是最完美的方法, 因为它只在目标函数的一个或多个极疏区的相邻区取样。MCMC方法已广泛应用于从一个复杂的分布函数中取样, 特别是当不知道函数的精确形式的时候。

在油藏模型研究中, MCMC技术用于研究贝叶斯推断的后验分布。最后分布只不过产生于一套案例中油藏模型所用的样本。首先, 一个随机模型取样于一个先验分布的样本, 这是马尔可夫链的出发点。其次, 第二个模型被随机选择但受限且关联于马尔可夫链中已存在的前一个模型。当这个马尔可夫链足够长后, 可以用此链中的模型产生后验分布。生成未知分布的另一种相关方法叫作遗传算法 (GAs) , 它有多个相关的应用软件。遗传算法是基于控制自然界遗传基因规则、类别涉及面广的一类最优化算法。

在遗传算法中的一个称为“繁殖”的过程中, 特种油藏模型的子模型产生于先前运行的模型混合参数值。最后应用所有的子模型生成一种分布。较之GAs, MCMC方法在统计上更严格。MCMC也被当作是遗传算法的一个类型, 因为下一个模型依赖于其前一个模型, 这在GAs就是其母模型。

1.2 实时数据与综合最小方差线性递推估算滤波器

综合最小方差线性递推估算滤波器 (卡尔曼滤波器) 技术广泛应用于天气预报, 近来应用于石油工程以进行概率预测。卡尔曼滤波器是一种蒙特卡洛手段, 是一种通过连续模型校正而进行不确定性量化的技术。较之传统的历史拟合中进行同步数据同化, EnKF (综合卡尔曼滤波器) 要顺序进行。EnKF从受静态数据 (如岩心、测井及构造资料) 约束的整体油藏模型开始运行。当数据可用时, EnKF用来自于综合模型的静态信息和模型预测数据更新每个模型的实现。因此, EnKF生成一套受生产历史约束的拟模型实现, 而且从理论上讲, 它应该尊重先前的静态或地质信息。然而, 由于受依赖于样本数量或模型实现的所有统计测量中固有的杂波影响, 对于小的整体规模, EnKF更新能导致地质不协调实现。因此, 如果最终的模型实现尊重历史生产数据, 那么模型可能就不符合先前的地质信息。在EnKF过程中使用大的整体规模可以减少运行的困难, 但对于油田规模计算费用可能比较昂贵。EnKF仍是为达到历史拟合而进行积极研究的主题, 并且该技术正致力于解决运行中的难题。

2 连续模拟过程

用连续MCMC方法进行历史拟合以及生成概率预测需要几个部分相结合:首先, 不确定性油藏的参数以及先验分布必须是确定的。其次, 既然将来的油藏动态的不确定性通常要通过一套油藏模型的模拟动态来评估, 那么必须有在参数空间取样以及生成油藏模型的方法。MCMC方法在这里被用来探究不确定性参数空间, 并自动生成运行模拟的编码。针对每一示例模型, 随着历史拟合的进行, 预测也一并进行, 同时运行结果也被存储起来。在连续的历史拟合过程中, 油田现场的新数据也一直被加入目标函数中。被更新的目标函数致力于后来的模拟运行。最后, 各个运行结果被合并入概率预测。

2.1 参数化

在作任何模拟之前都必须首先决定考虑哪些不确定性参数。一旦确定了相关的参数, 就要指定先验分布 (这里通常是连续的) 以在数量上描述参数的不确定性。正确的分布类型通常是建立在油藏特征数据上的。后期的分布就在贝叶斯框架内定义, 而先前的分布根据观测的动态数据作出修正, 表示如下:

式中dobs从油田现场上观测的动态数据, 如含水率;

m不确定性参数;

P (m) 不确定性参数的先验概率分布;

P (dobs|m) 与观测数据有关的似然函数;

P (m|dobs) 后期的分布函数。

如果假设先验模型和数据误差服从正态分布, 那么这里的后验分布P (m|dobs) 就成为如下形式 (Howard, 2005) :

式中g (m) 符合不确定性参数的先验概率分布m的模拟油藏响应 (曲线) ;

Cm参数协方差;

CD资料协方差。

2.2 目标函数

目标函数是从数量上评估一个独立模型再现油田现场观测数据的满意程度的函数。这一项被定义为如下后验分布函数的一部分:

在 (3) 式中, (m-μ) TCm-1 (m-μ) 是先验项;[g (m) -dobs]TCD-1[g (m) -dobs]是似然项;μ表示先验平均值。目标函数的结果就是先验信息与观测信息的结合是一个在贝叶斯框架内后验分布结构的推论。

2.3 Metropolis-Hasting MCMC算法

该MCMC方法的主要目标是构造一个平稳分布与后验分布相匹配的马尔科夫链。这种后验分布典型地定义在多维参数空间上并常有多路系统式。这里所使用的Metropolis-Hasting MCMC方法 (Hastings, 1970) 常被用来从复杂后验分布中采样, 其随机步进M-H采样过程如下:

(1) 从先验分布中随机采样一套参数, 表示为mt1;

(2) 从状态ti时刻到ti+1, mti+1=mti+σε, 其中, ε是一个标准的正态随机变量;σ是一个比例因子;

(3)

(4) 从0～1间的均匀分布中随机取一个数y。如果yR, 则接受链中的mti+1;否则又将mti放入链中;

(5) 回到第 (2) 步。

2.4 连续数据的同化作用

在连续模拟过程中的不同时间点上, 油田现场新数据就显得很有用了。由于通常会设想油田现场的信息越丰富则导致预测和不确定性评估的效果就越好, 因此尽快在模拟过程中引入新数据是有利的。当更多的数据加入, 则 (3) 式中将包括更多的观测及模拟数据点。这样, 观测数据的限定不符合目标函数中相关项的状况就将随着时间的推移而改变。尽管持续变化的目标函数在统计上是不严格的, 但可以假设在这种违规情况下能产生合理的概率预测。

2.5 概率预测

模拟过程的最后一步是将为单独的马尔可夫链样本所做的生产预测合并入概率预测中。在连续的MCMC过程中, 通过使用足够多的最新采样模型可在任何时候生成概率预测。要么用少量的模型以使改变目标函数的影响最小, 要么用大量的模型以使获取可代表后验分布的平稳系列分布的机会最大, 这二者之间需要寻求一个可以接受的平衡。

3 PUNQ-S3油藏实例研究

PUNQ-S3人造油藏作为在生产预测中量化不确定性的测试实例已应用了很多次。PUNQ-S3油藏模型是一个基于实际油田的5层3相人造油藏 (Bos, 1999) 。该模型包括13 300 (192825) 个网格块, 其中1761个网格块起作用。油藏被东部和南部的断层所限制并与西和北边相当厚的含水层有关系。带井位的结构如图1所示, 图中有一个小气顶在中央并以红色显示, 6口生产井表示为黑色的圆点。石油行业和学术理论上的几个合作伙伴用不同的方法来检验历史拟合与不确定性量化技术。这个PUNQ项目的目标就是确立一种方法可以通过适当方式将油藏模型、油藏参数、油井观察结果的不确定性与生产预测的不确定性结合起来。

用来确定不确定性参数及指定先验分布的过程与传统上模拟研究中评估输入不确定性时所做的工作是一致的。在PUNQ-S3模型中, 假设孔隙度呈正态分布, 而渗透率呈对数正态分布。另外, 研究中的不确定性参数不直接使用孔隙度和渗透率值而是孔隙度和渗透率的乘积, 在运行模拟过程中, 这些乘子应用于孔隙度和渗透率底图。这样做的效果与直接使用孔隙度和渗透率值的效果相同, 但是这种方法简化了实施过程。对每一层使用统一的属性值绘制底图。基于Barker 2001年提供的平均孔隙度常量值, 垂向和纵向渗透率通过 (4) 、 (5) 式 (Gu and Oliver, 2004) 计算得到。表1列出了底图中使用的平均孔隙度和渗透率值。

式中, kv表示垂向渗透率;kh表示横向渗透率。

用每层6个同源区对PUNQ-S3模型进行参数化。另外, 通过分裂油藏的方式将各区视为相互独立的区域, 这一点是基于实时地质描述的, 指示油藏广泛分布着东南走向的优质夹层 (Imperial College, 2007) 。这样, 在先验分布中, 5层6区3个属性 (孔隙度、垂向渗透率和横向渗透率) 就产生总共90个未确定的乘子。孔隙度乘子被指定服从中值为1、标准偏差为0.3的正态分布 (Barker et al, 2001) 。基于 (4) 、 (5) 式, 垂向和横向渗透率乘子被指定服从中值为1、标准偏差为1.35的对数正态分布。为防止出现极端和不实际的渗透率值, 此乘子分布被冠以上不过4而下不过0的限制。同样, 孔隙度范围在0～2.28。

1 mD=1.0210-3μm2

3.1 目标函数

由于在研究中确定的不确定性参数是孔隙度和渗透率, 目标函数变为:

式中, X代表不确定的乘子。

在实验过程中, 在第5、6、7、8、9年初顺序加入历史资料。因此, 当加入新数据时, 目标函数中的观测资料值的数量成倍数增加。

3.2 参数空间研究

在连续历史模拟过程中, 模拟运行过程是与历时8.04年的观测资料进行拟合。进行的预测是16.5年的生产经历。PUNQ-S3模型连续模拟开始于第4.5年, 持续到第9年底, 并作出16.5年的预测。在第一个半年 (从第4.5年到第5年) , 用MCMC采样方法的模拟中有4 500个模型作为样品, 并且在其余的年份里每年样本模型数为9 000, 最后模拟样品总数是49 500。所有这些都假设9 000个模型的运行都发生在一年中的实际时间, 即每小时运行一个模型, 这不是今天石油行业中用于油田模拟的非典型的运行频率。图2显示模拟运行的累积次数与油藏生产时期的对比。包括在目标函数计算中的额外数据也导致目标函数更大且引起目标函数曲线中太多的数值起落跳跃 (图3) 。

3.3 预测

一个独立的预测是与每一次历史拟合一起运行。在第5、6、7、8、9年和第10年, 利用在链中可用的样品模型在相应的时间进行概率预测。通过过去一年 (或第5年中的半年) 的模拟运行可以生成这些概率预测, 这些概率预测的累积分布一起显示在图4中。

3.4 不确定性估计的校准

从图5中可以看出, 前3个预测没能反映真实情形, 在4.5～5年间偏离明显。因此, 在前几年里不确定性明显淡化了, 这可能归因于先验分布中对不确定性的低估, 或是对观测资料误差的低估。即使没有大量的可用动态资料, 人们也总是要去量化不确定性。为了增加不确定性, 可以在先验分布或观测资料或两者中使用更大的标准差为了克服这些缺点, 增大了先验乘子的标准差。渗透率乘子标准差从1.35增大到20, 而孔隙度乘子标准差从0.3增大到0.5。伴随着标准差的增大, 先前的后验分布被重新建立, 直到它们在本质上支持所有继发的后验分布。

4 结论

MCMC方法是一种应用于历史拟合与不确定性量化的强有力工具。通过对油藏的历史拟合及其整个生命周期的连续预测, 随着时间的推移和观测资料的增加, 预测不确定性范围逐渐变小。较之MCMC方法传统上应用于过去一时间上的模拟研究, 连续MCMC方法可用很少的模型在特别的时间点上进行合理的概率预测, 这是因为连续方法得益于最新观测资料同化前的模型运行。在油藏的长期生产中持续建模可考虑以下方面:在参数化过程中获得更多的不确定性参数, 更多的样本以及由此造成的参数空间大得多的片断取样。这应该最终导致更可靠的概率预测。

连续的模拟方法也提供一种从始至终的不确定性估计校准机制。如果观测到后验分布除了变窄以外始终发生的跳跃起伏, 那么可以放大先验分布中的标准差或者增加观测资料以增加不确定性。应该采取调整措施直到继发的后验分布与所有先前的分布相一致。