R语言DCC-GARCH模型

本人的数据形式如下。从Excel中复制粘贴到了txt中。

导入数据，并绘制三个样本序列的价格序列图。

p=read.table('C:/Users/dqxq767/Documents/R/p.txt',sep = '\t',header = T) time=as.Date(p$'date',"%Y/%m/%d") p1=ts(p[,2:4]) install.packages('zoo') library(zoo) par(mfrow=c(1,1),oma=c(0.2,0.2,0.2,0.2)) plot(zoo(p1,time),xlab="time",ylab="价格美元", plot.type = "single",col=c("red","black","yellow","green"),lty=1:3,main="样本的价格序列"); par(mfrow=c(1,3),oma=c(0.2,0.2,0.2,0.2)) plot(zoo(p1[,1],time),xlab="time",ylab="价格美元",main="比特币价格序列") plot(zoo(xt[,2],time),xlab="time",ylab="价格美元",main="以太坊价格序列") plot(zoo(xt[,3],time),xlab="time",ylab="价格美元",main="瑞波币价格序列") 

计算相关系数矩阵

p1cor=cor(p1) 相关系数矩阵 p1cor 

计算基本统计量最大值、最小值、中位数、偏度、峰度、极值

install.packages('DistributionUtils') library(DistributionUtils) #基本统计量计算 data_outline = function(x){ 
    m = mean(x) d=max(x) xd=min(x) me = median(x) s = sd(x) kur=kurtosis(x) ske=skewness(x) R = max(x)-min(x) data.frame( Mean=m, Median=me, max=d,min=xd,std_dev=s, Skewness=ske, Kurtosis=kur, R=R) } for (i in 1:3){ 
   print(data_outline(p1[,i]))} 

绘制价格序列分布直方图

par(mfrow=c(1,3)) hist(p1[,1],main="比特币价格分布直方图",col="yellow",xlab="价格") hist(p1[,2],main="以太坊价格分布直方图",col="yellow",xlab="价格") hist(p1[,3],main="瑞波币价格分布直方图",col="yellow",xlab="价格") 

通常对时间序列对数化、差分。后面将处理后的序列称之为“收益率序列”

p1ld=diff(log(p1)) par(mfrow=c(1,1)) plot(zoo(p1ld,time), plot.type = "single",xlab="Time",ylab="value",col=c("red","black","yellow"),lty=1:3,main="样本的对数收益率线形图"); plot(zoo(p1ld[,1],time), plot.type = "single",xlab="Time",ylab="value",col=c("red"),main="比特币的对数收益率线形图"); plot(zoo(p1ld[,2],time), plot.type = "single",xlab="Time",ylab="value",col=c("black"),main="以太坊的对数收益率线形图"); plot(zoo(p1ld[,3],time), plot.type = "single",xlab="Time",ylab="value",col=c("yellow"),main="瑞波币的对数收益率线形图"); 

对收益率序列做频数分布直方图、正态图，充实一下内容

 par(mfrow=c(1,2)) #正态检验作图 hist(p1ld[,1],main="比特币对数差分收益率密度图",col="yellow",xlab="",ylim=c(0,10),probability=T) lines(density(p1ld[,1]),lwd=2); data_outline(p1ld[,1]) norm(p1ld[,1],); line(p1ld[,1]) hist(p1ld[,2],main="以太坊对数差分收益率密度图",col="yellow",xlab="",ylim=c(0,10),probability=T) lines(density(p1ld[,2]),lwd=2); data_outline(p1ld[,2]) norm(p1ld[,2],); line(p1ld[,2]) hist(p1ld[,3],main="瑞波币对数差分收益率密度图",col="yellow",xlab="",ylim=c(0,10),probability=T) lines(density(p1ld[,3]),lwd=2); data_outline(p1ld[,3]) norm(p1ld[,3],); line(p1ld[,3]) 

对收益率序列进行shapiro正态检验，p值越小，越拒绝原假设，表明序列不服从正态分布。原假设H0：数据服从正态分布

for (i in 1:3)print(shapiro.test(p1ld[,i])) 

为了防止可能造成的后续ARIMA模型伪回归的问题，对收益率序列进行ADF单位根检验，原假设序列存在单位根，当p值越小越拒绝原假设，可以得出该序列平稳。若不平稳，需在进行差分处理。再不行，额，我也不知道了

install.packages('tseries') library(tseries) for (i in 1:3)print(adf.test(p1ld[,i],alt="stationary")) 

