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)
评论7