##数据预处理
df<-read.csv(file="上证50股票数据(20200615-20201127).csv")         
dim(df)
head(df)
y<-ts(df[,1],start = 1,frequency = 5)
n<-length(y)
t<-1:n
##数据分组
group<-1:round(n*0.7)
y.train<-ts(y[group],start = 1,frequency = 5)
y.test<-y[-group]
n1<-length(y.train)
t1<-1:n1
n2<-length(y.test)

#平稳性检验
par(mfrow=c(1,2))
plot(t1,y.train,xlab="时间",ylab = "上证50指数",lty=1,lwd=2,type="l")
acf(y.train,main="")  #自相关图
#install.packages("fUnitRoots")
library(fUnitRoots) 
#单位根检验(原假设为非平稳序列)
for (i in 1:3) {print(adfTest(y.train,type = "nc",lag=i))}
for (i in 1:3) {print(adfTest(y.train,type = "c",lag=i))}
for (i in 1:3) {print(adfTest(y.train,type = "ct",lag=i))}

#纯随机性检验
for(i in 1:2){
  print(Box.test(y.train,lag=6*i,type ="Box-Pierce"))   #原假设为白噪声数据
}

#一阶差分使得序列平稳
y.traindiff1<-diff(y.train,1)
plot(y.traindiff1,xlab="时间",ylab = "上证50指数的一阶差分",lty=1,lwd=2,type="l")
acf(y.traindiff1,main="")
#纯随机性检验
for(i in 1:2){
  print(Box.test(y.traindiff1,lag=6*i,type ="Box-Pierce"))   #原假设为白噪声数据
}
 
#自动定阶建立模型
y.train1<-y.train[35:n1]
par(mfrow=c(1,2))
acf(y.train1,main="")  #自相关图
pacf(y.train1,ylab="PACF",main="")  #自相关图
#install.packages("forecast")
library(forecast)
auto.arima(y.train1)
y.fit<-arima(y.train1, order=c(1,0,0))
y.fit

#模型残差纯随机性检验
ArchTest<-function(rtn,m=10){
# 时间序列ARCH效应的拉格朗日乘子检验
# rtn: 时间序列
# m: 选择自回归的阶数
y=(rtn-mean(rtn))^2
T=length(rtn)
atsq=y[(m+1):T]
x=matrix(0,(T-m),m)
for (i in 1:m){
x[,i]=y[(m+1-i):(T-i)]
}
md=lm(atsq~x)
summary(md)
}

