R语言时间序列分析

最新推荐文章于 2024-03-18 09:11:50 发布

zhenglit

最新推荐文章于 2024-03-18 09:11:50 发布

阅读量1.8w

点赞数 36

文章标签：时间序列 R语言 ARIMA 数据分析

本文链接：https://blog.csdn.net/zhenglit/article/details/106265980

版权

R语言时间序列分析

时间序列	白噪声序列	终止分析，无信息可提取
	非白噪声序列	平稳序列	AR/MA/ARMA
		非平稳序列	ARIMA，实际应用中最常见

当拿到一个时间序列的时候，首先分析该时间序列的类型情况，不同类型的序列有不同的处理方式。本文包含以下几个部分：

1、时间序列数据准备

2、时间序列平稳性检验

3、拟合ARIMA模型

4、ARIMA模型的检验诊断

5、用ARIMA模型进行预测

一、时间序列数据准备

> stodata <- read.csv("Stk_Day.csv")
> head(stodata)
      代码       时间 开盘价 最高价 最低价 收盘价 成交量.股. 成交额.元.
1 SH600000 1999-11-10  29.50  29.80  27.00  27.75  174085000 4859102208
2 SH600000 1999-11-11  27.58  28.38  27.53  27.71   29403400  821582016
3 SH600000 1999-11-12  27.86  28.30  27.77  28.05   15007900  421591008
4 SH600000 1999-11-15  28.20  28.25  27.70  27.75   11921000  332952000
5 SH600000 1999-11-16  27.88  27.97  26.48  26.55   23223100  628908032
6 SH600000 1999-11-17  26.50  27.18  26.37  27.18   10052500  268995008
> library(xts)
> x1 <- stodata[,3]
> x2 <- stodata[,7]
> dt <- as.Date(stodata[,2]) #dt必须是时间格式
> tsdata1 <- xts(x1,dt)
> tsdata2 <- xts(x2,dt) #建立时间序列类型
> par(mfrow=c(2,1)) #画图时两行一列的排列
> plot(tsdata1,col="black")
> plot(tsdata2,col="red")

二、时间序列平稳性检验

（1）时序图主观检验

平稳时间序列的均值和方差都为常数。根据这一性质，平稳序列的时序图显示该序列值始终在一个常数附近随机波动，而且波动的范围有界；如果有明显的趋势性或周期性那它通常不是平稳序列。观察tsdata1的时间序列图，明显有趋势性，所以这个时间序列不是平稳序列（一般非平稳序列经过一定的差分计算会呈现出一定的平稳性）。

（2）绘制自相关函数图检验

> acf(tsdata1)

平稳序列具有短期相关性，这个性质表明对平稳序列而言通常只有近期的序列值对现时值的影响比较明显，间隔越远的过去值对现时值的影响越小。随着延迟期数k的增加，平稳序列的自相关系数ρ(k)会比较快的衰减趋向于零，并在零附近随机波动（平稳时序的自相关系数具有拖尾性，即模型的自相关系数ρ(k)呈指数的速度衰减，始终有非零取值，不会在k大于某个常数之后就恒等于零）。而非平稳序列的自相关系数衰减的速度比较慢，这就是利用自相关图进行平稳性检验的标准。

可以通过自相关和偏自相关图来判断时间序列是否平稳。平稳的序列的自相关图和偏相关图不是拖尾就是截尾。截尾是指在某阶之后，系数都为 0；拖尾是指有一个衰减的趋势，但是不都为0。

ACF图形状	序列模型
不衰减到0	不平稳
固定点上有大于0的值（周期性）	周期性
所有点都是0	完全随机（白噪声）
指数衰减到0（拖尾）	AR(p)
正负交替衰减到0（拖尾）	AR(p)
前几个大于0，后面都为0（截尾）	MA(q)

