##数据预处理
df<-read.csv(file="上证50股票数据(20200615-20201213).csv")         
dim(df)
head(df)
x1<-df$"中国石油"
x2<-df$"洛阳钼业"

residplot <- function(fit, nbreaks=10) {
  z <- fit
  hist(z, breaks=nbreaks, freq=FALSE,
       xlab="Stock Price",
       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( "Normal Curve", "Kernel Density Curve"),
         lty=1:2, col=c("blue","red"), cex=.7,bty='n')
}
par(mfrow=c(1,2))
residplot(x1)
mtext("股票价格", side=1,line=2)
mtext("密度", side=2,line=2)  
residplot(x2)  
mtext("股票价格", side=1,line=2)
mtext("密度", side=2,line=2)  
#分位数估计置信区间
alpha<-0.05
n<-length(x1)
p<-0.5
x1.qu<-quantile(x1,p) #分位数估计
x1.qu
x<-sort(x1)
r<-qbinom(alpha/2,n,p)
r
s<-qbinom(1-alpha/2,n,p)
s
print(paste0(p,"分位数的置信区间为","[",x[r],",",x[s],"]"))
#修正的分位数估计置信区间
r1<-floor(n*p+0.5+qnorm(alpha/2)*sqrt(n*p*(1-p)))
r1
s1<-ceiling(n*p-0.5+qnorm(1-alpha/2)*sqrt(n*p*(1-p)))
s1
print(paste0(p,"分位数的置信区间为","[",x[r1],",",x[s1],"]"))

##两样本均值比较
RI<-0
n1<-length(x1)
n2<-length(x2)
for(i in 1:n1){
  RI[i]<-length(which(x1[i]<x2))
}
Un.12<-sum(RI)/(n1*n2)
w<-0.5+qnorm(1-alpha)*sqrt((n1+n2)/(12*n1*n2))
if(Un.12>w){print("拒绝原假设,即认为x2比x1收盘价高")}else{print("不能拒绝原假设,即认为x1比x2收盘价高")}

##非参数检验
wilcox.test(x1,x2,alternative = "greater")             #符号秩和检验
wilcox.test(x1,x2,alternative = "greater",paired=TRUE) #秩和检验
mood.test(x1,x2)
cor.test(x1,x2,method = "pearson")  #原假设不相关
cor.test(x1,x2,method = "kendall")  #原假设不相关

##非参数回归
PZH<-read.csv(file="PZH.csv")
PZH<-data.frame(PZH)
YNBY<-read.csv(file="YNBY.csv")
YNBY<-data.frame(YNBY)
PZH.price<-PZH$X600436.ss.close
YNBY.price<-YNBY$X000538.sz.close
library(ggplot2)
dat<-data.frame(PZH.price,YNBY.price)
ggplot(data = dat,aes(x=PZH.price,y=YNBY.price))+
  geom_point(cex=2)

##分组
#[0,1]区间化
X<-PZH.price
y<-YNBY.price
n<-length(y)
set.seed(361)
group<-sample(1:n,round(n*0.6))
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~X.train)
summary(ls.fit)
ls.beta<-coef(ls.fit)
yhat.train<-Z.train%*%ls.beta
dat.ls<-cbind(X.train,yhat.train)
mean(abs(y.train-Z.train%*%ls.beta))
mean(abs(y.test-Z.test%*%ls.beta))

##N-W估计
#拇指法则选择带宽
Rhat<-quantile(X.train,0.75)-quantile(X.train,0.25)
h<-1.06*min(sd(X.train),Rhat/1.34)*n1^(-1/5)
NW.pre<-ksmooth(X.train,y.train,kernel = "normal",bandwidth=h/0.25,range.x = range(X))
NW.pre
NW.test<-ksmooth(X.train,y.train,kernel = "normal",bandwidth=h/0.25,range.x = range(X),x.points = X.test)
mean(abs(y.test-NW.test$y))

