時間序列分析在R中實踐
阿新 • • 發佈：2019-01-31
先安裝時間序列相關包，直接從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)