猪肉收益率ARIMA-GARCH_arima_garch_R语言_猪肉价格_ARIMA-GARCH_

setwd("C:/Users/lenovo/Desktop/毕业论文") #更改路径 data<-read.table("生猪出厂价.csv",sep=',',header=T) #画时序图 library(xts) library(plotrix) library(tseries) library(forecast) library(zoo) library(FinTS) library(timeDate) library(timeSeries) library(fBasics) library(fGarch) library(rugarch) library(quantmod) data$Week=as.Date(data$Week,"%Y-%m-%d") ### #### pre<-xts(data$Pigexfactoryprc,data$Week) rate1=vector() for (j in 2:length(data$Pigexfactoryprc)) { rate1[j]=(data$Pigexfactoryprc[j]-data$Pigexfactoryprc[j-1])/data$Pigexfactoryprc[j-1] rate1[j]=rate1[j]*100 } rate1[1]=rate1[2] pre=xts(rate1,data$Week) #绘时序图 plot(pre,col=4,pch=8,lwd=2) #LB检验 p6<-Box.test(pre,type="Ljung-Box",lag=6) p12<- Box.test(pre,type="Ljung-Box",lag=12) #自相关、偏相关检验 li<-lm(data$Pigexfactoryprc~data$Week) acf(rstudent(li),main="acf自相关系数") pacf(rstudent(li),main="pacf自相关系数") #单位根检验 adf.test(pre) #差分,确定差分的阶 pre1=diff(pre,differences = 1) pre1=pre1[-1] p<-adf.test(pre1) i=2 while (p$p.value>0.05) { pre1=diff(pre,differences = i) pre1=pre1[-(1:i)] p<-adf.test(pre1) i=i+1 } i=i-1#阶 plot(pre1,col=4,pch=8,lwd=2) p<-adf.test(pre1) #自相关、偏相关检验 acf(pre1,main="差分后acf自相关系数") pacf(pre1,main="差分后pacf自相关系数") #系统自动定阶 au=auto.arima(pre) #拟合ARIMA模型 x.fit=arima(pre,order = c(3,1,2)) x.fit #残差白噪音检查 for (j in 1:5) { print(Box.test(x.fit$residual,lag=6*j)) } #预测 x.fore=forecast(x.fit,h=12,level = c(99.9)) plot(x.fore) lines(x.fore$fitted,col="green",lty=2) lines(pre,col="green",lty=2) #portmanteau for (j in 1:5) { print(Box.test(x.fit$residual^2,type="Ljung-Box",lag=6*i)) } #ARCH-LM ArchTest(x.fit$residual,demean = FALSE) #GARCH(1,1) r.fit<-garchFit(~1+garch(1,1),data=x.fit$residual,trace=F) #拟合GARCH（1,1）模型 summary(r.fit) ### r.fore=predict(r.fit) plot(r.fore$meanForecast) #条件异方差置信区间和方差齐性置信区间比较图示 plot(pre) lines(r.fore[,1],col=4,lty=2) lines(r.fore[,2],col=4,lty=2) abline(h=1.96*sd(pre),col=4,lty=2) abline(h=-1.96*sd(pre),col=4,lty=2) #IGARCH int_igarch.sec= ugarchspec(variance.model = list(model="iGARCH"), mean.model = list(armaOrder=c(0,0),include.mean=T)) ir.fit=ugarchfit(spec=int_igarch.sec, data=x.fit$residual)#igarch i.usd=ugarchforecast(ir.fit,n.ahead = 10)#预测 #GARCH-M int_garchm.sec= ugarchspec( variance.model = list(model="sGARCH", garchOrder=c(1,1)), mean.model = list(armaOrder=c(0,0),include.mean=T,archm=T,archpow=2)) ugarchfit(spec=int_garchm.sec, data=x.fit$residual*100) #TGARCH int_igarch.sec=ugarchspec(variance.model = list(model="gjrGARCH",garchOrder=c(1,1)), mean.model =list(armaOrder=c(0,0))) ugarchfit(spec=int_igarch.sec, data=x.fit$residual) ##### ##### library(ggplot2) df<- as.data.frame(r.fit) ggplot( df, aes(x = t, y = mean)) + geom_line() + geom_ribbon( aes(x = t, ymin = mean - 2 * mean.se, ymax = mean + 2 * mean.se), color = "grey", alpha = 0.5) + geom_hline(color = "blue", yintercept = 0) + coord_cartesian(ylim = c(-0.5, 0.5)) ##### resi<-residuals(r.fit,standardize=T)#残差走势图 time<-c(1:3000)/240+2006 plot(time,resi,xlab="year",ylab="stand-resi,",type="l",lty=1,col="blue",main="标准化残差时序图") acf(resi,lag=40) pacf(resi^2,lag=40) a<-predict(r.fit,100)#预测 names(a) a1<-a$meanError#预测波动率走势图 date<-c(1:100)/240+2018.3583 plot(date,a1,xlab="time",ylab="预测的波动率",type="l",lty=3,col="blue",main="预测的波动率走势图") pred.model <- predict(r.fit, n.ahead = 10, trace =FALSE, mse = 'cond', plot=FALSE) View(pred.model)