时间序列分析：对非平稳时间序列进行建模

Posted 2021-05-02 R语言中文社区

时间序列分析：对非平稳时间序列进行建模

在这篇博客中，我将会简单的介绍一下时间序列分析及其应用。这里，我们将使用匹兹堡大学的教授David Stoffer所开发的R包astsa进行时间序列分析。而与之相关的课本，可以在 这里找到它的电子书，而且是免费的，它对于任何对这个话题感兴趣的读者都有一定的帮助。

时间序列就是一串基于具体时间区间的观察值。它在经济预测这块用有广泛的应用，而在预测未来一段时间的天气方面也有很广泛的应用。时间序列分析的本质就是利用一个具体的过往的观测值来预测未来的观测值。

首先，我们下载和加载一下astsa包。

 
   
   
 

  
    
    
  
install.packages("astsa")
  
    
    
  
require(astsa)
 
   
   
 

我们现在利用数据集gtemp数据集进行检验，它通过预测1880-2009年的气温变化，来预测1951-1980年的平均气温。

 
   
   
 

  
    
    
  
plot(gtemp, xlab="Year", ylab="Deviation", main="Global Temperature Deviations, 1880-2009")
 
   
   
 

得到gtemp图：

在建模之前，我们要检验一下这个时间序列是否平稳。如果一个时间序列是平稳的，它要满足三个条件：
1.常数均值稳定在t。
2.常数方差稳定在t。
3.两个观察值之间的自协方差仅由两个观测值之间的距离来决定，它可用log（h)来表示。
现实生活中，一般的数据通常都不完全满足上面的三个条件，除非这个数据集用与白噪声测量。从我们的肉眼来观察，gtemp里的时间序列是非平稳的。其均值是波动的，而且这个是很明显的上升趋势。不过，其方差就比较平稳了。
我们可以使用acf2()函数来进一步的检测它。如果它是平稳的，ACF/PACF图所显示的将会是在最前面的一小部分点之间的相关性小幅度波动以后就开始急剧的下降。

 
   
   
 

  
    
    
  
acf2(gtemp)
  
    
    
  

  
    
    
  
 ACF  PACF
  
    
    
  
 [1,] 0.88  0.88
  
    
    
  
 [2,] 0.82  0.23
  
    
    
  
 [3,] 0.79  0.14
  
    
    
  
 [4,] 0.78  0.17
  
    
    
  
 [5,] 0.72 -0.13
  
    
    
  
 [6,] 0.70  0.08
  
    
    
  
 [7,] 0.67  0.00
  
    
    
  
 [8,] 0.63 -0.07
  
    
    
  
 [9,] 0.58 -0.06
  
    
    
  
[10,] 0.58  0.14
  
    
    
  
[11,] 0.55 -0.07
  
    
    
  
[12,] 0.49 -0.13
  
    
    
  
[13,] 0.46  0.04
  
    
    
  
[14,] 0.46  0.11
  
    
    
  
[15,] 0.44 -0.02
  
    
    
  
[16,] 0.41  0.01
  
    
    
  
[17,] 0.40  0.09
  
    
    
  
[18,] 0.40  0.03
  
    
    
  
[19,] 0.36 -0.13
  
    
    
  
[20,] 0.34  0.00
  
    
    
  
[21,] 0.32 -0.03
  
    
    
  
[22,] 0.30 -0.03
 
   
   
 

这个蓝色的虚线描述了每个时间间隔的极值。显而易见的，这个时间序列是非平稳的，尤其是那些数值较大的过往观测值与未来的值是相关的。尽管回归方法允许给这个数据集的时间序列拟合一条光滑的曲线，时间序列所关注的就是除去尽可能多的趋势来确认回归线所抓取不到的信息的潜在因子。
对于所有的t值，我们将从Xt中减去Xt-1以实现数据序列从非平稳转为平稳。这种方法叫差分化，并且可通过diff()函数来实现。

 
   
   
 

  
    
    
  
 plot(diff(gtemp))
  
    
    
  
     acf2(diff(gtemp))
  
    
    
  

  
    
    
  
        ACF  PACF
  
    
    
  
 [1,] -0.29 -0.29
  
    
    
  
 [2,] -0.16 -0.28
  
    
    
  
 [3,] -0.12 -0.32
  
    
    
  
 [4,]  0.22  0.00
  
    
    
  
 [5,] -0.15 -0.20
  
    
    
  
 [6,]  0.02 -0.10
  
    
    
  
 [7,]  0.03 -0.04
  
    
    
  
 [8,]  0.11  0.05
  
    
    
  
 [9,] -0.20 -0.13
  
    
    
  
