扩散参数（共4篇）
扩散参数 篇1
MRI是目前诊断颈椎管内肿瘤等颈髓疾病的最佳影像学方法,但临床上有很多颈髓病变通过常规MR的T1加权像和T2加权像观察缺乏敏感性。扩散张量成像(diffusion tensor imaging,DTI)可以为脑部病变提供组织结构评估及生理学参数分析[1]。脊髓DTI已从大鼠体内外研究逐步进入人体的初级应用阶段,为下一步针对颈髓疾病的研究建立了一条基线。本研究通过观察健康志愿者颈髓DTI参数的正常值,探讨其分布特点、变异及影响因素。
1 资料与方法
1.1 研究对象
2014年2—12月招募36例健康成年志愿者,其中男21例,女15例;年龄18~77岁,平均(43.6±17.6)岁,年龄分布见图1。所有受检者均无神经系统症状和阳性体征,并排除既往有神经系统疾病史、MRI检查禁忌及MRI扫描T2像在颈髓内出现高信号者。
1.2仪器与方法
常规MR扫描采用Siemens 3.0T MR TRIOTIM扫描仪,通过快速自旋回波(FSE)序列对所有受检者行T1WI矢状位及T2WI轴位、冠状位、 矢状位扫描。T1WI扫描参数:TR 735 ms,TE11 ms。T2WI扫描参数:TR 3130 ms,TE 91 ms。 矢状位扫描视野(FOV)300 mm, 层厚3.0 mm, 层间距0.3 mm。轴位T2WI扫描第1颈椎至第7颈椎(C1~7)水平,FOV 160 mm,层厚3.6 mm,层间距0.4 mm。采集矩阵448 × 314。T2WI矢状位3D容积扫描参数:TR 3000 ms,TE 502 ms,层厚1.0 mm,层间距0,层数128层,FOV 250 mm × 250 mm,采集矩阵256 × 256。DTI采用敏感编码回波平面成像(SE-EPI),在3D扫描基础上进行矢状位tensor序列扫描。扫描参数:FOV 240 mm × 240 mm,TR 5600 ms,TE 93 ms,层厚3.5 mm, 层间距0, 采集矩阵128 × 128, 共38 层。另扫描一幅b =0 的无扩散加权图像。激励次数1 次,b =1000 s/mm2;扩散梯度方向数(NDGDs)12 ;相位编码方向:头- 脚方向。DTI扫描时间约为1 min 31 s。
1.3图像后处理
应用NUMARIS/4软件syngo MRB15版本后处理工作站,进行数据采取和影像重建,由2名影像科主治医师采用双盲法完成,所有DTI参数的最终测量值取两者实际测量值的平均值。
1.3.1各向异性分数(FA)和表观扩散系数(ADC)值测定
在Neuro 3D环境下启动tensor序列,于b值为0 的FA图上在C1~7水平勾画出感兴趣区(ROI)。每个ROI直径为7 mm,位置选择在正对每个节段颈椎椎体中央。此时软件可自行运算得出各个颈髓水平的FA值和ADC值。
1.3.2扩散张量纤维束示踪(diffusion tensor tracking,DTT)的生成
在扫描生成的38 层矢状位图像中选取正中层面生成DTT。设置纤维示踪停止条件:轨道角35°,FA =0.18。纤维束生成后与颈椎管3D矢状位图像对应,观察其走行及完整性。
1.4统计学方法
采用SPSS 21.0 软件,两组间计量资料比较采用成组资料t检验,多组间计量资料比较采用单因素方差分析,若差异有显著性,则进一步以Bonfferoni法行多组间均数的多重比较;两变量间的相关性采用Pearson相关分析,P<0.05表示差异有统计学意义。
2 结果
2.1 FA、ADC图效果及DTT效果
分别见图2、3。DTT图像可清晰显示颈髓内白质纤维束结构,未出现图像扭曲和失真现象。纤维束呈现头尾方向,具有3D效果,可以通过任意角度旋转观察。T2WI证实为正常颈髓的纤维束呈现出饱满、顺滑且完整的结构(图3)。
2.2 FA值、ADC值与年龄的关系
全颈髓FA值、ADC值与年龄的Pearson相关分析显示,FA值与年龄呈显著负相关(r=-0.801,P<0.001);ADC值与年龄呈正相关(r=0.426,P<0.05)。见图4。
2.3 FA值、ADC值与性别的关系
不同性别受检者全颈髓FA值和ADC值比较,FA、ADC值在男、女两组间差异无统计学意义(t=0.906、0.362,P>0.05),见表1。
2.4 FA值、ADC值与颈髓水平的关系
C1~7颈髓水平FA值比较,最高FA值出现在C2水平,最低在C6水平,各组间差异无统计学意义(P>0.05)。C1~7颈髓水平ADC值比较,最高ADC值出现在C6水平,最低在C2水平,各组间差异无统计学意义(P>0.05),见表1。
3 讨论
3.1 DTI用于定量分析的原理
扩散加权成像能够探测、衡量组织中水分子的移动。水分子发生位移现象后会对扩散加权扫描时的信号产生衰减效应。由于其特性所致,中枢神经系统中白质轴突结构会促使水分子主要以平行于轴突纤维的方向扩散,而不是垂直于轴突纤维[2]。垂直于纤维方向扩散受限的原因除髓鞘的因素,更多的是源于细胞膜[2]。DTI利用这种与方向有关的扩散特性,可以推断轴突纤维的定位和划定解剖界限。DTI采用张量的框架来表征三维空间内的分子运动。通过计算沿3个主轴方向的扩散系数,即可产生DTI参数,故所有DTI参数都能通过数学公式计算得到。尽管其计算公式复杂,直接进行手工计算困难,但是软件的应用使这项工作变得容易。与DTI有关的参数有ADC、FA、RA、VR及DTT。其中FA值介于0~1之间,1代表假想下最大各向异性的扩散,其数学计算公式为[3]:
其中,λ1、λ2和λ3分别为椭球体最长径、前后径和左右径的扩散强度,λ=(λ1+λ2+λ3)/3。
3.2颈髓DTI参数特点
目前关于颈髓DTI正常值的文献报道少见,且报道的数值差异较大[4,5,6]。研究样本量的大小及DTI采集参数如ROI所取面积的大小(体素的含量)和NDGDs等均会对FA值造成影响,其颈髓FA值为0.6~0.7,标准差约为平均值的10%,即0.06~0.07[4,5,6]。本研究在C1~7水平测得的平均FA值为0.678~0.721,标准差为0.056~0.071 ;较上述报道结果平均值略高,标准差略低,推测变异较小的原因可能与所取样本量较大、受试者年龄跨度较大有关。FA值、ADC值在各颈髓水平间的分布呈从头侧向尾侧FA值逐渐下降,ADC值逐渐上升[4],该分布特点的可能原因是颈膨大使低位颈髓中灰质成分增多,或臂丛神经根进出低位颈髓造成此处的FA值降低[6];但Chang等[7]报道无此差异。本组资料FA值在颈髓的这种变化趋势不明显,尽管其最高值出现在C2水平,最低值在C6水平,但FA值在不同颈髓水平的总体差异无统计学意义,其原因可能是b值取1000 s/mm2,造成低位颈髓信噪比较低,FA值被高估。同样,这种差异体现在ADC值上也不明显,但是ADC值与FA值呈显著负相关。各颈髓水平的ADC平均值为1.072~1.205,标准差为0.178~0.237,变异大于FA值,提示FA值是DTI更可靠的定量分析指标。在DTI检测脊髓损伤[8,9,10]和颈椎病[11,12]微结构变化方面,FA值是敏感指标。
目前关于FA、ADC值与年龄关系的报道局限于某些特定的颈髓水平。Mamata等[13]在研究一组颈椎病患者时发现,在C3~4正常脊髓区域,FA值与年龄呈负相关,ADC值与年龄呈弱正相关;Agosta等[14]针对上颈髓的研究表明FA值与年龄呈负相关;但Brander等[4]报道DTI参数值与年龄无关。本研究结果显示,DTI参数值与年龄有明显的相关性,其中FA值与年龄呈显著负相关,ADC值与年龄呈正相关,进一步证实FA值对与年龄相关的脊髓结构改变敏感,同时提示对脊髓病变行DTI研究时需要设立年龄匹配的对照组。
3.3颈髓DTI成像质量的影响因素
DTT是利用扩散张量所含的信息,运用各种示踪技术重建出三维的白质纤维束。目前DTT大致分为线形扩展技术和能量最小技术,后者又可分为快速行进法和模拟退火法。本研究采用线性扩展技术,FA阈值设为0.18,轨道角35°。示踪时沿具有最大本征值的本征向量方向扩展,当FA值小于所设阈值、两个夹角>35°时终止进程。低位颈髓易受部分容积效应的影响造成空间分辨率降低。本研究受益于3.0T MR扫描仪的应用,减少上述不利影响。3.0T扫描仪可生成具有更高分辨率和更高信噪比的图像。尽管不会直接影响DTI参数值,但是高场强可以提高参数值的精确度,增强灰白质区的对比度,并促使纤维束示踪进程更加准确、有效[15]。Holder等[16]报道,在轴位和矢状位上测量的FA、ADC值无明显差别。因此,本研究仅采用矢状位上的DTI参数值测量。NDGDs会对DTI图像质量造成影响。随着NDGDs增大,尽管可获得更高信噪比的图像,但扫描时间也随之延长,患者可能不能忍受而易产生伪影,造成图像质量下降。b值取值越高,产生的信号衰减越大,信噪比越小,得到的图像信息越少;b值取值过低,又会使对比度下降。激励次数设高后,信噪比提高,但同时扫描时间也成倍延长。扫描的层厚也会影响图像质量。层厚越厚,尽管信噪比越高,但部分容积效应也被放大。因此,具体实施时应经多方面权衡后选择合适的扫描参数组合。本研究在正式实施前用多种参数组合进行预实验,最后采用NDGDs =12、b值取1000 s/mm2、激励次数1、层厚3.5 mm的组合,认为由此产生的图像质量最好,一次扫描成功率高达97.2%。后处理过程中ROI的选取也可影响成像质量。ROI取值越小,容积效应越小,理应获得更准确的测量值,但也丧失了部分扩散信息。本研究的做法是在尽可能避免周围非脊髓结构的基础上选取最大面积脊髓的ROI,由此ROI获得的测量值比较稳定。
3.4颈髓DTI的局限性
由于椎管空间狭窄的特性、磁敏感性伪影的存在及呼吸循环脏器的运动,使得DTI技术在脊髓中的应用受到限制。获得足够的空间分辨率仍然是一个难题,并且也很难成像出白质中的每一条纤维束。扫描时间对于急性脊髓损伤患者尤其是一项限制,这些患者因疼痛等原因更难忍受附加的扫描时间。此外,信噪比在整个颈髓中呈现不均匀现象,越往尾端信噪比越低。低信噪比可导致各向异性测定值高估,特别是在低各向异性组织如灰质结构。虽然3.0T MR扫描仪的使用可提高信噪比,但目前还不能普及应用。此外,有必要开发出标准化的软件来处理张量图像,使之在临床常规使用成为可行。
总之,通过测定FA和ADC值,使颈髓MRI能够用于定量分析,加上示踪技术的应用可直观显示白质纤维束的结构,显著增强了MRI在颈髓疾病诊断中的作用,并使其用于评估颈髓超微结构所发生的细微变化成为可能。健康人群脊髓DTI的研究为下一步针对颈髓疾病的DTI研究建立了基线。全颈髓FA、ADC值与性别和颈髓水平无关,但与年龄有关,提示年龄是混杂因素,在研究颈髓疾病DTI时需设立年龄匹配的对照组进行分析。
扩散参数 篇2
当前对大气扩散参数的测定方法有2种, 即由湍流量测量风速脉动量确定和利用气象站常规仪器观测进行分类参数化的方法[1,2,3,4]。该文利用固定观测点的湍流脉动量直接计算大气扩散参数, 这样就能较准确地反映出兰州市西固区的扩散特征。
1 观测仪器和资料
采用的仪器为三轴风速仪, 仪器的2个水平探头分别指向N、W方向, 第3个探头指向天顶。测量时间段为1 h, 每2s测量1次数据[5]。由于观测资料处理、时间原因, 该文只对稳定和不稳定条件下的大气扩散参数进行分析。
2 兰州市西固区地形及气候特征
兰州市西固区位于东经103°19′~104°41′, 北纬35°38′~36°13′。总体地势西北高、东南低, 为半封闭的哑铃状河谷盆地, 东西长而南北窄, 南北两山对峙环绕, 山体走向为东北西南向, 由北逐渐向南倾斜;以黄土梁和黄土河谷为主。
兰州城区年平均风速为0.7 m/s, 且呈单峰型分布。其中春、夏2季风速比较大, 秋季次之, 冬季最小;城区年平均静风率为74.8%, 冬季静风率高达87%以上;兰州城区一年四季都有逆温层存在, 尤其以冬季的贴地逆温最为突出;气候属温带半干旱气候, 山区垂直性气温变化较显著。年平均气温8.5℃, 最高36.1℃, 最低-23.4℃, 年均降雨量330 mm。
3 观测原理
根据Taylor湍流统计理论, 利用风速脉动资料计算大气扩散参数方法[6]。
3.1 三维速度平均值
假设三轴风速仪的三维风速序列分别为ui0、vi0、wi0, 三维平均速度和此时的水平平均风速计算公式分别为:
方向角θ0 (以正Y0向正X0方向旋转的角度为正) 为:
在经过坐标变化, ui、vi、wi在新坐标 (x、y、z) 中表示为:
3.2 三维速度方差 (σ)
计算公式分别如下:
3.3 风向方差
计算公式分别如下:
式中:σθ和σΦ的标准偏差分别表征水平和垂直方向上风向脉动的涨落状况。
污染物的扩散是一种拉格朗日观测, 污染气体的小单位或无惯性粒子严格地随气流运动。用固定风速仪测得的速度 (欧氏) 变化快于测量随着气流运动的气块速度 (拉氏) 。而由于大气扩散参数必须用欧拉观测得到的数据计算, 因而要求建立欧拉和拉格朗日湍流之间的关系[1,2]。Hay和Pasquill已将拉氏相关系数R (ξ) 和R (t) 之间建立了联系, 拉氏相关系数的直接观测由欧氏的时间相关系数的观测来代替, 即关系式为:
3.4 大气扩散参数方程
由三轴风速仪计算拟合大气扩散参数方程, 经量纲分析, Taylor公式可以用下述函数来拟合:
式中:x为下风向距离, f (x) 是普适函数。
3.5 f (x) 函数值的确定
式中:x为扩散距离, L为奥布霍夫长度。
4 个例计算与分析 (以兰州市西固区为例)
以2005年1月兰州市西固区的大气观测资料进行计算。
4.1 不同稳定度下的实测各特征量的平滑平均值
不同稳定度下的实测各特征量σθ、σΦ和V的平滑平均值如表1所示。
4.2 不同稳定度下的扩散参数方程
稳定条件下, Hay和Pasquill方法确定的扩散参数方程如下:
不稳定条件下, Hay和Pasquill方法确定的扩散参数方程如下:
4.3 各个稳定度扩散曲线与P-G、Briggs扩散曲线的比较
P-G扩散曲线法是根据风速、云与太阳辐射状况对大气的稀释扩散能力进行判断, 以确定其稳定度扩散级别 (A、B、C、D、E、F, 共6类) [7]。Briggs是在一系列扩散试验的基础上, 分析曲线在不同距离的特性, 将其结合在一起, 得到一系列扩散曲线, 即Briggs扩散曲线[8]。基于大量的扩散试验资料和理论分析, 得出扩散参数随离源下风向距离 (x) 的变化曲线, 具体如图1、2所示。不论是在稳定条件下还是在不稳定条件下, 用脉动风速资料计算的水平扩散参数 (σy) 和垂直扩散参数 (σz) 均比P-G法和Briggs曲线法大。主要原因是P-G法和Briggs曲线法使用于平原地区的扩散参数, 而兰州市西固区的地形比较复杂, 用脉动资料确定大气扩散参数更能准确地反映当地大气扩散稀释能力。
5 小结
研究表明, 用P-G扩散、Briggs扩散计算复杂下垫面的扩散参数时, 结果偏小。对于地形复杂的地区 (如兰州市西固区) , 研究其扩散能力时, 应尽可能利用脉动观测资料确定大气扩散参数, 以便更为准确地反映当地大气的稀释扩散能力, 为扩散模式提供精确的参数, 真实地反映当地的环境质量状况[9,10], 从而为准确地预测当地天气变化提供科学的理论基础。
摘要：分析2005年1月兰州市西固区脉动风速的测量数据, 通过计算得到水平扩散参数 (σy) 和垂直扩散参数 (σz) , 并与P-G、Briggs扩散曲线得到的扩散参数进行比较。结果表明:在稳定条件或不稳定条件下, 用P-G模式及Briggs方法计算的结果都小于用风速脉动资料确定的结果。
关键词：三轴风速仪,扩散参数,湍流,甘肃兰州
参考文献
[1]牛生杰.沙尘气溶胶微结构及其对降水影响的观测和数值模拟研究[D].南京:南京气象学院, 2004.
[2]王冠.兰州市街道尘埃磁性特征研究[D].兰州:兰州大学, 2008.
[3]陈义珍, 柴发合.欧拉观测的大气扩散参数计算方法探析及应用[J].上海环境科学, 1997, 16 (11) :17-19.
[4]林金彪.用三轴风速仪确定大气扩散参数的方法初探[C]//湖北省气象学会.2007年湖北省气象学会学术年[1]会暨全省青年气象科技工作者学术交流会议学术论文详细文摘汇集.武汉:[出版者不详], 2007.
[5]桑建人, 孙继明, 牛生杰.沙漠地区大气扩散参数特征分析[J].中国沙漠, 2000, 20 (3) :301-304.
[6]文术, 杨柏生.利用湍流特征量的测定值确定扩散参数[J].四川气象, 1995, 15 (3) :26-27.
[7]张镭, 王继红, 任健美, 等.用风脉动资料确定低山丘陵地区大气扩散参数的研究[J].兰州大学学报:自然科学版, 1994, 30 (4) :152-156.
[8]杨洪斌, 张云海, 马雁军.朝阳地区大气扩散参数计算分析[J].环境保护科学, 2006, 32 (1) :5-7.
[9]兰州市清洁能源行动办公室..兰州市清洁能源行动规划[R].兰州, 2002.
扩散参数 篇3
Sh ibor是央行稳步推进利率市场化改革, 提高金融机构自主定价能力, 指导货币市场产品定价, 完善货币政策传导机制而培育的中国货币市场基准利率体系。自2007年央行推出shibor利率后, 在央行进8年的积极培育下, 已基本确立了shibor利率在我国货币市场的基准利率地位, 并逐渐成为传导货币政策、反映市场利率变动的重要指标, 并在市场化利率形成机制中发挥重要作用。同时, 市场上发展和创新了一大批以Shibor为基准金融产品。大量债券交易、存款、贷款、贴现和理财类产品定价逐步与Shibor挂钩。目前, 金融机构市场成员交易层面已有部分金融产品定价与Shibor相挂钩:衍生品市场利率互换、远期利率协议;货币市场同业借款、同业存款、货币互换、理财产品等;债券市场以Shibor为基准浮动利率金融债券、企业债券、企业短期融资券, 以Shibor为基准的票据转贴现、票据回购等;以Shibor为定价及交易基准的金融产品越来越多, 如债券买卖、票据贴现、个人及机构理财产品, 最终, 银行存贷款利率定价也与Shibor相挂钩。
因此, 随着越来越多地的产品与Shibor挂钩, 我们应积极培育我国商业银行金融衍生产品的自我定价能力, 让Shibor真正在我国金融自由化和利率市场化改革进程中起到货币政策利率传导的主导与核心作用。以Shibor为标的的金融衍生产品定价问题成为越来越多学者的研究对象。要做到这一点, 我们首先要研究Sbibor及其本身的利率期限结构, 对其进行利率建模, 从而能够有效刻画出其利率期限结构动态变化的特征, 进而对利率的未来变动进行科学的预测, 这样不管是对以Shibor为标的的金融衍生产品定价, 还是利率市场化条件下以Shibor作为货币市场基准利率风向标的金融风险管理都起着极其重要的作用。
二、shibor利率模型的选择
由于金融市场兴起较早而且比较完善, 国外在利率方面的研究已经发展到了相当高的水平。在利率动态模型方面, 单因素利率模型主要有Merton (1973) 、Vasicek (1997) 、CIR (Cox, Ingersoll和R oss) (1985) 、和CKLS (Ch an, Kar oly i, Longstaff和Sander s) (1992) 等。其中, CKLS模型为不同的利率期限结构建立了一个共同框架利率模型, 单其利率模型只描述了利率漂移项的均值回复特性与利率扩散项的水平效应。Merton (1976) 在标准布朗运动中加入跳跃因素来对金融市场中存在的不连续动态变化进行描述。随后也有国外学者, 如Das (2002) , Chen&Scott (2002) 等人为在模型中加入随机波动率更贴近现实, 包含了更多的利率资产波动率时变特征。
回顾传统的利率模型, 我们知道它们都假设利率服从连续扩散过程。但是大量的实证研究表明, 现实情况中利率经常会出现突然的跳跃和巨幅的震荡, 与传统假设相违背。Das把Vasieck模型扩展到跳跃扩散模型, 并通过实证表明跳跃过程能捕捉连续模型无法解释的利率特征, 认为加入随机波动了更贴近现实。Jarrow等发现跳跃行为能很好地解释“波动率微笑”之谜。Ahn和Thompson认为利率的样本轨迹存在非连续性特征, 由此提出了利率期限结构的跳跃扩散模型。Johannes提出了一个关于跳跃行为的统计检验量, 并发现三个月的短期国债利率存在显著的跳跃行为。因此, 我们可以这样认为:在shibor利率模型中加入跳跃项能够更好地刻画和解释实际中的利率波动情况。
本文借鉴潘璐 (2011) 所用的CKLS-SV-Jump模型:
其中, corr[d Z1, t, d Z2, t]=0, J1, t~N (ψ1, γ12) , J2, t~N (ψ2, γ22) , q1, t, q2, t分别服从频率为λ1*和λ2*的泊松分布, 其中, Pr{q1, t=1}=λ1Δt, Pr{q1, t=0}=1-λ1Δt, Pr{q2, t=1}=λ2Δt, Pr{q1, t=0}=1-λ2Δt。该模型在扩散模型的基础上加入了跳跃项, 其实证研究证明该模型能很好地刻画利率的均值回复、水平效应、跳跃的特征。
三、MCMC参数模拟方法
(一) 马尔可夫链模拟
在介绍马尔可夫模拟之前, 先回顾马尔可夫过程。考虑一个随机过程{Xt}, 假定每个Xt都在空间Θ上取值, 如果h>t当时, 随机过程{Xt}满足P{Xh|Xs, s≤t}=P (Xh|Xs) , 则称过程{Xt}为马尔可夫过程。即马尔可夫过程满足这样的性质:给定Xt的值, Xh (h>0) 的值不依懒于Xs (s<t) 的取值。
接下来, 考虑参数向量为θ和数据为X的推断问题, 其中θ∈Θ。为了做出推断, 我们要知道分布P (θ|X) 。马尔可夫链模拟的思想是在Θ上模拟一个马尔可夫过程, 这个过程收敛于平稳转移分布P (θ|X) 。模拟的关键是构造一个具有指定的平稳转移分布P (θ|X) 的马尔可夫过程, 并且充分长地运行这个模拟, 使得过程当前值的分布与平稳转移分布足够接近。对给定的P (θ|X) , 可以证明能够构造许多具有所需性质的马尔可夫链。这种利用马尔可夫链模拟来得到分布P (θ|X) 的方法称为MCMC方法。
MCMC方法构建的基础是Cliffor d-H ammer sley理论, 该理论认为基于参数的后验信息P (θ|x, y) 和状态变量的后验信息P (x|θ, y) , 即可唯一确定联合后验分布。而MCMC方法要求从联合后验分布中中抽取样本, 因此我们根据Clifford-Hammersley理论, 可以分别从参数的后验信息p (Θ|x, y) 和状态变量的后验信息p (x|Θ, y) 中抽取样本, 通常参数和状态变量其各自的后验分布都是比较容易得到的。由贝叶斯原理, 参数的后验分布可以写成参数的先验分布、似然函数和某个常数的乘积形式。MCMC方法的核心就在于将高维度的联合分布p (Θ, x|y) 分解为低维度的p (x|Θ, y) 和p (Θ, x|y) , 然后再从这两个分布中分别取样。
由于运用MCMC参数估计方法前需要对模型进行离散化, 因此, 我们接下来对欧拉离散法进行了介绍。
(二) 欧拉离散法
接下来对CKLS-SV-Jump模型进行离散化, 可以得到:
进一步可以得到下式:
最终我们整理成下式:
因此, 我们可以将rt+Δ基于rt的条件分布和θt+Δ基于θt的条件分布表示为:
其中, {εt}~iid.N (0, 1) ;{Jt}~iid.N (0, 1) 。
(三) Gibbs抽样
Gibbs抽样方法是MCMC方法中应用最广泛的一种抽样方法, 是用来获取一系列近似等于指定的多维概率分布的观察样本的方法。Gibbs抽样是单元素Metropolis-Hastings算法, 在Gibbs抽样中, 可能难以抽取, 而Metropolis方法具有很大的灵活性, 它可取为易于抽取的分布。
在足够的燃烧期后, Gibbs序列才能收敛到一个独立于初始值的平稳分布, 能否得到这个平稳分布是Gibbs抽样的关键。由于燃烧期内产生的参数值是不平稳的, 因此要放弃燃烧期内产生的样本。如果我们进行了n次抽样, 燃烧m次, 那么就将序列X (m+1) , X (m+2) , …, X (n) 作为马尔可夫链的实现值, 进行后续的参数估计和统计描述工作。
在Gibbs抽样中, 一个比较困难的事情是判断上述迭代过程是否收敛到均衡状态。不少作者都建议用的某些特征来判断。最常用的是观察遍历均值是否收敛来判断。比如在由Gibbs抽样得到的链中, 每隔一段距离计算一次参数的遍历均值, 未使用来计算平均值的变量近似独立, 通常可每隔一段取一个样本, 当这样得到的均值稳定后, 可认为Gibbs抽样收敛。
(四) Metropolis-Hastings算法
1953年, Metr opolis等人提出了一种构造转移核的方法, 该方法此后被Hastings加以推广, 形成了Metropolis-Hastings方法。具体思路如下:
MH算法的思想是构造一个以目标分布π (x) 为极限分布的Mar k ov链。任选一个不可约转移核概率p (·, ·) 以及一个函数a (·, ·) , 0<a (·, ·) ≤1, 对任一组合 (x, x') (x≠x') , 定义
则p (x, x') 形成一个转移核。
假设卡尔可夫链在时刻t处于状态x, 即X (t) =x, 则首先由q (·|x) 产生一个潜在的转移x→x', 然后根据概率a (x, x') 决定是否转移。即在潜在转移点x'找到后, 以概率a (x, x') 接受x'作为卡尔可夫链在t+1时刻的状态值, 以1-a (x, x') 拒绝转移到x', 即马尔可夫链在下一时刻的状态值仍为x。也就是说, 如果置X (t+1) =x', 我们就称之为“接受建议”;否则就称为“拒绝建议”。在实际计算中以概率a (x, x') 接受x', 可以这样实现:从U (0, 1) 产生一个随机数u, 如果u≤a (x, x') , 则置X (t+1) =x';否则置X (t+1) =x。
因此, 在有了x'后, 我们可以从[0, 1]上均匀分布抽一个随机数u, 有
这里, 分布q (·|x) 称为建议分布。由于我们的目标是使后验分布成为平稳分布, 所以在有了q (·|·) 后, 我们要选择一个a (·, ·) 使相应的p (x, x') 以π (x) 为其平稳分布, 我们通常选择
因此也有
该式产生的马尔可夫链是可逆的, 即π (x) p (x, x') =π (x') p (x', x) , 且π (x') 是该式确定的马尔可夫链的平稳分布。
(五) MH随机游走采样法
如果M-H算法中的建议函数q (x, x') 满足1、对称性q (x, x') =q (x', x) ;2、q (x, x') 仅与x'-x有关, 即q (x, x') =q (|x'-x|) , 那么算法就演变为通常意义下的随机游动采样法。最常见的一种随机游动采样法以正态分布, 即q (x'-x) =φ (x'-x) 为建议函数。
这里是任意一种多维正态分布的密度函数, 也就是说x'~N (x, Σ) 。其中, x为当前状态值, Σ是任意的正定矩阵。在实际计算中, 由于要求Σ是正定矩阵, 所以常取为Σ=σI, σ是一个参数。这时候, 如何选取建议函数的方差σ是影响算法效率的主要因素。对于一维情况, 就是要选择一个合适的σ值, 而对于多维情况, 我们需要确定给一个合适的正定方差矩阵Σ。一般认为在模拟多维正态随机数时, 如果调整σ的值使得在整个模拟过程中建议被接受的比例为20%左右, 则将得到比较好的模拟效果。
随机游动采样法是M-H算法中较为常用的一种形式, 们已经看到对于随机游动采样法而言, 如何选择建议函数的方差成为影响算法效率的主要因素, 对于一维情况, 就是要选择一个合适的σ值;而对于d>1的多维情况, 我们需要一个合适的正定方差矩阵Σ。在选择σ的过程中, 根据关于采用多维正态作为提议函数的随机游动采样法的一些理论结果, 一般认为如果目标分布是d维的且各分量之间独立, 则调整接受概率到0.234左右将最大限度地优化随机游动采样法的效率, 并且接受概率0.234与目标分布的具体形式无关。
基于这样的想法, 我们尝试构造一种自适应算法来寻找合适的提议函数的方差。大致的思想是让接受比率落入一个可以接受的范围内, 比如区间[a-ε, a+ε]。以一维的情况为例, 先给定一个初值σ0, 用这个值产生一条链, 来估计接受比率P0。若P0在给定的区间[a-ε, a+ε]中, 则符合要求, 算法停止;反之, 则以为中心做一次随机游动, 得到一个σ1, 再计算P1是否满足要求, 以此类推, 直至找到合适的σ值。具体实现步骤如下:
1) 给定σ一个初值σ0, 置n=0;
2) 以N (Xi-1, σn) 为建议函数, 用随机游动采样法产生一条长度为m的马尔可夫链;
3) 计算第2步中产生的马尔可夫链接受新状态的比率Pn;
5) 如果n>2且|Pn-a|<|Pn-1-a|, 则置σ=σn-1;
6) 从N (σn, σ) 中产生一个随机数, 记为σn+1, 置n=n+1, 转入第2步。
其中, 步骤5的判断是为了保证接受比率Pn要优于Pn-1, 即σn至少不会比σn-1坏;第6步中σ的选取没有固定的规律可循, 主要是看对初始σ0的估计是否接近于合适的σ值。若估计σ0与σ较为接近, 则σ可选择得小一些;若估计σ0与σ相差较远, 则σ可选择得大一些;对于d>1的多维情况, 如果要确定一个任意的正定矩阵∑, 则要确定d2个参量, 显得有些麻烦。更多情形下, 我们取∑=σI。这样一来要确定的参数就只有一个, 同时上述算法也可以用于多维情形, 只需在6步中“以N (Xi-1, σn) 为建议函数”改为“以N (Xi-1, fnJ) 为建议函数”即可。
四、模型参数估计结果
(一) 模型的先验分布
贝叶斯学派认为, 在进行观察以获得样本之前, 人们对θ也会有一些知识。因为是在试验观察之前, 故称之为先验知识。因此, 贝叶斯派认为, 应该把θ看作是随机变量。θ的分布函数记为H (θ) , θ的密度函数记为h (θ) , 分别称为先验分布函数和先验密度函数, 两者合称为先验分布。在MCMC方法中, 待估参数的先验分布的设定将会对该方法运行效率产生很大影响。本文为了方便起见, 参照了Siddhartha Chib, Federico Nardari和Neil Shephard (2002) 和潘璐 (2011) 设定的参数先验分布。根据共轭理论, 参数有相同的后验分布和先验分布。
(二) 抽样方法
波动率状态变量的完全联合后验分布可以表示为:
其中,
所以, 我们得到:
因为上述分布是不可识别的, 所以我们无法使用Gibbs抽样方法, 而需要使用Metropolis抽样方法进行抽样。而且由于它的复杂形式, 在没有采取近似处理的情况下, 用但不Metropolis方法是无法直接进行抽样的。因此, 根据Jacquier, Polson和Rossi (1994) 的研究成果, 我们采用随机游走Metropolis算法。因为我们可以对满条件分布p (r|rt-Δ, rt+Δ, Θ, r) , 尤其是它的尾部进行很好的近似。在上式中, 第一项是对数正态分布, 第二项是倒Gamma分布, 所以满条件分布p (r|rt-Δ, rt+Δ, Θ, r) 的建议分布可以很好地由某个合适的倒Gamma分布来近似。我们设这个建议分布为q (θt) , 真正的条件密度函数为π (rt) =p (r|rt-Δ, rt+Δ, Θ, r) 。
具体地, 我们将CKLS-SV-Jump模型的MCMC算法分解成以下步骤:
……
第十二步:利用Metropolis-Hastings自适应算法从p (r|rt-Δ, rt+Δ, Θ, r) 中抽取θt。
其中第十二步又可以分为以下几步:
(1) 给θ设定一个初始值θ0, 置t=0;
(2) 以N (rt-Δ, θt) 为建议分布, 用随机移动M-H算法产生一条长度为的马尔可夫链;
到此为止, 完成第一轮迭代过程。接下来, 利用新的参数作为初始值, 重复以上步骤, 得到更新参数, 重复前面的步骤次, 直至各参数和波动率潜在变量θt的模拟轨迹收敛为止, 得到一系列的随机抽取。我们可以认为 (αm*, βm*, ρm, ψ1, m, ψ2, m, γ21m, γ22m, λ1, m, λ2, m, υm, φm, τm2, θm) 渐进等价于参数估计值, 然后把收敛时的模拟值作为对各参数和潜在变量θt的估计值。
(三) 数据选择
扩散参数 篇4
短期利率不仅是刻画利率期限结构动态过程的重要状态变量, 同时也是固定收益证券与衍生品定价的重要依据, 因此有许多研究分别提出不同的连续时间模型来刻画短期利率的动态过程。短期利率rt的单变量连续时间模型一般设定为: drt=μ (rt) dt+σ (rt) dWt (1)
其中μ (rt) 和σ (rt) 分别是漂移函数和扩散函数, Wt是布朗运动 (Brownian motion) , t∈[0, T]。各种单变量连续时间模型之间的差异主要在于漂移函数和扩散函数的设定不同, 主要的模型有Vasicek[16], Richard[13], Brennan和Schwartz[4], Cox, Ingersoll和Ross[7], Chan, Karolyi, Longstaff和Sauders[5], Duffie和Kan[8]等等。Aït-Sahalia[2]给出了一个更加一般化的模型, 即:
drt= (α0+α1rt+α2rundefined+α3rundefined) dt+ (β0+β1rt+β2rundefined) -1/2dWt (2)
上述各种模型都是通过对参数α0, α1, α2, α3, β0, β1, β2和β3进行不同设定而得到的[2]通过对比在不同参数设定下利率数据的边际密度和没有参数设定下原始数据的边际密度, 发现所有的参数模型都与实际数据不符, 因而否定了参数模型正确性。Backus、Foresi和Zin[3]的研究则表明利率模型的错误设定会导致严重的定价误差。
对利率模型的参数施加较少约束的非参数模型近来有了较多的发展和应用。在国外的研究中, Aït-Sahalia[1]使用7天欧元存款利率的日观测数据, 先用参数方法估计漂移函数, 再结合非参数密度估计量和漂移函数的参数估计量从而获得扩散函数的非参数估计量。研究结果表明扩散函数的非参数估计量与参数估计量相比具有较大的波动性, 并且随着利率水平的提高其波动幅度也相应增大, 因而在整体上是非线性的。Jiang、Knight[12]在Florens-Zmirou[9]等研究的基础上, 提出了漂移函数和扩散函数的非参数核估计量, 并将该方法应用于加拿大3个月国债利率的日观测数据。研究结果表明漂移函数在利率水平较高或较低时均具有均值回复特性, 而在中等利率水平处则比较平缓, 不具有均值回复特性。由于扩散函数的波动幅度较大, 整体上也是非线性的且不具有均值回复特性。Stanton[15]使用泰勒级数展开式和无穷小生成元 (infinitesimal generator) , 得到基于不同阶数的漂移函数和扩散函数的近似估计, 并且在引入核函数后得到两者的非参数估计量。通过对美国3个月国债收益率日观测数据的研究, 发现漂移函数在较低和中等的利率水平处比较平缓, 而当利率水平较高时则具有较强的均值回复特性, 因而整体上是非线性的。扩散函数则随着利率水平的上升而逐渐增加, 具有非线性特性。
在国内的研究中, 李和金、郑兴山和李湛[19]使用Florens-Zmirou的方法对国债回购利率进行非参数估计, 研究结果表明漂移函数在中等利率水平处比较平缓, 而在利率水平较低和较高时均表现出较强的均值回复特性。扩散函数则随着利率水平的提高而增加, 而且是非线性的。宋永安和陆立强[20]将Florens-Zmirou的非参数方法, 应用于上海证券交易所1个月回购利率R028每日收盘数据, 其非参数估计结果表明漂移函数整体上是非线性的, 在利率水平较高时具有较强的均值回复特性, 而在其它利率水平下则比较平稳。扩散函数的波动幅度随着利率水平的提高而增加, 也呈现出非线性的特性。周荣喜等[21]使用Florens-Zmirou的非参数方法对上海证券交易所回购国债GC001数据进行了研究, 并且在估计中分别使用了抛物线核函数与高斯核函数。研究结果表明高斯核函数的估计效果较好, 其估计结果表明漂移函数整体上十分平缓, 在利率水平较高时略微具有波动性;扩散函数则随着利率水平的提高而增加, 而且是非线性的。胡瑾瑾和陈淼垚[17]使用Stanton[15]和Jiang[11]等研究的方法, 计算漂移函数和扩散函数的非参数估计量, 通过对上海证券交易所债券市场7天回购利率R007的研究, 发现漂移函数在利率水平较低和中等利率水平处十分平缓, 而在较高的利率水平处呈现较强的均值回复特性, 因而在整体上是非线性的。由于扩散函数则在利率水平较低和中等利率水平处估计值较低, 而在利率水平较高时的估计值较高, 因而在整体上也是非线性的。
Aït-Sahalia[2]的非参数检验认为参数模型之所以被拒绝, 主要是因为它们对漂移函数的错误设定, 即应该将漂移函数设定为非线性的而不是线性的。这一观点可以用前面的文献综述加以佐证, 即绝大部分结论都认为漂移函数是非线性的。但是这一结论也受到了其它研究的质疑, 比如Chapman和Pearson[6]认为当真实的漂移函数是线性的时, Aït-Sahalia[2]和Stanton[15]的非参数估计量也会倾向于得到非线性的漂移函数, 因而漂移函数具有非线性特征的结论是不稳健的。有鉴于此, Sam和Jiang[14]借鉴Hjort和Glad[10]的研究思路, 通过将参数试点估计量 (pilot estimator) 与非参数“修正因子”相结合来改进函数的估计。为了估计美国3个月国债利率日观测数据的动态模型, 还使用了期限分别为6个月、1年的国债数据和期限分别为3、5、10年的中期国债数据。使用新方法所得到的结果与使用Stanton方法所得到的结果相比, 两者的漂移函数在较低和中等的利率水平处均比较平缓, 并且在利率水平较高时具有均值回复特性, 但是前者的回复程度要明显减弱许多, 因而得出漂移函数在利率水平较高时只是略微具有均值回复特性的结论。本文采用Stanton[15]与Sam和Jiang[14]两种非参数估计方法, 对SHIBOR市场利率扩散过程的漂移函数进行估计, 通过对不同结果的比较来研究漂移函数的特性以及不同期限利率之间的相互影响。
二、利率模型的非参数估计
为了避免对漂移函数和扩散函数的错误设定, Stanton[15]提出了一种通过使用离散的观测数据对两个函数进行近似估计的有效方法。考虑式 (1) 中的扩散过程rt, 将条件期望Et[f (rt+Δ, t, Δ) ] 用泰勒级数展开为:undefined
其中A是{rt}的无穷小生成元[15]选取fμ (r, t) ≡r, fσ (r, t) ≡ (r-rt) 2, 分别用于估计两个函数, 所得到的一阶近似估计分别为:
undefined
利用核估计方法, 得到式 (5) 的非参数估计为:
undefined
其中K (z) =2 (π) -1/2e- (1/2) z2为正态核函数, h为窗宽 (window width) 。
Sam和Jiang[14]的非参数方法与现有参数方法和其它非参数方法的关键区别, 在于它使用的是由多个期限的利率序列构成的面板数据, 其中第一个利率的漂移函数是需要估计的, 而其它期限的利率序列则用作辅助序列参与估计。该方法借鉴了Hjort和Glad[10]的研究思路, 使用参数试点估计量和非参数“修正因子”来改进函数的估计。其优点在于, 当试点估计量能够较好地近似条件均值时, “修正因子”是一个平滑的函数, 因而更容易用非参数方法进行估计。该方法分两步实施:第一步, 同时使用各个期限的利率用以估计漂移函数的非参数混合估计量, 由于不同期限利率的漂移函数一般来说是不同的, 因而该估计量不是最优的;第二步, 使用非参数“修正因子”来纠正混合估计量中存在的偏误。
定义μp (r) =∑undefinedμj (r) wj (r) , 其中wj (r) = (pj (r) ) / (∑undefinedpj (r) ) 。pj (r) 是第j个利率序列的密度函数, j=1, …, J。μp (r) 就是要用各个期限的利率数据来估计的漂移函数的混合估计量, 利用它可以得到漂移函数的估计, 即:undefined
在式 (8) 中, c (r) 为“修正因子”, 且μp (r) ≠0。c (r) 的非参数估计undefined (r) 可表示为:
undefined
在实际应用中, 由于μp (r) 未知, 因此需要用它的一致估计量来代替, 即:
undefined
结合式 (9) 与式 (10) 就可以得到漂移函数μ1 (r) 的非参数估计量undefined1 (r) 。当使用日观测数据时, 式 (9) 和式 (10) 中的δ取值为1。
三、实证分析
本文使用上海银行间拆放利率 (SHIBOR) 数据作为研究对象, 该数据在国内的研究中已逐渐得到了广泛的应用, 同时也具有重要的参考价值。数据的样本期间从2006年10月8日到2011年7月29日, 总共有1 206个日观测值, 选择的利率期限分别为1周、2周、1个月、3个月、6个月、9个月和1年。考虑到国内外还没有相关研究使用期限为隔夜的利率数据, 因而本文选择估计期限为1周的利率, 而其它期限更长的利率则作为辅助利率, 期限为1周的利率及其一阶差分分别由图1和图2给出:
表1则给出了各个期限利率的描述统计量:
漂移函数的非参数估计分别使用Stanton与Sam&Jiang两种方法, 最优窗宽按照公式undefined进行计算, 其中s为利率序列的标准差, T为样本观测个数。相应的估计结果分别由图3与图4给出:
从图3中的结果可以看到漂移函数在中等利率水平以下的范围内十分平稳, 几乎为0。在利率水平较高时则具有均值回复特性, 此后又趋于平稳。该结果与国内外相关研究的发现较为接近。图4中漂移函数的非参数估计结果是按照Sam&Jiang方法估计的, 可以看出与图3有明显的差异。在利率水平较低时, 漂移函数十分平稳, 而在中等利率水平处则表现出波动特征, 且波动幅度较大。在利率水平较高时具有均值回复特性, 并且均值回复与平稳交替出现。
如前文所述, Sam&Jiang方法所构造的非参数估计量通过引入辅助利率来纠正以往方法中存在的偏误, 而图4中的结果却与该方法的初衷相违。考虑到其它相关研究中并没有发现类似的结果, 因而认为图中的结果可能与辅助利率的选择有关。在SHIBOR数据中, 期限为1个月及以下的利率为短端利率, 期限为3个月至1年的利率为长端利率。短端利率的波动幅度较大, 表现为波动趋势;长端利率的波动幅度较小, 表现为漂移趋势。当使用长端利率作为短端利率的辅助序列时, 两者在趋势上的差异对短端利率漂移函数的估计可能会造成一定的影响。长端利率对1周利率的漂移函数估计的影响可以通过比较图4与图5做初步判断:
图5与图4的区别在于辅助利率的选择不同, 即只选择期限为2周和1个月的短端利率作为辅助利率, 而没有将长端利率包括进来。从图5可以看出在中等利率水平处, 虽然漂移函数的波动幅度有所减小, 但是仍然具有明显的波动特征。可见长端利率只是加剧了波动幅度, 并不是造成估计结果之间差异的主要原因③。此外, 从表1所列出的描述统计量中可以发现, 即使短端利率之间也在偏度、峰度和最大值等描述统计量上存在明显差异。短端利率之间的差异或许是导致漂移函数波动特征的原因之一, 并且其影响主要体现在利率水平中等和较高时, 然而这一判断的正确性还有待进一步的研究与分析④。
为了进一步考察短端利率与长端利率趋势上的差异对估计结果的影响, 本文还对期限为1个月的利率分别按照两种非参数方法进行了估计, 其它期限更长的利率则作为辅助利率。从表1可知该期限的利率是所有短端利率中与长端利率在统计特性上最为接近的, 并且从走势图 (本文未给出) 上来看, 该利率与其它短端利率相同的是它具有波动趋势, 不同的是它与长端利率具有相似的线性趋势。估计结果分别在图6与图7中给出, 图7中的结果与图6的总体上差异不大, 只是在中等利率水平处略微具有波动特征, 但是波动幅度不大, 同时均值回复强度要大一些。因此, 长端利率对1个月利率的漂移函数估计结果的影响与图4中的相比要明显减弱许多, 这也跟该利率与长端利率在统计特性上最为接近的观察相一致。可见长端利率对短端利率的漂移函数估计结果的影响取决于两者在统计特性上的差异程度。当统计特性差异较大时长端利率的影响就大, 并且主要体现在中等利率水平处;当统计特性差异较小时长端利率的影响要小很多。
四、结论
本文使用SHIBOR市场利率数据, 结合Stanton与Sam&Jiang两种非参数方法, 对利率扩散过程的漂移函数进行了估计。对于期限为1周的利率, 在两种方法下所得到的漂移函数的非参数估计结果存在显著差异, 其中Stanton方法所得到的结果与以往研究较为接近, 而Sam&Jiang方法所得到的结果表明在中等利率水平处漂移函数具有较强的波动特征, 且在利率水平较高时的均值回复特性上也存在差异。这样的结果不但有违Sam&Jiang方法的初衷, 而且在以往研究中也没有与之相似的情形, 因而认为导致这一现象的原因在于辅助利率的选择。通过去掉辅助序列中的长端利率, 发现长端利率只是加剧了中等利率水平处漂移函数的波动幅度, 而短端利率之间在统计特性上的差异或许是导致这一现象的主要原因。为了进一步分析长端利率对短端利率漂移函数估计值的影响, 本文还估计了与长端利率在统计特性上最为接近的1个月利率的漂移函数, 发现两种方法所得到的估计结果差异很小, 表明长端利率对短端利率漂移函数估计值的影响取决于两者在统计特性上的差异程度, 统计特性差异越大影响就越大, 反之亦然。