图中的水平虚线是相关系数是否显著的分界线，在虚线之上的是显著的。显著的自相关是考虑应用自回归移动平均模型ARIMA来对时间序列建模的一个证据。可以看出，tsdata1不是平稳序列，但是是自相关的。

> d_tsdata1 <- diff(tsdata1) #diff函数计算差分，输出还是一个时间序列格式
> acf(d_tsdata1[-1]) #-1去掉差分后的第一项NA

tsdata1差分后的时间序列d_tsdata1明显具有“拖尾性”，可见tsdata1的差分序列是平稳的。这也是ARIMA模型原理的关键所在之处。

（3）检验时间序列的自相关

以上两种“看图”的方式都是有人为主观参与的检验方式，其实可以用严格的显著性检验来得到结论的。Box.test(ts)函数输出一个p值，一般情况下，如果p值小于0.05，说明时间序列有显著的自相关；而p值大于0.05，则没有证据证明自相关存在。本例中，p值<<0.05证明时间序列tsdata1是自相关的。

> Box.test(tsdata1)

	Box-Pierce test

data:  tsdata1
X-squared = 3670.1, df = 1, p-value < 2.2e-16

最后总结：自相关的序列一定是 [平稳序列] 或者 [非平稳序列]，肯定不是 [白噪声序列] ，通过自相关检验的序列一定可以用ARIMA模型拟合。

三、拟合ARIMA模型

> library(forecast)
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
> auto.arima(tsdata1) #ARIMA模型自动拟合
Series: tsdata1 
ARIMA(5,1,1) 

Coefficients:
         ar1     ar2     ar3     ar4      ar5      ma1
      0.5378  0.0076  0.0145  0.0415  -0.0980  -0.5958
s.e.  0.0830  0.0192  0.0188  0.0187   0.0165   0.0823

sigma^2 estimated as 0.3692:  log likelihood=-3390.94
AIC=6795.88   AICc=6795.91   BIC=6839.37
> m <- arima(tsdata1,order=c(5,1,1)) #如果提前知道最优阶数，也可以用直接用ARIMA函数，两种方式等价
> confint(m) #展示ARIMA模型系数的置信区间
           2.5 %      97.5 %
