第二十四章时间序列模型

第二十四 时间序列模的概率分布与时间t无关，则称该序列为严格的（狭义的）平稳时间序列。如果序列的一、二阶矩存在，而且对任意时刻t满足：协方差为时间间隔的函数。 通常用Tt表示长期趋势项，St表示季节变动趋势项，Ct表示循环变动趋势项，Rt

TtStCtTtStCtTtStStTtCtyE(R0E(R22 tM(1)t

N(ytyt

LytN11(

1(L )

)M(1) (y

t

t

t

t ytN M(1)

1(

L

)，tN,N1,L,T t

tN(yˆy(yˆy2 tN T最近N期序列值的平均值作为未来各期的预测结果。一般N取值范围：5N200。当历史序列的基本趋势变化不大且序列中随动成分较多时，N的N的取值应小一些。在有确定的季节变动周期的资料中，移动N值的一个有效方法是，比较若干模型的预测误11月～1111企业销售收月份123456月份789解：分别取N4N5yˆyˆ

ytyt1yt2yt3，t4,5,L,11ytyt1yt2yt3yt4，t5( y2tt( y2tt1S 1

11

) yˆ yˆy2t 11销售收入为993.6。 for%由于n的取值不同，yhat的长度不一致，下面使用了细胞数组s(i)=sqrt(mean((y(n(i)+1:m)-yhat{i}(1:end-2.2移动平均在简单移动平均中，每期数据在求平均时的作用是等同的。但是，每期数据设时间序列为y1,y2,L,yt,L；移动平均 w1ytw2y2LwNytN1，t w1w2L即以第t期移动平均数作为第t1期的预测值

例2 我国1979～1988年原煤产量如表2所示试用移动平均法预测1989年表2我国原煤产量统计数据及移动平均预测值年三年移动平均预测相对误差w13,w22,w31

3yt2yt1yt32计算三年移动平均预测值，其结果列于表2中。1989年我国原煤产量的预测

39.829.288.94619826.666.235

y100%y

)100%1

y=[6.35 fori=1:m-n+1在移动平均法中，wt的选择，同样具有一定的经验性。一般的原则是：近期 1(Mt(1)

yt1LytN1M(2)1(M(1)LM )M(2)1(M(1)M(1) tN t tyˆtmatbtm，m atyt1ytbtyt2yt2bt…ytN1yt(NM(1)ytyt1LytN1yt(ytbt)L[yt(N1)bt Nyt[12L(N1)]btyyM(1)N1

N1

M(1)N1 y M(1)M(1) t t M(1)M(2)N1 a2M(1)M 2

(

31965～1985319861987年的发3我国发电量及一、二次移动平均值t123456789解由散点图1可以看出，发电总量基本呈直线上升趋势，可用趋势移动平均法1原始数据散点M(1)3461.2，M(2) （12 2M(1)M(2) 2 (M(1)M(2)) 6 于是，得t21ˆ22ˆ23load yhat1(i)=sum(y(i:i+n-fori=1:m2-n+1yhat2(i)=sum(yhat1(i:i+n- 1 ；而NN1 N设时间序列为y1,y2,L,yt,L， 系数，01，一次指数平 S(1)y(1)S(1)S(1)(yS(1) t t tM(1)M(1)ytytN 以 以 t

t

yM 1 Mt1 t1 t1 N令SM(1)，即得式 S(1)y(1)S S(1)yt(1)[yt1(1)St(1)]L(1)jyt t（14）式表明S(1)是全部历史数据的加权平均，加权系数分别为t,(1),(1)2,L；显然(1)j

1(1)tt即

S也就是以第t期指数平滑值作为t1

在进行指数平滑时，系数的选择是很重要的。由式（15）可以看出，的大小规定了在新预测值中新数据和原预测值所占的。值越大，新数据所占的就ˆt 体现了修正的幅度，值愈大，修正幅度愈大；值愈小，修正幅度也愈小。虑任何新信息；若选取1，则yˆt1yt，即下期预测值就等于本期观测值，完全不相信过去的信息。这两种情况很难的预测。因此，值应根据时间序列的0～1之间选择。具体如何选择一般可遵循下列原则：①如果时间序列波动不大，比较平稳，则应取小一点，如（0.1～0.5。以减少修正幅度，使预测模型能（0.6～0.800例 解采用指数平滑法，并分别取0.2,0.5和0.8S(1)y1y2 即yˆS(1)

yˆt1yt(1)4某种电器销售额及指数平滑预测值计算表（单位：万元ttttt234567894可以看出，0.20.5和0.8时，预测值是很不相同的。究竟取何值为好，5预测的标准误S为yˆ198851.1754。loaddianqi.txt yt=dianqi;n=length(yt);alpha=[0.20.50.8];m=length(alpha);for一次指数平滑法虽然克服了移动平均法的缺点。但当时间序列的变动出现直线趋为S(1)y(1)S S(2)S(1)(1)S( t式中S(1)为一次指数的平滑值；S(2)为二次指数的平滑值。当时间序列{y}，从某 yˆtmatbtm，m a2S(1)S

531965～1985年的发电总量资料为例，试用二次指数平滑法预测1986年和1987年的发电总量。表6我国发电总量及一、二次指数平滑值计算 tyt16789解取0.3S(1和S(2S(1S(2)676 S(1S(2)，列于6 S(1)3523.1，S(2) （19， 2S(1)S(2)4013.7 (S(1)S(2)) 1 于是，得t21ˆ22ˆ23yˆ(2S(1)S(2))

(S(1)S(2)即

1 yˆ1 S(1) S(

1

1 令t1,2,L,20，由（20）可求出各期的模拟值。计算结果见表6。计算的程序如下：loadfadian.txt yt=fadian;n=length(yt);fori=2:n S(1)y(1)S t tS(2)S(1)(1) tS(3)S(2)(1)S t式中St(3)为三次指数平 abmCm2，m a3S(1)3S(2)S

( bt2(1)2[(65 2(54 (43)St ct2(1)2 St例 7某省全民所有制单位固定资产投资总额及一、二、三次指数平滑值计算表（单位：亿元tyt的估123456789 2某省固定资产投资总额趋势(0)S(0)S(0) y1y2y321.94S

S(1)151.77，S(2)101.28，S(3) （23，a11219.91，b1138.38，c11于是，得t11yˆ1989yˆ12yˆ111a11b11c11 （22并令m1，则得33 3 ( (1 (1)2 (1)2 令t0,1,2,L,10，（24）可求出各期的模拟值，见表7。计算的程序如下：loadtouzi.txt yt=touzi;n=length(yt);for指数平滑预测模型是以时刻t为起点，综合历史序列的信息，对未来进行预测的。选择合适的系数是提高预测精度的关键环节。根据实践经验，的取值范围一般以0.1～0.3为宜。值愈大，系数序列衰度，所以实际上取值大小起着控制参加平均的历史数据的个数的作用。值愈大意味着采用的数据愈少。因此，可以得到选择值的一些基本准则。如果预测目标的基本趋势已发生系统的变化，则值应取得大一些。这样，另外，由于指数平滑是递推计算，所以必须确定初始值S(1),S(2),S(3) §4当时间序列呈直线增加时，可运用一阶差分指数平滑模型来预测。其如下：ytytyˆt1yt(1

其中的为差分记号。式（25）表示对呈现直线增加的序列作一阶差分，构成一个平实际值迭加，作为变量下一期的预测值。对于这个的数学意义可作如下的解释。yt1yt1ytytyt1 的yt1也要改为预测值，亦即成为式（27。的平均数（指数平滑值）加上当前值的实际数进行预测，比一次指数平滑法只用变例7某工业企业1977～1986年锅炉消耗量资料如表8所示，试预测1987年表8某企业锅炉消耗量的差分指数平滑法计算表（0.4）（单位：百吨t22435273843解由资料可以看出，消耗量，除个别年份外，逐期增长量大体在200吨左右，即呈直线增长，因此可用一阶差分指数平滑模型来预测。我们取0.4，初始值为差分序列首项值，计算结果列于表8中。预测1987年消耗量为yˆ19872.494446.49（百吨 ytyt2ytyt

其中2表示二阶差分。yt1yt1ytytyt1yt(yt1yt)ytyt2yt1yt这种组合模型来取代。但是，对于指数平滑法存在的系数的选择问题，以及只§5 NNyˆt1w1ytw2yt1LwNytN1wiyt

ti1期的观测值，N为权数的个数。其调整权数的为w'w2k i1t 式中，i1,2,LN，tNN1,Ln，n为序列数据的个数，wi为调整前的第i个权数，w'为调整后的第i个权数，k为学习常数，e为第t1 差、原观测值和学习常数等三个因素。学习常数k的大小决定权数调整的速度。9某时间序列数据123456789开始，当t2（33，（34，w'w2k

i1tw'w2key 3w'w2key 3利用所得到的权数，计算第t14个观测值y1，增加一个新的观测值y3。ˆt1ˆ4w'yw'y1 21w'0.55420.90.130.312w'0.52720.90.130.22这样进行到t10yˆyˆw' w' 1 2但由于没有t11的观测值y11et1e11y11w'2.0，w' yˆw' w'y 1 2计算的程序如下：N=2;Terr=10000;forj=N+1:m-

w,在开始调整权数时，首先要确定权数个数N和学习常数k。一般说来，当时间序列的观测值呈季节变动时，N应取季节性长度值。如序列以一年为周期进行季节变动N12N4。如果时间序列无小的k值。初始权数的确定也很重要，如无其它依据，也可用1/Nw1(i1,2,L,N 趋势外推法是根据事物的历史和现时资料，寻求事物发展规律，从而推测出事物推；（e）（f）yy0其中系数y0和K值由历史数据利用回归方法求得。对式（35）取对数可得lnylny0

令Ylny，Aln则YAK0a00b1tabt0yˆtK

y1,y2,L,ym；第二部分ym1,ym2,L,y2m；第三部分：y2m1,y2m2,L,y3mmS1y,S

tt且

yt,S3t

t S1

(Kabt)mKab(1bb2Lbm1mm

(Kabt)mKabm1(1bb2Lbm1

(Kabt)mKab2m1(1bb2(1bb2Lbm1)(b1)bm则根据39)SmKabbm

bSmKabm1bmS3mK

b2m1bmb

Sb S 1a

S b

b(bm 1 ab(bm1)KS1 b 据进行检验，检验方法是看给定数据的逐期增长量的比率是否接近某一常数b。即yt1ytbytyt1

例 11某厂收音机解经计算可知yt1ytytyt

5n15m51969年作为开始年份t1。（38，得S1262.5，S2377.8，S3b0.9608，a143.2063，Ky179.7162143.2063y1986179.7162143.20630.960818110（万部functionchanliangglobalab

a=(s2-s1)*(b-1)/(b*(b^m-%定义预测函数functiony=yuce(t)globalabkCompertz曲线tyˆKabt，K0,0a1,0b t量的比率是否接近某一常数b。即lnyt1lnyt lnytlnCompertz曲线用于描述这样一类现象：初期增长缓慢，以后逐渐加快。当达到一ln记

logK(ln

得mS1y't,S2

y't,S3

y'tlnyt。则类似式（42），

Sb S 1a'

S b

b(bm 1 a'b(bm1)K'S1 b 9（8）10Gompertz曲线方程，求出各年收音机销售量的趋势值，并预测1986年的销售量。（48S119.7558，S221.6094，S3（49ba'1.2588，aK'4.8929，Kyˆ1986108.3143（万部）计算的程序如下：yuce=@(t,a,b,k)k*a.^(b.^t);%定义预测的函数yt=log(xsh);n=length(yt);m=n/3;a=(s2-s1)*(b-1)/(b*(b^m-Logistic曲线的一般数学模型是dyry(1y y 1式中c为常y ，K0，a0，0b K

近某一常数b。即1/yt11/yt1/yt1/

y't得ty'tKmS1y't,S2

y't,S3

（42

1S21 Sb S 1 1a12

S bb(bm 55 1

KS1 b 10（8）根据表10Logisticn15m5，根据式（54），得S10.0971，S20.0666，S3（55b0.8493，a0.0174，K

0.00850.01740.8493tloadxsh.txtxsh.txt中yt=1./xsh;n=length(yt);m=n/3;a=(s2-s1)*(b-1)/(b*(b^m-1(1(yyˆ S

定义1 给定随机过程{Xt,tT}。固定t，Xt是一个随 量，设其均值为t，当t变动时，此均值是t的函数，记为tE(Xt固定tX的方差为2。当t变动时，这个方差也是t

2Var(X)E[(X)2 称为随机过程的方差函数。方差函数的平方根t 定义 对随机过程{Xt,tT}，取定t,sT，定义其自协方差函数t,sCov(Xt,Xs)E[(Xtt)(Xss 为刻画{Xt,tT}在时刻t与s之间的相关性，还可将t,s标准化，即定义t

t tt

因此，自相关函数t,s是标准化自协方差函数。定义3设随机序列{Xtt0,1,2,L满足1）EXt)常数；2）tk,tk（k0,1,2,L）与t定义 设平稳序列{t,t0,1,2,L}的自协方差函数kk

2k

2

kk其中k,0k01k10平稳白噪声序列的方差是常数2，因为k0(k0)，则t的任意两个不同时点检验序列平稳性的方法很多，在此介绍其中一种，即Daniel检验。Daniel检验方法建立在Spearman相关系数的基础上。Spearman相关系数是一种秩相关系数。先讲样本的秩的概念。设x1,x2,L,xn是从nx(1x(2Lx(nxix(k，则kxiRi，对每一个i1,2,LnRi是第i个秩统计量。R1R2L,Rn总称为秩统计量。例如，对样本数据-0.8,-3.1,1.1,-5.2,-5.2,-3.1,-0.8,1.1,对于二维总体X,Y)的样本观测数据(x1y1x2y2L,(xn,yn)X,Y的一元样本数据x1,x2L,xny1y2L,ynx1,x2L,xn的秩统计量R1,R2,L,Rny1,y2,L,ynS1,S2,L,Snn这两组秩统计量的相关系数，即Spearman相关系数是n(RiR)(SiSn(SS2in(SS2ini(RRi1 1其中R Ri，S qxy1n(n21)diRiSii1,2,LnH0:XY0，H1:XY可以证明，当(X,Y)是二元正态总体，且H0成立时，统t

n2的t分布t(n2)ttt2n2 H0对于时间序列的样本X1X2L,Xn，记Xt的秩为RtRXt)，考虑变量对(tRtt1,2,LnSpearmanqs qs1n(n21)(tRt

T n11DanielXt计算(tRtt1,2,Ln则 则

t2n2例11某城镇1964～1983每年发生的偷窃率（每万人平均发生的偷盗次数）12（时间t:1~201964年～1983年。问偷盗率数据是否可看表12某城镇年偷盗率数123456789 对显著水平0.05，算得qs0.5624，计算得t统计量的值T2.8857，上2分位点的值t2n2)2.1009，所Tt2n1)，x0=[1.372.961.913.102.082.544.073.622.913.964.192.713.423.023.542.664.114.25 x0=x0';x0=x0(:)'; %计算Qs的值 %计算t统计量的值 %计算上alpha/2分位点XtX1X2LXnˆ1 t tn

k通常有如下两

1nk(Xnt1

X

，0k

且ˆkˆkk的估

1nk(Xnkt1

X

，0k

ˆk，0k

AR模型，即自回归序列（AutoRegressiveModel；Model；ModelAR(p)序设{Xtt0,1,2,LXt1Xt12Xt2LnXtp 其中是零均值、方差是2Xp 为AR(p)序列，(1,2,L,p称为自回归参数向量，其分量j,j1,2,L,p称为自回归系数。

Xt1，BkXtX(B)11B2B2LpBp，(B)XtMA(q)设{Xtt0,1,2,LXtt1t12t2Lqtq， 其中是零均值、方差是2Xq的滑动平均序列，简 记为MA(q)序列(1,2,L,q称为滑动平均参数向量，其分量j,j1,2,Lq性叠加，那么这一输出服从MA模型。Btt1，

(B)11B2B2LXt(B)ARMApq)设{Xtt0,1,2,LXt1Xt1LpXtpt1t1Lqtq， 其中是零均值、方差是2的平稳白噪声，则称Xpq的自回归滑动平均 序列，简记为ARMApqq0时，它是ARpp0MA(q)应用算子多项式(B),(B，式（70）(B)Xt(B)对于一般的平稳序列{Xtt0,1,2,LEXt(Xt)1(Xt1)Lp(Xtp)t1t1L 其中是零均值、方差是2的平稳白噪声，利用后移算子(B),(B)，式（71）可 为(B)(Xt)(B)t关于算子多项式(B),(B)，通常还要作下列假定：1）(B和(B无公共因子，又p0，q0)(B)0设{t0,1,2,L}是零均值平稳白噪声，Var(

2。若{Gk1,2,L}t

k

，G0

XtGktk k可以通过以往白噪声ttj,L， 的运算，其输入是平稳白噪声t，其输出是平稳线性序列Xt。XtE(Xt)GkE(tk)k tk,tE(XtkXt)E[(Gjtkj)(Giti 式中，只有当tkjtijikE

tk

t

)jik0 2 G（k0tk,t iki从而Xt是平稳的。严格地说，在条件（72）下，级数GikGi是收敛的引进算子G(B)G(B)GBk，G k则XtGktkk

例 AR(1)：Xt1Xt1解(B)11B(B11

t，

1(B) t。1t 1B2B2Lk1

kXkBk k即

k

1t Gk（k0,1,2,L 1注意到(B)1B0的根是 。当1时 1

1，即(B)对于ARMA(p,q)序列(B)Xt(B)当(B)满足平稳性条件（(B)的根全在单位圆外）时，可证Gk是负指数下降的，即存在g1,g20，使得Ggeg1k（k0 设{Xtt0,1,2,L是ARMApq(B)XtXtXtI1Xt1I2Xt2LXtIkXtkI(B)Xtt kI(B)1I

IkBk（I1kk

k 则称ARMApqXt具有逆转形式，而{Ikk0,1,2,L 当(B1BLBq的根全在单位圆外，即(B满足可逆条件时，可证存在h1,h20， 2k heh1k2k例13 MA(1)：Xt(1 1解t Xt1Xt1 kIkk0,1,2,L 设{Xtt0,1,2,L是MA(qXt(B)tt1t12t2Lqtq其中是零均值，方差是2Xt

2(122L2

k2( L

1k

k

qk

k

L

k k 1k qkq 1k 122L k0，k

k这表明对于MA(qtsqXtXs不相关，这一性质称为自协方差函数与自相关函数的q步截尾性。AR(p)序列的自相关函数先证明ARMA(p,q)的一个重要E(Xtkt)0，k XtGjtj，XtkGjtk

E(Xtkt)E(Gjtkjt)GjE(tkjt)0，(k Xt是ARpXt1Xt12Xt2LpXtptVar()2，用 （k0）乘上式两端，求均值 E(XtXtk)1E(Xt1Xtk)2E(Xt2Xtk)LpE(XtpXtk)E(tXtk由式（80）E(tXtk0(k0)，故k1k12k2Lpkp，k0即(B)k0k0这是自协方差函数k满足的差分方程。自相关函数k也满足同样的差分方程(B)k0，k在式（81）k1,2,Lp01kk1121Lp

2

Lp

p1p12p2L或者写为1 p11 2 p22 M MM p

p 还可证明，对于AR(p)序列

2

pp0j

2E(2)E[(XXLX)2]2 p由于jiji

1t1 tp jj

j1

jij 22jj

j

)

2j

j

0j

ji

0

cec1k（c,c0 对于自相关函数k也有同样性质。这说明k或k不能在有限步之后截尾，而是按负ARMA(pq)设{Xtt0,1,2,L是ARMApqXt1Xt1LpXtpt1t1L(B)

k2 k

q由（84）kq1,q2,Lqpq11q2q1Lpq1

1q12qLpq2qp1qp12qp2Lpq1 q1p1Lq2qq2p2M MM qp

q

q p 解此方程组，可得自回归系数向量(,,L, t其次，令Yt(BXtXt1Xt1LpXtp，则{Ytt0,1,2,L是参数为(1,2,L,q)T的MA(q序列。{Yt0,1,2,L的自协方差函数t(Y)E(Y

Y)

X

)(X

)]

，（ppp tkt ppp

tk t

i,

jkj k将(Y)代入（78）式k 2(122L2 kijk

2

i, (k1k1Lqkq 1k解此方程组，得滑动平均系数向量(1,2,L,q)T及2对于重要的ARMA(1,1序列Xt1Xt1t1t1111。122, 1,1

k

(1)( 1 12k1 k11

1 (111)(11)k1，k 122 1例 Xt0.8Xt1t0.4t1，其中t~N解利用如下的程序：));% forj=2:10000x(j)=0.8*x(j-1)+elps(j)-0.4*elps(j- %产生样 for 由样本算得样本自相关函数如下（其中的一次运行结果：ˆ10.5284，ˆ20.4287，ˆ30.3442，ˆ4（8810.5231，20.4185，30.3348，4偏相关函数及AR(p)序列偏相关函数的截尾性设{Xtt0,1,2,L是零均值平稳序列，从时间序列预报的角度引出偏相关函数的定义。如果已知{Xt1Xt2LXtk的值，要求对Xt作出预报。此时，可以考虑由{Xt1,Xt2,L,Xtk对Xt的线性最小均方估计，即选择系数k,1,k,2,L,k,k，将

E[(Xtk,jXtj)2] 02k,jjk,jk,i令k,

0j1,2,Lk

j1kjk,iji0，j1,2,L,两端同除0并写成矩阵形式，可知k,j应满足下列线性方 k1k,1 1 k2k,22 MM M

k k

1k,k

k我们称{k,k,k1,2,L为Xt的偏相关函数。式（89）Toeplitz矩阵。这一线性方程组的解k,1,L,k,k推1,1

k

k , 1

k

)( k

k,j1jj11

k, 0）k1,jk,jk1,k1k,k1j j1,2,L,ARp序列的偏相关函数{k,kk1,2,Lpk,kk 1kk

k对于MA(q)和ARMA(p,q)序列，其偏相关函数{k,k,k1,2,L}具有拖尾性。例15 讨论例14所述时间序列的偏相关特性。ˆ1,10.5351，ˆ2,20.2075，ˆ3,30.0683，ˆ4,4计算的程序如下：));% forj=2:10000x(j)=0.8*x(j-1)+elps(j)-0.4*elps(j- %产生样 for for forj=1:k-1s1=s1-rho(k-j)*f(k-s2=s2-rho(j)*f(k- forj=1:k-1f(k,j)=f(k-1,j)-f(k,k)*f(k-1,k- %式（90）中第三式的计 %直接利 工具箱计算偏相关函 (B)Xt(B)t首先要进行模型的识别与定阶，即要判断是ARpMA(qARMApq模型的ARMA设Xt是AR(p)序列Xt1Xt1LpXtpp其中p已知(1,2,L,)T未知，需要估计。在式（82）中ˆj代替jpj1,2,L,p，得1,2,L,p)T的矩估计ˆˆ,ˆ,L,ˆp)T所满足的方程组 ˆ 2 p22 M MMˆ ˆp

1ˆp在式（83）中，以ˆ代替ˆ代替，得2 jjˆ还可以由递推计算得到，式（90）中ˆjj，以ˆk,j代替kj

ˆjˆp,j，j1,2,L,

j1,2,L,MA(q)Xt为MA(qq已知，1,2,L,q)T未知。在式（78）ˆ代替，可得(,,L,)T和2的矩估计(ˆ,ˆ,L,ˆ)T、ˆ2 ˆ2ˆ(1ˆ2Lˆ2ˆˆˆ

k1k1Lqkqˆ2 k1,2,L, 下面介绍MA(q序列参数估计的线性迭代算法。给定ˆ(ˆ,ˆ,L,ˆ)T与ˆ2 一组初值（例如ˆ0,ˆ2ˆ），代入式（93）右边，由此得到等式左边的值ˆ2(1ˆˆˆˆTˆˆ ˆˆˆˆTˆˆ(1(1),(1),L,(1))，称为2 右边得到ˆ2(2)和ˆ(2)。以此方式迭代下去，直到某步的ˆ2(m)和ˆ(m)与前一步 ˆ2(m1)和ˆ(m1)相差不大时停止迭代ˆ2ˆ2(m)，ˆ ARMA(p,q) ,L,ˆ)Ttt

由式（91）得，自回归系数(1,2,L,p)T的矩估计(ˆ,ˆ,L,ˆp)T

ˆ ˆ 2 p22M MMˆ 1ˆp

p令YtXtˆ1Xt1LˆpXtp，则Yt的自协方差函 ，（ tk ijk i,以ˆ代替，得(Y ijk i,把Yt近似看作MA(q)序列，Ytt1t1Lq再根据前面介绍的MA(q序列的参数估计方法求出2和(,,L,)T的矩估计ˆ和(ˆ,ˆ,L,ˆ)T

计算实践表明，当n充分大时，AR(p)序列的矩估计效果较好，但MA(q)ARMA(p,q)序列的矩估计则不太理想，ARMA序列参数ARMA(p,q)序列参因为MA(q序列和ARMApq序列参数的矩估计涉及到非线性运算，所以只解，因而计算简便，快速，易于实现，便于对高阶ARMA序列参数进行初估计。方程组求解或采用偏相关递推求出，将它们近似地当作相应ARMA(p,q)模型的p*个逆函数估计IˆIˆ,LIˆ。 (1BB2LBq)Iˆ jmax(p,q)1,max(p,q)2,L,p得到滑动平均系数向量(,,L,)T的初估计(ˆ,ˆ,L,ˆ)T jIj1Ij12Ij2Lj1I1其中当jq时，j0。得到1,2,L,p的初估计ˆ1,ˆ2,L,ˆp值得注意的是，由上述三步所得的ARMA(p,q)模型参数的初估计，并不能保统计性不变，只需不断地将ˆ(B1ˆBˆB2LˆBq0 逐步用相应根的倒数代替。那么换根后所产生的*(B既满足可逆性又能保持模型的ARMA(p,q)序列参数的最小二乘估XtAR(p)序列，则Xt1Xt12Xt2LpXtpX1X2LXnnpXt1Xt1LpXtpt，tp1,p2,L,n

p)T，定义函2S()(2t

Lp

tpS(达到极小的ˆ为参数的最小二乘估计，记作ˆˆL,L,ˆL)TLS()0，得

(ATA)X X X1 Xp1A

X X

，

p2

M

M 令ˆL1

np n i0,1,L,pj1,2,L,p

i

t

t

ˆL1ˆL1

1 ˆLˆLˆ2 2p02 M MMˆL ˆLˆL p 若用ˆL代替，则得tˆLXˆL LˆL

，tp1,p2,L, t

t从而方差2的最小二nnˆ2

(XˆL LˆL )2

S(ˆL当n充分

npttt

t n |i此时，最小二乘估计ˆL与矩估计ˆ是非常接近的MA(q)与ARMA(p,q)序列参数的最小二乘估ARMA(pq与MA(q这两种序列参数的最小二乘估计方法基本相同，以下只叙述ARMApq)序列参数的估计方法。设Xt是ARMApq)Xt的一个观测样本X1X2L,Xn。如7.2.3小节讨论偏相关函数时，由{Xt1Xt2L,Xtk}预报Xt，现取kt1，即由{Xt1,Xt2L,X1}的线性组合预报Xtt2,3,Ln，估t

t1,jXtj，j2,3,L,由式（89，知t1,1,t1,2,t1,t1满足下列方程 t2t1,1 1 t3 t1,2 2，t2,3,L, M M M

t

1t1,t1

t1因为在理论上，是(, )T，(, )T的函数，所以 t1,,的函数，从而预报的残差平方和也是,

XtjS(,)(

t1,2

LLL LLL估计（CLS）和最大似然估计（ML。设Xt是ARMA(p,q)序列Xt1Xt1LpXtpt1t1L其中是零均值，方差为2X

XtIjXtj

(11B2B2L

p)Xt(11B2B2LqBq)

t(1I1BI2B2L)X

（98 11BLBp(1BLBq)(1IBIB211 1I 331I22I1jj1Ij12Ij2Lj1I1Ijq时，令j0ip时，令i0。将式（100）I11I222

I3331I2Ijjj1Ij1L

给定1,2,L,p)T，1,2,L,q)T（101）可以递推算得逆函数{Ij,I(B)XtXtIjXtj需要注意，,的取值必须使得ARMApqXt,Xt ˆt2(XtIjXtj)2 Xt0，当t0（X1X2LXnIj由式（101）递推算得，因为Ij是,的函数，故残差平方和也是,的函数，即 S(,)(

IjXtj

注意到上述平方和是,的非线性函数，最小化式（102）需用非线性最小二乘法。即X的平稳可逆域中寻求ˆ*,ˆ*，使得ˆ ˆS(ˆ*,*)

是因为它假定过去未观测值等于零。又2的估 rpq

n 假设{Xtt0,1,2,L是零均值正态ARMApqX XX,LX f(x|,)(2)22exp(1xT2其中(,L,, )T，x(x,x,L,x)T，E(XXT)令M21

xTf(x|,2)(22)2M2exp( 2

2 xTL(,|x)lnf(x|,)

ln2

ln(2)2

ln2

2xxx,Lx)T 是rijE(XiXj)r|i设Xt的传递形式是Xt Gjtj， G |i

k

kk|i 由此即M21与2无关xTMx有关，与2无关 S()xT则称S()为平方和函数。进ˆtE(t|x1,x2,L,xn

即ˆt是tx1x2Lxnx1x2Lxn时关于t的估计，则nS()

s可见S()是一残差平方和的形式。sS()ˆs称为最小平方和估计。L(,2|x)1lnMnln(2)nln2S( 2注意M,S()均与2无关

2

S()2(2

ˆ2S()xT n L(,ˆ2|x)

lnM

nln(xTMx) 其中c [ln(2)1lnn]。略去常数c，2L()1ln2

mmmmˆ

)2的最大似然估计（ML估计（108）即得

对于AR(p)序(B)XtXt1Xt1LpXtp S()xTMxd iji

i0其中01dijxtijxtdjii,j0,1,L,pt1 d1p d01 dD 2p，d02 M Mdd

dp

dpp

d0pdS()TD2TdS()达到最小的

ULS3.工具箱相关命中GARCH工具箱可以实现时间序列建模的功能。工具箱中的模型ARMAX(R,M,Nx)的标 ytCiytitjtjkX(t,k k其中{i}是自回归系数，{j是滑动平均系数；X是解释回归矩阵，它的每一个列是一个时间序列，X(t,k)表示第tk列的数据。AR表示自回归中的R个系数向量；Regress表示回归系数k。标函数值，Innovations是残差向量，Sigmas是对应于Innovations的标准差。16Xt0.8Xt1t0.4t1，t~N解编写如下程序));% forj=2:10000x(j)=0.8*x(j-1)+elps(j)-0.4*elps(j- %产生样spec=garchset('R',1,'M',1,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec,x) 11例 Xt0.5Xt10.3Xt2t，t~N10000的数据，对其参数进行估计。解编写如下程序));% x(1)=0;x(2)=0; forj=3:10000x(j)=-0.5*x(j-1)+0.3*x(j- %产生样本spec=garchset('R',2,'M',0,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec,x) 间序列，若已知服从ARMApqpq的值。实际上在建模过程中，首先要解决定阶问题，才能进行参数估计。在定阶与参数估计后，对建立的模型要进行检验。其基本做法是检验模型误差 是否为白噪声。若检验认为 是白 g(x|,L f(x)g(x|0I(f(),g(|))f(x) f(x)|00K－Lf(xg(x| ˆ(m,L)T（kk未知）的最大似然估 kkkAIC(k)2ln(L(ˆ(m)))2k kXt是ARMApqkpq1个，包括自回归参数1,2,L,p)T，滑动平均参数,,L,q)T及2。在式（106）中， 略去lnM L'(,2|x)nln2S( 2由式（108）可得，2的最大似然估计ˆ21 （17L'(,2|x)

nln2 因此ARMA(p,q)序列AIC定阶准则为：选p,q，使 ˆ值，则认为序列是ARMA(pˆ,qˆ)。(B)(Xt)这时，未知参数个数为kpq2，AICpq2

若拟合模型的残差记为ˆt，它是t的估计。例如，对ARp序列，设未知参数的估计是ˆ1,ˆ2,L,ˆp，则残差ˆtXtˆ1Xt1LˆpXtp，t1,2,Lp（设X0X1LX1p0记nk ，k1,2,L,nttLjung－Box2检验统计量是2n(n k1n

H0k0kmH1：对某些kmk002检验法：给定显著性水平，查表得上分位数2(mr22(mr)

22(mr 例 Xt0.8Xt1t0.4t1，t~N));% forj=2:10000x(j)=0.8*x(j-1)+elps(j)-0.4*elps(j- %产生样forforspec=garchset('R',i,'M',j,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec,x); %computeAkaikeandBayesianInformationCriteria 计算结果显示应为ARMA(1,1)序列若实际序列是ARMApqpq值计算其AIC。因此，计算量较大。只有认真反复计算，才能求得真阶p和q。例 Xt0.8Xt1t0.4t1，t~N解编写程序如下：));% forj=2:10000x(j)=0.8*x(j-1)+elps(j)-0.4*elps(j- %产生样spec=garchset('R',1,'M',1);%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec,x)%拟合参数 fore(j)=x(j)-coeffX.AR*x(j-1)-coeffX.MA*e(j- %计算%下面进行chi2检验，是否服从0均值的正态分布，nparam表示估时间序列的m步预报，是根据{Xk,Xk1,L}的取值对未来km时刻的随量Xkm(m0)Xˆk(mXk,Xk1,L的线性组合。引进估k{x|x

Xk

2j,

现设Xt是零均值正态ARMApq)序列。所谓平稳线性最小均方预报Xˆk(m)是指Xˆk(mk，且使kE[(XkmXˆ(m))2]k在正态性的条件下Xˆk(m)Xkm对于{Xk,Xk1,L的条件期望Xˆk(m)E(Xkm|Xk,Xk1,L)E(Xkm|k

因为{XkXk1,LXˆk(m用{XkXk1,L的元素的线性表示不甚方便。考虑到对于ARMA(p,q)序列，总存在着传递形式，令{|d

,d2} k

jk

则可以证明kkXˆk(m可以用k中的元素来表示。设零均值ARMApq序列{Xtt0,1,2,LXtGjtj则Xˆk(m)GjmkjGmkGm1k1Gm2k2

若令ek(mXkmXˆk(m，ek(m)kmG1km1LGm1k 有Var(e(m))2(1G2G2LG2 m1时，由式（126，ek(1)Xk1Xˆk(1)说明k时刻的一步预报误差是k1。Xˆk(m)E(Xkm|k)E(Xkm|k

E(cjXkj|k)cjE(Xkj|kE( |)E(

|)Xˆk(m),mk k km,mXE(km|k)E(km|k)k

m,m如{Xtt0,1,2,L是零均值ARMApq Xkm1Xkm12Xkm2LpXkmp1km1LqkmqkmpXˆk(m)E(Xkm|k)1E(Xkm1|k)LpE(Xkmp|k1E(km1|k)LqE(kmq|k)E(km|k1Xˆk(m1)2Xˆk(m2)LpXˆk(mXt是MA(qXˆk(m)0，mXˆk1(m)Xˆk(m1)Gm[Xk1Xˆk这是因为由式（125）及k1Xk1Xˆk(1)，

Xˆk1(m)Gjmk1jGmk1Gmjk1 Gm(Xk1Xˆk(1))Gm1ikiXˆk(m1)Gm(Xk1Xˆk（130）可以这样理解：Xk1Xˆk(1)表示取得观测值Xk1后带1时刻的m步预报值等于k时刻的m1步预报值与“新息”的和，这里取“权”为Gm。设ARMApqXttXtIjXt则 (

k(m) Ij

Xk1I(1) I(m)

m1II(mi

i

,Xk

,LI(m。又式jjI(m)的递推计算式由逆函数{Ij,j1,2,L决定jARMAp,q)AR(p)序列的预报又Xˆk(m)Xkm，m这就给出AR(p)序列的预报递推Xˆk(1)1Xk2Xk1LpXkXˆk(2)1Xˆk(1)2XkLpXk ˆ k ˆ k

(p1)2Xˆ

(p2)L

(1)pˆk(m)1Xˆk(m1)2Xˆk(m2)LpXˆk(mp),m)Xk,Xk1L,Xkp1。这是AR(p)MA(q)与ARMA(pq)序列的关于MA(q序列{Xtt0,1,2,LXˆk(m)0，m Xˆ(q)(Xˆ(1),Xˆ(2),L,Xˆ （130 kGmmm1,2,LqXˆk1(2)2Xˆk(1)Xˆk(3)2Xˆk1(q1)q1Xˆk(1)Xˆk(q)Xˆk1(q)qXˆk(1)q 1

0 1 Xˆ(q) MXˆ(q)2 k M k

0

) 对于ARMA(p,q)序列，由式义预报向量（134。令*

j j1,2,L, j 0 G1 0 G 2 Xˆ(q) MXˆ(q)M k k 1 Gq1 * * Gqq L1 Gq j

kq1j 020随机产生下列时间序列的10000个观测值Xt0.65Xt1解利用模拟数据计算得到其1，2，3步预测（其中的一次运行结果）为- -以下进行理论验证。模拟数据中X100001.4925，故理论预报值为（由式（133）ˆ10000ˆ10000ˆ100003步预报理论标准差计算如下。因G110.65G220.42251G30.2746Std表示，由式 Std[e10000(1)](1G(1G(1GG2

(3)]

));% forj=2:10000x(j)=-0.65*x(j-1)+elps(j);%产生样本spec=garchset('R',1);%指定模型的结构[coeffX,errorsX,LLFXgarchfit(spec,x)拟合参数%下面计算3步预测值，其中第一个参数返回值是标准差，第二个参数返回值为预[sigmaForecast,x_Forecast]=garchpred(coeffX,x,3) x_theory(3)=-0.65*x_theory(2)%计算三步理论预测值 ARMA序列的建模与预报。在实际中遇到的时间序列往往有三个特性：趋势性、季节性与非平稳性。本节主要采用Box－ARMA序列，再用上面介绍的方法去研究。我们先看一个例子。考虑研究时间序列{Xtt0,1,2,LXt1.5Xt10.5Xt2 改写为(BXtt，其中(B11.5B0.5B2(B)0B1，B221B11在单位圆周上，并非在单位圆外，即原序列非平稳，因而不是AR(2)序列。(XtXt1)0.5(Xt1Xt2)令XtXtXt1Wt0.5Wt1

Wt是AR(1)序列。式（138）定义的运算11阶差分运算，原来的非平稳序列Xt转化为平稳序列Wt。XtXtXt1(1B)2XtXt2Xt1Xt2(1B)2XdXt(1B)dX其中dd d d d

d (1B)11B2

L 1)Bd 设{Xtt0,1,2,LddX Xt(B)dXt(B)

若dX0，则dX (B)(dXt)(B)t，t Xt为一般ARIMApdq序列，若未知，可用dXtx估计。XtX1X2LXn1n1个；2X1X2L,XdtXtX1X2L,Xn，计算样本，说明已服ARMA模型。若自相关函数与偏相关函数1个不是截尾的或拖尾的，说Xt1阶差分Xtt2,3,Ln，并求其样本自相关函数与样本偏相关函数，再用上述方法讨论。这样，直至判断dX是平稳序列为止。在实际计tX1,X2,L,XdWtdZt，td1,L,（1）dt tXtX1Wj1XkWjk tk （2）d

XtX2(t2)(X2X1)(tjtXk(tk)(XkXk1)(tjk1)Wjk，tk

Xt是ARIMApdqp0时，称为IMA(dqq0时，称为ARI(p,d)序列。设{Xtt0,1,2,L是ARIMApdqd1d2当d1XtWtˆk ˆk

kk

kk

Xˆk(m)2Xˆk(m1)Xˆk(m2)WˆkmXˆk(m)Xkm(XkXk1) )ˆk例 XtXt1t1tXtXt

解经1阶差分以后成为MA(1,1)序列。式（145）是可逆的(B)111时，(B0的根1/1I(B)Xt即XtI1Xt1I2Xt2L

（15(1B)Xt(11B)（147(1B)Xt(11B)I(B)X

1B(11B)I即I(B)1B11B(11)B1(1)(BB22B3即1得

1 I(1)j1，j XtX(1

t记WtXtXt1t1t1，则Wt是MA(1)序列。由式（143，kˆkkXk11[Wˆk(1)Xk1Xk](11)Xk11Xˆk （143Xˆk(m)Xˆk(1)，m式（149）的含义是：k1时刻Xt的1步预报是Xk1,Xˆk(1)的 别是11,1。式（148）取tk1，并两边对k）Xˆ(1)(1 ）

k1表明Xt在k时刻的1步预报是Xk,Xk1,L 线性组合，其权是11，(11)1 ,L的指数平滑。这是时间序列预报的指数平滑法的1 k适的。式（149（150）是时间序列IMA(1,1预报的特点。例 13化工生产过程浓度数ARIMA(3,1,3)模型。计算得ˆ11.1762，ˆ21.1158，ˆ3ˆ1

1.8129，ˆ1.8174，ˆ (11.1762B1.1158B20.1706B3)(1B)X (11.8129B1.8174B20.8412B3) 6789 a=nonzeros(a')';%按照原来数据的顺序去掉零元素 %计算偏相关函数 %计算1阶差分 %计算偏相关函数 fori=0:3forspec=garchset('R',i,'M',j,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec,da); %computeAkaikeandBayesianInformationCriteria r=input('输入阶数R＝');m=input('输入阶数spec2=garchset('R',r,'M',m,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec2,da) [sigmaForecast,w_Forecast]=garchpred(coeffX,da,10) %计算10步预报值 %计算原始数据的10步预测值负荷1小时为采样间隔，显然Xt将包24小时的周期性变化。一般，上午有一sXt(1Bs)XD(1Bs)Dt时间序列{Xtt0,1,2,L如果满足下列模型：(Bs)DX(Bs tB(Bs)1BsB2sL B (Bs)11Bs2B2sLQ式（151）中，Et一般不必是白噪声，而可设它是另一个ARIMA(p,d,q)序(B)dEt(B)由式

(Bs)dDX(Bs)d (B)(Bs)dDX(Bs)(B)dE(Bs) 令WtdDX t(B)(Bs)Wt(B)(Bs t式（152）称为乘积型季节性模型，其阶数常用(p,d,q)(P,D,Q)表示。式（152）中的WX经差分dDX的 列相邻时刻与相隔为周期s的时刻之间复杂变化的规律。例如，W(1B)(1Bs

W(1BBs 1 它是MA(s1)模型，但是其系数1,2,L,s102Ls10，s1，s11）(0d,10D,1)sW(1B)(1Bs 2）(0d,11D,1)s(1Bs)W(1B)(1Bs 3）(0,d,2)0,D,2)sW(1B

B2)(1BsB2s

例 测得某地区一口井7年的水埋深数据如表 所示试预报第年全年 例 789111263647510687解（1）12个月的季Wt12对Wt进行稠密系数ARMA模型拟合。用选取的pq的各种阶数形式进行试算，用计算的程序如下：loadwater.txt;%water.txt % fori=s+1:m1y(i-s)=x(i)-x(i-m2=length(y);%周期差分后数据的个数 m3=length(w);%计算最终差分后数据的个数fori=0:3forspec=garchset('R',i,'M',j,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec,w); %computeAkaikeandBayesianInformationCriteria savebdataxywnm1m2 789spec2=garchset('R',1,'M',13,'Disy','off');%指定模型的结构[coeffX,errorsX,LLFX]=garchfit(spec2,w); [sigmaForecast,w_Forecast]=garchpred(coeffX,w,n)%求w的预报值 %求y的预报值forx(m1+j)=yhat(j)+x(m1+j- 例 89解（1）这一时间序列的图形见图3。3数据变化趋而且波动的幅度越来越大。这个例子是Box－Jenkins建模的范例。该范例采用了数据变换的方法使数据平稳化。对Xt作对数变换Ytln设Yt为一个乘积型季节性序列。对Yt再作差WttW(11B2B2L13B13)tloadhang.txt; x=log(x0);%做对数变换y=lnx %周期s=12 fori=s+1:m1y(i-s)=x(i)-x(i-m2=length(y);%周期差分后数据的个数 m3=length(w);%计算最终差分后数据的个数fori=0:3forspec=garchset('