时间序列分析实验报告

Harbin Institute of Technology

课程名称：设计题目：院 系：班 级：设 计 者：学 号：指导教师：设计时间：实验报告

时间序列分析 非平稳时间序列建模 电信学院 冀振元 2010-05-07

一、绪论

稳序列的直观含义就是序列中不存在任何趋势性和周期性,其统计意义就是一阶矩为常数，二阶矩存在且为时间间隔t的函数。但是在实际问题中，我们常遇到的序列，特别是反映社会、经济现象的序列，大多数并不平稳，而是呈现出明显的趋势性或周期性。这时，我们就不能认为它是均值不变的平稳过程，需要用如下更一般的模型——Xt??t?Yt来描述。其中,?t表示Xt中随时间变化的均值,它往往可以用多项式、指数函数、正弦函数等描述，而Yt是Xt中剔除趋势性或周期性?t后余下的部分，往往可以认为是零均值的平稳过程，因而可以用ARMA模型来描述。具体的处理方法可分为两大类：一类是通过某些数学方法剔除掉Xt中所包含的趋势性或周期性（即?t），余下的Yt可按平稳过程进行分析与建模，最后再经反运算由Yt的结果得出Xt的有关结果。另一类方法是具体求出?t的拟

?t，然后对残差序列{Xt???t}进行分析，该残差序列可以认为是合形式，求出???Y????，最后综合可得X?t。如果我们对?t的平稳的。利用前述方法可以求出Yttt形式并不敢兴趣，则可以采取第一类方法，否则可以用第二类方法。需要再强调的一点是，时间序列非平稳性的表现是多种多样的，这里我们所能分析处理的仅是一些较为特殊的非平稳性。

二、建模原理

2.1平稳化方法

2.1.1差分

一般而言，若某序列具有线性的趋势，则可以通过对其进行一次差分而将线性趋势剔除掉，然后对差分后的序列拟合ARMA模型进行分析与预测，最后再通过差分的反运算得到Xt的有关结果。做一次差分可记为?Xt，则

?Xt?Xt?Xt?1

(1) 如果对一阶差分结果再进行差分，则称为高阶差分，差分的次数称为差分的阶，d阶差分记为?dXt。

2.2.2 季节差分

反映经济现象的序列，不少都具有周期性，例如，刚收获的小麦，由于供应充足，价格一般是较低的，然后随着供应量的减少，价格会逐渐上涨，直至下一个收获季节又重新开始这一周期。设Xt为一含有周期S的周期性波动序列，则Xt,Xt_?s,Xt?2s,…为各相应周期点的数值，它们则表现出非常相近或呈现某一趋

1

势的特征，如果把每一观察值同下一周期相应时刻的观察值相减，这就叫季节差分，它可以消除周期性的影响。季节差分常用?s表示,?sXt?Xt?Xt?s其中S为周期。

2.2.3对数变换与差分运算的结合运用

如果序列Xt含有指数趋势，则可以通过取对数将指数趋势转化为线性趋势，然后再进行差分以消除线性趋势。

2.2齐次非平稳

在除去局部水平或趋势以外，某些非平稳时间序列显示出一定的同质性，即序列的某一部分与任何其他部分极其相似。这样的序列往往经过若干次差分后可转化为平稳序列，这种非平稳性称为齐次非平稳性，差分的次数称为齐次性的阶。实际中较为常见的是一阶和二阶的齐次非平稳性，表现为两种情形：第一种是序列呈现出水平非平稳性，除了局部水平不同，序列是同质的，也就是说序列的一部分和其他部分非常相似，只是在垂直方向上位置不同。这样的序列经过一次差分后可转化为平稳序列。第二种是序列呈现出水平和斜率的非平稳性，序列既没有固定的水平，也没有固定斜率，除此之外，序列是同质的，这样的序列经过两次差分后可转化为平稳序列。

2.3 ARIMA模型

对于d阶齐次非平稳序列Xt而言，?dXt是一个平稳序列，设其适合ARMRA（p,q）模型，即

?(B)(?dXt)??(B)at （2） 也可写作：

?(B)(1?B)dXt??(B)at （3） 其中：

?(B)?1??1B??2B2?....??pBp （4） ?(B)?1??1B??2B2?...??qBq （5）