对前面d=1阶差分平稳后的收益率序列进行ARMA（p,q）模型的定阶。看偏自相关图pacf用来确定p(我总是会忘记acf和pacf到底对应的是谁，干脆就这样记忆，pacf有一个p,所以确定p),和自相关图acf用来确定q，几阶结尾，就确定为几阶。

install.packages('stats') library(stats) for (i in 1:3){ 
    par(mfrow=c(1,2)) acf(p1ld[,i]) pacf(p1ld[,i]) } 

其实看图不如直接调用forecast程序包中的auto.arima函数根据AIC准则自动拟合ARIMA模型来的直接。当然为了扩充内容还是可以的

 install.packages('forecast') library(forecast) armamodel1=auto.arima(p1ld[,1],ic = 'aic') armamodel2=auto.arima(p1ld[,2],ic = 'aic') armamodel2 armamodel3=auto.arima(p1ld[,3],ic = 'aic') armamodel3 

Box.test(residuals(armamodel1),lag=5 ,type = "Ljung-Box") Box.test(residuals(armamodel2),lag=5 , type = "Ljung-Box") Box.test(residuals(armamodel3), lag=5 ,type = "Ljung-Box") 

进一步用LB检验残差项的平方是否存在自相关性，检验p值很小拒绝原假设，表明存在自相关性，满足建立GARCH模型的前提条件。

Box.test(residuals(armamodel1)^2,lag=5) Box.test(residuals(armamodel2)^2,lag=5) Box.test(residuals(armamodel3)^2,lag=5) 

收益率序列拟合后的ARMA模型残差会表现出异方差性。此时需要用ARCH检验来判断扰动项的条件方差对它前期方差依赖的程度，确保该序列可以建立GARCH 模型。

install.packages('MTS') library(MTS) archTest(residuals(armamodel1),lag=5) archTest(residuals(armamodel2),lag=5) archTest(residuals(armamodel3),lag=5) 

install.packages('fGarch') library(fGarch); gfit1=garchFit(~garch(1,1),data=p1ld[,1],include.mean = FALSE,trace = F,cond.dist = "norm") gfit1 gfit2=garchFit(~garch(1,1),data=p1ld[,2],include.mean = FALSE,trace = F,cond.dist = "norm") gfit2 gfit3=garchFit(~garch(1,1),data=p1ld[,3],include.mean = FALSE,trace = F,cond.dist = "norm") gfit3 

install.packages('rmgarch') library(rmgarch) meanSpec=list(armaOrder=c(0,0),include.mean=FALSE,archpow=1) distSpec=c("mvnorm") varSpec=list(model="sGARCH",garchOrder=c(1,1)) spec1=ugarchspec(mean.model=meanSpec,variance.model=varSpec) mySpec=multispec(replicate(2,spec1)) mspec=dccspec(mySpec,VAR=F,robust = F,lag=1,lag.max=NULL,lag.criterion = c("AIC"),external.regressors = NULL,robust.control = list(gamma=0.25,delta=0.01,nc=10,ns=500),dccOrder = c(1,1),distribution = distSpec,start.pars = list(),fixed.pars = list()) 上述是方程参数估计准备过程 fdcc12=dccfit(data=p1ld[,c(1,2)],mspec,out.sample = 10,solver = "solnp",solver.control = list(),fit.control = list(eval.se=TRUE,stationary=TRUE,scale=FALSE),parallel=TRUE,parallel.control=list(pkg=c("multicore"),cores=2),fit=NULL,VAR.fit=NULL) show(fdcc12) plot(fdcc12) fdcc13=dccfit(data=p1ld[,c(1,3)],mspec,out.sample = 10,solver = "solnp",solver.control = list(),fit.control = list(eval.se=TRUE,stationary=TRUE,scale=FALSE),parallel=TRUE,parallel.control=list(pkg=c("multicore"),cores=2),fit=NULL,VAR.fit=NULL) show(fdcc13) plot(fdcc13) fdcc23=dccfit(data=p1ld[,c(2,3)],mspec,out.sample = 10,solver = "solnp",solver.control = list(),fit.control = list(eval.se=TRUE,stationary=TRUE,scale=FALSE),parallel=TRUE,parallel.control=list(pkg=c("multicore"),cores=2),fit=NULL,VAR.fit=NULL) show(fdcc23) plot(fdcc23)