##局部多项式回归
LPR.fit<-loess(y.train~.,data=data.frame(x=X.train),normalize = FALSE,degree = 1)
LPR.fit
dat.lpr<-cbind(X.train,LPR.fit$fitted)
mean(abs(y.test-predict(LPR.fit,newdata=data.frame(x=X.test))))

##样条光滑
library(splines2)
r<-3
Kn<-floor(n1^(1/(2*r+3)))
XX<-scale(bSpline(X,df=Kn+r),scale = FALSE,center = TRUE)
XX.train<-XX[group,]
XX.test<-XX[-group,]
bs.fit<-lm(y.train~XX.train)
summary(bs.fit)
bs.beta<-coef(bs.fit)
dat.bs<-cbind(X.train,cbind(1,XX.train)%*%bs.beta)
mean(abs(y.train-cbind(1,XX.train)%*%bs.beta))
mean(abs(y.test-cbind(1,XX.test)%*%bs.beta))


##函数拟合图
plot(X,y,pch=16,col=8)
points(dat.ls[order(X.train),1],dat.ls[order(X.train),2],type = "l",lwd=2,col=1,lty=1)
points(NW.pre$x,NW.pre$y,type = "l",lwd=2,col=2,lty=2)
points(dat.lpr[order(X.train),1],dat.lpr[order(X.train),2],type = "l",lwd=2,col=3,lty=4)
points(dat.bs[order(X.train),1],dat.bs[order(X.train),2],type = "l",lwd=2,col=4,lty=5)
legend("topleft",c("线性模型","N-W估计","局部线性估计","B样条估计"),col=c(1,2,3,4),
       lty=c(1,2,3,4),lwd=c(2,2,2,2),bty='n')

#测试集预测结果箱线图
dat<-cbind(abs(y.test-Z.test%*%ls.beta),abs(y.test-NW.test$y),
           abs(y.test-predict(LPR.fit,newdata=data.frame(x=X.test))),
           abs(y.test-cbind(1,XX.test)%*%bs.beta))
colnames(dat)<-c("线性模型","N-W估计","局部线性估计","B样条估计")
boxplot(dat,col=2:5)


##可加模型和部分线性可加模型
##数据转化为【0,1】
colunif=function(X){
  Y=X
  if(is.vector(Y)){return((Y-min(Y))/(max(Y)-min(Y)))}
  for (i in 1:dim(X)[2]) {
    Y[,i]<-(X[,i]-min(X[,i]))/(max(X[,i]-min(X[,i])))
  }
  return(Y)
}
##三次B样条基函数矩阵
Bbisc<-function(X,Kn,r){
  n<-nrow(X)
  p<-ncol(X)
  b0<-matrix(,n,p)
  B<-bSpline(X[,1],df=Kn+r)
  if(p==1){return(B)}
  for (i in 2:p) {
    b0<-bSpline(X[,i],df=Kn+r)
    B<-cbind(B,b0) 
  }
  return(B)
}
df<-read.csv(file="奢侈品概念股票数据(20200101-20201231).csv")   
dim(df)
head(df)
X<-as.matrix(df[,-1])
X<-colunif(X)
y<-df[,1]
n<-length(y)
p<-ncol(X)
##数据分组
set.seed(136)
group<-sample(1:n,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~X.train)
summary(ls.fit)
ls.beta<-coef(ls.fit)
mean(abs(y.train-cbind(1,X.train)%*%ls.beta))
mean(abs(y.test-cbind(1,X.test)%*%ls.beta))

#可加模型
r<-3
Kn<-floor(n1^(1/(2*r+3)))
BS1<-Bbisc(X,Kn,r)
BS<-scale(BS1,center = TRUE,scale = FALSE)
XX.train<-BS[group,]
XX.test<-BS[-group,]
am.fit<-lm(y.train~XX.train)
summary(am.fit)
am.beta<-coef(am.fit)
mean(abs(y.train-cbind(1,XX.train)%*%am.beta))
mean(abs(y.test-cbind(1,XX.test)%*%am.beta))
#可加模型中非参数估计结果
kk<-matrix(seq(1,(Kn+r)*p,by=1),ncol=p)
par(mfrow=c(3,5))
xx<-c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","x13","x14")
f<-c("f1","f2","f3","f4","f5","f6","f7","f8","f9","f10","f11","f12","f13","f14")
for (j in 1:p){
  plot(X.train[,j],XX.train[,kk[,j]]%*%am.beta[kk[,j]+1],xlab=xx[j],ylab=f[j],col=4,pch=4)
}

