用R语言实现奶牛月产奶量的时间序列分析.docx
《用R语言实现奶牛月产奶量的时间序列分析.docx》由会员分享,可在线阅读,更多相关《用R语言实现奶牛月产奶量的时间序列分析.docx(17页珍藏版)》请在冰豆网上搜索。
![用R语言实现奶牛月产奶量的时间序列分析.docx](https://file1.bdocx.com/fileroot1/2022-12/14/8e1e1490-e447-405e-9d92-e51f97b26f2a/8e1e1490-e447-405e-9d92-e51f97b26f2a1.gif)
用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")