毕业实用模型（三）——时间序列forecast包的使用 ——————————————...

zhangtt123

00引言
2 I# M: v! y0 S3 b) N0 c& }在毕业实用统计模型（一）——时间序列1中介绍了时间序列的模型的基本建模思路，在毕业设计实用模型（二）——时间序列之SARIMA2中说明了ARMA、ARIMA、SARIMA三大模型的关系与参数调整。但是在实际建模中往往会遇到更加实际的需求：模型的评估、参数检验、预测图的展示、模型参数的调整。对于很对不喜欢编程的人来说很是痛苦。本文将会重点从上述几个方面讲述forecast包里的主要的函数，并给出实例。从本文中你将会学会：

更加高效的的模型预测图。
. o/ k; n+ t/ I8 I, Z+ h' i模型的参数检验
, @0 |* k" o; `2 ^$ V模型定阶的函数
模型的评估函数
% h8 e! T7 T9 z, t1 f拟合线性的模型函数
! Z# w; G( E8 W* S8 V相关图绘画，以及误差的标注
) P6 j. h* S& z5 x3 w  K, ]输出模型的预测误差
, \2 v% {+ ]- _, J! b- R注：本文部分代码案例来自forecast包，大家有疑问可以自行去查找3.如果进不去可以运行下面代码会找到。
library(forecast)
) h) ]0 Z% T4 |! b# N; `- s2 shelp(package = forecast)
% l/ g' e4 `8 |( X# v" U+ [5 S
2 B% v7 C' {5 S, o* D0 z: x1、accuracy函数
描述：输入参数是模型，输出下面的结果：
: E* N$ c  @4 s3 \" T# X' o8 \
7 C, y4 j. l7 b8 d3 y

/ ^: v0 ?3 f6 D函数示例：

0 z/ h$ }! E$ x; @, W

. ~0 O7 ^4 O; o& C, f最后的图片：