#部分线性可加模型
#install.packages("corrplot")
library(corrplot)
#线性相关性可视化
corrplot(cor(df[group,]),type = "upper", tl.col = "black", tl.srt = 45)
m.fun<-rep(2,p)
m.fun[abs(cor(df[group,])[1,-1])>=0.8]<-1
if(m.fun[1]==2){ZZ<-BS[,kk[,1]]}else{ZZ<-X[,1]}
for(j in 2:p){
  if(m.fun[j]==2){ZZ<-cbind(ZZ,BS[,kk[,j]])}else{ZZ<-cbind(ZZ,X[,j])}
}
ZZ.train<-ZZ[group,]
ZZ.test<-ZZ[-group,]
plam.fit<-lm(y.train~ZZ.train)
summary(plam.fit)
plam.beta<-coef(plam.fit)
mean(abs(y.train-cbind(1,ZZ.train)%*%plam.beta))
mean(abs(y.test-cbind(1,ZZ.test)%*%plam.beta))

par(mfrow=c(3,5))
k<-matrix(seq(1,(Kn+r)*p,by=1),ncol=p)
plot(X.train[,1],ZZ.train[,k[,1]]%*%plam.beta[k[,1]+1],xlab=xx[1],ylab=f[1],col=4,pch=4)
for (j in 2:p){
  if(m.fun[j]==2){
    k[,j]<-(k[Kn+r,j-1]+1):(k[Kn+r,j-1]+Kn+r)
    plot(X.train[,j],ZZ.train[,k[,j]]%*%plam.beta[k[,j]+1],xlab=xx[j],ylab=f[j],col=4,pch=4)
  }else{
    k[,j]<-rep(k[Kn+r,j-1]+1,Kn+r)
    plot(X.train[,j],ZZ.train[,k[1,j]]*plam.beta[k[1,j]+1],xlab=xx[j],ylab=f[j],col=4,pch=4)
  }
}


#变量选择后的部分线性可加模型
grouped<-function(m.fun,Kn,r){
  if(m.fun[1]==2){grouped<-rep(1,Kn+r)}else{grouped<-1}
  for (j in 2:p) {
    if(m.fun[j]==2){grouped<-c(grouped,rep(j,Kn+r))}else{grouped<-c(grouped,j)}
  }
  return(grouped)
}
library(grpreg)
set.seed(321)
cv.fit<-cv.grpreg(ZZ.train,y.train,group=grouped(m.fun,Kn,r),penalty="grSCAD",gamma=3.7)
plot(cv.fit)

vplam.beta<-cv.fit$fit$beta[,which.min(cv.fit$cve)]
mean(abs(y.train-cbind(1,ZZ.train)%*%vplam.beta))
mean(abs(y.test-cbind(1,ZZ.test)%*%vplam.beta))

par(mfrow=c(3,5))
plot(X.train[,1],ZZ.train[,k[,1]]%*%vplam.beta[k[,1]+1],xlab=xx[1],ylab=f[1],col=4,pch=4)
for (j in 2:p){
  if(m.fun[j]==2){
    k[,j]<-(k[Kn+r,j-1]+1):(k[Kn+r,j-1]+Kn+r)
    plot(X.train[,j],ZZ.train[,k[,j]]%*%vplam.beta[k[,j]+1],xlab=xx[j],ylab=f[j],col=4,pch=4)
  }else{
    k[,j]<-rep(k[Kn+r,j-1]+1,Kn+r)
    plot(X.train[,j],ZZ.train[,k[1,j]]*vplam.beta[k[1,j]+1],xlab=xx[j],ylab=f[j],col=4,pch=4)
  }
}