arima模型r语言 xreg r语言ar模型预测

R时间序列分析|S&P500股指的ARIMA模型预测与残差ARCH效应分析

• 前言

• 一、数据及分析目的
• 二、数据探索
• 三、ARIMA模型构建
• 四、残差分析
• 五、模型预测

前言

由于R语言对新手并不够友好，网上的资料相对也偏少，导致博主在使用R进行时间序列分析的过程中非常痛苦，参考和大量的资料和教科书做法。因此在project完成后，将代码、分析过程及必要的解释分享如下，希望可以帮到有需要的人。

创作不易，转载请注明。如有错漏之处，还请各位大佬指正。

一、数据及分析目的

Project要求对SP500的数据进行探索和分析，建立三个拟合模型，使用预测精度对模型进行筛选，并对模型残差包含的信息进行分析，最后做出预测。

主要用到的R包如下:

# 如果没有安装，直接使用install.packages()即可
library(depmixS4)
library(tseries)
library(forecast)
library(readxl)
library(stats)
library(TSA)

二、数据探索

package里自带有1950.2至2012.1的月末数据，这里选择收盘价作为train data。并且，在雅虎下载了2012.2至2013.1的数据作为test data。

#导入数据
data("sp500")
#导入12.02-13.01的数据
sp500_12_13 <- read_excel("E:/^GSPC.xlsx")#雅虎金融可直接下载
#转化为时间序列
sp500_12_13 <- ts(sp500_12_13['Close'],frequency=12, start=c(2012,2))
close_price <- ts(sp500['Close'],frequency=12, start=c(1950,2))
print(close_price)

原序列及对数序列可视化

win.graph(width=3.25,height=2.5,pointsize=8)
plot(close_price,ylab='SP500 Close Price')
plot(log(close_price),ylab='Log(SP500 Close Price)')

首先考虑将序列化为平稳序列。观察原序列图像及对数处理图像可以发现，对数处理可以在一定程度上缓解方差波动的问题。但序列仍存在明显的增长性，可以通过对序列进行分解直观观察。

#打印出四种趋势图
dc<-decompose(log(close_price))
plot(dc)

为了消除强烈的增长性，使用差分处理，并使用Augmented Dickey-Fuller Test进行序列平稳性检验。

adf.test(diff(log(close_price)))
#pp.test(diff(log(close_price)))#也可以用此方法进行检验

Augmented Dickey-Fuller Test
data: diff(log(close_price))
Dickey-Fuller = -8.6284, Lag order = 9, p-value = 0.01
alternative hypothesis: stationary

可以看到拒绝原假设，序列平稳，可以进行ARIMA建模。

三、ARIMA模型构建

首先对取对数后一阶差分的序列进行分解。

#打印出关于季节性趋势的图表
dc<-decompose(log(close_price))
season<-dc$figure
plot(season,type = "b",xaxt="n",xlab = "Month",ylab = "Season Effect")
#或者这样
boxplot(log(close_price)~cycle(log(close_price)))#显示季节性特征

没有观察到十分强烈的季节性特征（波动范围较小），序列增长性也没有十分强烈。值得注意的是，分解后的trend序列似乎仍存在一些波动趋势，但并不能被直观观察到（后续进行分析）。接着判断序列的自相关性。

#查看ACF和PACF图
win.graph(width=3.25,height=2.5,pointsize=8)
tsdisplay(diff(log(close_price)))
dc<-decompose(diff(log(close_price)))
plot(dc)

ACF图和PACF图都没有明显的截尾和拖尾特征，对AR()和MA()的之后阶数不好判断，但1~5期滞后都表现出较强的特征，可以考虑逐个建模比较AIC值（越小越好）。

但R提供了简洁的方法，使用函数自动寻找最优模型。这里分别使用AIC和BIC准则进行寻找。

#默认使用aic判断
auto.arima(log(close_price))
#使用bic判断
auto.arima(log(close_price),ic='bic')

#ARIMA(1,1,1)(0,0,1)[12] with drift
#ARIMA(0,1,0) with drift

给出了两个模型，一个是带漂移项的季节性ARIMA(1,1,1)(0,0,1)₁₂模型，另一个是带漂移项的随机游走ARIMA(0,1,0)模型。

漂移项$\mu$是一非零常数，是因为一阶差分为：
$\Delta\text{X}_t=\text{X} -\text{X}_{t-1}= \mu + \epsilon_t$这表明时间序列$\text{X}_t$向上或向下漂移，取决于$\mu$的符号是正还是负。通常ARIMA模型可通过再次差分消除漂移项。但本文依然包括了漂移项进行建模。

分别对两个模型进行建模。

注：arima()函数带漂移项建模之后，在预测阶段频频出现问题，博主能力有限无法解决，因此使用Arima()函数进行建模，两个函数几乎一样，只是该函数对带漂移项建模更加友好。

#ARIMA(1,1,1)(0,0,1)[12] with drift
model2=Arima(log(close_price),order=c(1,1,1),seasonal=list(order=c(0,0,1),period=12),include.drift = TRUE)
summary(model2)#aic = -2588.36

#ARIMA(0,1,0) with drift
model3=Arima(log(close_price),order=c(0,1,0),include.drift = TRUE)
summary(model3)#bic = -2576.33

对模型的summary整理如下。

