拓端tecdat|R语言代写ARMA，ARIMA（Box-Jenkins），SARIMA和ARIMAX模型用于预测时间序列数据

 原文链接:http://tecdat.cn/?p=5919

 

在本文中，我将介绍ARMA，ARIMA（Box-Jenkins），SARIMA和ARIMAX模型如何用于预测给定的时间序列数据。

 

使用后移运算符计算滞后差异

我们可以使用backshift运算符来执行计算。例如，后轴运算符可用于计算的时间序列值的滞后差异ÿy经由yi−Bk(yi),∀i∈k+1,…,tyi−Bk(yi),∀i∈k+1,…,t其中kk表示的差异滞后。对于k=1k=1，我们获得普通的成对差异，而对于k=2k=2我们获得相对于前任先前的成对差异。让我们考虑R中的一个例子。

使用R，我们可以使用diff函数计算滞后差异。函数的第二个参数表示所需的滞后kk，默认设置为k=1k=1。例如：

 
By <- diff(y) <span style="color:#888888"># y_i - B y_i</span>
B3y <- diff(y, <span style="color:#880000">3</span>) <span style="color:#888888"># y_i - B^3 y_i</span>
message(paste0(<span style="color:#880000">"y is: "</span>, paste(y, collapse = <span style="color: 
## y is: 1,3,5,10,20
## By is: 2,2,5,10
## B^3y is: 9,17

 自相关函数

 

要计算自相关，我们可以使用以下R函数：

<span style="color:#000000"><span style="color:#000000"><code>get_autocor <- <strong>function</strong>(x, lag) {
    x.left <- x[<span style="color:#880000">1</span>:(length(x) - lag)]
    x.right <- x[(<span style="color:#880000">1</span>+lag):(length(x))]
    autocor <- cor(x.left, x.right)
    <strong>return</strong>(autocor)
}</code></span></span>

 

<span style="color:#000000"><span style="color:#000000"><code><span style="color:#888888"># correlation of measurements 1 time point apart (lag 1)</span>
get_autocor(y, <span style="color:#880000">1</span>) </code></span></span>
## [1] 0.9944627
<span style="color:#000000"><span style="color:#000000"><code><span style="color:#888888"># correlation of measurements 2 time points apart (lag 2)</span>
get_autocor(y, <span style="color:#880000">2</span>)</code></span></span>
## [1] 0.9819805

数据的高度自相关表明数据具有明确的时间趋势。

部分自相关

由于观察到较大滞后的自相关可以是较低滞后的相关结果，因此通常值得考虑部分自相关函数（pACF）。pACF的想法是计算部分相关性，这种相关性决定了对变量的最近观察的相关性。pACF定义为：

 

φkk:=Corr(yt,yt−k|yt−1,⋯,yt−k+1)k=0,1,2,⋯φkk:=Corr⁡(yt,yt−k|yt−1,⋯,yt−k+1)k=0,1,2,⋯

 

使用pACF可以识别是否存在实际滞后的自相关或这些自相关是否是由其他测量引起的。

计算和绘制ACF和pACF的最简单方法是分别使用acf和pacf函数：

<span style="color:#000000"><span style="color:#000000"><code>par(mfrow = c(<span style="color:#880000">1</span>,<span style="color:#880000">2</span>))
acf(y) <span style="color:#888888"># conventional ACF</span>
pacf(y) <span style="color:#888888"># pACF</span></code></span></span>

​

在ACF可视化中，ACF或pACF被绘制为滞后的函数。指示的水平蓝色虚线表示自相关显着的水平。

分解时间序列数据

  •  StSt
  • TtTt
  • ϵtϵt

执行分解的方式取决于时间序列数据是加法还是乘法。

加法和乘法时间序列数据

加法模型假设数据可以分解为

 

yt = St + Tt + ϵt.yt = St + Tt + ϵt.

 

另一方面，乘法模型假设数据可以被分解为

 

 

  • 添加剂：每个时期的季节效应放大器相似。
  • 乘法：季节性趋势随时间序列的变化而变化。

AirPassengers数据集提供了乘法时间序列的示例。

<span style="color:#000000"><span style="color:#000000"><code>data(AirPassengers)
plot(AirPassengers)</code></span></span>

​

log(StTtϵt)=log(St)+log(Tt)+log(ϵt)log⁡(StTtϵt)=log⁡(St)+log⁡(Tt)+log⁡(ϵt)AirPassengers 数据集：

<span style="color:#000000"><span style="color:#000000"><code>plot(log(AirPassengers))</code></span></span>

​

正如我们所看到的，采用对数已经使季节性成分的幅度沿时间均衡。请注意，总体增长趋势没有改变。

在R中分解时间序列数据

要分解R中的时间序列数据，我们可以使用该decompose函数。请注意，我们应该通过type参数提供时间序列是加法的还是乘法的。

示例1：AirPassengers数据集

对于AirPassengers数据集，我们指定数据是乘法的并获得以下分解：

<span style="color:#000000"><span style="color:#000000"><code>plot(decompose(AirPassengers, type = <span style="color:#880000">"multiplicative"</span>))</code></span></span>

​

分解表明，多年来航空公司乘客总数在增加。此外，我们已经观察到的季节性影响已被清楚地捕捉到。

示例2：EuStockMarkets数据集

让我们考虑可以为EuStockMarkets数据集找到的分解：

<span style=" /span>] <span style="color:#888888"># DAX data</span>
<span style="color:#888888"># data do not seem to be multiplicative, use additive decomposition</span>
decomposed <- decompose(daxData, type = <

​

该图显示了1992年至1998年的DAX数据中的以下内容：

  • 整体价值稳步上升。
  • 季节性趋势强烈：每年年初，股价相对较低，并在夏季结束时达到相对最大值。
  • 除1997年和1998年之间的最终测量外，随机噪声的贡献可以忽略不计。

固定与非固定过程

生成时间序列数据的过程可以是静止的也可以是非静止的。 例如，数据EuStockMarkets和AirPassengers数据都是非平稳的，因为数据有增加的趋势。为了更好地区分固定和非固定过程，请考虑以下示例：

<span style="color "># climate data </span>
<strong>library</strong>(tseries)
data(nino)

​

左图显示了一个静止过程，其中数据在所有测量中表现相似。右图显示了一个非平稳过程，其中平均值随着时间的推移而增加。

介绍了与时间序列数据分析相关的最重要概念后，我们现在可以开始研究预测模型。

ARMA模型

ARMA代表自回归移动平均线。ARMA模型仅适用于固定过程，并具有两个参数：

  • p：自回归（AR）模型的顺序
  • q：移动平均（MA）模型的顺序

ARMA模型可以指定为

 

ÿ^Ť= c + εŤ+ Σi = 1pφ一世ÿt - 我- Σj = 1qθĴεt - j。y^t=c+ϵt+∑i=1pϕiyt−i−∑j=1qθjϵt−j.

 

使用以下变量：

  • cc
  • ϵtϵtttϵt∼N(0,σ2)ϵt∼N(0,σ2)
  • ϕ∈Rpϕ∈Rp
  • ytyttt
  • θ∈Rqθ∈Rq
  • ϵtϵttt

使用backshift运算符制定ARMA模型

使用backshift运算符，我们可以通过以下方式制定ARMA模型：

 

（ 1 - Σi = 1pφ一世乙一世）yŤ= （ 1 - Σj = 1qθĴ乙Ĵ）εĴ(1−∑i=1pϕiBi)yt=(1−∑j=1qθjBj)ϵj

 

ϕp(B)=1−∑pi=1ϕiBiϕp(B)=1−∑i=1pϕiBiθq(B)=1−∑qj=1θjBjθq(B)=1−∑j=1qθjBj

 

ϕp(B)yt=θq(B)ϵt.ϕp(B)yt=θq(B)ϵt.

 

ARIMA模型

dd

总之，ARIMA模型具有以下三个参数：

  • p：自回归（AR）模型的顺序
  • d：差异程度
  • q：移动平均（MA）模型的顺序

在ARIMA模型中，通过将替换差异，将结果转换为差异ytyt

 

(1−B)dyt.(1−B)dyt.

 

然后通过指定模型

 

ϕp(B)(1−B)dyt=θq(B)ϵt.ϕp(B)(1−B)dyt=θq(B)ϵt.

 

d=0d=0(1−B)0yt=yt(1−B)0yt=ytdd

 

(1−B)1yt(1−B)2yt=yt−yt−1=(1−2B+B2)yt=yt−2yt−1+yt−2(1−B)1yt=yt−yt−1(1−B)2yt=(1−2B+B2)yt=yt−2yt−1+yt−2

 

在下文中，让我们考虑ARIMA模型的三个参数的解释。

pp

p∈N0p∈N0d=0d=0Byt=yt−1Byt=yt−1ϕ1ϕ1yt−1yt−1yt−2yt−2ϕ1ϕ1ϕ2ϕ2

p=1p=1d=0d=0q=0q=0

 

y^t=μϵt+ϕ1yt−1y^t=μϵt+ϕ1yt−1

 

自回归的影响

我们可以使用该arima.sim函数模拟自回归过程。通过该功能，可以通过提供要使用的MA和AR项的系数来指定模型。在下文中，我们将绘制自相关图，因为它最适合于发现自回归的影响。

<span styl 0" /span>)
par(mfrow = c(<span style="color:#880000">2</span>, <span style="color:#880000">2</span>))
<span style="color:#888888">#  880000">"ARIMA(1,0,0)"</span>)
<span style="color:#888888"># plot partial acf</span>
acf(x, type = <span style="color:#880000">"partial"</span>, main = <span style="color:#880000">"Partial aut "color:#880000">0.65</span>, <span style="color:#880000">0.3</span>)), 
        n = <span style="color:#880000">1000</span>)
plot(x, main = <span style="color:#880000">"ARIMA(2,0,0)"</span>)
acf(x, type = <span style="color:#880000">"partial"< span></span>

​

第一个例子表明，对于ARIMA（1,0,0）过程，订单1的pACF非常高，而对于ARIMA（2,0,0）过程，订单1和订单2自相关都很重要。因此，可以根据pACF显着的最大滞后来选择AR项的顺序。

dd

 ARIMA（0,1,0）模型简化为随机游走模型

 

y^t=μ+ϵt+yt−1.y^t=μ+ϵt+yt−1.

 

 差异的影响

以下示例演示了差异对AirPassengers数据集的影响：

<span style="color:#0 "><span style /span>), main = <span style="color:#880000">"After differencing"</span>)</code></span></span>

​

虽然第一个图表显示数据显然是非静止的，但第二个图表明差异时间序列是相当静止的。

qq

 

 

其中当前估计值取决于先前测量值的残差。

移动平均线的影响

可以通过绘制自回归函数来研究移动平均线的影响：

<span style="color:#00 # Example for ARIMA(0,0,1)</span>
x <- arima.sim(list(ma = <span style="color:#880000">0.75</span>),
                n = <span style="color:#880000">1000</span>)
plot(x, main = <span style="color:#880000">"ARIMA(0,0,1)"</span>)
acf(x, main = <span style="color:#88000 n"</span>)</code></span></span>

​

请注意，对于自回归图，我们需要注意第一个x轴位置表示滞后为0（即标识向量）。在第一个图中，只有第一个滞后的自相关是显着的，而第二个图表明前两个滞后的自相关是显着的。为了找到MA术语的数量，适用与AR术语类似的规则：MA术语的顺序对应于自相关显着的最大滞后。

在AR和MA术语之间进行选择

为了确定哪个更合适，AR或MA术语，我们需要考虑ACF（自相关函数）和PACF（部分ACF）。使用这些图我们可以区分两个签名：

  • pp
  • rr

AR和MA术语的影响

AR和MA术语的组合导致以下时间序列数据：

<span style="color:# "># ARIMA(1,0,1)</span>
x <- arima.sim(list(order = c(<span style="color:#880000">1</span>,<span style="color:#880000">0</span>,<span style="color:#880000">1</span>), ar = <span style="color:#880000">0.8</span>, ma = <span style="color:#880000">0.8</span>), n = <span style="color:#880000">1000</span>)
plot(x, main = <span style="color:#880000">"ARIMA(1,0,1)"</span>)
acf(x, main = <span   ARIMA(2,0,2)</span>
x <- arima.sim( )
plot(x, main = <span style="color:#880000">"ARIMA(2,0,2)"</span>)
acf(x, main = <span style="color:#880000">"ACF"</span>)
pacf(x, main = <span style="color:#880000">"pACF"</span>)</code></span></span>

​

 SARIMA模型

  •  P：季节性自回归（SAR）项的数量
  • D：季节差异程度
  • 问：季节性移动平均线（SMA）的数量

 

ARIMAX模型

 

R中的预测

auto.arimaforecastppddqqPPDDQQstepwiseapproximationFALSE

SARIMA模型用于固定过程

我们将使用包中的nino数据展示ARMA的使用，该数据tseries给出了Nino Region 3.4指数的海面温度。让我们验证数据是否静止：

<span style="color:#000000"><span style="color:#000000"><code>plot(nino3.4)</code></span></span>

​

d=0d=0

为了验证是否存在任何季节性趋势，我们将分解数据：

<span style="color:#000000"><span style="color:#000000"><code>nino.components <- decompose(nino3.4) an>

​

没有整体趋势，这是固定过程的典型趋势。但是，数据存在强烈的季节性因素。因此，我们肯定希望包含对季节性影响进行建模的参数。

季节性模型

(P,D,Q)S(P,D,Q)SD=0D=0ninoS=12S=12

<span style="color:#000000"><span style="color:#000000"><code>nino.season <- nino.components$seasonal ode></span></span>

​

P=2P=2Q=0Q=0

非季节性模型

ppqq

<span style="color:#000000">< tyle="color:#000000"><code>par(mfrow = c(<span style="color:#880000">1</span>,<span style="color:#880000">2</span>))
acfp <- acf(n #888888"># transform lag from years to months</span>
acfp$lag <- acfp$lag * <span style="color:#880000">12</span>
plot(acfp, main = <span style="color:#880000"> months</span>
acfpl$lag <- acfpl$lag * <span style="color:#880000">12</span>
plot(acfpl, main = <span style="color:#880000">"pACF"</span>)</code></span></span>

​

 我们可以使用包中的Arima函数来拟合模型forecast。

<span style=" r:#888888"># non-seasonal model: (p,d,q)</span>
order.non.seasonal <- c(<span style="color:#880000">2</span>,<span style="color:#880000">0</span> r:#880000">2</span>,<span style="color:#880000">0</span>,<span style="color:#880000">0</span>)
A <- Arima(nino3.4, order = order.non.seasonal,
            seasonal = order.seasonal)</code> </span>

我们现在可以使用该模型来预测未来Nino 3.4地区的气温如何变化。有两种方法可以从预测模型中获得预测。第一种方法依赖于predict函数，而第二种方法使用包中的forecast函数forecast。使用该predict功能，我们可以通过以下方式预测和可视化结果：

<span sty #000000"><span style="color:#000000"><code><span style="color:#888888"># to construct a custom plot, we can us n>
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:forecast':
## 
##     autolayer
<span style="color:#00000 yle="color:#880000">0</span>), cbind(fortify(fore df$y + plot.df$sd * <span style="color:#880000">1.96</span>
plot.df$lower <- plot.df$y - plot.df$sd * <span style="color:#880000">1.96</span>
ggplot(plot.df, aes(x = x ,y = y)) +  
        ylab(<span style="color:#880000">"Temperature"</span>) + xlab(<span style="color:#880000">"Year"</span>)</code></span></span>

​

如果我们不需要自定义绘图，我们可以使用以下forecast函数更轻松地获取预测和相应的可视化：

# use the forecast function to use the built-in plotting function:
forecast <- forecast(A, h = 60) # predict 5 years into the future
plot(forecast)

用于非平稳数据的ARIMA模型

为了演示ARIMA模型对非平稳数据的使用，我们将使用包中的gtemp数据集astsa。该数据集提供全球平均陆地 - 海洋温度偏差的年度测量值。

<span style="color:#000000"><span style=" 0"><code><strong>library</strong>(astsa)
data(gtemp)

​

d=1d=1

<span style="color:#000000"><span style="

​

现在，数据似乎是静止的。

 

<span style="color:#000 0"><span style="color:#000000"><code>par(mfrow = c(<span style="color:#880000">1</span>,<span s code></span></span>

​

p=0p=0q=1q=1

<span style="color:#000000"><span style="

我们现在可以预测未来几年平均陆地 - 海洋温度偏差将如何变化：

<span style="color:#000000"> tyle="color:#000000"><cod 8"># predict 30 years into the future</span> pan></span>

​

该模型表明，未来几年平均陆地 - 海洋温度偏差将进一步增加。

关于空气质量数据集的ARIMAX

为了展示ARIMAX模型的使用，我们将使用臭氧数据集 。

让我们加载臭氧数据集并将其划分为测试和训练集。请注意，我们已确保训练和测试数据包含连续的时间测量。

<span style="color:#000000"><span style="color:#000000"><code>data(airquality)
ozone <- subset(na.omit(airquality))
set.seed(<span style="color:#880000">123</span>)
N.train <- ceiling(<span style="color:#880000">0.7</span> * nrow(ozone))
N.test <- nrow(ozone) - N.train
<span style="color:#888888"># ensure to take only subsequent measurements for time-series analysis:</span>
trainset <- seq_len(nrow(ozone))[<span style="color:#880000">1</span>:N.train]
testset <- setdiff(seq_len(nrow(ozone)), trainset)</code></span></span>

由于数据集未指示相对时间点，我们将手动创建此类注释：

为此，我们将在臭氧数据集中创建一个新列，该列反映了相对时间点：

<span :#000000"><span st style="color:#880000">"%j"</span>))
max.date <- as.numeric(format(max(dates), <span style="color:#880000">"%j"</span>))
ozone.ts <- ts(ozone$Ozone, start = min.date, end = max.date, frequency = <span style="color:#880000">1</span>)
ozone.ts <- window(ozone.ts, <span style="color:#880000" 21</span>, <span style="color:#880000">231</span>) <span style="color:#888888"># deal with repetition due to missing time values</span>
ozone$t <- seq(start(ozone.ts t</span></code></span></span>

现在我们有了时间维度，我们可以绘制臭氧水平的纵向行为：

<span style="color:#000000"><span style="color:#000000"><code><strong>library</strong>(ggplot2)
ggplot(ozone, aes(x = t, y = Ozone)) + geom_line() +
        geom_point()</code></span></span>

​

时间序列数据似乎是固定的。让我们考虑ACF和pACF图，看看我们应该考虑哪些AR和MA术语

<span style="color:#000000">< :#880000">"partial"</span>)</code></span></span>

​

自相关图非常不清楚，这表明数据中实际上没有时间趋势。因此，我们会选择ARIMA（0,0,0）模型。由于具有参数（0,0,0）的ARIMAX模型没有传统线性回归模型的优势，我们可以得出结论，臭氧数据的时间趋势不足以改善臭氧水平的预测。让我们验证一下：

<span style="color eviously developed weighted negative binomial model</span>
<strong>library</strong>(MASS)
get.weights <- <strong>function</strong>(ozone) {
    z.scores <- (ozone$Ozone - mean(ozone$Ozone)) / sd(ozone$Ozone)
    weights <- exp(z.scores)
    weights <- l$pred, ozone[testset, <span style="color:#880000">"Ozone"</span>])^<span style="color:#880000">2</span>
print(Rsquared.linear)</code></span></span>
## [1] 0.7676977
<span style="color:#000000"><span style="color:#000000"><code>print(Rsquared.temporal)</code></span></span>
## [1] 0.7569718

我们可以看到具有负二项式可能性的线性模型优于ARIMAX模型。

关于空气质量数据集的ARIMAX

要在更合适的数据集上演示ARIMAX模型，让我们Icecream从Ecdat包中加载数据集：

<span style="color:#000000"><

该Icecream数据集包含以下变量：

  • 缺点：人均品脱的冰淇淋消费量。
  • 收入::美元平均每周家庭收入。
  • 价格：每品脱冰淇淋的价格。
  • temp：华氏温度的平均温度。

测量结果是从1951-03-18到1953-07-11的四周观测。

我们将模拟缺点，冰淇淋消费作为时间序列，并使用收入，价格和平均值作为外生变量。在开始建模之前，我们将从数据框中创建一个时间序列对象。

<span style="col  style="color:#880000">"1951-03-18"</span>), as.Date(<span style="color:#880000">"1953-07-11"</span>)))
months <- c(seq(<span style="color:#880000">3</span>,<span style="color:#880000">12</span>), se or:#880000">1</span>, <span style="color:#880000">52</span>, <span style="color:#880000">4</span>))
ice.ts <- ts(Icecream$cons, start = c(<span style="color:#880000">1951</span>, <span style="color:#880000">3</span>), end = c(<span style="color:#880000">1953</span>, <span style="color:#880000">6</span>), frequency = <span style="color:#880000">52</span>/<span style="color:#880000">4</span>)</code></span></span>

我们现在调查数据：

<span style="color:#000000"><span style="color:#000000"><code>plot(decompose(ice.ts))</code></span></span>

​

因此，数据有两种趋势：

  1. 总体而言，1951年至1953年间，冰淇淋的消费量大幅增加。
  2. 冰淇淋销售在夏季达到顶峰。

ppqq

<span style="color:#000000"><span style 
acf(ice.ts, type = <span style="color:#880000">"partial"</span>)</code></span></span>

​

由于季节性趋势，我们可能适合ARIMA（1,0,0）（1,0,0）模型。但是，由于我们知道温度和外生变量的收入，因此它们可以解释数据的趋势：

​

<span style="color:#000000"><span style="color:#000000"><code>plot(Icecream$temp) <span style="color:#888888"># explains the seasonal trend</span></code><

​

由于income解释了整体趋势，我们不需要漂移术语。此外，由于temp解释了季节性趋势，我们不需要季节性模型。因此，我们应该使用ARIMAX（1,0,0）模型进行预测。为了研究这些假设是否成立，我们将使用以下代码将ARIMAX（1,0,0）模型与ARIMA（1,0,0）（1,0,0）模型进行比较

 yle="color:#880000">"income"</span>, <span style="color:#880000">"temp"</span>)],
        order = c(<span style="color:#880000">1</span>,<span style="color:#880000">0</span>,<span style="color:#880000">0</span>))
preds <- forecast(A, xreg = Icecream[test, c(<span style="color:#880000">"income"</span>, <span style="color:#880000">"temp"</span>)]) 
plot(preds)
lines(window(ice.ts, c(<span style="color:#880000">1951</span>, <span style="color:#880000">22</span>), c(<span style="color:#880000">1951</span>, <span  ="color:#880000">24</span>)
lines(x = as.numeric(rownames(as.data.frame(preds))), y = as.data.frame(preds)[,<span style="color:#880000">2</span>], lty = <span style="color:#880000">2</span>)</code></span></span>

​

ARIMAX（1,0,0）模型的预测显示为蓝色，而ARIMA（1,0,0）（1,0,0）模型的预测显示为虚线。实际观察值显示为黑线。结果表明，ARIMAX（1,0,0）明显比ARIMA（1,0,0）（1,0,0）模型更准确。

但请注意，ARIMAX模型在某种程度上不像纯ARIMA模型那样有用于预测。这是因为，ARIMAX模型需要对应该预测的任何新数据点进行外部测量。例如，对于冰淇淋数据集，我们没有超出1953-07-11的外生数据。因此，我们无法使用ARIMAX模型预测超出此时间点，而ARIMA模型可以实现：

<span style="color:#000000"><span style="color:#000000"><code>preds <- forecast(A.season, h = <span style="color:#880000">60</span>)
plot(preds)</code></span></span>

​

如果您有任何疑问，请在下面发表评论。 

posted @ 2019-05-06 17:55  拓端tecdat  阅读(1420)  评论(0)  编辑  收藏  举报