8 M3 {, v2 e* j2 U: C' Q
9 Y) v$ R( M. W( z" r2、Acf、Pacf、taperedacf、taperedpacf
1 E; m6 [6 D/ X+ W9 p这四个函数分别是自相关函数、偏自相关函数、带有误差的自相关函数、带有误差的偏自相关函数。前两个很熟悉和内置的acf、pacf一样，下面只给出后两个的示例：



- F$ k5 V7 y2 h0 h2 d3、arfima
可以建立长期记忆的时间序列模型。
" ?7 [, ^4 X8 `3 ?, S直接给例子了哦：
+ W0 K/ k' _" z
+ p8 H% Y2 W/ x  W' G3 t
输出模型效果：
, n3 M; L2 z  o1 o

. U7 c3 T1 |9 `$ }画出残差信息：
tsdisplay(residuals(fit))
! W% B7 J( n# Y4 ^; ?5 d. H
4 s7 ?& p1 W/ c7 w2 W
+ b9 u4 x2 Z: c: D4 t4、Arima函数
函数介绍：这个函数可以拟合平稳或者非平稳的且已经知道参数的时间序列模型。也可以带有季节因素。
6 e1 X+ P5 t& I- k* J* _+ Z0 M6 R
3 [8 y4 \% _' l( N$ o; [n = 50set.seed(0)x = rnorm(n, 1,3)  # 生成数据library(ggplot2)  # 载入画图包x %>%  Arima(order=c(3,1,1)) %>%  forecast(h=5) %>%  autoplot

$ W7 e4 p( [4 \# G

上述代码用了管道函数，ggplot2包。对模型从数据到预测到建模一部到位。给大家看看Arima的参数。
" H  O' P- }3 e1 t7 u% s! [$ a
- b1 t. o* ~$ ]! e. Hfunction (y, order = c(0, 0, 0), seasonal = c(0, 0, 0), xreg = NULL,
    include.mean = TRUE, include.drift = FALSE, include.constant,
+ I1 F  v' G  W2 a0 j2 D    lambda = model$lambda, biasadj = FALSE, method = c("CSS-ML",
0 z% ^/ N" L' V- _2 b/ X. u" Z        "ML", "CSS"), model = NULL, x = y, ...)
. c! v' |. v# r5 a! K
/ T) S/ W/ M: k9 s5、arima.errors函数
该函数使用简单，用于输出模型预测每一期的误差。注意和residuals函数进行区分。
# 生成数据建模n = 50set.seed(0)x = rnorm(n, 1,3)fit <- arfima(x)
2 `) o, u; J; o  J! e# 计算结果
2 \, W) D* K, r2 O1 \; ?; ?4 o0 r> arima.errors(fit)  #  预测误差
9 Z1 Q; V1 [1 {Deprecated, use residuals.Arima(object, type='regression') instead
3 {; v; V$ B0 d4 q' JTime Series:
Start = 1
2 ~) s7 d) r! z$ |9 YEnd = 50
Frequency = 1
[1]  4.78886285  0.02129992  4.98939779  4.81728796  2.24392430 -3.61985013
[7] -1.78570110  0.11583866  0.98269848  8.21396017  3.29078038 -1.39702775
[13] -2.44297103  0.13161528  0.10235465 -0.23453250  1.75667034 -1.67576338
# R/ X0 s/ \5 W: v7 @( d[19]  2.30704990 -2.71261527  0.32719634  2.13218694  1.40000908  3.41256853
[25]  0.82867968  2.51082392  4.25730809 -1.07286152 -2.85379806  1.14017852
[31]  0.29288033 -0.62866477 -0.29993095 -0.94841494  3.18025224  4.45573526
2 |" s' a4 T7 n[37]  3.97648110 -0.28853933  4.71491230  0.16196115  6.27370927  2.68223827
9 B- I4 q4 F2 R  z[43] -0.35835192 -1.49612989 -2.49971164 -2.19677174 -3.69134615  4.46961099
: d4 M( s6 }) F4 A$ V' f) k[49]  3.49614139  0.31801393
) G, }, x/ X& }+ v> residuals(fit)  # 残差
Time Series:
& v; p3 y9 B- i. }) R/ B0 IStart = 1
End = 50
9 \7 t* O9 [, Q! b4 NFrequency = 1
1 F& w* u; D& O% X8 H/ y/ s [1]  3.71706994 -1.28784846  3.87358515  3.45503095  0.78349981 -4.98056775
, C- Y" R) X/ l1 X. Q/ P0 d [7] -2.74303846 -0.77184435  0.03669199  7.20508533  1.79698739 -2.80377609
& B% n0 Y% f, g7 P/ V0 |. M[13] -3.54880551 -0.77827016 -0.86325716 -1.20139553  0.81773766 -2.72392791
[19]  1.43279254 -3.76501026 -0.47528735  1.24438107  0.37471968  2.37151650
5 M3 t6 Y% L) q  E2 @' O[25] -0.35743144  1.41677218  3.08432283 -2.39421040 -3.91173996  0.30272726
[31] -0.68356344 -1.59285606 -1.19943417 -1.83611913  2.34677217  3.38917930
[37]  2.72730075 -1.60715574  3.61473632 -1.17228009  5.12889059  1.21871171
3 p# ?8 ~4 b8 k) d6 ^7 D1 \; w[43] -1.73495615 -2.66242875 -3.50034594 -3.04300026 -4.46595111  3.84607121
[49]  2.44020165 -0.85475604

