时间序列分析实验报告

HHaarrbbiinnIInnssttiittuutteeooffTTeecchhnnoollooggyy实实验验报报告告课程名称：时间序列分析设计题目：非平稳时间序列建模院系：电信学院班级：设计者：学号：指导教师：冀振元设计时间：2010-05-071一、绪论稳序列的直观含义就是序列中不存在任何趋势性和周期性,其统计意义就是一阶矩为常数，二阶矩存在且为时间间隔t的函数。但是在实际问题中，我们常遇到的序列，特别是反映社会、经济现象的序列，大多数并不平稳，而是呈现出明显的趋势性或周期性。这时，我们就不能认为它是均值不变的平稳过程，需要用如下更一般的模型——tttXY来描述。其中,t表示tX中随时间变化的均值,它往往可以用多项式、指数函数、正弦函数等描述，而tY是tX中剔除趋势性或周期性t后余下的部分，往往可以认为是零均值的平稳过程，因而可以用ARMA模型来描述。具体的处理方法可分为两大类：一类是通过某些数学方法剔除掉tX中所包含的趋势性或周期性（即t），余下的tY可按平稳过程进行分析与建模，最后再经反运算由tY的结果得出tX的有关结果。另一类方法是具体求出t的拟合形式，求出tˆ，然后对残差序列{ttXˆ}进行分析，该残差序列可以认为是平稳的。利用前述方法可以求出tYˆ，最后综合可得tttYXˆˆˆ。如果我们对t的形式并不敢兴趣，则可以采取第一类方法，否则可以用第二类方法。需要再强调的一点是，时间序列非平稳性的表现是多种多样的，这里我们所能分析处理的仅是一些较为特殊的非平稳性。二、建模原理2.1平稳化方法2.1.1差分一般而言，若某序列具有线性的趋势，则可以通过对其进行一次差分而将线性趋势剔除掉，然后对差分后的序列拟合ARMA模型进行分析与预测，最后再通过差分的反运算得到tX的有关结果。做一次差分可记为tX，则1tttXXX(1)如果对一阶差分结果再进行差分，则称为高阶差分，差分的次数称为差分的阶，d阶差分记为tdX。2.2.2季节差分反映经济现象的序列，不少都具有周期性，例如，刚收获的小麦，由于供应充足，价格一般是较低的，然后随着供应量的减少，价格会逐渐上涨，直至下一个收获季节又重新开始这一周期。设tX为一含有周期S的周期性波动序列，则,2_,,ststtXXX…为各相应周期点的数值，它们则表现出非常相近或呈现某一趋2势的特征，如果把每一观察值同下一周期相应时刻的观察值相减，这就叫季节差分，它可以消除周期性的影响。季节差分常用s表示,stttsXXX其中S为周期。2.2.3对数变换与差分运算的结合运用如果序列tX含有指数趋势，则可以通过取对数将指数趋势转化为线性趋势，然后再进行差分以消除线性趋势。2.2齐次非平稳在除去局部水平或趋势以外，某些非平稳时间序列显示出一定的同质性，即序列的某一部分与任何其他部分极其相似。这样的序列往往经过若干次差分后可转化为平稳序列，这种非平稳性称为齐次非平稳性，差分的次数称为齐次性的阶。实际中较为常见的是一阶和二阶的齐次非平稳性，表现为两种情形：第一种是序列呈现出水平非平稳性，除了局部水平不同，序列是同质的，也就是说序列的一部分和其他部分非常相似，只是在垂直方向上位置不同。这样的序列经过一次差分后可转化为平稳序列。第二种是序列呈现出水平和斜率的非平稳性，序列既没有固定的水平，也没有固定斜率，除此之外，序列是同质的，这样的序列经过两次差分后可转化为平稳序列。2.3ARIMA模型对于d阶齐次非平稳序列tX而言，tdX是一个平稳序列，设其适合ARMRA（p,q）模型，即ttdaBXB)())((（2）也可写作：ttdaBXBB)()1)((（3）其中：ppBBBB....1)(221（4）qqBBBB...1)(221（5）称此模型为求和自回归滑动平均模型（IntegreatedAutoregressiveMovingAverageModel），简记为ARIMA（p,d,q），其中p,d,q分别表示自回归阶数、差分阶数和移动平均阶数。之所以称之为求和自回归滑动平均模型,是因为差分的反运算即位求和运算。常见的ARIMA模型有以下几种：1.ARIMA（0,1,1）ttaBXB)1()1(1（6）也可写作：111ttttaaXX（7）这就是说，tX是1阶齐次非平稳序列，一次差分后适合MA（1）模型，值得3注意的是，不能认为tX是平稳ARMA（1,1）序列，因为其特征根r=1，不在单位圆内。2.ARIMA（0,2,2）ttaBBXB)1()1(2212（8）即序列tX两次差分后适合MA(2)模型。3.ARIMA（1,1,1）ttaBXBB)1(])1)[(1(11（9）即tX经一次差分适合ARMA（1,1）模型。三、仿真试验如图3-1所示,为某市1985年-1993年各月工业生产总值(数据见附录1)。可以看出Xt具有明显的周期性，做一次差分Yt=Xt-Xt-12，剔除掉周期性。这样就可以按照平稳序列线性模型的知识来进行模式识别，参数估计等。1.求出差分后的数据的均值，并使序列零均值化，也就是将Yt-μ,11NjjNY得到的序列为零均值的平稳随机序列tW,如图3-2所示。2.求Wt的样本自相关函数kˆ和样本偏相关函数kkˆ，本例中选取的k=0，1，2，…，24,以保证k相对于n不能取太大。1122...ˆ;k0,1,2,,24kknknk（10）0ˆˆˆ/;k0,1,2,,24kkrr（11）kkkkkkkˆˆˆ1ˆˆˆˆˆˆˆ1ˆˆˆˆ1ˆˆˆ2111211221112121（12）图3-1:某市1985年-1993年各月工业生产总值Xt4图3-2:零均值化后的平稳序列Wt3.根据样本自相关函数kˆ和样本偏相关函数kkˆ确定模型类别和定阶。如图3-3和图3-4可以看出，当k2时，有kkˆ204.096/2,并且kˆ呈现拖尾现象，故可判定此模型为AR（2）模型。图3-3:样本自相关函数kˆ图3-4:样本偏相关函数kkˆ4.参数估计。)ˆ1/()ˆ1(ˆˆ21211，)ˆ1/(ˆˆˆ212122，)ˆˆˆˆ1(ˆˆ221102a，带入数据可得，3720.0ˆ1，1314.0ˆ2，50086.1ˆ2a。这样就得到了Wt的随机线性模型。tttta211314.03720.0。在进行预报时，可以先对tW进行预报，然后加上均值得到ttWYˆˆ的预报值,然后在反差分得到原始序列的预报值12ˆˆtttXYX。6附录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、样本自相关函数kˆ的值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、样本偏相关函数kkˆ的值1234560.42820.13140.0291-0.09300.0855-8.1917e-04789101112-0.03830.04470.0852-0.0834-0.0865-0.187813141516171870.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);fori=1:18X2((1+6*(i-1)):6*i)=X(i,:);endplot(M,X2)%X2为原始序列M2=linspace(1986,1994,96);Y=zeros(1,96);fori=1:968Y(i)=X2(i+12)-X2(i);endfigureplot(M2,Y)%Y为季节差分后的序列saveX2saveYloadYs=sum(Y);m=s/length(Y);%求序列均值Y2=zeros(1,length(Y));fori=1:length(Y)Y2(i)=Y(i)-m;endplot(M2,Y2);%Y2为零均值化后的序列K=24