用R语言实现奶牛月产奶量的时间序列分析.docx
《用R语言实现奶牛月产奶量的时间序列分析.docx》由会员分享，可在线阅读，更多相关《用R语言实现奶牛月产奶量的时间序列分析.docx（17页珍藏版）》请在冰豆网上搜索。
用R语言实现奶牛月产奶量的时间序列分析
奶牛月产奶量的时间序列分析
本文应用R软件对奶牛月产奶量建立时间序列模型并进行预测。
文章主要从以下几个方面进行：
1.描述性统计
2.模型识别
3.参数估计
4.模型诊断
5.预测
6.其他建模方法及效果对比
7.结论
最终通过多方面对比，我们选择了ARIMA（0,1,1）×（0,1,1）12模型用于以后数据的预测。
一、描述性统计
1.1数据的选取
本文引用的是DataMarket中的时间序列数据“Monthlymilkproduction:
poundspercow.Jan62–Dec75”，包括从1962年1月到1975年12月共168个月度数据，单位为pounds/month。
数据如下：
从中我们将62-74年,共156条数据作为训练集，75年的12个月数据作为测试集，用于最后评价模型预测效果的参考。
1.2数据的描述性统计
变量统计表1-1
数据类型
最小值
下四分数
中位数
均值
上四分数
最大值
数值型数据
553.0
677.8
761.0
754.7
824.5
969.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的乘积季节模型可表为：
at代表独立干扰项或随机误差项，s的值是一个季节循环中观测的个数，
表示同一周期内不同周期点的相关关系，
则描述了不同周期中对应时点上的相关关系，二者结合起来便同时刻画了2个因数的作用。
二、模型识别
2.1序列平稳化
我们首先画出奶牛月产奶量的时序图。
奶牛月产奶量时序图2-1
由上面奶牛月产奶量的时序图可见，序列呈现明显的趋势。
为消除趋势，我们首先对数据做一阶差分，一阶差分后的时序图和自相关分析图如下：
一阶差分时序图2-2
一阶差分ACF图2-3
如图所示，序列的趋势基本消除，但当k=12、24时，序列的样本自相关系数显著不为0，表明季节性存在，且季节周期为12。
因此我们再对一阶差分后的序列做12步季节差分，得到序列的折线图。
一阶12步差分折线图2-4
此时可以发现序列已经消除了趋势和周期性，基本处于平稳状态，我们可以进一步通过ADF单位根检验序列的平稳性，得到：
由于p值小于0.05，所以拒绝原假设，选择备择假设，即此时序列是平稳的。
2.2模型识别
一个序列如果是纯随机的，就意味着它每一次新的变化都无迹可寻，我们可以从中挖掘不出对预测有益的信息。
一旦我们发现某个序列是纯随机序列，说明这个序列已经没有什么可以挖掘的有用信息，因而可以停止对它分析。
通常我们可以用Ljung-Box检验，简称LB检验，LB检验的原假设是所检验的序列是纯随机序列。
我们在这对平稳后的序列做BL检验，检查序列是否为白噪声序列。
由对一次差分和季节差分后的序列做LB检验，发现p值小于0.05，说明该序列不是白噪声序列，可以进一步提取序列中的信息。
此时画出序列的自相关与偏自相关分析图，如下：
差分后序列ACF、PACF图2-5
首先考虑季节自相关特征，考察延迟12阶、24阶等以周期长度为单位的自相关系数和偏自相关系数，自相关图显示延迟12阶自相关系数显著非零，但延迟24阶自相关系数落入2倍标准差范围。
而偏自相关图显示延迟12阶和延迟24阶的偏自相关系数都显著非零。
所以可以认为季节自相关特征是自相关系数截尾，偏自相关系数拖尾，这时以12步为周期的ARMA（0,1）12模型提取差分后序列的季节自相关信息。
再考虑序列12阶以内的自相关系数和偏自相关系数的特征，以确定短期相关模型。
相关系数图和偏自相关系数图显示12阶以内的自相关系数和偏自相关系数并没有表现出明显的趋势，因此无法判断ARIMA模型中p和q。
我们可以用扩展的自相关法（EACF）法来帮助识别模型参数。
从图中可以看出，从0-10阶非零三角区显示p=0或1，q=1的模型比较适合。
下面我们建立这两种参数的模型，进行拟合效果比较。
通过建立模型，我们得到的两种模型的AIC如下：
低阶乘法季节模型AIC值表2-1
模型
mod1
mod2
（p,d,q）*（0,1,1）12
0,1,1
1,1,1
AIC值
986.17
988,08
由表可以发现AIC值最小的乘法季节模型是ARIMA（0,1,1）×（0,1,1）12
三、参数估计
对上述乘法季节模型给出最大似然估计，及其标准误差。
牛奶模型的参数估计：
牛奶产量模型的参数估计
系数估计值
标准误差
MA
（1）
-0.2579
0.0784
MA（12）
-0.6116
0.0664
=54.09:
对数似然函数=-491.09,AIC=986.17
四、模型诊断
4.1异常点检验
我们首先画出模型ARIMA（0,1,1）×（0,1,1）12的标准参差图，如下图所示：
从标准残差图中我们可以发现在1971年附近可能存在异常值，需要对模型做进一步的异常值检验。
在时间序列中，可识别的异常值有两种，即可加异常值和新息异常值，通常记为AO与IO。
我们对模型进行检验，结果如下：
结果表明时间序列没有存在的可加异常值，但109个数是一个较为明显的新息异常值。
此时我们考虑异常值的影响重新修正模型，在进行异常点检测
所得结果如下：
修正牛奶产量模型的参数估计
系数估计值
标准误差
MA
（1）
-0.2486
0.0723
MA（12）
IO-109
-0.5877
32.3011
0.0650
6.7619
=47.2:
对数似然函数=-481.07,AIC=968.15
我们发现，考虑异常值后，模型各项的系数并没有太大的改变，但标准误差都有所减小，并且对数似然函数值变大，AIC值也变小了，并且IO效应是非常显著的，总之，模型的效果得到了改善。
此时预测模型的表达式为：
Yt=Yt-1+Yt-12-Yt-13+et+0.2486et-1+0.5877et-12+0.2486*0.5877et-13
4.2系数显著性检验
我们对改进后的模型ARIMA（0,1,1）×（0,1,1）12的系数显著性进行检验，首先运用confint（）函数计算模型中系数的置信区间。
通过检验，每一项的系数都不包含0，因此在5%的置信水平下，每一项的系数都是显著的。
4.3白噪声检验
如果残差序列为白噪声序列，则说明我们的模型已经充分提取了序列的
信息，我们无法再通过调整模型从模型中获取更多的信息，因而模型的建立是成功的，此时残差的表现应该大概像那些独立、同分布的具有零均值和相同标准差的变量。
如果残差序列是非白噪声序列，则说明我们的模型是不完善的，需要对其进行修正。
对残差序列的纯随机性的检验的一个快速有效的方法是tsdiag（）函数，然后观察标准化的残差、残差的ACF和LB的检验的p值。
残差检验图4-1
从第一张图可以看出，标准参差基本在±3之内，除一点外，基本没有偏差值。
在ACF图中，所有阶数基本为0，因此我们可以认为残差相互独立。
第三张图我们可以看出LB检验42阶之内的p值均大于0.05，进一步说明模型的残差时白噪声序列。
进一步我们可以做42阶的LB检验，来进一步来确定残差序列互不相关。
我们对残差序列进行LB检验，可以发现p值大于0.05，接受原假设，即残差序列是白噪声过程。
此时，我们可以判断模型的残差均值为0，且相互独立同分布。
4.3正态性检验
我们首先画出残差的Q-Q图来粗略判断是否服从正态分布，若服从正态分布，散点应基本位于斜线附近。
我们可以发现右上角有几个点落在了直线之外，因此我们进一步进行Shapiro-Wilk正态性检验：
我们发现p值大于0.05，因此接受原假设，即残差符合正态分布。
综上所述，我们建立的乘法季节模型ARIMA（0,1,1）×（0,1,1）12是合理的。
五、预测
利用上述时间序列模型我们给出12期的预测值，并将预测结果与测试集数据对比。
预测公式：
Yt=Yt-1+Yt-12-Yt-13+et+0.2486et-1+0.5877et-12+0.2486*0.5877et-13
预测图5-1
其中红线是真实值，蓝色线是预测值，灰色区域是置信度为99.5%的置信区间。
通过图我们可以发现预测值和真实值非常相近，因此模型的拟合效果不错。
六、其他建模方法及效果对比
除上述建模方法，我们还采用了其他两种建模方法，如自动建模和霍兰德建模，这两种建模方法比较简单，但准确性不高，所以所建模型并未采用，但我们仍将方法介绍如下。
6.1自动建模
我们可以用auto.arima（）函数进行自动建模，它的估计依据是AIC/BIC标准。
我们可以发现这种方法建立的模型没有进行一次差分，而只进行了季节差分。
我们画出只进行季节差分的序列并进行单位根检验，如下：
季节差分时序图6-1
我们发现只进行季节差分后的序列并不平稳，且单位根检验的p值大于0.05，接受原假设，即序列不平稳。
因此在此基础上建立的模型并不可靠。
6.2霍兰德建模
霍特季节性指数平滑（Holt-WintersExponentialSmoothing）
如果时间序列满足增量模型，存在升降趋势，并且是季节性的数据。
这时，可用霍特季节性指数平滑法来做短期预测。
霍特季节性指数平滑法，包含三个参数α、β和γ。
α平滑常数，β是趋势常数，γ是季节常数。
三个参数的取值范围都是0-1。
在R中的实现，可以使用HoltWinters（）方法来实现。
alpha，beta和gamma的值，分别是0.546、0、0.751。
alpha的值比较小，表明该时间序列的某一时间点的水平预测值，是基于近期观测值和远期观测值。
beta为0，表明时间序列趋势部分值不随时间变化而改变的，也就是所有时间点上，趋势的预测值都是初始值。
gamma为0.956，这说明季节性部分的预测值，对近期观测值所取得权重很大。
我们查看该模型预测结果，如下
模型拟合图6-2
从图上，可以看出，霍特季节性指数平滑方法在做这个时间序列的季节性预测时很成功，预测结果与原始数据很相近。
下面对预测结果做检验，看预测模型是否需要改进。
同样的方法计算相关性和做Ljung-Box检验。
过程及结果如下：
残差ACF图6-3
我们可以发现LB检验的p值小于0.05，说明在滞后1-20期时，残差之间是有相关性的，从残差的ACF图我们可以看出，滞后1期、12、13期18、19期都存在相关性，因此该模型并不适用。
七、结论
本文结合R软件，经过模型识别、参数估计、模型检验等步骤，对奶牛月产奶量建立了乘法季节模型ARIMA（0,1,1）×（0,1,1）12，模型的预测公式如下：
Yt=Yt-1+Yt-12-Yt-13+et+0.2486et-1+0.5877et-12+0.2486*0.5877et-13
该模型能很好的拟合原始数据，并能准确预测未来数据。
相关代码：
#jiangqiaoS2015060192016.5.4
#timeseriesannlysis
#Monthlymilkproduction
rm（list=ls（））
library（tseries）
library（zoo）
library（forecast）
library（"TSA"）
#################读入数据############################################
setwd（"E:
\\研究生\\时间序列"）
milk1<-read.csv（"monthly-milk-production-pounds-p.csv",header=F,sep=","）
milk<-ts（milk1[1:
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,differences=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-Box"）
par（mfrow=c（2,1））
acf（as.vector（train_d_s）,lag.max=40,ci.type="ma"）#确定p,q,P,Q
pacf（as.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="ML"）;
#######################模型拟合#######################################
#极大似然估计
milk_ARIMA1<-arima（train,order=c（0,1,1）,seasonal=list（order=c（0,1,1）,period=12）,
method="ML"）
milk_ARIMA1
###################诊断性检验########################################
par（mfrow=c（1,1））
plot（window（rstandard（milk_ARIMA1）,start=c（62,1））,ylab="StandardizedResiduals"）;
abline（h=0）
detectAO（milk_ARIMA1）;
detectIO（milk_ARIMA1）
milk_ARIMA2<-arima（train,order=c（0,1,1）,seasonal=list（order=c（0,1,1）,period=12）,
io=c（109）,method="ML"）
milk_ARIMA2
confint（milk_ARIMA3）#系数显著性检验
tsdiag（milk_ARIMA2,gof.lag=42）#残差检验
Box.test（window（rstandard（milk_ARIMA2）,start=c（62,1））,lag=42,type="Ljung-Box"）
#正态性检验
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））））
##########################预测#######################################
par（mfrow=c（1,1））
milk_ARIMA1_f<-forecast.Arima（milk_ARIMA1,h=12,level=c（99.5））
milk_ARIMA1_f
plot.forecast（milk_ARIMA1_f）;lines（test,col="red"）
#######################################################
#自动拟合法
milk_ARIMA2<-auto.arima（train,stationary=F,seasonal=T,ic="aic"）
milk_ARIMA2
train_s<-diff（train,12）
plot（train_s）
adf.test（train_s,alt="stationary"）
######################################
#霍特季节性指数平滑方法
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="Ljung-Box"）
升级会员 