用R语言实现奶牛月产奶量的时间序列分析

1、奶牛月产奶量的时间序列分析本文应用R软件对奶牛月产奶量建立时间序列模型并进行预测。文章主要从 以下几个方面进行：.描述性统计.模型识别.参数估计.模型诊断.预测.其他建模方法及效果对比.结论最终通过多方面对比，我们选择了 ARIMA(0,1,1) 0,1,1)12模型用于以后数 据的预测。描述性统计数据的选取本文引用的是 Data Market中的时间序列数据 “Monthly milk production: pounds per cow. Jan 62 - Dec 75 ，包括从 1962 年 1 月至U 1975 年 12 月共 168 个 月度数据，单位为pounds/month。数据

2、如下：Jan62636455566768697071727374了5600628S04826821Feb 561 566 Jan62636455566768697071727374了5600628S04826821Feb 561 566 61E 622 635 667 696 690 707 756 799 F 775 782Mar 640 65370g 736 762 775 7B5 607 sea 890 6BJ 859 692254654806232 5 g g O 2 7 0 g o O 77778889699May 727 74 2 770 F?811 837 858 871 6S6

3、 94 2 961 957 969 966Jun657716736756 7g8817 26&4 585991393592494 793710 0 8口 4671J 6 6 618 68 0 93 9 8Aug 5gg 517 6Mg 双 647r 722 740 764 785 834 855 837 的7 S58Sep 568 583 604 615 瓯 681 701 72 5 740 790 809 764815 817oct 577 5B7 611际 687 706 723747 800810 791 812 8Z7NOV 553 565 594 5。2 64 5 660 677 6

4、90 711 763 766 760 773 797C2845S8141O5 233689 33891350 0 014D556666777 8 8 8B8从中我们将62-74年，共156条数据作为训练集75年的12个月数据作为测试集，用于最后评价模型预测效果的参考。数据的描述性统计变量统计表1-1数据类型最小值卜四分数中位数均值上四分数最大值数值型数据553.0677.8761.0754.7824.5969.0时间序列的分布图和时间序列的分解如下:62646870727476time时间序列分解图1-1由图可以看出，时间序列含有明显的季节性和上升趋势，且没有波动集群现象，可以考虑季节模型，最

5、常用的是 ARIMA真型。乘法季节模型乘法季节模型是随机季节模型与 ARIMA模型的结合。统计学上纯RIMA (p, d, q )模型记作：O 。其中t代表时间，Xt表示响应序列，B是后移算子，R=1-B, p、d、q分别表示自回归阶数、差分阶数和移动 平均阶数；(B)表示自回归算子；(B)表示滑动平均算子。一个阶数为(P, d, q ) X (P, D, Q ) s的乘积季节模型可表为：oat代表独立干扰项或随机误差项，s的值是一个季节循环中观测的个数， 表示同一周期内不同周期点的相关关系，则描述了不同周期中对应时点上的相关关系，二者结合起来便同时刻画了 2个因数的作用。模型识别模型识别序列

6、平稳化我们首先画出奶牛月产奶量的时序图Time奶牛月产奶量时序图 2-1由上面奶牛月产奶量的时序图可见， 序列呈现明显的趋势。为消除趋势，我 们首先对数据做一阶差分，一阶差分后的时序图和自相关分析图如下：PP出Time阶差分时序图2-2阶差分时序图2-2LL一阶差分ACF图2-3如图所示，序列的趋势基本消除，但当 k=12、24时，序列的样本自相 关系数显著不为0,表明季节性存在，且季节周期为12。因此我们再对一阶 差分后的序列做12步季节差分，得到序列的折线图。ppI Time一阶12步差分折线图2-4此时可以发现序列已经消除了趋势和周期性，基本处于平稳状态，我们 可以进一步通过ADF单位根

7、检验序列的平稳性，得到：Augmented Dickey-Ful1 er restdata: tra1n_d_sDickey-FulHer - -5.901, Lag order - 5, p-value - 0.01 alternative hypothesis: stationary由于p值小于0.05,所以拒绝原假设，选择备择假设，即此时序列是平稳的。模型识别一个序列如果是纯随机的，就意味着它每一次新的变化都无迹可寻，我们可以从中挖掘不出对预测有益的信息。一旦我们发现某个序列是纯随机序列，说明 这个序列已经没有什么可以挖掘的有用信息，因而可以停止对它分析。通常我们 可以用Ljung-Bo

8、x检验，简称LB检验，LB检验的原假设是所检验的序列是纯 随机序列。我们在这对平稳后的序列做 BL检验，检查序列是否为白噪声序列。Box-Ljung testdata: as.vector(train_d_s)X-squared = 71.139, df = 40, p-value = 0_00176由对一次差分和季节差分后的序列做 LB检验，发现p值小于0.05,说明该 序列不是白噪声序列，可以进一步提取序列中的信息。此时画出序列的自相关与偏自相关分析图，如下：(XI o dsrecrAO det ect10(mi 1k_ARiMA13aind 109.0000001ambdal 4.8B4

9、691结果表明时间序列没有存在的可加异常值，但109个数是一个较为明显的新息异常值。此时我们考虑异常值的影响重新修正模型，在进行异常点检测所得结果如下：Call;ar imafx = train, order = c. (0 n 1T 工)，= list (or ccr = c (0, 1,1)period = 12) T method - ml, fo .工109)coefficients;皿 10naimcg0,240b7.587731.9GC3e.0.07230.D65O6.9111ta5 47.2: loo 1 ikeHhooc = -461.07 . sic = 06S. 15修正牛

10、奶产量模型的参数估计MA(1)-0.24860.0723MA(12)-0.58770.0650IO-10932.30116.7619=47.2:对数似然函数=-481.07,AIC=968.15我们发现，考虑异常值后，模型各项的系数并没有太大的改变， 但标准误差 都有所减小，并且对数似然函数值变大， AIC值也变小了，并且IO效应是非常 显著的，总之，模型的效果得到了改善。此时预测模型的表达式为:Yt二Yt-i+Yt-i2-Yt-i3+et+O,24860-i+0.58770_i2+O,2486 0.58770-134.2系数显著性检验我们对改进后的模型ARIMA(0,1,1) X(0,1,1)12的系数显著性进行检验，首先 运用coint()函数计算模型中系数的置信区间。3 5 %97.5%mal-0.3902517-0.1069572smal-0.715186S-0.460215910-109IB.42084864 5, 5118083通过检验，每一项的系数都不包含 0,因此在5%的置信水平下，每一项的 系数都是显著的。4.3白噪声检验如果残差序列为白噪声序列，则说明我们的模型已经充分提取了序列的 信息，我们无法再通过调整模型从模型中获取更多的信息，因而模型的建