6、arimaorder
函数功能：输出函数的阶数，用于时间序列的自动定阶。下面用自动定阶函数auto.arima和管道函数结合输出参数。给出示例。
0 N9 ?& v, B2 c2 {$ w6 s( M2 [
6 w( X% [$ Y& Sn = 50set.seed(100)x = rnorm(n, 3, 5) + 1:50  # 生成数据x %>% auto.arima %>% arimaorderp d q 3 1 0
可以配合arima函数进行模型的参数的确定。

* s- ^. f& g+ ^3 ?1 v2 L/ Y7、auto.arima
& c- u* u0 h' t3 c, p函数功能：可以使用函数进行自动定阶。下面时该函数的参数。这里会进行讲解并给出具体的实例：
auto.arima(
- w- Y9 b; R* [" T  y, # 数据
  d = NA,  # 原始数据差分的阶数
  D = NA, # 季节因素的差分阶数
, v9 s5 W$ i8 Y$ K4 t  max.p = 5,  # 遍历达到的最大的AR模型的阶数。
$ g: g1 t* {/ P- {: n  max.q = 5,  # 遍历达到的最大的MA模型的阶数。
  max.P = 2,  # 季节因素遍历达到的最大的AR模型的阶数。
  max.Q = 2,  # 季节因素遍历达到的最大的MA模型的阶数。
# N' }& F  W( _( M' h9 N# k  max.order = 5,
  max.d = 2,
  max.D = 1,
  start.p = 2,  # 开始遍历的AR模型的阶数。
  start.q = 2,  # 开始遍历的MA模型的阶数。
* h/ }) z% s4 V2 K; o  start.P = 1,
  start.Q = 1,
* ?( z! A2 K  V7 \) I  stationary = FALSE,  #
  seasonal = TRUE,
! |7 M* C. ?* [0 z5 z* l* D  ic = c("aicc", "aic", "bic"),  # 常用的选择模型的标准
  stepwise = TRUE,
* U# f& s. O$ ]+ |  nmodels = 94,
" ~, T6 |! b4 O% s& U7 A+ u  trace = FALSE,# 是否跟踪模型的选择
5 w8 _7 \5 j3 Z1 U9 r, h! z  approximation = (length(x) > 150 | frequency(x) > 12),
  method = NULL,
  D' x9 m* z9 C, \" t  truncate = NULL,
. U8 |6 e7 D0 H* X9 p  xreg = NULL,
, s/ g2 f1 x: }) g! z  test = c("kpss", "adf", "pp"),  # 平稳性检验的方式
7 s! C5 t# Q. b$ s7 Z8 y  test.args = list(),
& n- W! M0 b/ x7 s- {1 v# f  seasonal.test = c("seas", "ocsb", "hegy", "ch"),  # 季节性分析
" l1 }. u$ I  \& \( M  seasonal.test.args = list(),
  allowdrift = TRUE,  # 可以去掉漂移项
  allowmean = TRUE,  # 可以去掉均值项
) F6 q8 `# X( B3 z( X) ~3 r  lambda = NULL,
  biasadj = FALSE,
  parallel = FALSE,
$ v- l' `3 p- t, Y6 k# E" j& ?/ c  num.cores = 2,
2 e* u" ^( Z, \( Z  x = y,
5 A2 w' c9 L$ _8 O" H$ G  ...
" {' A! ], N7 O( ]% ~% R)

) H1 j' u8 c6 T' bauto.arima函数的参数众多，我把上述函数的参数给了一定的备注。大家可以根据这个注释和数据的其他的检验依据自己去逐一的验证调试。最终得出一个比较满意的模型。
下面给出部分调参的实例：
5 c3 N  H  B7 M, a# N. K7 T- v# 画时序图plot(ts(x))
- K- K  o9 ]+ ~! v$ V( i+ G

! r! v7 u2 T6 b1 A0 |" k下面展示追踪模型的参数：


$ R1 |& ?+ l0 f3 M7 Z6 \1 u5 b& B更改准则AICc成AIC继续追踪模型的参数：
" ]1 q" K" ?6 n7 }; I; H$ z
' p! w+ c" e4 O& Y
  s- P0 X4 Z6 t. d( p. G0 i* X这个函数的参数就介绍到这里，大家用到自己继续探索。
& l# t# H% F: P
; W. j; v$ w. q5 V3 g8、ggplot2中的时间序列相关系数
这里把函数内置的实例给出来。

1 J0 C; {) a' ]6 Y9 Z5 n. N2 f# 载入画图的ggplot2包
7 H- S3 C7 L6 L) U9 ~/ w- Vlibrary(ggplot2)
1 G. V( H2 M: v- U; L* _4 b& c
# v$ M% E2 K; G8 i下面的例子使用winneind数据。是个时间序列数据。

8 D8 i  i5 c& [$ h6 I1 ggAcf(wineind)
) E6 q; t; v3 C& `; |" _
/ h8 v2 P+ v6 l. A) Y5 w1 g

1 wineind %>% taperedacf(plot=FALSE) %>% autoplot

: @6 j( h* l# c. V6 K; x( Q( K8 i
. n8 y7 i7 P+ z; {
  J+ _6 u8 V+ p- M1 ggCcf(mdeaths, fdeaths)

1 u$ R  F6 C* S: N3 b
/ B) J  c" ?0 ^- U+ v
' X* _, `6 Z8 b6 {% W) V9、tslm
+ d- g! M, r- t  E9 T函数功能：用时间序列分量拟合线性模型
为了我方便第一次见到的同学学习，先贴出模型的参数：

7 Q) R$ S! j6 H; x& z1 > tslm
2 function (formula, data, subset, lambda = NULL, biasadj = FALSE,
3    ...)
0 K/ O3 F" h" D' z: \, m; C( `0 o
formula参数是公式、data是数据。下面使用函数例子展示函数功能。

1  > y <- ts(rnorm(120,0,3) + 20*sin(2*pi*(1:120)/12), frequency=12) # 构造数据
2  > fit1 <- tslm(y ~ trend + season)  # 趋势+季节
" \8 e; u3 v  n% [. x+ a3  > fit1
- ?$ z0 }$ \" l- m& I7 S- F4  Call:
5  tslm(formula = y ~ trend + season)
/ ^0 ?3 U! I0 d( `6 D" |4 A# I! G6  Coefficients:
9 T: m) i( X' C# j7 O3 n7  (Intercept)        trend      season2      season3      season4  
8     10.20906      0.01175      7.31274      8.82928      7.04245  
9      season5      season6      season7      season8      season9  
10   -2.04022    -10.24111    -20.57768    -27.61414    -30.16768  
11   season10     season11     season12  
12  -27.56840    -21.53062    -10.35272  
5 ?5 L( w0 k% i+ T& T: h- h
% Y7 t, D% i2 Z& }1 plot(forecast(fit1, h=20))  # 画出模型1的预测图


- L. T$ p) S7 z% s& _1  > fit2 <- tslm(y ~ season)
& ?. U# l; m( {2  > fit2
3  Call:
4  tslm(formula = y ~ season)
5  Coefficients:
' ^1 e  u; \# y( f. K: S6  (Intercept)      season2      season3      season4      season5  
7     10.856        7.324        8.853        7.078       -1.993  
8    season6      season7      season8      season9     season10  
9    -10.182      -20.507      -27.532      -30.074      -27.463  
10   season11     season12  
/ z% x9 y) o: a/ k9 _0 A- F& X11    -21.413      -10.223
- V5 w) L: o. ~7 J5 F
1  plot(forecast(fit2, h=20))  # 画出模型2的预测图
: o! l- |# d4 P" a' P9 Y: `" d
# T2 D: j7 _8 V7 }5 b1 o
10、CV交叉验证
& d  Z; K9 b5 T5 E) L4 d: QCV函数显示模型的CV AIC AICc BIC AdjR2值，用上述模型直接给出例子。

1  > CV(fit1)
0 o9 h* s9 ]4 M9 U2          CV        AIC       AICc        BIC      AdjR2
- S' e/ J0 t; Z/ B) @1 ?, N3   12.851627 306.718526 310.718526 345.743411   0.945285
6 E7 X2 W5 J1 y: j2 p9 D4  > CV(fit2)
5           CV         AIC        AICc         BIC       AdjR2
6   12.7985986 306.6337582 310.0677205 342.8711509   0.9449195
" M6 R: S$ z+ q& ^- J4 j% W0 q& T
11、forecast
2 Z% T9 F0 k9 [% D# Y这个函数可以给出模型的预测值用于画图和分析。
' d/ i4 }4 `7 B' ?6 I+ |下面给出上述模型一的预测值。
8 k* z% ~2 p. |0 C8 e
4 d9 o$ V$ `' ]8 Q
0 l( M1 R  B/ _& a. A结果说明：第一列是时间、第二列是预测值、后面4列是预测值的80%、95%的置信下限和置信下限。

12、geom_forecast
1 X" H; u  Z" d) T" L+ o这个图是ggplot2的预测图，下面贴出代码给出效果。这个函数也需要结合ggplot2包进行。
: v# u' M  s5 d7 S1 F. C* ^
" R* `* D6 q6 U1 w1  library(ggplot2)
2 v. _) l- Y1 P7 m2  autoplot(USAccDeaths) + geom_forecast()  # 图一
  K8 q) X1 w7 i9 ^3 R3  lungDeaths <- cbind(mdeaths, fdeaths)
4  autoplot(lungDeaths) + geom_forecast()  # 图二
' K2 i# w7 S0 c- K! t


13、总结
" ^4 r% M% B: E8 ]forecast包里的函数众多，上面只是介绍很少的一部分，更多函数的使用方式留着给大家探索。
- @; S% \# O/ \* ?% [8 b
- E5 W& G5 ]6 c) [. x7 o1 r" b14、参考文献
https://blog.csdn.net/weixin_46111814/article/details/105348265 ↩︎

https://blog.csdn.net/weixin_46111814/article/details/105370507 ↩︎
! o7 n5 A( B4 s: w% X8 |- H" A
http://127.0.0.1:17627/library/forecast/html/00Index.html ↩︎
————————————————
( y* a0 s* o  y( W  ]# N( A版权声明：本文为CSDN博主「逆天者顺A」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
+ n0 e7 c$ }2 k0 f原文链接：https://blog.csdn.net/weixin_46111814/article/details/105583080


7 P* B% R$ i. p7 T) w* D% ^

星辰的微光

点个赞，感谢分享
3 e  j: m4 A/ N+ Z% d: C5 k* X* v  t