目录
00引言
上面这个场景大家应该都见过,但是对于很多毕业生来确实不是不想写,就是无从下手。根源来说就是源于自己对毕业设计的要求太高。看下图是不是很真实:
拿到数据面对海量的方法文献无从下手,想到方法看不懂生涩的理论也具体实现不了,跑不出模型结果。本文介绍时间序列的建模思路与实现方法。希望能为各位想用此模型的兄弟姐妹提供整体的思路,以及完整的案例。
本文的案例是用R语言实现。注重代码的实现。若是熟悉其他软件的姐妹们仅供思路参考。
通过读本文你可以获得:
1、王燕老师写的时间序列教材的内容导读。
2、平稳时间序列的案例。
3、完整的平稳时间序列的建模思路。
4、平稳时间序列的R语言代码讲解。
1、课本导读
1.1教材的选择
系统的学习模型教材肯定是比文献要更加全面的,对于时间序列的教材,推荐以下10本书,鉴于版权原因就不给大家贴链接了,需要的单独索要:QQ: 1758714024
选择着本书的原因仅仅是学的时候用的这本书1。理论部分推荐大家看王燕老师编著的第四版。那本基于sas,理论相对这本有所修正。
在系统的介绍之前,先给大家介绍一下这本书,简单导读一下。方便没有用过的童鞋能快速熟悉脉络。
1.2导读
本书共分为六章:
第一章是从定义和发展出发介绍一下时间序列模型以及基础的R语言基础。熟悉过程和R语言的同学可以直接跳过或者简单一了解。
第二章则是时间序列的预处理。主要是介绍了平稳性和纯随机性的定义和检验工具。
第三章是本书学习重点,从理论、建模、估参、检验、优化、预测六部分介绍了平稳性模型ARIMA。多读这一章有利于理解模型脉络和原理。
第四章第五章是非平稳的处理手段:第四章是确定性因素介绍了重要的分解定理;第五章是随机性分析,主要从模型假设出发对ARIMA进行了一些必要推广。
第六章介绍了多元时间序列以及两个重要的检验,单位根和协整。
2、建模思路流程图
3、数据准备
哎呀呀,实在不想找数据,咋么办呢随机生成把,em…只能这个样子了。
# 数据模拟
set.seed(0) # 设置随机种子、保证结果的可重复性
# 注意这个函数只能生成平稳随机序列,生成时注意ar参数的设置
da <- arima.sim(n = 100,list(ar = c(0.8, -0.5), ma = c(0.2, 0.5)))
er <- rnorm(100,0,0.1)
data = da + er
> data
Time Series:
Start = 1
End = 100
Frequency = 1
[1] 0.12389653 1.49945239 1.20745886 1.31747199 1.57587201 0.36540926
[7] -1.44492485 -1.78292312 -1.48271702 -1.03006832 -0.54900628 -1.11676718
[13] -0.05195605 1.28991098 2.63158995 1.76024674 1.66887178 0.13920201
[19] 1.70160454 1.96635121 1.25980331 -0.69688467 -2.77764077 -3.63842647
[25] -3.98345147 -1.05019197 1.46007723 2.03231032 1.79755660 -0.32303663
[31] 1.51065917 1.05375105 0.98906708 0.12789197 0.22674246 0.30845434
[37] -1.90705909 -3.48193446 -2.77681165 -1.14464315 0.03006204 0.21778821
[43] -1.07278976 -0.85279759 -2.00234148 -0.80575079 -0.16724405 0.53839246
[49] 0.92228923 0.46849446 -0.53525609 -0.86716155 -0.17052567 1.57628817
[55] 1.79469323 1.24031480 -0.88594329 -3.11117116 -3.51114947 -0.98442478
[61] 1.56417591 2.42945315 1.02397262 -1.01183100 -0.57115941 -0.32765595
[67] 1.76948693 2.42886603 3.01508993 0.98020403 -0.05131142 -2.07124064
[73] -0.95345380 0.13839126 0.42572337 1.90721663 1.59843191 2.07847050
[79] 1.87443764 0.30767669 -0.92514337 -1.44593647 -2.02971665 -0.40349496
[85] -1.16088431 1.11945315 0.57036474 2.28298753 1.34368371 -0.37764917
[91] -1.41256753 -1.33152368 -1.90564256 -0.14162506 -1.07420129 0.03475444
[97] 0.57130049 0.61898680 0.67894434 0.22853495
上述就是生成的数据,在生成残差时,方差特意搞得小了点,意思下就性了,别太难为人家模型了。arima.sim
这个函数估计大家也用不着就不讲解用法了,用到时候保证ar
参数最后平稳就行。
plot(data) # 虽说平稳,但是还是让大家看看效果
4、模型建立
现在完全有种数据在手,模型我有的感觉。下面进行第一步预处理。
4.1预处理
预处理主要是对数据进行插补,填时间,看平稳性,当然我的数据是平稳的就不需要做平稳性检验。插值方法主要用样条插值和线性插值;在这里总结一下平稳性检验的方法。
1、时序图。很多时序图一看就会有明显的趋势,可以看出不平稳。
2、自相关图和偏自相关图。平稳的序列具有相关系数截尾性,可以通过这两个图来判断时序的平稳性。
3、当然经常用的还会有统计量的检验方式,ad
和adf
检验,也称为单位根检验。2R语言中做单位的函数有:fUnitRoots包中的UnitrootTests
和adfTest
。2、用tseries包中的adf.test
和pp.test
。具体的做法见函数的hlep例子和参数解释。
4、白噪声检验是检验纯随机性的。
4.2 acf、pacf
下面是定阶依据1,来源于文献[1]74页。
自相关系数 | 偏自相关系数 | 阶数 |
---|---|---|
拖尾 | p阶截尾 | AR(p) |
q阶截尾 | 拖尾 | MA(q) |
拖尾 | 拖尾 | ARMA(p,q) |
接下来是自相关图和偏自相关图。
从上面两个图,可以看出这个是ARMA的混合模型,无法达到定阶的目的。画这两个图需要注意的是:
1、需要序列检验平稳,不然画这个图无意义。
2、有时候这个图无法达到定阶的目的,就像ARMA这种模型。
3、看这两个图比较主观,没有固定标准,一般不用其确定模型形式。
4.3参数估计
参数估计一般用以下三种方法:ML:极大似然估计,CSS:条件最小二乘估计,CSS-ML:二者的混合方法分别在估计参数时设置不同的参数。到目前为止我们还没有定参数。用auto.arima
函数定下参数,该函数的参数设置灵活,我们贴出他的函数参数:
function (y, d = NA, D = NA, max.p = 5, max.q = 5, max.P = 2,
max.Q = 2, max.order = 5, max.d = 2, max.D = 1, start.p = 2,
start.q = 2, start.P = 1, start.Q = 1, stationary = FALSE,
seasonal = TRUE, ic = c("aicc", "aic", "bic"),
stepwise = TRUE, nmodels = 94, trace = FALSE, approximation = (length(x) >
150 | frequency(x) > 12), method = NULL, truncate = NULL,
xreg = NULL, test = c("kpss", "adf", "pp"),
test.args = list(), seasonal.test = c("seas", "ocsb",
"hegy", "ch"), seasonal.test.args = list(),
allowdrift = TRUE, allowmean = TRUE, lambda = NULL, biasadj = FALSE,
parallel = FALSE, num.cores = 2, x = y, ...)
输入参数有:原始序列、差分阶数、遍历的p、q范围等等。他可以输出参数估计的结果,我们也可以用arima
函数来估计参数:
function (x, order = c(0L, 0L, 0L), seasonal = list(order = c(0L,
0L, 0L), period = NA), xreg = NULL, include.mean = TRUE,
transform.pars = TRUE, fixed = NULL, init = NULL, method = c("CSS-ML",
"ML", "CSS"), n.cond, SSinit = c("Gardner1980",
"Rossignol2011"), optim.method = "BFGS",
optim.control = list(), kappa = 1e+06)
其他参数就不介绍了,method
时用来设置估计方式的。
贴出两种方式的建模结果。
> (fit <- auto.arima(data))
Series: data
ARIMA(2,0,3) with zero mean
Coefficients:
ar1 ar2 ma1 ma2 ma3
1.1219 -0.4597 -0.2205 0.3663 -0.6751
s.e. 0.1247 0.1004 0.1180 0.0921 0.1169
sigma^2 estimated as 0.7361: log likelihood=-125.76
AIC=263.53 AICc=264.43 BIC=279.16
> # 参数估计,ML极大似然估计,CSS条件最小二乘估计,CSS-ML二者的混合方法
> (x.fit <- arima(data,order = c(2,0,3),method = "ML"))
Call:
arima(x = data, order = c(2, 0, 3), method = "ML")
Coefficients:
ar1 ar2 ma1 ma2 ma3 intercept
1.1219 -0.4597 -0.2205 0.3663 -0.6751 -0.0001
s.e. 0.1248 0.1004 0.1180 0.0921 0.1169 0.1197
sigma^2 estimated as 0.6992: log likelihood = -125.76, aic = 265.53
两个函数的结果一致。可以看出其参数真实结果和真实参数相差不大。模型效果还可以。
4.4模型评价
模型评价主要的几个指标是:
1、计算平均相对误差。
2、残差是否为白噪声序列,代码看附录的函数说明。
3、系数的一阶范数。
4、AIC信息准则。
5、模型的解释性。
6、模型的复杂度。
7、参数的显著性检验。
4.5预测
最后进行预测:
> x.fore <- forecast(x.fit, h = 10)
> x.fore
> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
101 -0.102760845 -1.174409 0.968887 -1.741705 1.536183
102 0.006071458 -1.436700 1.448843 -2.200457 2.212600
103 -0.056892787 -1.803108 1.689322 -2.727498 2.613713
104 -0.066670782 -1.814055 1.680713 -2.739064 2.605723
105 -0.048694801 -1.872936 1.775546 -2.838631 2.741241
106 -0.024032019 -1.931831 1.883767 -2.941759 2.893695
107 -0.004626320 -1.951014 1.941761 -2.981369 2.972117
108 0.005807276 -1.948516 1.960131 -2.983073 2.994687
109 0.008591691 -1.945835 1.963019 -2.980447 2.997630
110 0.006919018 -1.948378 1.962216 -2.983450 2.997288
> plot(x.fore)
最终画出预测图:
上面是往后预测10期的预测图。
5、附录
除了本文用到的函数,下面还有一些时间序列里用到的函数并给出说明:
函数 | 含义 |
---|---|
ts | 构造时间序列 |
abline | 添加垂直水平线 |
Box.test | 白噪声检验 |
filter | 拟合AR或MA模型 |
ar | 拟合ar模型的参数 |