用R语言实现奶牛月产奶量的时间序列分析Word文档格式.docx

1、553.0677.8761.0754.7824.5969.0时间序列的分布图和时间序列的分解如下：时间序列分解图1-1由图可以看出，时间序列含有明显的季节性和上升趋势，且没有波动集群现象，可以考虑季节模型，最常用的是ARIMA模型。1.3乘法季节模型乘法季节模型是随机季节模型与 ARIMA 模型的结合。统计学上纯 RIMA （p， d， q ） 模型记作：。其中 t 代表时间，Xt 表示响应序列，B是后移算子， R=1-B，p、 d、 q 分别表示自回归阶数、差分阶数和移动平均阶数；（B）表示自回归算子； （B）表示滑动平均算子。 一个阶数为 （P， d， q ） （P， D， Q ） s 的

2、乘积季节模型可表为： at 代表独立干扰项或随机误差项， s 的值是一个季节循环中观测的个数，表示同一周期内不同周期点的相关关系，则描述了不同周期中对应时点上的相关关系，二者结合起来便同时刻画了 2 个因数的作用。二、 模型识别2.1序列平稳化我们首先画出奶牛月产奶量的时序图。奶牛月产奶量时序图2-1由上面奶牛月产奶量的时序图可见，序列呈现明显的趋势。为消除趋势，我们首先对数据做一阶差分，一阶差分后的时序图和自相关分析图如下：一阶差分时序图2-2一阶差分ACF图2-3如图所示，序列的趋势基本消除，但当k=12、24时，序列的样本自相关系数显著不为0，表明季节性存在，且季节周期为12。因此我们再

3、对一阶差分后的序列做12步季节差分，得到序列的折线图。一阶12步差分折线图2-4此时可以发现序列已经消除了趋势和周期性，基本处于平稳状态，我们可以进一步通过ADF单位根检验序列的平稳性，得到：由于p值小于0.05，所以拒绝原假设，选择备择假设，即此时序列是平稳的。2.2模型识别一个序列如果是纯随机的，就意味着它每一次新的变化都无迹可寻，我们可以从中挖掘不出对预测有益的信息。一旦我们发现某个序列是纯随机序列，说明这个序列已经没有什么可以挖掘的有用信息，因而可以停止对它分析。通常我们可以用Ljung-Box检验，简称LB检验，LB检验的原假设是所检验的序列是纯随机序列。我们在这对平稳后的序列做BL

4、检验，检查序列是否为白噪声序列。 由对一次差分和季节差分后的序列做LB检验，发现p值小于0.05，说明该序列不是白噪声序列，可以进一步提取序列中的信息。此时画出序列的自相关与偏自相关分析图，如下： 差分后序列ACF、PACF图2-5首先考虑季节自相关特征，考察延迟12阶、24阶等以周期长度为单位的自相关系数和偏自相关系数，自相关图显示延迟12阶自相关系数显著非零，但延迟24阶自相关系数落入2倍标准差范围。而偏自相关图显示延迟12阶和延迟24阶的偏自相关系数都显著非零。所以可以认为季节自相关特征是自相关系数截尾，偏自相关系数拖尾，这时以12步为周期的ARMA（0,1）12模型提取差分后序列的季节

5、自相关信息。再考虑序列12阶以内的自相关系数和偏自相关系数的特征，以确定短期相关模型。相关系数图和偏自相关系数图显示12阶以内的自相关系数和偏自相关系数并没有表现出明显的趋势，因此无法判断ARIMA模型中p和q。我们可以用扩展的自相关法（EACF）法来帮助识别模型参数。从图中可以看出，从0-10阶非零三角区显示p=0或1，q=1的模型比较适合。下面我们建立这两种参数的模型，进行拟合效果比较。通过建立模型，我们得到的两种模型的AIC如下：低阶乘法季节模型AIC值 表2-1模型mod1mod2（p,d,q）*（0,1,1）120,1,11,1,1AIC值986.17988,08由表可以发现AIC值

6、最小的乘法季节模型是ARIMA（0,1,1）（0,1,1）12三、 参数估计对上述乘法季节模型给出最大似然估计，及其标准误差。牛奶模型的参数估计：牛奶产量模型的参数估计系数估计值标准误差MA（1）-0.25790.0784MA（12）-0.61160.0664=54.09: 对数似然函数=-491.09 ,AIC=986.17四、 模型诊断4.1异常点检验我们首先画出模型ARIMA（0,1,1）（0,1,1）12的标准参差图，如下图所示：从标准残差图中我们可以发现在1971年附近可能存在异常值，需要对模型做进一步的异常值检验。在时间序列中，可识别的异常值有两种，即可加异常值和新息异常值，通常记