for(i in 1:2){
  print(Box.test(y.fit$residuals,lag=6*i,type ="Box-Pierce"))   #原假设为白噪声数据
}
for(i in 1:6){
  print(Box.test(y.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(y.fit$residuals,i))
}

#模型系数的检验
t.value<-y.fit$coef/c(0.0943,15.9765)    #分母为系数标准差
pt(t.value,df=length(y.train1)-length(t.value),lower.tail = FALSE)

#模型的预测
y.predict<-forecast(y.fit, h=n2)
plot(y.predict)
mean(abs(y.test-y.predict$mean))
plot(ts(y,start = 1),xlab="时间",ylab = "上证50指数",lty=1,lwd=2)
points(t[-(1:n1)],y.predict$mean,lwd=2,lty=1,type = 'l',col=3)

##趋势拟合
#install.packages("splines2")
library(splines2)
fun.t<-scale(bSpline(t,df=4,degree = 2),center = TRUE,scale = FALSE)      #保证模型的可识别性
fun.t1<-fun.t[1:n1,]
ls.fit<-lm(y.train~fun.t1)
ls.beta<-ls.fit$coef
mean(abs(ls.fit$residuals))
for(i in 1:6){
  print(Box.test(ls.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(ls.fit$residuals,i))
}
mean(abs(y.test-cbind(1,fun.t[-(1:n1),])%*%ls.beta))
#par(mfrow=c(1,2))
plot(ts(y,start = 1),xlab="时间",ylab = "上证50指数",lty=1,lwd=2)
points(cbind(1,fun.t1)%*%ls.beta,lty=1,lwd=2,type="l",col=2)
points(t[-(1:n1)],cbind(1,fun.t[-(1:n1),])%*%ls.beta,lwd=2,lty=1,type = 'l',col=3)
abline(v=max(group),lty=2)
legend(80,3050,c("真实曲线","拟合曲线","预测曲线"),col=1:3,
       lty=c(1,1,1),lwd=c(2,2,2),bty='n')
##平滑法
#移动平均
#install.packages("TTR")
library(TTR)
order0<-10
ma.fit<-SMA(y.train,n=order0)
ma.predict<-function(order0,y,h=1){
  kk<-1
  ma.fit<-SMA(y,n=order0)
  n<-length(ma.fit)
  ma.pr<-0
  ma.pr[1]<-ma.fit[n]
  while(kk<h){
    ma.fit<-SMA(c(y,ma.pr),n=order0)
    n<-length(ma.fit)
    kk<-kk+1
    ma.pr[kk]<-ma.fit[n]
  }
  return(ma.pr)
}
ma.pre<-ma.predict(order0,y.train,h=n2)
mean(abs(y.test-ma.pre))
plot(ts(y,start = 1),xlab="时间",ylab = "上证50指数",lty=1,lwd=2)
points(t1,ma.fit,lty=1,lwd=2,type="l",col=2,pch=16)
points(t[-(1:n1)],ma.pre,lwd=2,lty=1,type = 'l',col=3)
abline(v=max(group),lty=2)
legend(80,3050,c("真实曲线","拟合曲线","预测曲线"),col=1:3,
       lty=c(1,1,1),lwd=c(2,2,2),bty='n')




##时间序列STL分解
y.additive<-stl(y.train,s.window="period")
plot(y.additive)
mtext("                  时间", side=1,line=3)
mtext("无规则部分           趋势部分           季节部分              数据", side=2,line=3)
y.multiplicative<-stl(log(y.train),s.window="period")
plot(exp(y.multiplicative$time.series),main="")
mtext("                  时间", side=1,line=3)
mtext("无规则部分                 趋势部分                季节部分", side=4,line=1)
#单指数模型
library(forecast)
one.fit<-ses(y.train,h=n2)
one.fit$model
for(i in 1:6){
  print(Box.test(one.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(one.fit$residuals,i))
}
mean(abs(y.test-one.fit$mean))
plot(one.fit)
#双指数模型
two.fit<-holt(y.train,h=n2,damped = TRUE)
two.fit$model
for(i in 1:6){
  print(Box.test(two.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(two.fit$residuals,i))
}
mean(abs(y.test-two.fit$mean))
plot(two.fit)
#三指数模型
three.fit<-hw(y.train, seasonal ="additive",h=n2)
three.fit$model 
for(i in 1:6){
  print(Box.test(three.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(three.fit$residuals,i))
}
mean(abs(y.test-three.fit$mean))
plot(three.fit)
#乘法模型
three.fit1<-hw(y.train, seasonal ="multiplicative",h=n2)
three.fit1$model 
for(i in 1:6){
  print(Box.test(three.fit1$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(three.fit1$residuals,lag=i))
}
mean(abs(y.test-three.fit1$mean))
plot(three.fit1)

#条件异方差
#install.packages("FinTS")
#install.packages("tseries")
library(FinTS)
library(tseries)
garch.fit<-garch(ls.fit$residuals,order = c(0,1))
summary(garch.fit)
plot(garch.fit)

garch.fit$coef[1]+ls.fit$residuals*garch.fit$coef[2]
predict(garch.fit)
plot(ts(garch.fit$fitted.values),xlab="时间",ylab="-------------------------",main="")
for(i in 1:6){
  print(Box.test(ls.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(ls.fit$residuals,i))
}


#门限自回归模型
#install.packages("TSA")
library(TSA)
setar.fit<-tar(y.train,p1=8,p2=8,d=1,print = TRUE)
setar.pre<-predict(setar.fit,n.ahead = n2)
mean(abs(y.test-setar.pre$fit))
for(i in 1:6){
  print(Box.test(setar.fit$residuals^2,lag=i,type ="Ljung-Box"))   #原假设为不存在条件异方差
  print(ArchTest(setar.fit$residuals,i))
}