时间序列分析随机模拟

一、随机模拟实验

1.实验题目

？）和（，是否适用于很大)2(AR 1AR ),1()(

Y

2.实验目的和意义

（1）实验目的：检验公式是否适用于AR(1)和AR （2）的预测估计。

（2）实验意义：若题目成立，则对于所有的AR （1）和AR （2）模型，其预测会趋向于一条水平之直线，

3.简述实验方法和步骤

（1）首先模拟一个AR （1）序列，生成K 个数列，将n 个数搁置起来，预测搁置的n 个值， 参数估计，是否符合模型。最后在估计的序列均值上画一条水平线。

（2）首先模拟一个AR （2）序列，生成K 个数列，将n 个数搁置起来，预测搁置的n 个值， 参数估计，是否符合模型。最后在估计的序列均值上画一条水平线。

4.具体实施过程

（1）AR （1）过程

首先模拟一个过程的，)1(1008.0AR 。模拟48个值，将最后八个值搁置起来，与预测值比较。

（a ）验证 和 的极大似然估计：

图表4.1极大似然估计

从图表4.1可以看出，该模型符合AR （1）模型，所以我们继续下一步。

（b ）预测接下来的8个值，并画出带这8个预测值的序列，在估计的序列均值上画一条水平线。画出预测及其95%预测极限。

图表4.2预测及估计均值水平线

从图4.2中可以看出，预测值落在预测区间内，并且趋向于一条水平直线。此时仅仅是 很小的时候趋势已经很明显了，所以当 越大，)(

Y 越趋向于一个均值 。

（2）AR （2）过程

首先模拟一个过程的，，)2(10075.0-5.121AR 。模拟52个值，将最后12个值搁置起来，与预测值比较。 （a ）验证 和 的极大似然估计：

图表4.3极大似然估计

从图表4.3可以看出，该模型符合AR （2）模型，所以我们继续下一步。

（b ）预测接下来的12个值，并画出带这12个预测值的序列，在估计的序列均值上画一条

水平线。画出预测及其95%预测极限。

图表4.4预测及估计均值水平线

从图4.4中可以看出，预测值落在预测区间内，并且趋向于一条水平直线。和AR(1)一样仅仅是 很小的时候趋势已经很明显了，所以当 越大，)(

Y 越趋向于一个均值 。所以AR(2) 也满足公式。

5.实验结果分析和讨论

我们很好地模拟了AR （1）和AR （2）模型，其预测值也很好的落在预测区间内，两个模

型的预测均趋向于一个均值，所以 ).2(AR 1AR ),1()(）和（，适用于很大

Y （注：程序在附录）

附录：

set.seed(132456)

series=arima.sim(n=48,list(ar=0.8))+100 future=window(series,start=41) series=window(series,end=40) #(a)

model=arima(series,order=c(1,0,0)) model #（b ）

plot(model,n.ahead=8,ylab='Series&Forecasts',col=NULL,pch=19) abline(h=coef(model)[names(coef(model))=='intercept'])

plot(model,n.ahead=8,ylab='Series,Forecasts,Actuals&Limits',pch=19) points(x=(41:48),y=future,pch=3)

abline(h=coef(model)[names(coef(model))=='intercept'])

AR(2)模型

library(TSA)

set.seed(132456)

series=arima.sim(n=52,list(ar=c(1.5,-0.75)))+100

actual=window(series,start=41)

series=window(series,end=40)

#(a)

model=arima(series,order=c(2,0,0))

model

#(b)

result=plot(model,n.ahead=12,ylab='Series&Forecasts',col=NULL,pch=19) abline(h=coef(model)[names(coef(model))=='intercept'])

forecast=result$pred

cbind(actual,forecast)

plot(model,n1=25,n.ahead=12,ylab='Series,Forecasts,Actuals&Limits',pch=19) points(x=(41:52),y=actual,pch=3)

Abline

(h=coef(model)[names(coef(model))=='intercept'])

二、案例分析

1.问题题目和意义

服装消费是与人类生活密不可分的生活方式，服装消费增长给我国经济的增长带来巨大贡献。本文以2002-2016年各季度我国服装销售量为研究对象,运用ARIMA模型做时间序列分析,并预测2017-2018年各季度的服装销售量,对服装销售量预测分析提供理论基础。

2.数据来源

查阅《中华人民共和国国家统计局》网站季度数据，给出2002-2016年各季度服装销售量如下表：

3.简述采用的方法步骤

（1）方法：建立时间序列模型

（2）步骤：根据所给的数据，分别进行模型识别、模型拟合、模型诊断以及模型预测

4.具体实施过程

(1)模型识别

为了合理地应用ARMA模型，先观察时间序列图，图表2.1左图中的上升趋势将会建立一个非平稳模型。故我们需要对数据进行对数变换处理.

图表2.1服装销售总量时间序列图

图表2.2 服装销售量的样本ACF

图表2.2是服装销售量的ACF，我们需要进行一阶差分，使它有一定趋势。

图2.3 服装销售量的一阶差分

图表2.3显示的是经一次差分后服装销售量的时间序列图。由图可以知道，经过该处理数据表现出平稳性，但是仍不平稳，因此还需进一步对数据进行处理。

图表2.4 服装销售量的一次差分序列的样本ACF

图表2.4表示，我们应用季节差分法所得的序列来建立模型。

图表2.5显示的是服装销售量数据经过一次差分和季节差分后的时间序列图。此时大部分的季节性已经消失了。

图表2.5 工业生产总值经一次和季节差分后的时间序列图

图表2.6印证了经两次差分后的时间序列已经几乎不再具有自相关性。此图

也说明建立自相关的简单模型就行了。

图表2.6 服装销售量经一次和季节差分后的样本ACF

考虑识别乘法季节4)1,1,0()1,1,0( ARIMA 此模型满足上述诸多要求。模型构建的诊断阶段还要修正。

（2）模型拟合

对于工业生产总值的 4)1,1,0()1,1,0( ARIMA 模型 ，图表2.7给出了极大似然 估计及其标准误差。

图表2.7 服装销售量模型的参数估计

（3）诊断性检验

为了对估计后的 4)1,1,0()1,1,0( ARIMA 模型进行检验，首先图表2.8显示的是残差直方图。

图表2.8 4)1,1,0()1,1,0( ARIMA 模型的残差

图表2.9显示的是残差的QQ 正态图。图表展示了残差的分位数-分位数图，这些点看起来非常接近一条直接—特别是中间部分。该图使我们不能拒绝模型误差项是正态的假设

图表2.9残差QQ 图：4)1,1,0()1,1,0( ARIMA 模型

（4）预测

预测是时间序列分析的意义所在，图表2.10给出了服装销售量序列及其前置8期的预测，以及上下95%的预测极限。图中也显示出了1997-2016年的观测数据。预测很好地模仿出了序列的随机周期性，而且预测极限也显示出预测之精度令人满意。

5.结果分析和讨论

通过平稳性检验，阶数识别，模型诊断等过程，对我国12002-2016年的国内

服装销售量构建了ARIMA 模型，从拟合的效果来看，预测很好地模仿出了序列的随机周期性，说明模型拟合的效果非常好，具有一定的可信度，但是还有待于做进一步完善。

（注：程序见附录）

附录：

图表2.1

data<-read.csv("C:\\Users\\lenovo\\Desktop\\cc.csv",header = T)

data=ts(data[,2],frequency=4,start=c(2002,1))

library(TSA)

win.graph(width=6.5,height=3,pointsize=8); oldpar=par; par(mfrow=c(1,2)) plot(data,ylab='costume',type='o');

图表2.2

acf(as.vector(data))

图表2.3

plot(diff(data),ylab='First Difference of costume',type='o')

图表2.4

acf(as.vector(diff(data)))

图表2.5

plot(diff(diff(data),lag=4),ylab='First and Seasonal Difference of costume')图表2.6

acf(as.vector(diff(diff(data),lag=4)))

图表2.7

m1.cos=arima(data,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4)) m1.cos

图表2.8

win.graph(width=3,height=3,pointsize=8);

hist(window(rstandard(m1.cos)),xlab='Standardized Residuals')

图表2.9

win.graph(width=2.5,height=3,pointsize=8);

qqnorm(window(rstandard(m1.cos))); qqline(window(rstandard(m1.cos)));图表2.10

win.graph(width=6.5,height=3,pointsize=8)

plot(m1.cos,n1=c(2002,1),n.ahead=8,pch=19,ylab='costume')