ar1  0.375048689  0.70050206
ar2 -0.030063447  0.04527976
ar3 -0.022297674  0.05135013
ar4  0.004946544  0.07812118
ar5 -0.130321386 -0.06562303
ma1 -0.757130275 -0.43451464
> m <- arima(tsdata1,order=c(5,1,1),fixed=c(NA,0,0,NA,NA,NA)) #用fixed参数提出模型中不显著的系数
Warning message:
In arima(tsdata1, order = c(5, 1, 1), fixed = c(NA, 0, 0, NA, NA,  :
  一些AR参数是固定的：把transform.pars设成FALSE

模型的阶数由3个整数来表示，(p,d,q)，这里的p是自回归系数的个数，d是差分的阶数，q是移动平均系数的个数。当建立arima模型时，一般不知道什么是合适的阶数。一般会用函数auto.arima()，让它来为我们选择合适的模型阶数，而不是手动繁琐地寻找p，d，q的最优组合。

ARIMA模型有一个令人头痛的问题：并不是所有的系数都是显著，也就是回归的系数可能等于0。当存在不显著的系数的时候就需要进行剔除，函数arima()含有一个向量参数fixed，可以对ARIMA模型中不显著的系数进行剔除。在本例中，ar2，ar3的置信区间包含0值，因此可以剔除。这里函数输出的警告信息不用管。

四、ARIMA模型的检验诊断

> tsdiag(m)

函数tsdiag()给出一组（3个）图形。通过对这三个图形进行分析可对ARIMA模型进行诊断：

（1）标准化残差图。一个好的ARIMA模型的标准化残差没有波动聚集性，本例中具有一定的聚集性，不过仍可用于预测分析。

（2）残差自相关图。一个好的ARIMA模型的自相关函数没有显著的自相关性，本例中的ACF图表现很好，没有显著的自相关性。

（3）Ljung-Box统计量p值图。如果Ljung-Box统计量的p值都比较大，表明残差没有任何模式，模型已经提取了数据中的所有信息，剩余的仅仅是“噪声”。本例是一个表现很好的例子。

五、用ARIMA模型进行预测

> p <- predict(m,n.ahead=50)
> p
$pred
Time Series:
Start = 3688 
End = 3737 
Frequency = 1 
 [1] 16.47158 16.63245 16.68664 16.60102 16.67366 16.71499 16.72320 16.71790
 [9] 16.72714 16.72682 16.72301 16.71997 16.71938 16.71816 16.71737 16.71718
[17] 16.71736 16.71745 16.71757 16.71771 16.71780 16.71784 16.71785 16.71786
[25] 16.71785 16.71784 16.71783 16.71782 16.71782 16.71782 16.71782 16.71782
[33] 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782
[41] 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782
[49] 16.71782 16.71782

$se
Time Series:
Start = 3688 
End = 3737 
Frequency = 1 
 [1] 0.6072344 0.8367099 1.0064440 1.1473611 1.2840543 1.3885973 1.4779744
 [8] 1.5589136 1.6348650 1.7055611 1.7734991 1.8394046 1.9034873 1.9657107
[15] 2.0262971 2.0852772 2.1426909 2.1986054 2.2531249 2.3063352 2.3583241
[22] 2.4091757 2.4589676 2.5077667 2.5556329 2.6026193 2.6487736 2.6941385
[29] 2.7387531 2.7826530 2.8258714 2.8684387 2.9103836 2.9517324 2.9925098
[36] 3.0327389 3.0724413 3.1116372 3.1503454 3.1885837 3.2263689 3.2637167
[43] 3.3006419 3.3371585 3.3732799 3.4090186 3.4443864 3.4793948 3.5140544
[50] 3.5483755

函数predict()会根据模型计算未来的观测值和他的标准误差。返回值是一个有两个元素的列表，一个元素pred是预测值，另一个元素se是标准误差。这里的pred和se都编码为时间序列格式，因此R会打印预测的数值，同时还有他们的时间序列属性，如果需要预测未来多个观测值，需要设置参数n.ahead，本例中向后预测50个观测值。需要注意的是，每向后预测一个，标准误差会变大，这是因为越往后预测就会变得越不确定。

本例的预测效果不是太好，主要表达使用过程。

tseries

上面的时间序列分析的一些原理以及xts包，其实做时间序列分析还有许多包，下面记一下tseries包。

一、数据准备与查看

> library(tseries)

    ‘tseries’ version: 0.10-47

    ‘tseries’ is a package for time series analysis and
    computational finance.

    See ‘library(help="tseries")’ for details.
> library(forecast)
> air <- AirPassengers #关于ts数据格式的构建，可查询资料学习，本例直接应用一个给出的ts格式
> air
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
1952 171 180 193 181 183 218 230 242 209 191 172 194
1953 196 196 236 235 229 243 264 272 237 211 180 201
1954 204 188 235 227 234 264 302 293 259 229 203 229
1955 242 233 267 269 270 315 364 347 312 274 237 278
1956 284 277 317 313 318 374 413 405 355 306 271 306
1957 315 301 356 348 355 422 465 467 404 347 305 336
1958 340 318 362 348 363 435 491 505 404 359 310 337
1959 360 342 406 396 420 472 548 559 463 407 362 405
1960 417 391 419 461 472 535 622 606 508 461 390 432
> plot(air)
> tsdisplay(air) #直接使用tsdisplay来观察，它包含了时序图，以及acf、pacf两个相关图
> air_train <- ts(as.vector(air[1:132]),frequency=12,start=c(1949,1)) #训练集
> air_test <- ts(as.vector(air[133:144]),frequency=12,start=c(1960,1)) #测试集

二、数据处理与检验

> tsdisplay(air_train)
#air_train存在一个向上的趋势，采用一阶差分消除
> air_train1 <- diff(air_train,1)
> tsdisplay(air_train1)
#差分后的时间序列基本围绕在0波动，基本平稳
> adf.test(air_train1) #单位根检验，p<0.05，序列air_train1为平稳序列

        Augmented Dickey-Fuller Test

data:  air_train1
Dickey-Fuller = -6.679, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(air_train1) : p-value smaller than printed p-value 
#差分后的序列acf图中显示存在一个12阶的季节性影响因素，然后通过做12阶差分图看一下是否可消除
> tsdisplay(diff(air_train1,12))

三、模型拟合与诊断

> mo <- auto.arima(air_train)
> mo
Series: air_train 
ARIMA(1,1,0)(0,1,0)[12] 

Coefficients:
          ar1
      -0.2431
s.e.   0.0894

sigma^2 estimated as 109.8:  log likelihood=-447.95
AIC=899.9   AICc=900.01   BIC=905.46
> confint(mo) #查看模型系数，检验是否有需要剔除的系数
         2.5 %      97.5 %
ar1 -0.4184353 -0.06782791
> tsdiag(mo) #进行模型诊断，三个图表明是个很好的模型
> qqnorm(mo$residuals)
> qqline(mo$residuals) #用QQ图来检验残差的正态性，从而确定模型质量。如果所有点都落在45度线上表明模型表达非常好。
> summary(mo)
Series: air_train 
ARIMA(1,1,0)(0,1,0)[12] 

Coefficients:
          ar1
      -0.2431
s.e.   0.0894

sigma^2 estimated as 109.8:  log likelihood=-447.95
AIC=899.9   AICc=900.01   BIC=905.46

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE
Training set 0.579486 9.907267 7.483159 0.1187348 2.880429 0.2457523
                   ACF1
Training set 0.01227544

四、模型预测与评估

> pred <- forecast(mo,h=12,level=c(99.5))
> pred
         Point Forecast  Lo 99.5  Hi 99.5
Jan 1960       424.1099 394.6963 453.5234
Feb 1960       407.0557 370.1672 443.9442
Mar 1960       470.8257 426.8166 514.8349
Apr 1960       460.8817 410.9544 510.8090
May 1960       484.8681 429.6094 540.1268
Jun 1960       536.8714 476.7621 596.9807
Jul 1960       612.8706 548.2716 677.4695
Aug 1960       623.8708 555.0751 692.6664
Sep 1960       527.8707 455.1199 600.6215
Oct 1960       471.8707 395.3690 548.3725
Nov 1960       426.8707 346.7936 506.9479
Dec 1960       469.8707 386.3711 553.3704
> plot(pred)
> lines(air_test,col="red") #将预测和实际值画在一个坐标图中，图中红线是实际值，蓝线是预测值

推荐几篇很有用的文章，同时也是本文主要参考的文章：

1.https://blog.csdn.net/gdyflxw/article/details/55509656

2.https://zhuanlan.zhihu.com/p/27919755

3.https://blog.csdn.net/LiuYuan_BJTU/article/details/67068404

4.https://zhuanlan.zhihu.com/p/27453841

zhenglit

关注

36
点赞
踩
361

收藏

觉得还不错? 一键收藏
2
评论
R语言时间序列分析

时间序列白噪声序列终止分析，无信息可提取非白噪声序列平稳序列 AR/MA/ARMA 非平稳序列 ARIMA,实际中最常见当拿到一个时间序列的时候，首先分析该时间序列的类型，不同类型的序列有不同的处理方式。本文包含以下几个部分：1、时间序列数据准备 2、时间序列平稳性检验 3、拟合ARIMA模型 4、ARIMA模型的检验诊断 5、用ARIMA模型进行预测 6、预测效果的评估一、时间序列数据准备> stodata &l
复制链接

扫一扫