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>

​

如果您有任何疑問，請在下面發表評論。