时间序列分析在R中实践

先安装时间序列相关包，直接从CRAN上安装，会自动安装依赖的包，再导入包

#####################包准备#####################
install.packages("lmtest")
install.packages("tseries")
install.packages("fGarch")
install.packages("FinTS")
install.packages("forecast")
install.packages("MSBVAR")

library(xts)
library(lmtest)
library(tseries)
library(fGarch)
library(FinTS)
library(forecast)
library(MSBVAR)

本文选择2005年1月4日至2007年10月17日共673个数据，上证指数和深圳指数，收盘价和收益率数据部分如下：

date	sz_prc	sh_prc	sz_rat	sh_rat
2005/1/5	315.25	1251.94	0.0064	0.0032
2005/1/6	311.98	1239.43	-0.0045	-0.0044
2005/1/7	312.61	1244.75	0.0009	0.0019
2005/1/10	315.85	1252.4	0.0045	0.0027
2005/1/11	316.42	1257.46	0.0008	0.0018
…	…	…	…	…

sz_prc，sh_prc分别表示深证指数和上证指数的收盘价，sz_rat，sh_rat分别表示深证指数和上证指数的收盘率：

 

  

#####################数据准备#####################
d <- read.csv("E:/Work/R模型/数据表.csv",header=T)
dt <- as.Date(d[,"date"])
rh <- d[,"sh_rat"]
rz <- d[,"sz_rat"]
rt.rh <- xts(rh,dt)
rt.rz <- xts(rz,dt)

开始分别对RH和RZ进行自相关分析

#判断自相关阶数
pacf(rt.rh)
ar_rh <- auto.arima(rt.rh,max.p = 20)
l_rh <- ar_rh$arma[1]

根据滞后阶数定义新的数据序列

rt.rhl <- lag(rt.rh,k=-l_rh,na.pad=FALSE)
rt.rhsl <- rt.rh[1:length(rt.rhl)]

对序列直接做线性回归

#自相关滞后回归
fit_rh <- lm(rt.rhsl ~ rt.rhl)
summary(fit_rh)

r_rh <- resid(fit_rh)
#是否存在ARCH效应
ArchTest(as.ts(r_rh))

#GARCH(1,1)模型
l_rh   #6
garchFit(~arma(6,0)+garch(1,1),data=rt.rh,algorithm="nlminb+nm",trace=F,include.mean=F)

#判断自相关阶数
pacf(rt.rz)
ar_rz <- auto.arima(rt.rz,max.p = 20)
l_rz <- ar_rz$arma[1]

根据滞后阶数定义新的数据序列

rt.rzl <- lag(rt.rz,k=-l_rz,na.pad=FALSE)
rt.rzsl <- rt.rz[1:length(rt.rzl)]

对序列直接做线性回归

#自相关滞后回归
fit_rz <- lm(rt.rzsl ~ rt.rzl)
summary(fit_rz)

r_rz <- resid(fit_rz)
#是否存在ARCH效应
ArchTest(as.ts(r_rz))

#GARCH(1,1)模型
l_rz   #7
garchFit(~arma(7,0)+garch(1,1),data=rt.rz,algorithm="nlminb+nm",trace=F,include.mean=F)

分别对序列做单位根检验

#####################单位根检验#####################
adf.test(rt.rh)
adf.test(rt.rz)

#####################协整检验#####################
fit_xz <- lm(rt.rh ~ rt.rz)
r_xz <- residuals(fit_xz)
adf.test(as.ts(r_xz))

#####################Granger因果检验#####################
rt.gg <- cbind(rt.rh,rt.rz)
colnames(rt.gg) <- c("rh","rz")
granger.test(rt.gg,p=1)

从上图可以看rh是导致rz变化的granger原因。

#####################建立误差修正模型#####################
dy_rh <- diff(rt.rh,na.pad=FALSE)
dy_rz <- diff(rt.rz,na.pad=FALSE)
r_lag<- lag(r_xz,k=+1,na.pad=FALSE)
fit_ECM <- lm(dy_rh ~ r_lag + dy_rz)
summary(fit_ECM)