##数据预处理
df<-read.csv(file="上证50股票数据(20200615-20201127).csv")         
dim(df)
head(df)#查询数据集中前六行的数据
X<-as.matrix(df[,-1])
Z<-cbind(1,X)
y<-df[,1]
n<-length(y)
p<-ncol(X)
##数据分组
group<-1:round(n*0.7)
X.train<-X[group,];Z.train<-cbind(1,X[group,])
y.train<-y[group]
n1<-length(y.train)
X.test<-X[-group,];Z.test<-cbind(1,X[-group,])
y.test<-y[-group]

ls.fit<-lm(y.train~.,data=data.frame(X.train))
ls.beta<-coef(ls.fit)
summary(ls.fit)        
mean(abs(y.train-predict(ls.fit)))
mean(abs(y.test-Z.test%*%ls.beta))
#指数追踪与预测
plot(1:n,y,type="o",col="black",xlab="时间",ylab = "上证50指数",lty=1,cex=0.8,pch=16)
points(group,predict(ls.fit),type="o",col="blue",lty=2,cex=0.8,pch=0)
points((1:n)[-group],Z.test%*%ls.beta,type="o", col="red",lty=2,cex=0.8,pch=2)
abline(v=max(group),lty=2)
legend(90,3050,c("真实值","拟合值","预测值"),col=c('black','blue','red'),
       lty=c(1,2,2),pch=c(16,0,2),bty='n')

##传统的变量选择方法
#install.packages("leaps")
library(leaps)
model<-regsubsets(X.train,y.train,nbest=1,nvmax=p,method= "backward",intercept=TRUE)
par(mfrow=c(2,2))
plot(model, scale="bic")        #BIC
plot(model, scale="Cp")         #Cp准则
plot(model, scale="adjr2")      #调整的R2 
plot(model, scale="r2")         #R2 

regsu.beta<-rep(0,p+1)
names(regsu.beta)<-c("常数项",colnames(X.train))
number<-which.min(summary(model)$bic)  #采用BIC准则确定模型
regsu.beta[names(coef(model,number)[-1])]<-coef(model,number)[-1]
regsu.beta[1]<-coef(model,number)[1]
mean(abs(y.train-Z.train%*%regsu.beta))
mean(abs(y.test-Z.test%*%regsu.beta))

model1<-regsubsets(X.train,y.train,nbest=1,nvmax=10,method= "backward",intercept=TRUE)
plot(model1, scale="bic")        #BIC
regsu.beta1<-rep(0,p+1)
names(regsu.beta1)<-c("常数项",colnames(X.train))
number1<-which.min(summary(model1)$bic)  #采用BIC准则确定模型
regsu.beta1[names(coef(model1,number1)[-1])]<-coef(model1,number1)[-1]
regsu.beta1[1]<-coef(model1,number1)[1]
mean(abs(y.train-Z.train%*%regsu.beta1))
mean(abs(y.test-Z.test%*%regsu.beta1))

names(coef(model1,number1)[-1])

#逐步回归
step.fit<-step(ls.fit) #AIC
summary(step.fit)
step.beta<-rep(0,p+1)
names(step.beta)<-c("常数项",colnames(X.train))
step.beta[names(coef(step.fit)[-1])]<-coef(step.fit)[-1]
step.beta[1]<-coef(step.fit)[1]
mean(abs(y.train-Z.train%*%step.beta))
mean(abs(y.test-Z.test%*%step.beta))


##现代变量选择方法
#install.packages("lars")
#install.packages("glmnet")
#install.packages("ncvreg")
#install.packages("grpreg")
#install.packages("nnlasso")
library(lars)
library(glmnet)
library(ncvreg)
library(grpreg)
library(nnlasso)
#LASSO
lm.lasso<-lars(X.train,y.train,normalize = FALSE, intercept =TRUE,type = "lasso")
plot(lm.lasso)
mtext("选出系数比例", side=1,line=2)
mtext("标准化系数", side=2,line=2)

lassocp.beta<-c(mean(y.train-X.train%*%coef(lm.lasso)[which.min(lm.lasso$Cp),]),
                coef(lm.lasso)[which.min(lm.lasso$Cp),])
which(lassocp.beta==0)
mean(abs(y.train-Z.train%*%lassocp.beta))
mean(abs(y.test-Z.test%*%lassocp.beta))

#LASSO
set.seed(361)
cv.lasso<-cv.glmnet(X.train,y.train,alpha = 1,intercept = TRUE,standardize=FALSE)
plot(cv.lasso)
lamda0<-min(cv.lasso$lambda[which(cv.lasso$nzero==10)])
lasso.fit<-glmnet(X.train,y.train,alpha = 1,lambda=lamda0,intercept = TRUE,standardize=FALSE)
lasso.beta<-coef(lasso.fit)
which(lasso.beta!=0)
mean(abs(y.train-Z.train%*%lasso.beta))
mean(abs(y.test-Z.test%*%lasso.beta))

colnames(X)[which(lasso.beta!=0)[-1]]

#SCAD
set.seed(361)
cv.scad<-cv.grpreg(X.train,y.train,type='grSCAD',gamma=3.7)
plot(cv.scad)
cv.scad$lambda
scad.fit<-grpreg(X.train,y.train,lambda=cv.scad$lambda[36],type='grSCAD',gamma=3.7)
scad.beta<-scad.fit$beta
which(scad.beta!=0)
mean(abs(y.train-Z.train%*%scad.beta))
mean(abs(y.test-Z.test%*%scad.beta))

colnames(X)[which(scad.beta!=0)[-1]]

par(mfrow=c(1,2))
plot(cv.lasso,xlab="log(λ)",ylab="均方误差(LASSO)")
plot(cv.scad,xlab="log(λ)",ylab="均方误差(SCAD)",main="")
#EN
set.seed(361)
alpha<-seq(0.1,1,0.1)
lambda.al<-cv.al<-0
for(i in 1:length(alpha)){
  cv.en<-cv.glmnet(X.train,y.train,alpha = alpha[i],intercept = TRUE,standardize = FALSE)
  la.value<-cv.en$lambda[which(cv.en$nzero==10)]
  la.cv<-rep(Inf,length(la.value))
  for (j in 1:length(la.value)) {
    en.fit<-glmnet(X.train,y.train,alpha = alpha[i],lambda=la.value[j],intercept = TRUE,standardize = FALSE)
    en.beta<-coef(en.fit)
    if(length(which(en.beta!=0))==11){
      la.cv[j]<-cv.en$cvm[which(la.value[j]==cv.en$lambda)]
    }
  }
  lambda.al[i]<-la.value[which.min(la.cv)]
  cv.al[i]<-la.cv[which.min(la.cv)]
}
alpha0<-alpha[which.min(cv.al)]
lambda0<-lambda.al[which.min(cv.al)]
en.fit<-glmnet(X.train,y.train,alpha = alpha0,lambda=lambda0,intercept = TRUE,standardize = FALSE)
en.beta<-coef(en.fit)
which(en.beta!=0)
mean(abs(y.train-Z.train%*%en.beta))
mean(abs(y.test-Z.test%*%en.beta))

colnames(X)[which(en.beta!=0)[-1]]