称此模型为求和自回归滑动平均模型（Integreated Autoregressive Moving Average Model），简记为ARIMA（p,d,q），其中p,d,q分别表示自回归阶数、差分阶数和移动平均阶数。之所以称之为求和自回归滑动平均模型,是因为差分的反运算即位求和运算。

常见的ARIMA模型有以下几种： 1.ARIMA（0,1,1）

(1?B)Xt?(1??1B)at （6） 也可写作：

Xt?Xt?1?at??1at?1 （7）

这就是说，Xt是1阶齐次非平稳序列，一次差分后适合MA（1）模型，值得

2

注意的是，不能认为Xt是平稳ARMA（1,1）序列，因为其特征根r=1，不在单位圆内。

2.ARIMA（0,2,2）

(1?B)2Xt?(1??1B??2B2)at （8） 即序列Xt两次差分后适合MA(2)模型。

3.ARIMA（1,1,1）

(1??1B)[( 1?B)Xt]?(1??1B)at （9）即Xt经一次差分适合ARMA（1,1）模型。

三、仿真试验

如图3-1所示,为某市1985年-1993年各月工业生产总值(数据见附录1)。可以看出Xt具有明显的周期性，做一次差分Yt=Xt-Xt-12，剔除掉周期性。这样就可以按照平稳序列线性模型的知识来进行模式识别，参数估计等。

1.求出差分后的数据的均值，并使序列零均值化，也就是将Yt-μ,??1N?j?1Yj得到的序列为零均值的平稳随机序列Wt,如图3-2所示。

N?，本例中选取的k=0，1，?k和样本偏相关函数?2.求Wt的样本自相关函数?kk2，…，24,以保证k相对于n不能取太大。

WW?W2W2?k?...?Wn?kWn?k?11?k;k?0,1,2,?,24 （10） r

n?k?r?k/r?0;k?0,1,2,?,24 （11） ? ?1???k1???????1?????2 ?k2???????????????kk?????k?1?2???k?1????1???1??2???????????1??2??1?1????1??1?????1??2???（12） ???

???????k?

图3-1:某市1985年-1993年各月工业生产总值Xt

3

图3-2:零均值化后的平稳序列Wt

?确定模型类别和定阶。如图?k和样本偏相关函数?3.根据样本自相关函数?kk?<2/96?0.204,并且??k呈现拖尾现3-3和图3-4可以看出，当k>2时，有?kk象，故可判定此模型为AR（2）模型。

?k 图3-3:样本自相关函数?

? 图3-4:样本偏相关函数?kk?2???2???12/(1???12)， ?1???1(1???2)/(1???12)，?4.参数估计。??1?0.3720，??2?0.1314，?0(1???1??1???2??2)，带入数据可得，??a2???4

?a2?1.0086。 ?这样就得到了Wt的随机线性模型。Wt?0.3720Wt?1?0.1314Wt?2?at。在进行

??W???的预报值,然后在预报时，可以先对Wt进行预报，然后加上均值得到Ytt??Y??X。 反差分得到原始序列的预报值Xttt?12

5

附录1

1、 某市1985-1993年各月工业生产总值(单位:万元) 1985.01 10.93 9.34 11.00 10.98 1985.07 10.62 10.90 12.77 12.15 1986.01 9.91 10.24 10.41 10.47 1986.07 11.32 11.73 12.61 13.04 1987.01 10.85 10.30 12.74 12.73 1987.07 13.18 13.75 14.42 13.95 1988.01 12.94 11.43 14.36 14.57 1988.07 15.18 15.94 16.54 16.90 1989.01 13.70 10.88 15.79 16.36 1989.07 16.62 16.96 17.69 16.40 1990.01 13.73 12.85 15.68 16.79 1990.07 16.80 17.27 20.83 19.18 1991.01 15.73 13.14 17.24 17.93 1991.07 17.70 19.87 21.17 21.44 1992.01 17.88 16.00 20.29 21.03 1992.07 21.55 22.01 22.68 23.02 1993.01 19.61 17.15 22.46 24.19 1993.07 22.91 24.03 23.94 24.12 11.29 12.24 11.51 13.14 13.08 14.53 14.25 16.88 17.22 17.51 17.59 21.40 18.82 22.14 21.78 24.55 23.40 25.87 11.84 12.30 12.45 14.15 14.27 14.91 15.86 18.10 17.75 19.73 18.51 23.76 19.12 22.45 22.51 24.67 26.26 28.25

