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) # y_i - B y_i
B3y <- diff(y, 3) # y_i - B^3 y_i
message(paste0("y is: ", paste(y, collapse =
## y is: 1,3,5,10,20
## By is: 2,2,5,10
## B^3y is: 9,17
自相关函数
要计算自相关,我们可以使用以下R函数:
get_autocor <- function(x, lag) {
x.left <- x[1:(length(x) - lag)]
x.right <- x[(1+lag):(length(x))]
autocor <- cor(x.left, x.right)
return(autocor)
}
get_autocor(y, 1)
## [1] 0.9944627
# correlation of measurements 2 time points apart (lag 2)
get_autocor(y, 2)
## [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
函数:
par(mfrow = c(1,2))
acf(y) # conventional ACF
pacf(y) # pACF
在ACF可视化中,ACF或pACF被绘制为滞后的函数。指示的水平蓝色虚线表示自相关显着的水平。
分解时间序列数据
StSt
TtTt
ϵtϵt
执行分解的方式取决于时间序列数据是加法还是乘法。
加法和乘法时间序列数据
加法模型假设数据可以分解为
yt = St + Tt + ϵt.yt = St + Tt + ϵt.
另一方面,乘法模型假设数据可以被分解为
添加剂:每个时期的季节效应放大器相似。
乘法:季节性趋势随时间序列的变化而变化。
AirPassengers
数据集提供了乘法时间序列的示例。
data(AirPassengers)
plot(AirPassengers)
log(StTtϵt)=log(St)+log(Tt)+log(ϵt)log(StTtϵt)=log(St)+log(Tt)+log(ϵt)
AirPassengers
数据集:
plot(log(AirPassengers))
正如我们所看到的,采用对数已经使季节性成分的幅度沿时间均衡。请注意,总体增长趋势没有改变。
在R中分解时间序列数据
要分解R中的时间序列数据,我们可以使用该decompose
函数。请注意,我们应该通过type
参数提供时间序列是加法的还是乘法的。
示例1:AirPassengers数据集
对于AirPassengers
数据集,我们指定数据是乘法的并获得以下分解:
plot(decompose(AirPassengers, type = "multiplicative"))
分解表明,多年来航空公司乘客总数在增加。此外,我们已经观察到的季节性影响已被清楚地捕捉到。
示例2:EuStockMarkets数据集
让我们考虑可以为EuStockMarkets
数据集找到的分解:
daxData <- EuStockMarkets[, 1] # DAX data# data do not seem to be multiplicative, use additive decomposition
decomposed <- decompose(daxData, type = "additive")
plot(decomposed)
该图显示了1992年至1998年的DAX数据中的以下内容:
整体价值稳步上升。
季节性趋势强烈:每年年初,股价相对较低,并在夏季结束时达到相对最大值。
除1997年和1998年之间的最终测量外,随机噪声的贡献可以忽略不计。
平稳与非平稳过程
生成时间序列数据的过程可以是平稳的也可以是非平稳的。例如,数据EuStockMarkets
和AirPassengers
数据都是非平稳的,因为数据有增加的趋势。为了更好地区分平稳和非平稳过程,请考虑以下示例:
par(mfrow = c(1,2))
data(nino)
x <- nino3.4
plot(x, main = "Stationary process")
plot(AirPassengers, main = "Non-stationary process")
左图显示了一个平稳过程,其中数据在所有测量中表现相似。右图显示了一个非平稳过程,其中平均值随着时间的推移而增加。
介绍了与时间序列数据分析相关的最重要概念后,我们现在可以开始研究预测模型。
ARMA模型
ARMA代表自回归移动平均线。ARMA模型仅适用于平稳过程,并具有两个参数:
p:自回归(AR)模型的顺序
q:移动平均(MA)模型的顺序
ARMA模型可以指定为
使用以下变量:
cc
ϵtϵtttϵt∼N(0,σ2)ϵt∼N(0,σ2)
ϕ∈Rpϕ∈Rp
ytyttt
θ∈Rqθ∈Rq
ϵtϵttt
使用backshift运算符制定ARMA模型
使用backshift运算符,我们可以通过以下方式制定ARMA模型:
ϕ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模型
总之,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模型的三个参数的解释。
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项的系数来指定模型。在下文中,我们将绘制自相关图,因为它最适合于发现自回归的影响。
第一个例子表明,对于ARIMA(1,0,0)过程,订单1的pACF非常高,而对于ARIMA(2,0,0)过程,订单1和订单2自相关都很重要。因此,可以根据pACF显着的最大滞后来选择AR项的顺序。
ARIMA(0,1,0)模型简化为随机游走模型
y^t=μ+ϵt+yt−1.y^t=μ+ϵt+yt−1.
差分的影响
以下示例演示了差异对AirPassengers
数据集的影响:
虽然第一个图表显示数据显然是非静止的,但第二个图表明差异时间序列是相当静止的。
其中当前估计值取决于先前测量值的残差。
移动平均线的影响
可以通过绘制自回归函数来研究移动平均线的影响:
请注意,对于自回归图,我们需要注意第一个x轴位置表示滞后为0(即标识向量)。在第一个图中,只有第一个滞后的自相关是显着的,而第二个图表明前两个滞后的自相关是显着的。为了找到MA术语的数量,适用与AR术语类似的规则:MA术语的顺序对应于自相关显着的最大滞后。
在AR和MA之间进行选择
为了确定哪个更合适,AR或MA术语,我们需要考虑ACF(自相关函数)和PACF(部分ACF)。使用这些图我们可以区分两个签名:
pp
rr
AR和MA术语的影响
AR和MA术语的组合导致以下时间序列数据:
SARIMA模型
P:季节性自回归(SAR)项的数量
D:季节差异程度
问:季节性移动平均线(SMA)的数量
ARIMAX模型
R中的预测
auto.arima
forecast
ppddqqPPDDQQstepwise
approximation
FALSE
SARIMA模型用于平稳过程
我们将使用数据展示ARMA的使用,该数据tseries
给出了Nino Region 3.4指数的海面温度。让我们验证数据是否静止:
d=0d=0
为了验证是否存在任何季节性趋势,我们将分解数据:
没有整体趋势,这是平稳过程的典型趋势。但是,数据存在强烈的季节性因素。因此,我们肯定希望包含对季节性影响进行建模的参数。
季节性模型
(P,D,Q)S(P,D,Q)SD=0D=0
nino
S=12S=12
P=2P=2Q=0Q=0
非季节性模型
我们可以使用包中的Arima
函数来拟合模型forecast
。
我们现在可以使用该模型来预测未来Nino 3.4地区的气温如何变化。有两种方法可以从预测模型中获得预测。第一种方法依赖于predict
函数,而第二种方法使用包中的forecast
函数forecast
。使用该predict
功能,我们可以通过以下方式预测和可视化结果:
#
# Attaching package: 'ggplot2'## The following object is masked from 'package:forecast':
#
# autolayer
如果我们不需要自定义绘图,我们可以使用以下forecast
函数更轻松地获取预测和相应的可视化:
forecast <- forecast(A, h = 60) # 预测未来5年
plot(forecast)
用于非平稳数据的ARIMA模型
为了演示ARIMA模型对非平稳数据的使用,我们将使用数据集astsa
。该数据集提供全球平均陆地 - 海洋温度偏差的年度测量值。
d=1d=1
p=0p=0q=1q=1
我们现在可以预测未来几年平均陆地 - 海洋温度偏差将如何变化:
该模型表明,未来几年平均陆地 - 海洋温度偏差将进一步增加。
关于空气质量数据集的ARIMAX
为了展示ARIMAX模型的使用,我们将使用臭氧数据集 。
让我们加载臭氧数据集并将其划分为测试和训练集。请注意,我们已确保训练和测试数据包含连续的时间测量。
由于数据集未指示相对时间点,我们将手动创建此类注释:
为此,我们将在臭氧数据集中创建一个新列,该列反映了相对时间点:
现在我们有了时间维度,我们可以绘制臭氧水平的纵向行为:
时间序列数据似乎是平稳的。让我们考虑ACF和pACF图,看看我们应该考虑哪些AR和MA术语
自相关图非常不清楚,这表明数据中实际上没有时间趋势。因此,我们会选择ARIMA(0,0,0)模型。由于具有参数(0,0,0)的ARIMAX模型没有传统线性回归模型的优势,我们可以得出结论,臭氧数据的时间趋势不足以改善臭氧水平的预测。让我们验证一下:
print(Rsquared.linear)
## [1] 0.7676977print(Rsquared.temporal)
## [1] 0.7569718
我们可以看到具有负二项式可能性的线性模型优于ARIMAX模型。
关于空气质量数据集的ARIMAX
要在更合适的数据集上演示ARIMAX模型,让我们加载数据集:
该Icecream
数据集包含以下变量:
缺点:人均品脱的冰淇淋消费量。
收入::美元平均每周家庭收入。
价格:每品脱冰淇淋的价格。
temp:华氏温度的平均温度。
测量结果是从1951-03-18到1953-07-11的四周观测。
我们将模拟缺点,冰淇淋消费作为时间序列,并使用收入,价格和平均值作为外生变量。在开始建模之前,我们将从数据框中创建一个时间序列对象。
我们现在调查数据:
因此,数据有两种趋势:
总体而言,1951年至1953年间,冰淇淋的消费量大幅增加。
冰淇淋销售在夏季达到顶峰。
由于季节性趋势,我们可能适合ARIMA(1,0,0)(1,0,0)模型。但是,由于我们知道温度和外生变量的收入,因此它们可以解释数据的趋势:
由于income
解释了整体趋势,我们不需要漂移术语。此外,由于temp
解释了季节性趋势,我们不需要季节性模型。因此,我们应该使用ARIMAX(1,0,0)模型进行预测。为了研究这些假设是否成立,我们将使用以下代码将ARIMAX(1,0,0)模型与ARIMA(1,0,0)(1,0,0)模型进行比较
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模型可以实现:
如果您有任何疑问,请在下面发表评论。
点击标题查阅往期内容
更多内容,请点击左下角“阅读原文”查看
案例精选、技术干货 第一时间与您分享
长按二维码加关注
更多内容,请点击左下角“阅读原文”查看