7、为AO与IO。我们对模型进行检验，结果如下：结果表明时间序列没有存在的可加异常值，但109个数是一个较为明显的新息异常值。此时我们考虑异常值的影响重新修正模型，在进行异常点检测所得结果如下：修正牛奶产量模型的参数估计-0.24860.0723IO-109-0.587732.30110.06506.7619=47.2: 对数似然函数=-481.07 ,AIC=968.15我们发现，考虑异常值后，模型各项的系数并没有太大的改变，但标准误差都有所减小，并且对数似然函数值变大，AIC值也变小了，并且IO效应是非常显著的，总之，模型的效果得到了改善。此时预测模型的表达式为：Yt=Yt-1+Yt-12-Y

8、t-13+et+0.2486et-1+0.5877et-12+0.2486*0.5877et-134.2系数显著性检验我们对改进后的模型ARIMA（0,1,1）（0,1,1）12的系数显著性进行检验，首先运用confint（）函数计算模型中系数的置信区间。 通过检验，每一项的系数都不包含0，因此在5%的置信水平下，每一项的系数都是显著的。4.3白噪声检验如果残差序列为白噪声序列，则说明我们的模型已经充分提取了序列的信息，我们无法再通过调整模型从模型中获取更多的信息，因而模型的建立是成功的，此时残差的表现应该大概像那些独立、同分布的具有零均值和相同标准差的变量。如果残差序列是非白噪声序列，则说明

9、我们的模型是不完善的，需要对其进行修正。对残差序列的纯随机性的检验的一个快速有效的方法是tsdiag（ ）函数，然后观察标准化的残差、残差的ACF和LB的检验的p值。残差检验图4-1从第一张图可以看出，标准参差基本在3之内，除一点外，基本没有偏差值。在ACF图中，所有阶数基本为0，因此我们可以认为残差相互独立。第三张图我们可以看出LB检验42阶之内的p值均大于0.05，进一步说明模型的残差时白噪声序列。进一步我们可以做42阶的LB检验，来进一步来确定残差序列互不相关。我们对残差序列进行LB检验，可以发现p值大于0.05，接受原假设，即残差序列是白噪声过程。此时，我们可以判断模型的残差均值为0，

10、且相互独立同分布。4.3正态性检验 我们首先画出残差的Q-Q图来粗略判断是否服从正态分布，若服从正态分布，散点应基本位于斜线附近。我们可以发现右上角有几个点落在了直线之外，因此我们进一步进行Shapiro-Wilk正态性检验：我们发现p值大于0.05，因此接受原假设，即残差符合正态分布。综上所述，我们建立的乘法季节模型ARIMA（0,1,1）（0,1,1）12是合理的。五、 预测利用上述时间序列模型我们给出12期的预测值，并将预测结果与测试集数据对比。预测公式：预测图5-1其中红线是真实值，蓝色线是预测值，灰色区域是置信度为99.5%的置信区间。通过图我们可以发现预测值和真实值非常相近，因此模

11、型的拟合效果不错。六、 其他建模方法及效果对比除上述建模方法，我们还采用了其他两种建模方法，如自动建模和霍兰德建模，这两种建模方法比较简单，但准确性不高，所以所建模型并未采用，但我们仍将方法介绍如下。6.1自动建模我们可以用auto.arima（）函数进行自动建模，它的估计依据是AIC/BIC标准。我们可以发现这种方法建立的模型没有进行一次差分，而只进行了季节差分。我们画出只进行季节差分的序列并进行单位根检验，如下：季节差分时序图6-1我们发现只进行季节差分后的序列并不平稳，且单位根检验的p值大于0.05，接受原假设，即序列不平稳。因此在此基础上建立的模型并不可靠。6.2霍兰德建模霍特季节性指

12、数平滑（Holt-Winters Exponential Smoothing）如果时间序列满足增量模型，存在升降趋势，并且是季节性的数据。这时，可用霍特季节性指数平滑法来做短期预测。霍特季节性指数平滑法，包含三个参数、和。平滑常数，是趋势常数，是季节常数。三个参数的取值范围都是0-1。在R中的实现，可以使用HoltWinters（）方法来实现。alpha，beta和gamma的值，分别是0.546、0、0.751。alpha的值比较小，表明该时间序列的某一时间点的水平预测值，是基于近期观测值和远期观测值。beta为0，表明时间序列趋势部分值不随时间变化而改变的，也就是所有时间点上，趋势的预测值