[10,]  0.15  0.09
  
    
    
  
[11,]  0.04  0.10
  
    
    
  
[12,] -0.07 -0.02
  
    
    
  
[13,] -0.17 -0.12
  
    
    
  
[14,]  0.15 -0.03
  
    
    
  
[15,]  0.06  0.01
  
    
    
  
[16,] -0.08 -0.07
  
    
    
  
[17,]  0.00  0.03
  
    
    
  
[18,]  0.14  0.10
  
    
    
  
[19,] -0.14 -0.06
  
    
    
  
[20,]  0.04  0.09
  
    
    
  
[21,]  0.00  0.03
  
    
    
  
[22,]  0.11  0.08
 
   
   
 

这看起来需要点技巧，这时，我们在1个时间间隔后面出去所有显著相关性。是时候使用sarima()函数来建立时间序列模型了。serima()函数有3个基本参数（q,d,p），它们分别表示自回归序列、差分度、移动平均序列。如果你对这些术语不熟悉，我建议你快速浏览这篇文章：。
ACF/PACF图告诉我们我们使用了多少参数度。如果ACF图是光滑的、几何衰减的而且PACF在log(p)处中止，我们应当使用纯AR(p)模型。给定ACF的模式没有展示出它是光滑衰减的，而衰减度也在1个区间以后就变得不再显著。于是，我们建立MA(1)模型来拟合这些数据。我们选择d=1作为我们的差分度，其原因则源于diff()函数。

 
   
   
 

  
    
    
  
sarima(gtemp,0,1,1)
  
    
    
  
initial  value -2.213853 
  
    
    
  
iter   2 value -2.283711
  
    
    
  
iter   3 value -2.285917
  
    
    
  
iter   4 value -2.299585
  
    
    
  
iter   5 value -2.304774
  
    
    
  
iter   6 value -2.320202
  
    
    
  
iter   7 value -2.323954
  
    
    
  
iter   8 value -2.327090
  
    
    
  
iter   9 value -2.328063
  
    
    
  
iter  10 value -2.328249
  
    
    
  
iter  11 value -2.328250
  
    
    
  
iter  12 value -2.328250
  
    
    
  
iter  13 value -2.328251
  
    
    
  
iter  14 value -2.328251
  
    
    
  
iter  14 value -2.328251
  
    
    
  
iter  14 value -2.328251
  
    
    
  
final  value -2.328251 
  
    
    
  
converged
  
    
    
  
initial  value -2.326014 
  
    
    
  
iter   2 value -2.326045
  
    
    
  
iter   3 value -2.326046
  
    
    
  
iter   4 value -2.326047
  
    
    
  
iter   5 value -2.326048
  
    
    
  
iter   5 value -2.326048
  
    
    
  
iter   5 value -2.326048
  
    
    
  
final  value -2.326048 
  
    
    
  
converged
  
    
    
  
$fit
  
    
    
  

  
    
    
  
Call:
  
    
    
  
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
  
    
    
  
    Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
  
    
    
  
    reltol = tol))
  
    
    
  

  
    
    
  
Coefficients:
  
    
    
  
          ma1  constant
  
    
    
  
      -0.6550    0.0063
  
    
    
  
s.e.   0.0792    0.0030
  
    
    
  

  
    
    
  
sigma^2 estimated as 0.0095:  log likelihood = 117.02,  aic = -228.03
  
    
    
  

  
    
    
  
$AIC
  
    
    
  
[1] -3.625671
  
    
    
  

  
    
    
  
$AICc
  
    
    
  
[1] -3.608821
  
    
    
  

  
    
    
  
$BIC
  
    
    
  
[1] -4.581555
 
   
   
 

从我们建立的模型，可以看到，这种拟合效果很不错。其残差在0左右徘徊而QQ图里，除了末尾各别数据点距离较大以外，整体的拟合效果还不错。我们现在对我们的ARIMA模型来预测2010-2015年这5年的大致趋势。蓝线标记了其中的标准误差。

 
   
   
 

  
    
    
  
sarima.for(gtemp,5,0,1,1)
  
    
    
  
$pred
  
    
    
  
Time Series:
  
    
    
  
Start = 2010 
  
    
    
  
End = 2014 
  
    
    
  
Frequency = 1 
  
    
    
  
[1] 0.5467731 0.5531159 0.5594588 0.5658016 0.5721445
  
    
    
  

  
    
    
  
$se
  
    
    
  
Time Series:
  
    
    
  
Start = 2010 
  
    
    
  
End = 2014 
  
    
    
  
Frequency = 1 
  
    
    
  
[1] 0.09746909 0.10310668 0.10845160 0.11354521 0.11841992
 
   
   
 

模型的拟合效果相当不错。
这个模型的状态转换给原始数据集gtemp里的一些趋势进行了中和，研究者们则通过分析先前没注意到的一个成分来增强模型的预测能力。

严禁修改，可以转载，请注明出自和原文链接。