《用R语言实现奶牛月产奶量的时间序列分析》由会员分享，可在线阅读，更多相关《用R语言实现奶牛月产奶量的时间序列分析（19页珍藏版）》请在金锄头文库上搜索。
1、奶牛月产奶量的时间序列分析本文应用R软件对奶牛月产奶量建立时间序列模型并进行预测。文章主要从 以下几个方面进行：1. 描述性统计2. 模型识别3. 参数估计4. 模型诊断5. 预测6. 其他建模方法及效果对比7. 结论最终通过多方面对比，我们选择了 ARIMA(0,1,1)x(0,1,1)12模型用于以后数 据的预测。一、描述性统计1.1数据的选取本文引用的是 Data Market中的时间序列数据“Monthly milk production: pounds per cow. Jan 62 - Dec 75”，包括从 1962 年 1 月到 1975 年 12 月共 168 个 月度数据，
2、单位为pounds/month。数据如下：8 7 3 7 45 7 1136 6 7 7 72 5 4 fi- 5 2 5 Do 9 o 7 7 7 7 no2 5 7 6 02 3 6 9 96 6 6 6 60 12 3 4 57 7 7 7 7 70 4 6 18 45O222B7 8 no B 6 82 5 7 3 L O 3 6 Bn- 77 7 7 89 9 4 L 8 616 9 8 0 98 8 8 8 9 8969669一I69615 621 602 635 661 667 645 6BS 6S1 6S7 660 硒g 701 706 677 711 725 723 690
3、734 740 747 711 751 790 S00 763 B00 S09 810 766 B05 7S4 791 760 B02： 315 012 773 BIB 317 S27 797 S4BJanFebMarAprMayJunJU1AugSepOCTNOV DCC62589561640656727697640599568577553 5S26360056665?7427166606175835B7565 55B6462861B7057707.3667S639604611594 634从中我们将62-74年,共156条数据作为训练集，75年的12个月数据作为测 试集，用于最后评价模型预
4、测效果的参考。1.2数据的描述性统计变量统计表1-1数据类型最小值下四分数中位数均值上四分数最大值数值型数据553.0677.8761.0754.7824.5969.0时间序列的分布图和时间序列的分解如下:62&4666870727476time时间序列分解图1-1由图可以看出，时间序列含有明显的季节性和上升趋势，且没有波动集群现 象，可以考虑季节模型，最常用的是ARIMA模型。1.3乘法季节模型乘法季节模型是随机季节模型与ARIMA模型的结合。统计学上纯RIMA（p， d，q）模型记作：0O 。其中t代表时间，Xt表示响应序列，B是后移算子，R=1-B，p、d、q分别表示自回归阶数、差分阶数
5、和移动 平均阶数；0 （B）表示自回归算子；0 （B）表示滑动平均算子。一个阶数为（P， d，q ） X （P， D， Q ） s的乘积季节模型可表为：00at代表独立干扰项或随机误差项，s的值是一个季节循环中观测的个数， 0表示同一周期内不同周期点的相关关系，则描述了不同周期中对应时点上的相关关系，二者结合起来便同时刻画了 2个因数的作用。二、模型识别2.1序列平稳化我们首先画出奶牛月产奶量的时序图。Ti me奶牛月产奶量时序图2-1由上面奶牛月产奶量的时序图可见，序列呈现明显的趋势。为消除趋势，我 们首先对数据做一阶差分，一阶差分后的时序图和自相关分析图如下：Time一阶差分时序图2-2一
6、阶差分ACF图2-3如图所示，序列的趋势基本消除，但当k=12、24时，序列的样本自相 关系数显著不为0,表明季节性存在，且季节周期为12。因此我们再对一阶 差分后的序列做12步季节差分，得到序列的折线图。646668707274Ti me一阶12步差分折线图2-4此时可以发现序列已经消除了趋势和周期性，基本处于平稳状态，我们 可以进一步通过ADF单位根检验序列的平稳性，得到：Augmented Dnckey-Ful1 er Testdata: rraln_d_sDickey-Pul1 er = -5.9091, Lag order = 5, p-valug = . DI alternativ
7、e hypothtsi s: st atio n ar y由于p值小于0.05,所以拒绝原假设，选择备择假设，即此时序列是平稳的。2.2模型识别一个序列如果是纯随机的，就意味着它每一次新的变化都无迹可寻，我们可 以从中挖掘不出对预测有益的信息。一旦我们发现某个序列是纯随机序列，说明 这个序列已经没有什么可以挖掘的有用信息，因而可以停止对它分析。通常我们 可以用Ljung-Box检验，简称LB检验，LB检验的原假设是所检验的序列是纯 随机序列。我们在这对平稳后的序列做BL检验，检查序列是否为白噪声序列。Box-Ljung testdata:. vect r (trai n_d_s-)X-5qua
8、red = 71.139, df = 40, p-value = .00176由对一次差分和季节差分后的序列做LB检验，发现p值小于0.05,说明该 序列不是白噪声序列，可以进一步提取序列中的信息。此时画出序列的自相关与偏自相关分析图，如下：LLoqtuPJoiCL00010203040Lag差分后序列ACF、PACF图2-5首先考虑季节自相关特征，考察延迟12阶、24阶等以周期长度为单位的自 相关系数和偏自相关系数，自相关图显示延迟12阶自相关系数显著非零，但延 迟24阶自相关系数落入2倍标准差范围。而偏自相关图显示延迟12阶和延迟 24阶的偏自相关系数都显著非零。所以可以认为季节自相关特征
9、是自相关系数 截尾，偏自相关系数拖尾，这时以12步为周期的ARMA(0,1)12模型提取差分后 序列的季节自相关信息。再考虑序列12阶以内的自相关系数和偏自相关系数的特征，以确定短期相 关模型。相关系数图和偏自相关系数图显示12阶以内的自相关系数和偏自相关 系数并没有表现出明显的趋势，因此无法判断ARIMA模型中p和q。我们可以 用扩展的自相关法(EACF)法来帮助识别模型参数。AR/MA1234567S91011121.3XX1XX2MMmMMMMM4XXXXX5XXXXX6XXXXXX7XXXXXXX从图中可以看出，从0-10阶非零三角区显示p=0或1,q=1的模型比较适合。 下面我们建立
10、这两种参数的模型，进行拟合效果比较。通过建立模型，我们得到的两种模型的AIC如下：低阶乘法季节模型AIC值 表2-1模型mod1mod2(p,d,q) *(0,1,1)120,1,11,1,1AIC值986.17988,08由表可以发现AIC值最小的乘法季节模型是ARIMA(0,1,1)X(0,1,1)12三、参数估计对上述乘法季节模型给出最大似然估计，及其标准误差。牛奶模型的参数估计：call :jrimA= tr n, order =1.3 1)，= 1 dst(Drder =1 , 1) , peri od = IN)method - rMLrCoeff I ci errts : mal
11、 smal -0.2579 -0.6116 s.e. 0. D7 840.0664slgmaz estimated as 54.CJ9： log Hkelihood = -+91.09, a1 c = 9B6.17牛奶产量模型的参数估计系数估计值标准误差MA(1)-0.25790.0784MA(12)-0.61160.0664=54.09:对数似然函数=-491.09 ,AIC=986.17四、模型诊断4.1异常点检验我们首先画出模型ARIMA(0,1,1) X (0,1,1)J的标准参差图，如下图所示:stonp 一另工盘N_pepcEffi寸二 LO从标准残差图中我们可以发现在1971年附
12、近可能存在异常值，需要对模型 做进一步的异常值检验。在时间序列中，可识别的异常值有两种，即可加异常值 和新息异常值，通常记为AO与IO。我们对模型进行检验，结果如下： detectAO(ml1k_ARiMA13;1 no ao derected detectzo(ml1k_ARiMA13J-ind 109.000000lambdal 4.BB4691结果表明时间序列没有存在的可加异常值，但109个数是一个较为明显的新 息异常值。此时我们考虑异常值的影响重新修正模型，在进行异常点检测所得结果如下：Call;arimafx = trai na order =匚(0 , la 1), seasona
13、l = list (oirdeir = c(0a 1 , 1) 3 period = 12) a method = ,rML,rs i = cC109)Coeffici ent5:mat smal IO-10Q-0.2486 -0.5877 31.9663s.e. 0.07230.06 506. SillsiamaA? esti nated as 47. 2:1 i kielihood = -481.07. aic = 968.15修正牛奶产量模型的参数估计系数估计值标准误差MA(1)-0.24860.0723MA(12)-0.58770.0650IO-10932.30116.7619=47.2: 对数似然函数=-481.07 ,AIC=968.15我们发现，考虑异常值后，模型各项的系数并没有太大的改变，但标准误差 都有所减小，并且对数似然函数值变大，AIC值也变小了，并且IO效应是非常 显著的，总之，模型的效果得到了改善。此时预测模型的表达式为:Y=Y +Y -Y +e+0.2486e +0.5877e +0.2486*0.5877et t-1t-12t-13 tt-1t-12t-134.2系数显著性检验我们对改进后的模型ARIMA(0,1,1)
