[论文]时间序列分析-非平稳时间序列建模实验报告.doc

实验报告课程名称： 时间序列分析 设计题目： 非平稳时间序列建模 院 系： 电信学院 班 级： 设 计 者： 学 号： 指导教师： 设计时间： 2010-05-07 一、绪论稳序列的直观含义就是序列中不存在任何趋势性和周期性,其统计意义就是一阶矩为常数，二阶矩存在且为时间间隔t的函数。但是在实际问题中，我们常遇到的序列，特别是反映社会、经济现象的序列，大多数并不平稳，而是呈现出明显的趋势性或周期性。这时，我们就不能认为它是均值不变的平稳过程，需要用如下更一般的模型来描述。其中,表示中随时间变化的均值,它往往可以用多项式、指数函数、正弦函数等描述，而是中剔除趋势性或周期性后余下的部分，往往可以认为是零均值的平稳过程，因而可以用arma模型来描述。具体的处理方法可分为两大类：一类是通过某些数学方法剔除掉中所包含的趋势性或周期性（即），余下的可按平稳过程进行分析与建模，最后再经反运算由的结果得出的有关结果。另一类方法是具体求出的拟合形式，求出，然后对残差序列进行分析，该残差序列可以认为是平稳的。利用前述方法可以求出，最后综合可得。如果我们对的形式并不敢兴趣，则可以采取第一类方法，否则可以用第二类方法。需要再强调的一点是，时间序列非平稳性的表现是多种多样的，这里我们所能分析处理的仅是一些较为特殊的非平稳性。二、建模原理2.1平稳化方法2.1.1差分一般而言，若某序列具有线性的趋势，则可以通过对其进行一次差分而将线性趋势剔除掉，然后对差分后的序列拟合arma模型进行分析与预测，最后再通过差分的反运算得到的有关结果。做一次差分可记为，则 (1)如果对一阶差分结果再进行差分，则称为高阶差分，差分的次数称为差分的阶，d阶差分记为。2.2.2 季节差分反映经济现象的序列，不少都具有周期性，例如，刚收获的小麦，由于供应充足，价格一般是较低的，然后随着供应量的减少，价格会逐渐上涨，直至下一个收获季节又重新开始这一周期。设为一含有周期s的周期性波动序列，则为各相应周期点的数值，它们则表现出非常相近或呈现某一趋势的特征，如果把每一观察值同下一周期相应时刻的观察值相减，这就叫季节差分，它可以消除周期性的影响。季节差分常用表示,其中s为周期。2.2.3对数变换与差分运算的结合运用如果序列含有指数趋势，则可以通过取对数将指数趋势转化为线性趋势，然后再进行差分以消除线性趋势。2.2齐次非平稳在除去局部水平或趋势以外，某些非平稳时间序列显示出一定的同质性，即序列的某一部分与任何其他部分极其相似。这样的序列往往经过若干次差分后可转化为平稳序列，这种非平稳性称为齐次非平稳性，差分的次数称为齐次性的阶。实际中较为常见的是一阶和二阶的齐次非平稳性，表现为两种情形：第一种是序列呈现出水平非平稳性，除了局部水平不同，序列是同质的，也就是说序列的一部分和其他部分非常相似，只是在垂直方向上位置不同。这样的序列经过一次差分后可转化为平稳序列。第二种是序列呈现出水平和斜率的非平稳性，序列既没有固定的水平，也没有固定斜率，除此之外，序列是同质的，这样的序列经过两次差分后可转化为平稳序列。2.3 arima模型对于d阶齐次非平稳序列而言，是一个平稳序列，设其适合armra（p,q）模型，即 （2）也可写作： （3）其中： （4） （5）称此模型为求和自回归滑动平均模型（integreated autoregressive moving average model），简记为arima（p,d,q），其中p,d,q分别表示自回归阶数、差分阶数和移动平均阶数。之所以称之为求和自回归滑动平均模型,是因为差分的反运算即位求和运算。常见的arima模型有以下几种：1.arima（0,1,1） （6）也可写作： （7）这就是说，是1阶齐次非平稳序列，一次差分后适合ma（1）模型，值得注意的是，不能认为是平稳arma（1,1）序列，因为其特征根r=1，不在单位圆内。2.arima（0,2,2） （8）即序列两次差分后适合ma(2)模型。3.arima（1,1,1） （9）即经一次差分适合arma（1,1）模型。三、仿真试验如图3-1所示,为某市1985年-1993年各月工业生产总值(数据见附录1)。可以看出xt具有明显的周期性，做一次差分yt=xt-xt-12，剔除掉周期性。这样就可以按照平稳序列线性模型的知识来进行模式识别，参数估计等。1.求出差分后的数据的均值，并使序列零均值化，也就是将yt-,得到的序列为零均值的平稳随机序列,如图3-2所示。2.求wt的样本自相关函数和样本偏相关函数，本例中选取的k=0，1，2，24,以保证k相对于n不能取太大。 （10） （11） （12）图3-1:某市1985年-1993年各月工业生产总值xt图3-2:零均值化后的平稳序列wt3.根据样本自相关函数和样本偏相关函数确定模型类别和定阶。如图3-3和图3-4可以看出，当k2时，有,并且呈现拖尾现象，故可判定此模型为ar（2）模型。图3-3:样本自相关函数图3-4:样本偏相关函数4.参数估计。，带入数据可得，。这样就得到了wt的随机线性模型。在进行预报时，可以先对进行预报，然后加上均值得到的预报值,然后在反差分得到原始序列的预报值。附录11、 某市1985-1993年各月工业生产总值(单位:万元)1985.0110.939.3411.0010.9811.2911.841985.0710.6210.9012.7712.1512.2412.301986.019.9110.2410.4110.4711.5112.451986.0711.3211.7312.6113.0413.1414.151987.0110.8510.3012.7412.7313.0814.271987.0713.1813.7514.4213.9514.5314.911988.0112.9411.4314.3614.5714.2515.861988.0715.1815.9416.5416.9016.8818.101989.0113.7010.8815.7916.3617.2217.751989.0716.6216.9617.6916.4017.5119.731990.0113.7312.8515.6816.7917.5918.511990.0716.8017.2720.8319.1821.4023.761991.0115.7313.1417.2417.9318.8219.121991.0717.7019.8721.1721.4422.1422.451992.0117.8816.0020.2921.0321.7822.511992.0721.5522.0122.6823.0224.5524.671993.0119.6117.1522.4624.1923.4026.261993.0722.9124.0323.9424.1225.8728.252、样本自相关函数的值1234560.42820.29070.18780.04210.08700.04787891011120.00230.04580.08490.0035-0.0483-0.1972131415161718-0.00698-0.0568-0.00630.15230.14110.11651920212223240.06540.0852-0.0595-0.0880-0.0112-0.0633、样本偏相关函数的值1234560.42820.13140.0291-0.09300.0855-8.1917e-04789101112-0.03830.04470.0852-0.0834-0.0865-0.18781314151617180.1346-0.00170.04830.17050.0707-0.0456192021222324-0.06620.1134-0.1357-0.11980.0966-0.0721附录2matlab源代码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,:);endplot(m,x2) %x2为原始序列m2=linspace(1986,1994,96);y=zeros(1,96);for i=1:96 y(i)=x2(i+12)-x2(i);endfigureplot(m2,y) %y为季节差分后的序列save x2save yload ys=sum(y);m=s/length(y); %求序列均值y2=zeros(1,length(y);for i=1:length(y) y2(i)=y(i)-m;endplot(m2,y2); %y2为零均值化后的序列k=24; %取k=24r0=(1/length(y2)*sum(y2.2); %求r0r=zeros(1,24); %求rkk=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;endp=r/r0; % p为自相关函数%求偏相关函数fai=zeros(1,24);fai(1)=p(1);for i=2:24; p2=zeros(1,i); p2(1)=1; p2(2:i)=p(1:(i-1); t=toeplitz(p2); t2=inv(t); an=t2*p(1:i); fai(i)=an(i);endfai2=zeros(1,25);fai2(1)=1;fai2(2:25)=fai; %fai2为加上fai00后的偏相关函数p3=zeros(1,25);p3(1)=1;p3(2:25)=p;