#相关系数临界值计算
Ralpha=function(df,alpha=0.05){
ta=qt(alpha/2,df)
Rs=ta^2/(df+ta^2)
sqrt(Rs)
}
Ralpha(1)
Ralpha(1,0.01)
Ralpha(2)
Ralpha(2,0.01)

#案例6.1
rm(list=ls())
read.csv("szzs.csv",header=T)->DA#文件的第一行是变量名
lm(DA$cf~DA$jr+DA$dc)->lm.sz#作线性回归，并存入变量lm.sz
summary(lm.sz)#显示模型的基本数据
matrix(c(DA[,2],DA[,23],DA[,24]),408,3)#显示本例分析数据部分
### 求线性模型系数的区间估计
beta.int<-function(fm,alpha=0.05){ A<-summary(fm)$coefficients; df<-fm$df.residual
left<-A[,1]-A[,2]*qt(1-alpha/2, df); right<-A[,1]+A[,2]*qt(1-alpha/2, df)
rowname<-dimnames(A)[[1]]; colname<-c("Estimate", "Left", "Right")
matrix(c(A[,1], left, right), ncol=3, dimnames = list(rowname, colname ))}
beta.int(lm.sz)

#案例6.2 沪深300
rm(list=ls())
DA=read.csv("hs300.csv",header=T)
lm.reg=lm(Y~0+.,data=DA)#无常数项去掉0+就含有常数项
summary(lm.reg)
D=scale(DA[,1:300])
R=cov(D)
eigen(R)$values
#计算残差平方和
x=as.matrix(DA[,1:300])
y=DA[,301]
coef<- coef(lm.reg)[2:301]
coef
r<- y-491.4-as.vector(coef%*%t(x))
Res<-t(r)
t(r)%*%r#计算残差平方和

#案例6.3 上证50
rm(list=ls())
DA=read.csv("sz50.csv",header=T)
DA=DA[3:53]
lm.reg=lm(Y~0+.,data=DA)
summary(lm.reg)
D=scale(DA[,1:50])
R=cov(D)
eigen(R)$values
#计算残差平方和
x=as.matrix(DA[,1:50])
y=DA[,51]
coef<- coef(lm.reg)[1:50]
coef
r<- y-as.vector(coef%*%t(x));
RSS=t(r)%*%r;#计算残差平方和
Res<-t(r)
shapiro.test(Res)#正态性检验

#指数跟踪图示
y.hat=as.vector(coef%*%t(x))
Res.n=r/sqrt((RSS/(nrow(DA)-ncol(DA))))
Time<-1:nrow(DA)
par(xpd=T)
plot(Time,y.hat,type="l",col="red",xlab="Time",ylab="上证50指数")
lines(Time,y,type="l",col="blue",lty=2)
legend("bottomright",expression("拟合值" , "实际值"),col=c("red","blue"),lty=c(1,2),box.lty=0)

#通过预测的指数跟踪效果
#前1000个数据用于拟合
sz50=read.csv("sz50.csv",header=T)
DA=sz50[1:1000,3:53]
lm.reg=lm(Y~0+.,data=DA)
summary(lm.reg)
D=scale(DA[,1:50])
R=cov(D)
eigen(R)$values
#计算残差平方和
x=as.matrix(DA[,1:50])
y=DA[,51]
coef<- coef(lm.reg)[1:50]
coef
r<- y-as.vector(coef%*%t(x));
RSS=t(r)%*%r;RSS#计算残差平方和
#残差图
Time=1:1000
plot(Time,r,type="l")
lines(c(-15,1000),c(0,0))

#预测后面的2168个指数数据
sz50.pre=sz50[1001:3168,3:53]
x.pre=sz50.pre[,1:50]
y.pre=sz50.pre[,51]
y.hat.pre=as.vector(coef%*%t(x.pre))
r.pre<- y.pre-y.hat.pre;
Res.pre<-t(r.pre);
RSS.pre=t(r.pre)%*%r.pre;RSS.pre#计算残差平方和
#残差图
Time.pre<-1:nrow(sz50.pre)
plot(Time.pre,Res.pre,type="l")
lines(c(-15,2168),c(0,0))
#预测指数跟踪图
par(xpd=T)
plot(Time.pre,y.hat.pre,type="l",col="red",xlab="Time",ylab="上证50指数")
lines(Time.pre,y.pre,type="l",col="blue",lty=2)
legend("bottomright",expression("预测值" , "实际值"),col=c("red","blue"),lty=c(1,2),box.lty=0)


#案例6.4 沪深300正回归拟合
#pls计算，参考沪深300最小市值股票/总市值=692636/1854621409=0.000373465，以此为阈值，利用nnls计算系数
rm(list=ls())
library(nnls)
DA=read.csv("hs300.csv",header=T)
x=as.matrix(DA[,1:300])
y=as.vector(DA[,301])
a=1:300
b=0.000373465
D=data.frame(a,b)
b0=as.vector(D[,2])
y1=y-x%*%b0
nnr<-nnls(x,y1)
nnr
coef<- coef(nnr)+0.000373465
coef
r<- y-as.vector(coef%*%t(x))
Res<-t(r)
t(r)%*%r#计算残差平方和
#图6.3 正回归残差图
Time<-1:480
plot(Time,Res,type="l")
lines(c(-15,515),c(0,0))

#图6.4
yv<-t(y1)
pv<-t(as.vector(coef%*%t(x)))
Time<-1:480
plot(Time,yv,type="l",col="red",xlab="Time",ylab="沪深300指数")
lines(Time,yv,type="l",col="blue",lty=2)
legend("bottomright",expression("预测值" , "实际值"),col=c("red","blue"),lty=c(1,2),box.lty=0)

#案例6.5 汇率跟踪一篮子货币
rm(list=ls())
DA=read.csv("rmbdata.csv",header=T)
DA=DA[2:12]
lm.rmb=lm(CNYUSD~.,data=DA)
summary(lm.rmb)

#逐年篮子预测
DA1=DA[23:283,]
lm.rmb1=lm(CNYUSD~.,data=DA1)
summary(lm.rmb1)
DA2=DA[284:545,]
lm.rmb2=lm(CNYUSD~.,data=DA2)
summary(lm.rmb2)
DA3=DA[546:805,]
lm.rmb3=lm(CNYUSD~.,data=DA3)
summary(lm.rmb3)
DA4=DA[806:1064,]
lm.rmb4=lm(CNYUSD~.,data=DA4)
summary(lm.rmb4)
DA5=DA[1065:1324,]
lm.rmb5=lm(CNYUSD~.,data=DA5)
summary(lm.rmb5)

#正回归
library(nnls)
x=as.matrix(DA[,1:10])
y=as.vector(DA[,11])
nnr<-nnls(x,y)
nnr
coef(nnr)

x1=as.matrix(DA1[,1:10]);y1=as.vector(DA1[,11]);nnr1<-nnls(x1,y1);coef(nnr1)
x2=as.matrix(DA2[,1:10]);y2=as.vector(DA2[,11]);nnr2<-nnls(x2,y2);coef(nnr2)
x3=as.matrix(DA3[,1:10]);y3=as.vector(DA3[,11]);nnr3<-nnls(x3,y3);coef(nnr3)
x4=as.matrix(DA4[,1:10]);y4=as.vector(DA4[,11]);nnr4<-nnls(x4,y4);coef(nnr4)
x5=as.matrix(DA5[,1:10]);y5=as.vector(DA5[,11]);nnr5<-nnls(x5,y5);coef(nnr5)