#a#数据预处理
df<-read.csv(file="上证50股票数据(20200615-20201127).csv")         
dim(df)
head(df)
X<-as.matrix(df[,-1])
y<-df[,1]
n<-length(y)
p<-ncol(X)
##协变量数据标准化与分组
group<-1:round(n*0.7)
X.train<-X[group,]
y.train<-y[group]
n1<-length(y.train)
X.test<-X[-group,]
y.test<-y[-group]


##建立多元线性模型
ls.fit<-lm(y.train~.+0,data=data.frame(X.train))
ls.beta<-coef(ls.fit)
summary(ls.fit)        
mean(abs(y.train-predict(ls.fit)))
mean(abs(y.test-X.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],X.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("car")
library(car)
crPlots(ls.fit)                               
##残差的正态性检验
par(mfrow=c(1,2))
qqPlot(ls.fit,xlab="t分位数",ylab="学生化残差",main = "QQ图")   #学生化残差QQ图
shapiro.test(rstudent(ls.fit))     #学生化残差正态性检验(原假设为正态分布)
residplot <- function(fit, nbreaks=10) {
  z <- rstudent(fit)
  hist(z, breaks=nbreaks, freq=FALSE,
       xlab="学生化残差",ylab="密度",
       main="误差分布图")
  rug(jitter(z), col="brown")
  curve(dnorm(x, mean=mean(z), sd=sd(z)),
        add=TRUE, col="blue", lwd=2)
  lines(density(z)$x, density(z)$y,
        col="red", lwd=2, lty=2)
  legend("topleft",
         legend = c( "正态曲线", "核密度曲线"),
         lty=1:2, col=c("blue","red"), cex=.7,bty='n')
}
residplot(ls.fit)                             #学生化残差密度曲线

##自相关性和异方差检验
ncvTest(ls.fit)                               #异方差得分检验(原假设同方差)
set.seed(123)
durbinWatsonTest(ls.fit)                      #自相关性的D-W检验(原假设为不相关)

#离群点检验
outlierTest(ls.fit)                           #Bonferroni离群点检验.
#高杠杆点判断
hatvalues(ls.fit)                             #杠杆值
which((n1-p)/p*(hatvalues(ls.fit)-1/n1)/(1-hatvalues(ls.fit))>qf(0.95,p,n1-p))  #判断高杠杆点

#b#数据预处理
df<-read.csv(file="上证50股票数据(20200615-20201127).csv")         
dim(df)
head(df)
X<-as.matrix(df[,-1])
y<-df[,1]
n<-length(y)
p<-ncol(X)
##协变量数据标准化与分组
group<-1:round(n*0.7)
X.train<-X[group,]
y.train<-y[group]
n1<-length(y.train)
X.test<-X[-group,]
y.test<-y[-group]

##建立多元线性模型
ls.fit<-lm(y.train~.+0,data=data.frame(X.train))

#install.packages("car")
library(car)

##异常点判读
par(mfrow=c(2,1))
rstandard(ls.fit)                             #标准化残差
plot(abs(rstandard(ls.fit)),pch=16,xlab = '指标',ylab = '标准化残差绝对值',ylim=c(0,3.5))
abline(h=sd(abs(rstandard(ls.fit)))*2,lwd=2,col=3)
abline(h=sd(abs(rstandard(ls.fit)))*3,lwd=2,col=2)
for (i in 1:max(group)) {
  lines(c(i,i),c(0,abs(rstandard(ls.fit)[i])),col=1)
  if(abs(rstandard(ls.fit)[i])>sd(abs(rstandard(ls.fit)))*3){text(i,abs(rstandard(ls.fit)[i])+0.3,i)}
}
rstudent(ls.fit)                              #学生化残差
plot(abs(rstudent(ls.fit) ),pch=16,xlab = '指标',ylab = '学生化残差绝对值',ylim=c(0,3.8))
abline(h=sd(abs(rstudent(ls.fit) ))*2,lwd=2,col=3)
abline(h=sd(abs(rstudent(ls.fit) ))*3,lwd=2,col=2)
for (i in 1:max(group)) {
  lines(c(i,i),c(0,abs(rstudent(ls.fit) [i])),col=1)
  if(abs(rstudent(ls.fit)[i])>sd(abs(rstudent(ls.fit)))*3){text(i,abs(rstudent(ls.fit)[i])+0.3,i)}
}

#删除异常点后的模型
influ<-which((n1-p-1)*rstandard(ls.fit)^2/(n1-p-rstandard(ls.fit)^2)>qf(0.95,1,n1-p-1)) #异常点判断
inls.fit<-lm(y.train[-influ]~.+0,data=data.frame(X.train[-influ,]))
inls.beta<-coef(inls.fit)
mean(abs(y.train[-influ]-predict(inls.fit)))
mean(abs(y.test-X.test%*%inls.beta))

influencePlot(ls.fit,xlab = '帽子值',ylab = '学生化残差')                         #回归影响图

par(mfrow=c(2,2))
dffits(ls.fit)                                #WK距离
plot(abs(dffits(ls.fit)),pch=16,xlab = '指标',ylab = 'WK距离',ylim = c(0,9))
abline(h=sqrt(p/n1)*2,lwd=2,col=2)
for (i in 1:max(group)) {
  lines(c(i,i),c(0,abs(dffits(ls.fit)[i])),col=1)
  if(abs(dffits(ls.fit)[i])>sqrt(p/n1)*2){text(i,abs(dffits(ls.fit)[i])+0.4,i)}
}
cooks.distance(ls.fit)                        #Cook距离
plot(cooks.distance(ls.fit),pch=16,xlab = '指标',ylab = 'Cook距离',ylim = c(0,1.5))
abline(h=4/(n1-p-1),lwd=2,col=2)
for (i in 1:max(group)) {
  lines(c(i,i),c(0,abs(cooks.distance(ls.fit)[i])),col=1)
  if(cooks.distance(ls.fit)[i]>4/(n1-p-1)){text(i,cooks.distance(ls.fit)[i]+0.07,i)}
}
#基于Cook距离删除异常点后的模型
influcook<-which(cooks.distance(ls.fit)>4/(n1-p-1))
inlscook.fit<-lm(y.train[-influcook]~.+0,data=data.frame(X.train[-influcook,]))
inlscook.beta<-coef(inlscook.fit)
mean(abs(y.train[-influcook]-predict(inlscook.fit)))
mean(abs(y.test-X.test%*%inlscook.beta))

covratio(ls.fit)                              #CovRatio准则
plot(covratio(ls.fit),pch=16,xlab = '指标',ylab = 'CovRatio准则',ylim = c(0,35))
abline(h=1,lwd=2,col=3)
abline(h=10,lwd=2,col=2)
for (i in 1:max(group)) {
  lines(c(i,i),c(0,covratio(ls.fit)[i]),col=1)
  if(covratio(ls.fit)[i]>10){text(i,covratio(ls.fit)[i]+2,i)}
}
#基于CovRatio准则删除异常点后的模型
influcovratio<-which(covratio(ls.fit)>10)
inlscovratio.fit<-lm(y.train[-influcovratio]~.+0,data=data.frame(X.train[-influcovratio,]))
inlscovratio.beta<-coef(inlscovratio.fit)
mean(abs(y.train[-influcovratio]-predict(inlscovratio.fit)))
mean(abs(y.test-X.test%*%inlscovratio.beta))

#平均拟合度量
MF.t<-function(t0,fit){
  sigam.hat<-summary(fit)$sigma
  beta<-coef(fit)
  sigam.hatI<-sqrt((n1-p-rstandard(fit)^2)/(n1-p-1)*sigam.hat^2)
  MF<-0
  for (i in 1:n1){
    beta.i<-beta-fit$residuals[i]/(1-hatvalues(fit)[i])*solve(t(X.train)%*%X.train)%*%X.train[i,]
    MF[i]<-mean(abs((y.train-X.train%*%beta.i)/(sigam.hatI[i]*sqrt(hatvalues(fit)[i])))^t0)^(1/t0)
  }
  return(MF)
}
MFt<-MF.t(2,ls.fit)
plot(MFt,pch=16,xlab = '指标',ylab = '平均拟合度量',ylim = c(0,1.5))
abline(h=mean(MFt)+sd(MFt),lwd=2,col=2)
for (i in 1:max(group)) {
  lines(c(i,i),c(0,MFt[i]),col=1)
  if(MFt[i]>mean(MFt)+sd(MFt)){text(i,MFt[i]+0.07,i)}
}
#基于平均拟合度量删除异常点后的模型
influmft<-which(MFt>mean(MFt)+sd(MFt))
inlsmft.fit<-lm(y.train[-influmft]~.+0,data=data.frame(X.train[-influmft,]))
inlsmft.beta<-coef(inlsmft.fit)
mean(abs(y.train[-influmft]-predict(inlsmft.fit)))
mean(abs(y.test-X.test%*%inlsmft.beta))

#分位数估计
#install.packages("quantreg")
library(quantreg)
qr.fit<-rq(y.train~X.train+0, tau=0.3)
qr.beta<-coef(qr.fit)
mean(abs(y.train-X.train%*%qr.beta))
mean(abs(y.test-X.test%*%qr.beta))

plot(1:n,y,type = 'o',xlab="时间",ylab = "上证50指数",lty=1,cex=0.8,pch=16)
points((1:n)[group],X.train%*%qr.beta,type="o", col="blue",lty=2,cex=0.8,pch=0)
points((1:n)[-group],X.test%*%qr.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')

##多重共线性的检验
vif(ls.fit)
colnames(X)[which(sqrt(vif(ls.fit))>10)]      
eigen(t(X.train)%*%X.train)$value
kappa(ls.fit)

#c# 有偏估计

d#数据预处理
df<-read.csv(file="上证50股票数据(20200615-20201127).csv")         
dim(df)
head(df)
X<-as.matrix(df[,-1])
y<-df[,1]
n<-length(y)
p<-ncol(X)
##协变量数据标准化与分组
group<-1:round(n*0.7)
X.train<-X[group,]
y.train<-y[group]
n1<-length(y.train)
X.test<-X[-group,]
y.test<-y[-group]

##建立多元线性模型
ls.fit<-lm(y.train~.+0,data=data.frame(X.train))

#(1)stein估计
sigam.hat<-summary(ls.fit)$sigma
tau.hat<-sigam.hat^2/sum(ls.beta^2)*sum(1/eigen(t(X.train)%*%X.train)$value)
c0<-ifelse(tau.hat>0.25,0,0.5+sqrt(0.25-tau.hat))
dmax<-2*(n1-p)/(n-p+2)*(min(eigen(t(X.train)%*%X.train)$value)*sum(1/eigen(t(X.train)%*%X.train)$value)-2)
1-dmax*sigam.hat^2/sum(predict(ls.fit)^2)
c0<-0.9999
stein.beta<-c0*ls.beta
mean(abs(y.train-X.train%*%stein.beta))
mean(abs(y.test-X.test%*%stein.beta))

#(2)岭估计
#install.packages("MASS")
library(MASS)                                     
model<-lm.ridge(y.train~.+0,data=data.frame(X.train),lambda =seq(0,1,length=50))
plot(model)
select(model)                                        ###模型的选择
lambda0<-model$kHKB
ridge.fit<-lm.ridge(y.train~.+0,data=data.frame(X.train),lambda =lambda0)
ridge.beta<-coef(ridge.fit)
mean(abs(y.train-X.train%*%ridge.beta))
mean(abs(y.test-X.test%*%ridge.beta))

#(3)Liu估计
#install.packages("bootstrap")
library(bootstrap)
beta.fit.liu<-function(X,y,x,d0){
  p<-ncol(X)
  beta<-solve(t(X)%*%X+diag(1,p))%*%(t(X)%*%y+d0*x)
  return(beta)
}
beta.predict<-function(x,X){
  return(X%*%x)
}
d<-seq(0,1,0.01)
turning<-0
for (la in 1:length(d)) {
  cv.liu<-crossval(X.train,y.train,beta.fit.liu,beta.predict,ls.beta,d[la],ngroup = 10)
  turning[la]<-mean(abs(y.train-cv.liu$cv.fit))
}
d0<-d[which.min(turning)]
liu.beta<-beta.fit.liu(X.train,y.train,ls.beta,d0)
mean(abs(y.train-X.train%*%liu.beta))
mean(abs(y.test-X.test%*%liu.beta))

#(4)主成分估计
pr<-princomp(X.train,scores=TRUE)                       ###主成分分析
screeplot(pr,type="lines",main="碎石图")              ###碎石图
Cumulative.var<-0
for (i in 1:p) {
  Cumulative.var[i]<-sum(pr$sdev[1:i])/sum(pr$sdev)
}
SVD<-svd(t(X.train)%*%X.train)
SVD$d
#r<-min(which(Cumulative.var>=0.99))   
r<-30
U<-SVD$u
Z.train<-X.train%*%U[,1:r]
pcr.beta<-U[,1:r]%*%matrix(lm(y.train~Z.train+0)$coef,ncol=1)
mean(abs(y.train-X.train%*%pcr.beta))
mean(abs(y.test-X.test%*%pcr.beta))

#(5)单参数主成分估计
SPPCR<-function(X,y,theta,adjust){
  #if(theta>adjust){print("theta必须小于adjust");break}
  SVD<-svd(t(X)%*%X)
  U<-SVD$u
  r<-max(which(SVD$d>=adjust))
  A_theta<-function(theta,r){
    A<-diag(SVD$d*theta)
    diag(A)[1:r]<-(SVD$d[1:r]-adjust+theta)/SVD$d[1:r]
    return(A)
  }
  A<-A_theta(theta,r)
  Z<-X%*%U
  alpha_hat<-lm(y~Z+0)$coef
  beta_bo<-U%*%A%*%matrix(alpha_hat,ncol=1)
  return(beta_bo)
}
adjust<-1
theta<-seq(SVD$d[p],adjust,0.001)
turning<-0
for (la in 1:length(theta)) {
  cv.single<-crossval(X.train,y.train,SPPCR,beta.predict,theta[la],adjust,ngroup = 10)
  turning[la]<-mean(abs(y.train-cv.single$cv.fit))
}
theta0<-theta[which.min(turning)]
singlepcr.beta<-SPPCR(X.train,y.train,theta0,adjust)
mean(abs(y.train-X.train%*%singlepcr.beta))
mean(abs(y.test-X.test%*%singlepcr.beta))


##非负约束估计
#install.packages("nnls")
library(nnls)
nls.fit<-nnls(X.train,y.train)
nls.beta<-coef(nls.fit)
mean(abs(y.train-X.train%*%nls.beta))
mean(abs(y.test-X.test%*%nls.beta))

plot(1:n,y,type = 'o',xlab="时间",ylab = "上证50指数",lty=1,cex=0.8,pch=16)
points((1:n)[group],X.train%*%nls.beta,type="o", col="blue",lty=2,cex=0.8,pch=0)
points((1:n)[-group],X.test%*%nls.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')