13、都是初始值。gamma为0.956，这说明季节性部分的预测值，对近期观测值所取得权重很大。我们查看该模型预测结果，如下模型拟合图6-2从图上，可以看出，霍特季节性指数平滑方法在做这个时间序列的季节性预测时很成功，预测结果与原始数据很相近。下面对预测结果做检验，看预测模型是否需要改进。同样的方法计算相关性和做Ljung-Box检验。过程及结果如下：残差ACF图6-3我们可以发现LB检验的p值小于0.05，说明在滞后1-20期时，残差之间是有相关性的，从残差的ACF图我们可以看出，滞后1期、12、13期18、19期都存在相关性，因此该模型并不适用。七、 结论本文结合R软件，经过模型识别、参数估计、

14、模型检验等步骤，对奶牛月产奶量建立了乘法季节模型ARIMA（0,1,1）（0,1,1）12，模型的预测公式如下：该模型能很好的拟合原始数据，并能准确预测未来数据。相关代码：#jiang qiao S201506019 2016.5.4#time series annlysis #Monthly milk productionrm（list = ls（）library（tseries）library（zoo）library（forecast）library（TSA）#读入数据#setwd（E:研究生时间序列milk1-read.csv（monthly-milk-production-pounds

15、-p.csv,header = F,sep = ,milk-ts（milk11:168,2,start = c（62）,frequency = 12）summary（milk） #描述性统计plot（stl（milk,12） #数据可视化分析par（mfrow=c（1,1） train-window（milk,start=62,end=74+11/12） #训练集test-window（milk,start=75） #测试集#模型识别#plot（train）acf（as.vector（train）,lag.max=36） #序列自相关函数train_d-diff（train,differenc

16、es=1） #一次差分plot（train_d）acf（as.vector（train_d）,lag.max=36）train_d_s-diff（diff（train）,lag=12） #一次差分和季节差分plot（train_d_s）adf.test（train_d_s,alt=stationary） #平稳性检验：平稳序列Box.test（as.vector（train_d_s）, lag=40, type=Ljung-Boxpar（mfrow=c（2,1） acf（as.vector（train_d_s）,lag.max = 40,ci.type=ma） #确定p,q,P,Qpacf（a

17、s.vector（train_d_s）,lag.max = 40）eacf（as.vector（train_d_s）#考虑识别乘法季节ARIMA（p,1,q）*（0,1,1）12模型milk_ARIMA1-arima（train,order=c（0,1,1）,seasonal = list（order=c（0,1,1）,period=12）,method = ML）;milk_ARIMA2-arima（train,order=c（1,1,1）,seasonal = list（order=c（0,1,1）,period=12）,method = #模型拟合#极大似然估计-arima（train,

18、order=c（0,1,1）,seasonal = list（order=c（0,1,1）,period=12）, method = milk_ARIMA1#诊断性检验#par（mfrow=c（1,1）plot（window（rstandard（milk_ARIMA1）,start=c（62,1）,ylab=Standardized Residualsabline（h=0）detectAO（milk_ARIMA1）;detectIO（milk_ARIMA1） io=c（109）,method = milk_ARIMA2confint（milk_ARIMA3） #系数显著性检验tsdiag（ m

19、ilk_ARIMA2,gof.lag=42） #残差检验Box.test（window（rstandard（milk_ARIMA2）,start=c（62,1）, lag=42, type=#正态性检验qqnorm（window（rstandard（milk_ARIMA2）,start=c（62,1）qqline（window（rstandard（milk_ARIMA2）,start=c（62,1）shapiro.test（window（rstandard（milk_ARIMA2）,start=c（62,1）#预测#milk_ARIMA1_f-forecast.Arima（milk_ARIMA

20、1,h=12,level=c（99.5）milk_ARIMA1_fplot.forecast（milk_ARIMA1_f）;lines（test,col=red#自动拟合法-auto.arima（train,stationary=F,seasonal=T,ic=aictrain_s-diff（train,12）plot（train_s）adf.test（train_s,alt=） #霍特季节性指数平滑方法milk_foreca-HoltWinters（train）milk_foreca#alpha是平滑指数，beta是趋势指数，gamma是季节指数#预测结果plot（milk_foreca）#模型检验acf（as.vector（milk_foreca2$residuals）,lag.max=20）Box.test（as.vector（milk_foreca2$residuals）, lag=40, type=