应用时间序列分析
大作业
姓 名：陈叮 学号： 5061214012 专业班级：应用统计18 院系：信息工程学院数学系 时间：2017/5/22
题目：对苏格兰异性结婚数据的时序分析
摘要：
本文以苏格兰1855年至2015年异性结婚数据为研究对象，首先运用R软件对1855-2010年的结婚数据绘制时序图、自相关图和做差分进行相关分析，得出一阶差分后的数据是趋于平稳的，然后根据主观确定拟合模型为MA(2)，并运用
.arima()函数进行模型的自动选择，得出ARIMA(0,1,2)模型即R软件里面的autoMA(2)模型最优，故我们所选择的拟合模型MA(2)是最优的，最后运用MA(2)模
型预测并进行预测残差检验，得出了苏格兰2011-2015年异性结婚数据的预测值（29200.45，28905.94，28905.94，28905.94，28905.94）与实际值（29135，30534，27547，28702，28020）相比，相差不大，这说明模型拟合较好，能反映数据的真实水平，而且残差检验也表明预测残差是平均值为0且方差为常数的正态分
布(服从零均值、方差不变的正态分布)，这进一步说明MA(2)模型是可以提供非常合适
预测的模型。
关键词：苏格兰；arima()函数；auto.arima()函数；R软件；预测
二、数据来源
本文的数据是1855-2015年苏格兰的结婚数据(Marriages, Scotland, 1855 to 2015 )，数据可以从网上(https://www.nrscotland.gov.uk/statistics-and-data/statistics/statistics-by-theme/vital-events/marriages-and-civil-partnerships/marriages-time-series-data)下载，数据见附件一。
三、模型的定阶与确定
3.1模型的定阶
3.1.1序列预处理[1]
首先，我们对苏格兰1855年至2010年的时间序列进行时序图和自相关分析，分析结果如图3.1.1.1和图3.1.1.2所示，程序见附录一。
1855-2010年苏格兰结婚数据的时序图50000结婚数据200003500018501900时间19502000
图3.1.1.1苏格兰1855年至2010年异性结婚数据的时序图
Series dataseriesACF-0.20.20.61.00510Lag1520
图3.1.1.2 苏格兰1855年至2010年异性结婚数据的自相关图
图3.1.1.1显示苏格兰的结婚数值的均值和方差变动很大，随着时间的增加，具有明显的上升趋势，是典型的非平稳序列。
图3.1.1.2显示该序列的自相关系数都超出了两倍标准误差，所以进一步证明了该序列是非平稳的。
综上所述，该序列是非平稳序列。 对于该非平稳时间序列，首先我们对数据进行1阶差分处理，以便消除其具
第 1 页 共 3 页
有的强烈的趋势性，来观察数据是否大致趋于平稳。因此得到的1阶差分时间序列图如下：
1855-2010年苏格兰一阶差分结婚数据的时序图一阶差分数据-100000500018501900时间19502000
图3.1.1.3苏格兰1855年至2010年异性结婚数据1阶差分后的时序图
从图3.1.1.3中我们可以看出，经过1阶差分后，该序列的平均值和方差是大致平稳的，所以我们使用ARMIA(p,1,d)模型是很合适的。通过一阶差分，我们去除了结婚数据的趋势部分，剩下了不规则部分。接下来我们可以检验不规则部分中邻项数数值是否具有相关性；如果有的话，可以帮助我们建立一个预测模型来预测苏格兰异性结婚数据的数值趋势。
3.1.2平稳性检验
由图3.1.1.3可以认为该序列一阶差分后，序列基本平稳，为了进一步判断其平稳性，考察差分序列的自相关图和偏自相关图，如图五和图六所示。
图3.1.2.1自相关图显示延迟2阶、3阶、4阶和5阶时的自相关值超出了2倍标准差范围，但是其他在延迟1-25阶的自相关系数都落入2倍标准差范围以内，从而判断该序列有很强的短期相关性，是2阶截尾，所以可以初步认为1阶差分后序列平稳。
图3.1.2.2偏自相关图显示,在延迟2阶和4阶时的偏自相关系数超出了2倍标准差范围，从lag4之后缩小至0，是4阶截尾，该序列趋于平稳。
综上所述，我们可以认为该序列的一阶差分序列自相关图2阶截尾和偏自相关图4阶截尾。
Series dataseriesdiff1ACF-0.20.20.61.00510Lag152025
图3.1.2.1该序列一阶差分后的自相关图
第 2 页 共 4 页
Series dataseriesdiff1Partial ACF-0.3-0.10.1510Lag152025
图3.1.2.2该序列一阶差分后的偏自相关图
3.1.3纯随机性检验
为了判断序列是否有分析价值，必须对序列进行纯随机性检验，即白噪声检验。如表3.1.3.1所示，P值远远小于0.05的临界值，因此，拒绝原假设，即可以认定1阶差分后的序列是平稳非白噪声序列，需要建立模型来拟合该序列的变化趋势。
表3.1.3.1 纯随机性检验
Box.test(dataseriesdiff1,type=\代码
Box-Ljung test
Data: Dataseriesdiff1 X-squared=83.411 Df=30 P-value=6.313e-07
3.2模型确定
3.2.1根据阶数确定模型
由该序列一阶差分的自相关图和偏自相关图，知道自相关值在滞后2阶之后
为0,且偏自相关值在滞后4阶之后缩小至0，那么意味着接下来的ARIMA模型对于一阶时间序列有如下性质：
1、ARMA(4,0)模型：即偏自相关值在滞后4阶之后缩小至0且自相关值缩小至0，则是一个阶层p=4自回归模型。
2、ARMA(0,2)模型：即自相关图在滞后2阶之后为0且偏自相关图缩小至0，则是一个阶数q=2的移动平均模型。
3、ARMA(p,q)模型：即自相关图和偏自相关都缩小至0，则是一个具有p和q大于0的混合模型。
接下来我们利用简单的原则来确定哪个模型是最好的：即我们认为具有最少参数
ARMA(4,0)有4个参数，ARMA(0,2)有2个参数，的模型是最好的。而ARMA(p,q) 第 3 页 共 5 页
2)模型至少有两个变量。因此，ARMA(0,2)模型被认为是最好的模型。ARMA(0，是二阶的移动平均模型，或者称作MA(2)。这个模型可以写作：
Xt??t??1?t?1??2?t?2 (3.2.1)
移动平均模型通常用于建模一个时间序列，此序列具有邻项观测值之间短期的相
关特征。直观地，可以很好理解MA模型可以用来描述苏格兰异性结婚数据中的不规则部分。
3.2.2运用auto.arima()函数[2]自动选择模型
表3.2.2.1 auto.arima()函数运行的结果
代码 Series: 最优模型
auto.arima(dataseries);
dataseries
ARIMA(0,1,2)
Coefficients:
Ma1 Ma2 0.1022 -0.4311 s.e 0.0763 0.0800
sigma^2 estimated as 4121992: log likelihood=-1399.62 AIC=2805.24 AICc=2805.4 BIC=2814.37 从表3.2.1.1中可以得出ARIMA(0,1,2)模型最适合该序列，这与我们前面通过主观确定的模型一样，这说明ARIMA(0,1,2)非常适合拟合该序列。
3.3模型的参数检验
对于ARIMA(0,1,2)模型的参数估计问题，我们运用arima()函数来估计，估计结果如下：
表3.3.1 ARIMA(0,1,2)模型的参数检验
代dataseriesarima=arima(dataseries,order=c(0,1,2));dataseriesarima 码
Call:
arima(x = dataseries, order = c(0, 1, 2))
Coefficients:
Ma1
第 4 页 共 6 页
Ma2
续表：
s.e.
0.1022 0.0763
-0.4311 0.0800
sigma^2 estimated as 4068802: log likelihood = -1399.62, aic = 2805.24
表3.3.1显示，?1?0.1022，?2??0.4311是比较显著地参数，所以模型的方程式确定为:
??Xt??t?0.1022?t?1??0.4311?t?2 3.3.1
3.4模型预测以及预测误差的检验
3.4.1 模型的预测
预测就是要利用已观测到的样本值对序列在未来某个时刻的取值进行估计。为了对随机序列未来发展进行预测，我们对原序列进行短期(h=5)预测，并与实际值进行对比，观测预测效果，预测结果如下表3.4.1.1所示。
表3.4.1.1运用ARIMA(0,1,2)模型预测2010-2015年的结婚数据
year 2011 2012 2013 2014 2015
Point Forecast 29200.45 28905.94 28905.94 28905.94 28905.94
实际观测值 29135 30534 27547 28702 28020
Lo 80 26615.40 25058.73 24685.64 24342.95 24024.26
Hi 80 31785.50 32753.16 33126.25 33468.94 33787.63
Lo 95 25246.95 23022.14 22451.54 21927.44 21440.05
Hi 95 33153.95 34789.75 35360.35 35884.44 36371.84
1，2)模型的拟合表3.4.1.1显示预测值与实际值十分接近，这说明ARIMA(0，效果非常好，很适合该时间序列的拟合。接下来，我们通过绘制预测图，直观的
看预测效果，预测图表明预测效果很好。
Forecasts from ARIMA(0,1,2)2000035000500001850190019502000
图3.4.1.1ARIMA(0,1,2)预测图
第 5 页 共 7 页
3.4.2预测误差的检验
在指数平滑模型下，观测ARIMA模型的预测误差是否是平均值为0且方差为常数的正态分布(服从零均值、方差不变的正态分布)是个好主意，同时也要观测连续预测误差是否自相关。
表3.4.2.1 预测误差的纯随机性检验
代码 Data:
Series dataforecast$residualsBox.test(dataforecast$residuals,type=\
Box-Ljung test
dataforecast$residuals
X-squared=42.036 Df=30 P-value=0.07107
ACF-0.20.20.61.00510Lag1520
图3.4.2.1预测误差的自相关图
相关图显示出在滞后1-20阶中样本自相关值都没有超出显著边界，而且
Ljung-Box检验的p值为0.7107，所以我们推断在滞后1-20阶中没有明显证据说明预测误差是非零自相关的。
为了调查预测误差是否是平均值为0且方差为常数的正态分布(服从零均值、方差不变的正态分布)，我们可以做预测误差的时间曲线图和直方图：
因此，把预测误差看作平均值为0方差为常数正态分布是合理。 既然依次连续的预测误差看起来不是相关，而且看起来是平均值为0方差为常数的正态分布，那么对于苏格兰异性结婚数据，ARIMA(0,1,2)模型看起来是可以提供非常合适预测的模型。
Histogram of forecasterrors0.000000.00020Density-15000-10000-50000forecasterrors50001000015000
图3.4.2.2预测误差的分布图
第 6 页 共 8 页
dataforecast$residuals-1000018500100001900Time19502000
图3.4.2.3预测误差的时序图
四、总结
时间序列是指同一种现象在不同时间上的相继连续的观察值排列而成的一组数字序列，时间序列分析相当重要，应用也相当广泛。
本文的研究对象是苏格兰异性之间的结婚数据，经过建立模型，对模型预测，选取最优模型。利用时间序列方法预测结婚数的未来变化做出改变。本文就是通过ARMIA模型对苏格兰异性之间的结婚数据作了时间序列分析的，并对其进行了预测。从而能客观的看出苏格兰异性之间的结婚数据，有利于政府做出决策。
五、参考文献
[1]王燕.应用时间序列分析（第四版）[M].北京：中国人民大学出版社，2005：65-102. [2]吴喜之，刘苗.时间序列分析R软件陪同[M].北京：机械工业出版社，2016：66-67.
六、附录
时间序列分析的R程序
data=read.csv(\datax=c(data$x);
dataseries=ts(datax,freq=1,start=c(1855)); plot(dataseries,xlab=\
时
间
\
结
婚
数
据
\年苏格兰结婚数据的时序图\acf(dataseries,lag.max=20,plot=TRUE); dataseriesdiff1=diff(dataseries); plot(dataseriesdiff1,xlab=\
时
间
\
一
阶
差
分
数
据
第 7 页 共 9 页
\年苏格兰一阶差分结婚数据的时序图\
acf(dataseriesdiff1,lag.max=25,plot=TRUE); pacf(dataseriesdiff1,lag.max=25,plot=TRUE); Box.test(dataseriesdiff1,type=\auto.arima(dataseries);
dataseriesarima=arima(dataseries,order=c(0,1,2)); dataforecast=forecast.Arima(dataseriesarima,h=5); plot(dataforecast);
plotForecastErrors<-function(forecasterrors)#检验预测误差是均值为零的正态分布 + {
+ #forecasterrors<-rainseriesforecasts2$residuals + mybinsize<-IQR(forecasterrors)/4 + mysd<-sd(forecasterrors)
+ mymin<-min(forecasterrors)-mysd*3 + mymax<-max(forecasterrors)+mysd*3 + mybins<-seq(mymin,mymax,mybinsize)
+ hist(forecasterrors,col=\+ mynorm<-rnorm(10000,mean=0,sd=mysd)
+ myhist<-hist(mynorm,plot=FALSE,breaks=mybins)
+ points(myhist$mids,myhist$density,type=\+ }
plotForecastErrors(dataforecast$residuals);
第 8 页 共 10 页