第20鞋第14期 6l20No.14 2008年7月 李卫国，等：基于MCMC模拟的相关系数平稳序列模型及其应用 jul,2008 态随机过程，则称{x),reT}为相关系数平稳正态过程。记 从而得到相关系数平稳AR()序列模型的似然函数为： 相关系数平稳正态过程的均值、方差函数分别为： L0=(无，，5x色b,n心) E(x()}=μ(a,) (2) =4o4,∑(-4ai-八TP,)(I1) var(x(t)]=o2(b.t) (3) 式中a=(a,4a,,b=(你，…b,)了均为待定参数。函数 式中的x:a,)表示均值为a,标准差为B的正态分布密 4(,σ()具体形式可用文献4提出的时间序列均值方差 度函数。设模型的各个参数的先验分布相互独立，则有： 函数的确定方法以及结合所研究问题的物理性质、经验公式 p()=p(a)p(b)p(n)p(oi) (12) 等求得。 根据贝叶斯公式便有：参数0的联合后验分布密度函数 12相关系数平稳序列 f(): f(a,b,n,o2,…,x)eL(0)p(0) (13) 定义3：，设{x(),teT}是相关系数平稳过程，记 x=(,,xy为它的一个样本，其中无=，)，=6+山， 2.2MCMC方法、Gibbs抽样及Metropolis-Hasting 。,△山分别是采样初始时刻和采样间隔。设y=(%,2,…y,)y 算法 是传统的AR(p以、MA(q)或ARMA(P,q)序列，若 MCMC方法5-刀的基本思想是构造一条具有指定平稳分 y=-4a】 (4) 布π(X)的马尔可夫链，即它的转移分布收敛到的后验分 1(b,4) 布。然后充分长地运行链，到链上取值的分布与其平稳分布 E{x)j=4(a,) (5) 足够接近时，把来自链上的样本作为来自(X)的样本， var(x(t))=l2(b,t) (6) 基于这些样本进行各种统计推断。例如求均值、方差等。不 上式中a=(a,a…a,),b=(,hb,,哈均为待定参 同的MCMC方法主要是构造Markov链的方法不同，常用 数，则称x=(:,,…x了为相关系数平稳ARp)序列、MA(q 来构造链的方法有Gibbs算法和Metropolis-Hasting算法。 序列或ARMA(p,q序列，此时也称(xt),t∈T]为相关系数 Gibbs抽样的思想是利用条件分布族来得到以()为 ARp、MA(q)或(ARMA P,q)过程。这样一般的相关系数 平稳分布的马尔可夫链。 平稳ARMA(p,q)序列模型为： 设π=π(8，…，8.) (14) 本a2-2n西-a】+4-2⑦ 记0=(g,,8"),8四=(g,…,8"): Ib,)台 1(b,-) 再记第1个至第(k-1)个分量固定为2，…，，同时将第 (7式中E{x(4)}=4(a,4),var{x)}=o2(b,), (k+1)至第m个分量固定为，…，的条件下，第k个分 ELELE:]=026x, 量8的条件分布为： 其中当k=i时6=1，当k≠i时，6=0，乃…np,…入为 π(88…，，，…，)，定义转移概率： 常数。当q=0时，（⑤）式表示的就是相关系数平稳AR(p)序列 Pagn兰p(0-→g2) 模型，当p-0时，(5)式是相关系数平稳MA(q)序列模型。 下文主要研究相关系数平稳AR(P)序列模型参数的贝叶斯 =1xaa,8m,8.,,8) (15) 估计。 可以证明{P}为一个转移矩阵，π()是它的平稳分 2相关系数平稳AR(1)序列模型的参数估计 布。这里并不需要知道π的具体表达式，只需知道固定其它 一切分量的条件下，余下一个分量的条件分布。 21模型的贝叶斯分析 如果条件分布不是来自一个标准分布族，则估计参数就 湍要用Metropolis-Hasting算法。该算法最早由Metropolis 设，e,是高斯过程，给定样本观测值(0，，，…x,)条 于1953年提出，1970年Hasting对它加以推广。由于分布 件下，记0=(a,b,n02y,a=(a,a…a,)b=(6,…b), 族未知，抽样无法直接从原来分布族进行，此时，从另一个 刀=(，…n,)厂，则相关系数平稳AR(p)序列模型中各个量 被称为建议分布的已知分布开始。假定建议分布为T(x,y), 的关系可表示为： 当前状态为x,其算法为： (xo0)-N(ao,bio2) (8) ①从建议分布T(x,y)产生y: (,,-p,)- (9) ②②由均匀分布U(0,)抽取u: N吃-以ai-Arb.D ③③计算比率(x,y),其中 (a,b,n.o2)-P(a,b,n.o2) (10) r(y)=min(l)T(. (16) π(x)T(xy) 式中x2表示给定元条件下的x,N(μ，o)表示均值为4， ④x+0= [y,if usr(x,y) (17) 方差为σ2的正态分布函数，P(表示关于0的分布函数。 x,other ·3649· 万方数据第20卷第14期 2008年7月 李卫国，等：基于McMc模拟的相关系数平稳序列模型及其应用 VU．20No．14 Jul．，2008 态随机过程，则称{工(f)，f∈r}为相关系数平稳正态过程。记 相关系数平稳正态过程的均值、方差函数分别为： E{工(f)}=∥(a，f) (2) var{工(f)l=仃2(b，f) (3) 式中a=‰，国…珥)7，b=(％，岛…仇)r均为待定参数。函数 ∥(·)，仃2(·)具体形式可用文献【4】提出的时间序列均值方差 函数的确定方法以及结合所研究问题的物理性质、经验公式 等求得。 1．2相关系数平稳序列 定义3：．设{工(f)，f∈丁)是相关系数平稳过程，记 工=“，死，…而)7为它的一个样本，其中玉=石(‘)，^=岛+池， f0，△f分别是采样初始时刻和采样间隔。设y=(为，弛，…儿)r 是传统的AR(p)、MA(q)或ARMA(p，q)序列，若 yj：譬娑掣 一 (4) ，(b，‘) ～ 耳z(f))=∥(a，f) (5) var{工(f)}=研，2(b，f) (6) 上式中a=(嘞，砚…珥)7’b=(％，籼…眈)r，D吾均为待定参 数，则称石=(礼花，…而)7为相关系数平稳AR(p)序列、MA(q) 序列或ARMA(p，q)序列，此时也称{工(f)，f∈7’)为相关系数 AR(p)、MA(q)或(ARMA p，q)过程。这样一般的相关系数 平稳ARMA(p，q)序列模型为： 铧：圭吼堕掣+&一窆概。(7) J(b，ft) 智” ，(b，“) ‘智一’ ～ (7)式中E{j【(“)}=∥(a，“)， vaf{工(“)}=盯2，2(b，“)， E【岛岛】=盯2最， 。 其中当七=f时况=l，当七≠i时，况=o，研…仉，丑…五为 常数。当q．0时，(5)式表示的就是相关系数平稳AR(p)序列 模型，当p=o时，(5)式是相关系数平稳MA(q)序列模型。 下文主要研究相关系数平稳AR(p)序列模型参数的贝叶斯 估计。 2相关系数平稳AR(1)序列模型的参数估计 2．1模型的贝叶斯分析 设而，局是高斯过程，给定样本观测值(知，再，恐，…而)条 件下，记0=(a’b，lI，盯2)7，a=(口o，口1…口，)7，b=(％，岛…阢)’， 玎=(砀，…仉)7，则相关系数平稳A尺(p)序列模型中各个量 的关系可表示为： (知fo)～Ⅳ(珊，蹭cr2) (8) (再I再一I'．一，再一p，o)一 (9) P Ⅳ(∑叩，(再一，一∥(a，f一_『))，盯2，2(b，f) 卢l (a，b，q，仃2)一P(a，b，q，盯2) (10) 从而得到相关系数平稳AR(p)序列模型的似然函数为： 以0)=以知，再，·-j矗khqD2) =奴如；匈，如)．n豇玉；∑仿(再叫一“a'l一脚，D2，2(b’f))(11) 式中的纵工；口，∥)表示均值为口，标准差为∥的正态分布密 度函数。设模型的各个参数的先验分布相互独立，则有： p(0)=P(a)p(b)p(11)p(口2) (12) 根据贝叶斯公式便有：参数O的联合后验分布密度函数 ，(·)： ，(a，b，'I，cr2 l葡，而…，五。)氍L(o)p(o) (13) 2．2 McMC方法、Gibbs抽样及MetmpoIis．H嬲ting 算法 MCMC方法‘"1的基本思想是构造一条具有指定平稳分 布石(臼Ix)的马尔可夫链，即它的转移分布收敛到的后验分 布。然后充分长地运行链，到链上取值的分布与其平稳分布 足够接近时，把来自链上的样本作为来自石(口Ix)的样本， 基于这些样本进行各种统计推断。例如求均值、方差等。不 同的McMC方法主要是构造Markov链的方法不同，常用 来构造链的方法有Gibbs算法和Me仃opolis．H越ting算法。 Gibbs抽样的思想是利用条件分布族来得到以万(口)为 平稳分布的马尔可夫链。 设石=万(岛，·一，以) (14) 记伊1)=(研”，…，碟’)，伊2)=(岛。’，…，簖’)； 再记第1个至第(七一1)个分量固定为研舶，…，碰2，，同时将第 (七+1)至第m个分量固定为绣¨l，…，创’的条件下，第k个分 量倪的条件分布为： 万(晚l纠”，…，联二，戗?1’．”，酬’)，定义转移概率： p掣‘w2)皇p(伊1’—》目‘2’) =兀万(岛I岛‘孙，…，B一1‘21，皖+1m，…，以‘1’) (15) 可以证明{办咐：，}为一个转移矩阵，石(∞是它的平稳分 布。这里并不需要知道石的具体表达式，只需知道固定其它 一切分量的条件下，余下一个分量的条件分布。 如果条件分布不是来自一个标准分布族，则估计参数就 需要用Me昀polis—Has“g算法。该算法最早由Metmpolis 于1953年提出，1970年Hasting对它加以推广。由于分布 族未知，抽样无法直接从原来分布族进行，此时，从另一个 被称为建议分布的已知分布开始。假定建议分布为r(工，y)， 当前状态为石“)，其算法为： ①从建议分布丁(五y)产生y； ②②由均匀分布u(o，1)抽取“； ③③计算比率，(工，y)，其中 啾朋刊n{1，嬲器)； (16) 雾塞￡辚黧篡篇茄霈’ ④∥枷镀乞∥∥”∥’ ∽， 方差为仃2的正态分布函数，P(·)表示关于0的分布函数。 【一”， Dt妇r 、’ ·3649· 万方数据