?k的值 2、样本自相关函数?1 0.4282 7 0.0023 13 -0.00698 19 0.0654 2 0.2907 8 0.0458 14 -0.0568 20 0.0852 3 0.1878 9 0.0849 15 -0.0063 21 -0.0595 4 0.0421 10 0.0035 16 0.1523 22 -0.0880 5 0.0870 11 -0.0483 17 0.1411 23 -0.0112 6 0.0478 12 -0.1972 18 0.1165 24 -0.063 ?的值 3、样本偏相关函数?kk1 0.4282 7 -0.0383 13 2 0.1314 8 0.0447 14 3 0.0291 9 0.0852 15 4 -0.0930 10 -0.0834 16 6

5 0.0855 11 -0.0865 17 6 -8.1917e-04 12 -0.1878 18

0.1346 19 -0.0662 -0.0017 20 0.1134 0.0483 21 -0.1357 0.1705 22 -0.1198 0.0707 23 0.0966 -0.0456 24 -0.0721 附录2

Matlab源代码

X=[10.93,9.34,11.00,10.98,11.29,11.84, 10.62,10.90,12.77,12.15,12.24,12.30, 9.91,10.24,10.41,10.47,11.51,12.45, 11.32,11.73,12.61,13.04,13.14,14.15, 10.85,10.30,12.74,12.73,13.08,14.27, 13.18,13.75,14.42,13.95,14.53,14.91, 12.94,11.43,14.36,14.57,14.25,15.86, 15.18,15.94,16.54,16.90,16.88,18.10, 13.70,10.88,15.79,16.36,17.22,17.75, 16.62,16.96,17.69,16.40,17.51,19.73, 13.73,12.85,15.68,16.79,17.59,18.51, 16.80,17.27,20.83,19.18,21.40,23.76, 15.73,13.14,17.24,17.93,18.82,19.12, 17.70,19.87,21.17,21.44,22.14,22.45, 17.88,16.00,20.29,21.03,21.78,22.51, 21.55,22.01,22.68,23.02,24.55,24.67, 19.61,17.15,22.46,23.19,23.40,26.26,

22.91,24.03,23.94,24.12,25.87,28.25]; %M=linspace(1985,1994,108); X2=zeros(1,108); for i=1:18

X2((1+6*(i-1)):6*i)=X(i,:); end

plot(M,X2) %X2为原始序列 M2=linspace(1986,1994,96); Y=zeros(1,96); for i=1:96

7

输入原始数据

Y(i)=X2(i+12)-X2(i); end figure

plot(M2,Y) %Y为季节差分后的序列 save X2 save Y load Y s=sum(Y);

m=s/length(Y); %求序列均值 Y2=zeros(1,length(Y)); for i=1:length(Y) Y2(i)=Y(i)-m; end

plot(M2,Y2); %Y2为零均值化后的序列 K=24; %取K=24

r0=(1/length(Y2))*sum(Y2.^2); %求r0 r=zeros(1,24); %求rk k=1:24; su=0; for i=1:24

for j=1:length(Y2)-i su=su+Y2(j)*Y2(j+i); end

r(i)=(1/length(Y2))*su; su=0; end

p=r/r0; % p为自相关函数

%求偏相关函数 fai=zeros(1,24); fai(1)=p(1); for i=2:24; p2=zeros(1,i); p2(1)=1;

8

p2(2:i)=p(1:(i-1)); t=toeplitz(p2); t2=inv(t); an=t2*p(1:i)'; fai(i)=an(i); end

fai2=zeros(1,25); fai2(1)=1;

fai2(2:25)=fai; úi2为加上fai00后的偏相关函数 p3=zeros(1,25); p3(1)=1;

p3(2:25)=p; %p3为加上p0后的自相关函数 n=[0:24]; figure plot(n,p3); grid figure plot(n,fai2); grid

%模型参数估计

rou=zeros(1,2); %取φ2j=φj，j=0，1 rou(1)=p(1); rou(2)=p(2); p4=zeros(1,2); p4(1)=1; p4(2)=p(1); T=toeplitz(p4);

F=inv(T)*rou'; %F为φ1和φ2组成的向量

tao=r0*(1-F(1)*p(1)-F(2)*p(2)); %tao为噪声方差 save F

9

本文来源：https://www.bwwdw.com/article/shq6.html