毕业实用统计模型(一)——时间序列

00引言

同学B 同学A 导师 同学A, 论文写怎么样了? 还没开始。 同学B,我也没开始。 这是我的论文终稿 老师,你好,还没开始写。 导师:同学 B顺利毕业. 同学B 同学A 导师

上面这个场景大家应该都见过,但是对于很多毕业生来确实不是不想写,就是无从下手。根源来说就是源于自己对毕业设计的要求太高。看下图是不是很真实:
在这里插入图片描述
拿到数据面对海量的方法文献无从下手,想到方法看不懂生涩的理论也具体实现不了,跑不出模型结果。本文介绍时间序列的建模思路与实现方法。希望能为各位想用此模型的兄弟姐妹提供整体的思路,以及完整的案例。
本文的案例是用R语言实现。注重代码的实现。若是熟悉其他软件的姐妹们仅供思路参考。
通过读本文你可以获得:
1、王燕老师写的时间序列教材的内容导读。
2、平稳时间序列的案例。
3、完整的平稳时间序列的建模思路。
4、平稳时间序列的R语言代码讲解。

1、课本导读

1.1教材的选择

系统的学习模型教材肯定是比文献要更加全面的,对于时间序列的教材,推荐以下10本书,鉴于版权原因就不给大家贴链接了,需要的单独索要:QQ: 1758714024
在这里插入图片描述
选择着本书的原因仅仅是学的时候用的这本书1。理论部分推荐大家看王燕老师编著的第四版。那本基于sas,理论相对这本有所修正。

在这里插入图片描述

在系统的介绍之前,先给大家介绍一下这本书,简单导读一下。方便没有用过的童鞋能快速熟悉脉络。

1.2导读

本书共分为六章:
第一章是从定义和发展出发介绍一下时间序列模型以及基础的R语言基础。熟悉过程和R语言的同学可以直接跳过或者简单一了解。
第二章则是时间序列的预处理。主要是介绍了平稳性和纯随机性的定义和检验工具。
第三章是本书学习重点,从理论、建模、估参、检验、优化、预测六部分介绍了平稳性模型ARIMA。多读这一章有利于理解模型脉络和原理。
第四章第五章是非平稳的处理手段:第四章是确定性因素介绍了重要的分解定理;第五章是随机性分析,主要从模型假设出发对ARIMA进行了一些必要推广。
第六章介绍了多元时间序列以及两个重要的检验,单位根和协整。

2、建模思路流程图

Created with Raphaël 2.2.0 获得数据 预处理是否平稳 计算自相关、偏自相关 建立模型形式 估计参数 模型参数检验 对模型进行优化 预测解释 yes no yes no yes no

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、当然经常用的还会有统计量的检验方式,adadf检验,也称为单位根检验2R语言中做单位的函数有:fUnitRoots包中的UnitrootTestsadfTest。2、用tseries包中的adf.testpp.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模型的参数

6、参考文献


  1. 王燕.时间序列分析:基于R.中国人民大学出版社.2015.3. ↩︎ ↩︎

  2. https://baike.baidu.com/item/单位根检验/5574482?fr=aladdin ↩︎

发布了14 篇原创文章 · 获赞 49 · 访问量 3万+

猜你喜欢

转载自blog.csdn.net/weixin_46111814/article/details/105348265