df<-read.csv(file="银行板块财务数据(2014-2019).csv")         
dim(df)
head(df)

X<-df[,-c(1:3,6,9)]
X1<-scale(X)
y<-df[,1]
n<-length(y)
p<-ncol(X)

##数据分组
group<-which(df[,3]!=2019)
X.train<-X1[group,];Z.train<-as.matrix(cbind(1,X1[group,]))
y.train<-y[group]
n1<-length(y.train)
X.test<-X1[-group,];Z.test<-as.matrix(cbind(1,X1[-group,]))
y.test<-y[-group]

#install.packages("geepack")
#install.packages("MESS")
#install.packages("car")
library(nlme)
library(geepack)
library(MESS)
library(MASS)
library(car)
#定义个体编号
id.names<-1:n1
for (i in 1:length(table(df[group,2]))) {
  id<-which(df[group,2]==names(table(df[group,2]))[i])
  id.names[id]<-rep(i,length(id))
}
id.names

##建立纵向数据的线性模型
#最小二乘法
ls.fit<-lm(y.train~.,data = data.frame(X.train))
summary(ls.fit)
ls.beta<-coef(ls.fit)
mean(abs(y.train-Z.train%*%ls.beta))
mean(abs(y.test-Z.test%*%ls.beta))

#加权最小二乘法或极大似然估计(等相关结构)
gls.ml<-gls(y.train~.,data = data.frame(X.train), correlation = corCompSymm(form = ~1|id.names), method="ML")
summary(gls.ml)
glsml.beta<-coef(gls.ml)
mean(abs(y.train-Z.train%*%glsml.beta))
mean(abs(y.test-Z.test%*%glsml.beta))

#约束极大似然估计(等相关结构)
gls.reml<-gls(y.train~.,data = data.frame(X.train),correlation = corCompSymm(form = ~1|id.names), method="REML")
summary(gls.reml)
glsreml.beta<-coef(gls.reml)
mean(abs(y.train-Z.train%*%glsreml.beta))
mean(abs(y.test-Z.test%*%glsreml.beta))

cbind(ls.beta,glsml.beta,glsreml.beta)



##边际模型
#判断自变量与因变量是否为线性关系
crPlots(ls.fit,main="",ylab="")  

#独立的工作矩阵
gee.ind<-geese(y.train~.,mean.link = "log",data =data.frame(X.train), id=id.names, family=gaussian, corstr = "independence")
summary(gee.ind)
geeind.beta<-gee.ind$beta
mean(abs(y.train-exp(Z.train%*%geeind.beta)))
mean(abs(y.test-exp(Z.test%*%geeind.beta)))

#等相关矩阵
gee.exc<-geese(y.train~.,mean.link = "log",data =data.frame(X.train), id=id.names, family=gaussian, corstr = "exchangeable")
summary(gee.exc)
geeexc.beta<-gee.exc$beta
mean(abs(y.train-exp(Z.train%*%geeexc.beta)))
mean(abs(y.test-exp(Z.test%*%geeexc.beta)))

#AR(1)
gee.ar1<-geese(y.train~.,mean.link = "log",data =data.frame(X.train), id=id.names, family=gaussian, corstr = "ar1")
summary(gee.ar1)
geear1.beta<-gee.ar1$beta
mean(abs(y.train-exp(Z.train%*%geear1.beta)))
mean(abs(y.test-exp(Z.test%*%geear1.beta)))

cbind(geeind.beta,geeexc.beta,geear1.beta)