QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 4309|回复: 1
打印 上一主题 下一主题

[个人总经验] 毕业实用模型(三)——时间序列forecast包的使用 ——————————————...

[复制链接]
字体大小: 正常 放大

326

主题

32

听众

1万

积分

  • TA的每日心情
    慵懒
    2020-7-12 09:52
  • 签到天数: 116 天

    [LV.6]常住居民II

    管理员

    群组2018教师培训(呼和浩

    群组2017-05-04 量化投资实

    群组2017“草原杯”夏令营

    群组2018美赛冲刺培训

    群组2017 田老师国赛冲刺课

    跳转到指定楼层
    1#
    发表于 2020-5-20 11:02 |只看该作者 |倒序浏览
    |招呼Ta 关注Ta
    00引言
    / z; H1 }4 g8 _0 \' \6 @在毕业实用统计模型(一)——时间序列1中介绍了时间序列的模型的基本建模思路,在毕业设计实用模型(二)——时间序列之SARIMA2中说明了ARMA、ARIMA、SARIMA三大模型的关系与参数调整。但是在实际建模中往往会遇到更加实际的需求:模型的评估、参数检验、预测图的展示、模型参数的调整。对于很对不喜欢编程的人来说很是痛苦。本文将会重点从上述几个方面讲述forecast包里的主要的函数,并给出实例。从本文中你将会学会:
    / a: {) [4 z4 o9 T
    . T6 k7 b- m) t0 y! f) ]  Y更加高效的的模型预测图。
    ; j/ h: H! ?; Y: ?. c模型的参数检验
    ; j8 K) }. Z  n+ w. Y) R4 n模型定阶的函数9 T* H, \' ^+ N
    模型的评估函数* @% X& \8 g* F3 P4 Z3 y% s
    拟合线性的模型函数
    0 ^* k' _% a! k/ b( p6 V相关图绘画,以及误差的标注
    3 A! P- S6 [- o2 S( K输出模型的预测误差& d, A& ^) N$ v# b. r! S( @8 m: R7 j1 [& U, K
    注:本文部分代码案例来自forecast包,大家有疑问可以自行去查找3.如果进不去可以运行下面代码会找到。
    7 m2 j! v- R' T" X7 flibrary(forecast)
    % D, Y# v2 n# J) y( J$ F! Rhelp(package = forecast)$ T5 Y" ]7 j  ~& A' R

    ) k5 x) Q3 _# N6 g7 S1、accuracy函数# o0 c8 l. C3 K/ h
    描述:输入参数是模型,输出下面的结果:' L  v( u$ W6 }4 T5 q( ~
    ! A9 C; r: B, r) _/ X! x3 V9 \3 T
    1.png
    " R: ^! B2 ^2 r; J( |& i3 F
    " y& }: m& N, s+ t+ {函数示例:
    6 I5 y- f( s% ^$ a- R" i5 w  {& e0 B; a6 r. z! |1 D
    2.png
    ( r8 ]  Q, X/ O1 S& I
    % G0 g, Z: s4 b, [4 Z: R/ n4 r最后的图片:
    + Y* K+ ?: F  n% l" y
    - V) C" G' s" j( b3 h& r% H$ |
    ; Z* Q5 q1 B  i7 ~  d# n2、Acf、Pacf、taperedacf、taperedpacf3 a6 K+ Y' ^: s+ n+ O
    这四个函数分别是自相关函数、偏自相关函数、带有误差的自相关函数、带有误差的偏自相关函数。前两个很熟悉和内置的acf、pacf一样,下面只给出后两个的示例:0 p: w$ S: J& ~  Z2 P$ `2 \8 m! p# T4 B
    3.png
    1 C9 q% x$ Z( ?9 {9 o: u
    ( R* m" t+ h. p$ a- B5 }9 f5 f% U" E( B7 d
    3、arfima7 R( @  l/ ^9 ?! w" f; E0 h
    可以建立长期记忆的时间序列模型。. p  j( v, \1 [% D5 c
    直接给例子了哦:( H5 p+ y: S( C$ B

    - b0 R8 h0 N! E( w8 q& D: y 4.png
    5 ?5 r. x/ u1 E9 T4 F9 W$ ?输出模型效果:
    ' g( g1 A' T9 j1 f 5.png # ?4 y: c: i$ ~

    7 T  _/ \# I  t/ Y/ @& z( H: P画出残差信息:
    + `7 I, K0 L: Q. Stsdisplay(residuals(fit))
    - F+ M: `# G9 d, W9 `0 q, r( _. h: i2 h8 `" c7 B5 e  H' P+ h
    6 P$ ^# ~* w1 ?2 a+ H% c
    4、Arima函数
    5 y/ P/ X6 h3 T: b函数介绍:这个函数可以拟合平稳或者非平稳的且已经知道参数的时间序列模型。也可以带有季节因素。/ j% S' I# i; N1 ?' d- |
    - s+ Z. U6 X5 B. S2 F
    n = 50set.seed(0)x = rnorm(n, 1,3)  # 生成数据library(ggplot2)  # 载入画图包x %>%  Arima(order=c(3,1,1)) %>%  forecast(h=5) %>%  autoplot

      ; Q! m( r3 F* T: o+ P1 W( \( J5 ]5 ]! d' ^
    上述代码用了管道函数,ggplot2包。对模型从数据到预测到建模一部到位。给大家看看Arima的参数。8 W5 K+ O3 f  R$ e' B* X

    / Z: U( }; A2 F3 I2 Pfunction (y, order = c(0, 0, 0), seasonal = c(0, 0, 0), xreg = NULL,
    3 @) N. B) j8 [6 b8 b9 y    include.mean = TRUE, include.drift = FALSE, include.constant, 4 Q& @, }/ K' m+ f8 D( @5 B
        lambda = model$lambda, biasadj = FALSE, method = c("CSS-ML",
    9 d8 j( n- J1 H9 z* |+ L        "ML", "CSS"), model = NULL, x = y, ...) ; a0 n4 o1 J+ U: n
    5 f- E' w& P3 w; `# q; U4 X5 p( k
    5、arima.errors函数
    $ \4 k; [: w" A) s# Y该函数使用简单,用于输出模型预测每一期的误差。注意和residuals函数进行区分。: m' m  P5 a% \. c- r& i* f
    # 生成数据建模n = 50set.seed(0)x = rnorm(n, 1,3)fit <- arfima(x), c. o& z/ d9 ~3 u8 o
    # 计算结果
    . d. ]1 E6 H: D/ U6 U2 H$ h! ]> arima.errors(fit)  #  预测误差
    1 e1 L: V1 ?; ^Deprecated, use residuals.Arima(object, type='regression') instead* v7 a; q( t) u+ \
    Time Series:
    / m, J# q7 Y' }0 pStart = 1 & v0 {* x4 Z8 N4 A% V
    End = 50   W: @* y$ [- X4 _. R9 L
    Frequency = 1
    $ h" N8 S5 L. n* g [1]  4.78886285  0.02129992  4.98939779  4.81728796  2.24392430 -3.61985013, U; @1 w6 z4 m9 Z
    [7] -1.78570110  0.11583866  0.98269848  8.21396017  3.29078038 -1.39702775  W- D: L& W. s0 I  ]
    [13] -2.44297103  0.13161528  0.10235465 -0.23453250  1.75667034 -1.67576338! a& a$ e& I' U( s  R
    [19]  2.30704990 -2.71261527  0.32719634  2.13218694  1.40000908  3.412568537 w2 c1 y5 G" Q; L0 U
    [25]  0.82867968  2.51082392  4.25730809 -1.07286152 -2.85379806  1.14017852
      ^7 \8 k( \, E[31]  0.29288033 -0.62866477 -0.29993095 -0.94841494  3.18025224  4.45573526: B& H0 I2 @; i! n- i
    [37]  3.97648110 -0.28853933  4.71491230  0.16196115  6.27370927  2.68223827
    6 i3 Y* m' M- _( X! X, T[43] -0.35835192 -1.49612989 -2.49971164 -2.19677174 -3.69134615  4.46961099
    ! n2 ^: N9 k6 S$ l, s[49]  3.49614139  0.31801393
    % H' Q1 m+ y- E: k> residuals(fit)  # 残差4 D7 Q2 u3 n6 w* n' G
    Time Series:
    0 [, j: @! {( M+ _* ]Start = 1
    ) B* l- m0 F8 Y! ~; x* J7 r, lEnd = 50
    - d5 r" p! q8 N" N2 ]9 C, c( PFrequency = 1
    9 q2 j6 z' y, b [1]  3.71706994 -1.28784846  3.87358515  3.45503095  0.78349981 -4.98056775
    * W$ @, S; U/ a) _. l [7] -2.74303846 -0.77184435  0.03669199  7.20508533  1.79698739 -2.803776093 G# m; k( o; [  |' N/ g
    [13] -3.54880551 -0.77827016 -0.86325716 -1.20139553  0.81773766 -2.72392791) ~2 Z0 `. m) j6 @3 U8 q
    [19]  1.43279254 -3.76501026 -0.47528735  1.24438107  0.37471968  2.37151650
    3 x9 |' t  p$ ]& \; ]0 C" Y[25] -0.35743144  1.41677218  3.08432283 -2.39421040 -3.91173996  0.30272726# [* E/ j( r, T; ~& {
    [31] -0.68356344 -1.59285606 -1.19943417 -1.83611913  2.34677217  3.389179301 A. R8 i4 {# z3 n
    [37]  2.72730075 -1.60715574  3.61473632 -1.17228009  5.12889059  1.21871171! w+ c3 C& G- m
    [43] -1.73495615 -2.66242875 -3.50034594 -3.04300026 -4.46595111  3.846071217 X9 v+ o) S% X8 N( \4 S
    [49]  2.44020165 -0.85475604. _- ~6 Z# b' d: c4 }9 k7 ?

    ! h- l. t5 }4 t' P7 \% p6、arimaorder& V* m% r. p( p# k- k
    函数功能:输出函数的阶数,用于时间序列的自动定阶。下面用自动定阶函数auto.arima和管道函数结合输出参数。给出示例。& r; \! n7 g1 S% H

    ' b$ {7 I7 P8 |2 y1 \! x$ Cn = 50set.seed(100)x = rnorm(n, 3, 5) + 1:50  # 生成数据x %>% auto.arima %>% arimaorderp d q 3 1 0 3 ?2 a) }0 v9 r8 M
    可以配合arima函数进行模型的参数的确定。& |8 j" Y& I9 ], D

    ( g1 t* o) Q# z2 ?# ~2 Z' x7、auto.arima$ i# B! H" [+ @; z, i1 l; u) V
    函数功能:可以使用函数进行自动定阶。下面时该函数的参数。这里会进行讲解并给出具体的实例:
    ( v% G' w& O, t9 vauto.arima(
    9 ~0 B, z0 P. m+ ?* G6 G& D: J  y, # 数据- B  x) h- ~/ X4 y6 ~) v
      d = NA,  # 原始数据差分的阶数6 y4 I4 w4 u* N9 e1 j8 Y/ h
      D = NA, # 季节因素的差分阶数1 }; P( Z- C/ H3 b
      max.p = 5,  # 遍历达到的最大的AR模型的阶数。
    4 Y) G$ }7 ~- l2 S5 Y/ U  max.q = 5,  # 遍历达到的最大的MA模型的阶数。
    2 D) B' r: b. r& X% x  max.P = 2,  # 季节因素遍历达到的最大的AR模型的阶数。2 p2 h& [, k) U% j" R' D
      max.Q = 2,  # 季节因素遍历达到的最大的MA模型的阶数。5 Y& X- G% D) Z3 _5 d+ K* S
      max.order = 5,
    8 ^5 R; v+ b4 |  H( ?: S  max.d = 2,
    / `7 N8 Q- E, K+ a) B: X. o  max.D = 1,
    - |' I( j5 \+ X# H' H  start.p = 2,  # 开始遍历的AR模型的阶数。* T. g! U2 y: |& G4 @: L
      start.q = 2,  # 开始遍历的MA模型的阶数。1 G7 u9 A3 |' x) J- |) }
      start.P = 1,' D" q! l9 j1 i$ I, b; ^& p% o
      start.Q = 1,
    ! @7 p7 a& \4 e9 h8 j  R  stationary = FALSE,  # 4 p% n# T/ b/ _8 M1 R, z$ o( ^
      seasonal = TRUE,
    3 f3 J  a. C3 E' j( f  ic = c("aicc", "aic", "bic"),  # 常用的选择模型的标准
    1 d% G6 I0 _! {& y$ b% K+ a  stepwise = TRUE,
    ; R# `! C+ Q- ]- T5 N5 x  nmodels = 94,
    : e6 h8 a* w) |' ^  trace = FALSE,# 是否跟踪模型的选择. j+ k: a) f" N* i
      approximation = (length(x) > 150 | frequency(x) > 12),
    " X+ H( [5 J' z. S7 M0 \  method = NULL,
    1 p( u1 X/ M/ g" w; q6 y  truncate = NULL,! A: a3 U0 `+ l& n6 X6 O/ z  f
      xreg = NULL,
    8 R0 M0 n% y% u( O; H  w2 S, z  test = c("kpss", "adf", "pp"),  # 平稳性检验的方式6 y$ [1 ?0 w1 @  ^) o
      test.args = list(),
    2 h3 @9 n+ A7 q# Q( r  seasonal.test = c("seas", "ocsb", "hegy", "ch"),  # 季节性分析
    7 Z' }! x4 J; X1 B. y% b  seasonal.test.args = list(),
    ( T) s1 @! a# w! Z  allowdrift = TRUE,  # 可以去掉漂移项. o' I' t% r/ u
      allowmean = TRUE,  # 可以去掉均值项3 V/ C  L, ?. D# b
      lambda = NULL,9 b9 D/ Y; s, j$ |
      biasadj = FALSE,
    ; O7 A% h1 q/ y  parallel = FALSE,
    0 v9 v* k+ j; [: P  num.cores = 2,: ?$ ~# ~4 h) ~) j
      x = y,
    6 u8 j  ?# i% l  ...
    - I$ ~. r( w' F6 a! m)6 Q$ T. |1 S& ~' w7 U  U

      M  i+ C0 H5 @) O5 Z" kauto.arima函数的参数众多,我把上述函数的参数给了一定的备注。大家可以根据这个注释和数据的其他的检验依据自己去逐一的验证调试。最终得出一个比较满意的模型。
    8 P; r! i- R' c4 e; D. b" n- M下面给出部分调参的实例:2 p7 H# E) U) M+ j& c- g) Y; {
    # 画时序图plot(ts(x))+ G1 L' A! o- ~+ N" q! H

    ( X4 n! M" q/ {- @/ k+ \* j% B+ s: j" c" l8 F
    下面展示追踪模型的参数:7 A  }" j! P) ?" C
    7.png
    7 E/ @+ `; s" w  S( {. S
    5 b+ b% ~* O9 O7 n更改准则AICc成AIC继续追踪模型的参数:
    6 `3 M% z' ?: z9 P- u4 g% ^) B 8.png 3 L7 x( M5 g9 o. j, t
    ' o* L5 i( q' }
    这个函数的参数就介绍到这里,大家用到自己继续探索。8 r! H4 G# X7 @

    ) s! h, a4 ~3 k$ w! r$ z, S8、ggplot2中的时间序列相关系数  ~) |6 j0 R: ?( e; A5 f
    这里把函数内置的实例给出来。
    9 u- T8 z% L2 D4 ]. `# T
    0 c- n# Z8 D8 s# 载入画图的ggplot2包
    3 ]  G5 R. y# X& Q8 ylibrary(ggplot2)8 O+ I+ E  W( x3 L3 b- k

    " G4 @3 `8 e' R" _; w( [5 x下面的例子使用winneind数据。是个时间序列数据。4 L5 i* x% S  b) j

    # _; d6 n) {8 m6 i& u- }3 U1 ggAcf(wineind)
    0 g1 W; D  D1 I  O  f: ]: `! w  J; D) y% r& R8 x; z

    0 p' T+ ^( Y: x. e; r  o& O1 G) e, T4 j, b2 ~* {6 l: L2 _" k7 b
    1 wineind %>% taperedacf(plot=FALSE) %>% autoplot. z/ T' ~0 P9 S+ E6 v

    2 o3 x2 b8 R3 V1 f+ E# q1 L$ b
    5 s0 E0 ]7 O7 P  I! _' ]8 ^' K( g+ i  `3 ~6 `
    1 ggCcf(mdeaths, fdeaths)
    / K. o' _" w$ {" P8 a# P1 a% y7 i
    " |! P. R' }2 y6 [) v3 J- k( U
    5 x8 g. C4 ]$ Y7 d6 V7 R, d
    3 O/ }5 X) F/ g: c5 q) Q/ |- b9、tslm3 W% x% Z. P) W
    函数功能:用时间序列分量拟合线性模型# Y1 h- N, x4 r# Y1 D0 d
    为了我方便第一次见到的同学学习,先贴出模型的参数:* S% C  N; |- V1 b) ?

    ) D, u" q5 S! ^1 > tslm# W7 ?8 M+ y# h9 q  F9 [- I9 _
    2 function (formula, data, subset, lambda = NULL, biasadj = FALSE, ; f5 V2 @2 `! r5 O, p
    3    ...)
    ' _: q9 ], [: d& Y/ }7 A# J0 N% y( j) m: S: p6 k, W5 `
    formula参数是公式、data是数据。下面使用函数例子展示函数功能。' h" f8 G* D: h1 e) P7 L

    ) I& r  Z% A6 d' [9 p1  > y <- ts(rnorm(120,0,3) + 20*sin(2*pi*(1:120)/12), frequency=12) # 构造数据
    , }' G1 ?* X, ]1 x5 @& i, h2  > fit1 <- tslm(y ~ trend + season)  # 趋势+季节8 u7 G7 i' \2 n
    3  > fit1
    $ j: h  h4 x( ^9 @9 X/ L; Z4  Call:
    4 h+ u) W- J) a5  tslm(formula = y ~ trend + season)2 I- c) A4 q* Y  H) z
    6  Coefficients:& A: ~% V6 r+ g. g/ J. d  e
    7  (Intercept)        trend      season2      season3      season4  1 p4 b% |3 ~' B
    8     10.20906      0.01175      7.31274      8.82928      7.04245  - \' }8 O8 t, s! a
    9      season5      season6      season7      season8      season9  4 G; n. K3 r. Q# Y( R# O0 o. E6 H
    10   -2.04022    -10.24111    -20.57768    -27.61414    -30.16768  
    ! Z  m  h) ?) y3 i11   season10     season11     season12  
    / `; {) |1 D( t* R! W  C$ i) g12  -27.56840    -21.53062    -10.35272  
    ; X. w1 {* C/ z& U" s- a1 _( q6 \0 \' W0 F7 y
    1 plot(forecast(fit1, h=20))  # 画出模型1的预测图
    # v" H- j3 U. _* ]
    1 x7 x( Y2 H  S5 @9 f( j& D' N6 l" w8 e" a
      c& U0 ?1 P6 m4 s; l# ]1  > fit2 <- tslm(y ~ season)& j% z: V' r4 T+ F8 L
    2  > fit2
    ) n9 f- ?4 L& c, y* y3  Call:  H  T; ~! B3 `1 }
    4  tslm(formula = y ~ season). |! I& x  @+ v. t
    5  Coefficients:* I( e/ M- l# Q5 l
    6  (Intercept)      season2      season3      season4      season5  : C' h- ?* z2 g2 C( }6 `0 R
    7     10.856        7.324        8.853        7.078       -1.993  5 j6 X0 q6 k7 n0 j/ n% A
    8    season6      season7      season8      season9     season10  7 G4 s, U; r5 m. [) B' f  Z
    9    -10.182      -20.507      -27.532      -30.074      -27.463  4 x: m2 p  N8 S$ d
    10   season11     season12  4 l. s+ S; U4 H2 a
    11    -21.413      -10.223 + V0 W" w- V" n9 Z

    + k- _7 T8 V0 i, ?' q/ S1  plot(forecast(fit2, h=20))  # 画出模型2的预测图
    ) t5 Q4 q3 H5 t% i! v1 j) G! y4 L+ p; j
    1 b8 t, d8 Z) j
    10、CV交叉验证
    ( l" ~) W* h2 }& |2 ?& jCV函数显示模型的CV AIC AICc BIC AdjR2值,用上述模型直接给出例子。
    9 w+ B' e5 i2 {7 X9 U' x3 f% o
    - B* [" i5 n) Q: ]% i1  > CV(fit1)
    ) v2 r& c) j7 G6 M5 f. @2          CV        AIC       AICc        BIC      AdjR2
    * v( [  R  ?1 W( r/ ]3   12.851627 306.718526 310.718526 345.743411   0.945285 . T( U0 i( u* O) z) x, K
    4  > CV(fit2)0 v( l9 X6 k$ k# B
    5           CV         AIC        AICc         BIC       AdjR2
    / ~5 p* s" R+ _; X' M6   12.7985986 306.6337582 310.0677205 342.8711509   0.9449195 * G* K& Q' m6 S2 ^+ ^

    * f6 O$ Q. o! P. A- b$ v11、forecast2 l# D% C6 n  W; y* |9 c1 _4 `) `
    这个函数可以给出模型的预测值用于画图和分析。
    9 l9 A7 S" t6 U. z$ E! j" _, k下面给出上述模型一的预测值。7 }( q) `: A3 B
    " c: C0 U3 m- ?$ l- q) E
    8.png
    0 r. R4 G. |- d% h+ h结果说明:第一列是时间、第二列是预测值、后面4列是预测值的80%、95%的置信下限和置信下限。
    ( l1 D1 d0 p' V) [$ v9 R8 y3 a) i8 B+ `% i0 ^
    12、geom_forecast8 M$ g8 ~+ E8 c  S
    这个图是ggplot2的预测图,下面贴出代码给出效果。这个函数也需要结合ggplot2包进行。
    . }$ \3 W! |- y$ @4 ~8 k, \. E9 q6 Z) t" r
    1  library(ggplot2)7 K* \0 x& f* _9 R0 T8 h
    2  autoplot(USAccDeaths) + geom_forecast()  # 图一
    ' ^6 A3 F; m, w3  lungDeaths <- cbind(mdeaths, fdeaths)3 ^# U. F' j: y
    4  autoplot(lungDeaths) + geom_forecast()  # 图二
    " T* |9 k5 ?; q* q3 @. }" c4 y5 K0 e4 a4 d/ U: J' O

    & U4 O3 X& q) ]. u3 s1 C7 P6 y& c$ O, q0 _' |, h/ X' ]
    13、总结+ M  r/ k" z) X' H
    forecast包里的函数众多,上面只是介绍很少的一部分,更多函数的使用方式留着给大家探索。* H5 d' U4 @0 U4 v. H
    ; O' L4 q- ?9 a0 q; n
    14、参考文献% c1 a9 t' ^/ M$ Z
    https://blog.csdn.net/weixin_46111814/article/details/105348265 ↩︎
    ) [( t" D5 l1 k3 {# M+ i+ I) @% Y9 |( Y: `. j. g: q9 r4 @
    https://blog.csdn.net/weixin_46111814/article/details/105370507 ↩︎
    7 @5 E. G4 ]8 @" z
    / O6 Q1 V! O4 i6 Thttp://127.0.0.1:17627/library/forecast/html/00Index.html ↩︎
    3 _+ b7 {$ K; }9 c————————————————6 A1 Z) h0 i5 v8 x( A+ P5 x
    版权声明:本文为CSDN博主「逆天者顺A」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    " s. h: v9 p& I3 O; [. b' G原文链接:https://blog.csdn.net/weixin_46111814/article/details/1055830803 n# F2 ?8 n. }3 K( A. j, F

      ~9 n6 }0 {4 i- d
    : w! ]  f( L7 O5 F- f9 E8 i
    zan
    转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信

    0

    主题

    1

    听众

    6

    积分

    升级  1.05%

  • TA的每日心情
    开心
    2020-5-21 09:14
  • 签到天数: 1 天

    [LV.1]初来乍到

    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-4-16 20:31 , Processed in 0.470696 second(s), 59 queries .

    回顶部