各种误差的意思如下：
预测误差=实际值 - 相对值
相对误差=（实际值 - 相对值）/实际值
预测精度=1 - 相对误差
ME: Mean Error=$\dfrac{1}{n}\Sigma(预测误差)$
RMSE: Root Mean Squared Error
MAE: Mean Absolute Error=$\dfrac{1}{n}\Sigma|预测误差|$
MPE: Mean Percentage Error=$\dfrac{1}{n}\Sigma|相对误差|$
MAPE: Mean Absolute Percentage Error=$\dfrac{1}{n}\Sigma\dfrac{|预测误差|}{实际值}$
MASE: Mean Absolute Scaled Error
ACF1: Autocorrelation of errors at lag 1.

接着是写出模型的表达式，季节性ARIMA模型的表达式比较复杂，其形式如下：

因此，ARIMA(1,1,1)(0,0,1)₁₂ with drift的表达式应为：

$(1-\phi_1B)(1-B)\ln\text{y}_t=(1+\theta_1B)(1+\Theta_1B^{12})\epsilon_t+\mu$其中$B$为滞后算子，$B^n$即滞后n期，$\mu$为漂移项。

根据建模的结果，model2的模型表达式为：

model3的模型表达式更加简单：

比较两个模型的预测精度，即误差越小，精度越高。下表可见model2的精度稍微更高一些。（model1为线性回归模型，本文不展示）

四、残差分析

#残差分析
#残差相关性质图
win.graph(width=6.5,height=6); tsdiag(model2)
#残差分布
win.graph(width=4,height=3,pointsize=8)
hist(residuals(model2),xlab='Residuals')
#残差Q-Q图
win.graph(width=6.5,height=6)
qqnorm(residuals(model2)); qqline(residuals(model2))
#残差正态性检验
shapiro.test(residuals(model2))

残差没有太多的异常值和自相关性。 但是，在某些年份（例如2000年和2008年）中似乎存在一些方差群集，这似乎是方差非齐，稍后进行检查。

Q-Q图显示残差序列最下部分似乎偏离了正态曲线，应该是个偏态分布。Shapiro-Wilk normality test 验证了我们的想法。该检验的原假设是序列为正态分布。

模型的残差似乎还蕴含了更多的信息，考虑进行ARCH效应检验。

#LM检验图
win.graph(width=6.5,height=6);McLeod.Li.test(y=residuals(model2))
#或者用这个进行LM检验（是否存在Arch效应），即检验残差方差齐性（H0）
library(FinTS)
ArchTest(residuals(model2))

上图显示，不管序列滞后多少期，LM检验都是显著的，检验表也是相同的结论，即模型的残差序列存在arch效应。

考虑建立Garch模型，使用eacf图对p,q定阶。

eacf(abs(residuals(model2)))
eacf(residuals(model2)^2)

绝对值数据和平方数据的示例EACF如下所示。

这部分不是特别懂，x表示显著，大概是看o围成的三角形。教科书中对十分相近的eacf图的描述原话是这样的：

Both of them convincingly suggest an ARMA(1,1) model, and therefore a GARCH(1,1) model for the original data.
他们两个都令人信服地建议使用ARMA（1,1）模型，因此建议使用原始数据的GARCH（1,1）模型。

建立Garch(1,1)模型。

r_garch=garch(x=residuals(model2), order=c(1,1) ,reltol=0.000001)
summary(r_garch)

所有系数均显著。Jarque Bera检验的结果表明Garch模型的残差是非正态的，而Box-Ljung检验的结果表明残差的平方是白噪声序列。

plot(residuals(r_garch),type='h',ylab='standard residuals')
acf(residuals(r_garch)^2,na.action=na.omit)
gBox(r_garch,method='squared')

模型的残差看起来很正常。

残差平方序列看起来没有自相关性。

Generalized Portmanteau Tests for GARCH Models也都未显著。

因此，对model2的残差建立Garch(1,1)模型进行拟合是合适的。

基于Garch(1,1)画出条件方差图。

plot((fitted(r_garch)[,1])^2,type='l',ylab='conditional variance',xlab='t')

上图显而易见，model2的残差的方差确实会随时间变化（方差非齐），并且条件方差在1975、1988和2009年出现三个峰值，正好符合历史上的股市崩盘。

model2的残差包含了标准普尔500指数方差集群信息，这对我们观察经济波动很有帮助。

五、模型预测

#ARIMA(1,1,1)(0,0,1)[12] with drift
model=Arima(log(close_price),order=c(1,1,1),seasonal=list(order=c(0,0,1),period=12),include.drift = TRUE)
#预测未来12期，99.5%置信区间
forecast<-forecast(model,h=12,level=c(99.5))
#可视化
win.graph(width=6.5,height=6)
#只显示2000至2013年
plot(exp(forecast$fitted),ylim=c(500,2300),xlim=c(2008,2013),col="green")
lines(close_price,col="red")
lines(exp(forecast$mean),col="green")#预测值
lines(sp500_12_13,col="red")#预测值对应的真实值
lines(exp(forecast$upper),lty=2)#预测的上界
lines(exp(forecast$lower),lty=2)#预测的下界
#预测精度
accuracy(exp(forecast$mean),sp500_12_13)

train set和test set的误差和预测精度如下。

模型的预测能力看起来还行，误差也在可接受的范围内，平均预测精度达到了96.87%，成功预测了SP500稍微向上的走势。

注：上图后半部分线断开的原因不是数值缺失，而是前后序列不同引起的。