#案例5.1
R=matrix(c(1,0.58,0.51,0.39,0.46,0.58,1,0.6,0.39,0.32,0.51,0.6,1,0.44,0.43,0.39,0.39,0.44,1,0.52,0.46,0.32,0.43,0.52,1),5,5)
factor=function(S, m){
       p<-nrow(S)
       diag_S<-diag(S)
       sum_rank<-sum(diag_S)
       rowname<-paste("X", 1:p, sep="")
       colname<-paste("Factor", 1:m, sep="")
       A<-matrix(0, nrow=p, ncol=m,
       dimnames=list(rowname, colname))
       eig<-eigen(S)
       for (i in 1:m)
       A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
       h<-diag(A%*%t(A))
       rowname<-c("SS loadings","Proportion Var","Cumulative Var")
       B<-matrix(0, nrow=3, ncol=m,
       dimnames=list(rowname, colname))
       for (i in 1:m){
       B[1,i]<-sum(A[,i]^2)
       B[2,i]<-B[1,i]/sum_rank
       B[3,i]<-sum(B[1,1:i])/sum_rank
       }
       list(loadings=A,
       var=cbind(common=h, spcific=diag_S-h), B=B)
       }
fa<-factor(R, m=1); fa
fa<-factor(R, m=2); fa
#求残差矩阵
A=loadings(fa)
R-A%*%t(A)-diag(c(0.34,0.19,0.31,0.27,0.22))

#也可用R内置的函数factanal()，但该函数采用的是极大似然法，不是主成分法，结果会有不同。
fa<-factanal(factors=1,covmat=R,rotation="none");fa
fa<-factanal(factors=2,covmat=R,rotation="none");fa

#案例5.2 上证股票周回升率数据
read.csv("zhsl.csv",header=T)->DA
R=cor(DA)
factor=function(S, m){
       p<-nrow(S)
       diag_S<-diag(S)
       sum_rank<-sum(diag_S)
       rowname<-paste("X", 1:p, sep="")
       colname<-paste("Factor", 1:m, sep="")
       A<-matrix(0, nrow=p, ncol=m,
       dimnames=list(rowname, colname))
       eig<-eigen(S)
       for (i in 1:m)
       A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
       h<-diag(A%*%t(A))
       rowname<-c("SS loadings","Proportion Var","Cumulative Var")
       B<-matrix(0, nrow=3, ncol=m,
       dimnames=list(rowname, colname))
       for (i in 1:m){
       B[1,i]<-sum(A[,i]^2)
       B[2,i]<-B[1,i]/sum_rank
       B[3,i]<-sum(B[1,1:i])/sum_rank
       }
       list(loadings=A,
       var=cbind(common=h, spcific=diag_S-h), B=B)
       }
fa<-factor(R, m=1); fa
fa<-factor(R, m=2); fa
fa<-factor(R, m=3); fa
#求残差矩阵
A=loadings(fa)
R-A%*%t(A)-diag(c(0.38,0.13,0.14,0.23,0.15,0.51,0.26,0.25,0.3))

#例5.2.1
R=matrix(c(1,-0.3333,0.6667,-0.3333,1,0,0.6667,0,1),3,3)
eigen(R)
A1=sqrt(1.745)*eigen(R)$vectors[,1];A1#求载荷矩阵
A2=eigen(R)$vectors[,2];A2
#作正交旋转
0.934^2;0.418^2+0.894^2;0.835^2+0.447^2
c(-0.418,0.894)/sqrt(0.974);c(0.835,0.447)/sqrt(0.897)
0.423^2-0.906^2;0.882^2-0.472^2
2*(-0.423)*0.906;2*0.882*0.472
1-0.642+0.555;0.833-0.766;1+0.642^2+0.555^2-0.766^2-0.833^2;2*(0.642*0.766+0.555*0.833)
(1.91-2*0.913*0.067/3)/(0.44-(0.913^2-0.067^2)/3)
C=0.25*atan(11.42)
cos(C);sin(C)
#旋转因子载荷矩阵
matrix(c(A1,A2),3,2)%*%matrix(c(0.932,0.362,-0.362,0.932),2,2)
((0.87^2/0.872)^2+(0.07^2/0.974)^2+(0.94^2/0.897)^2)/3-((0.87^2/0.872+0.07^2/0.974+0.94^2/0.897)/3)^2
((0.34^2/0.872)^2+(0.98^2/0.974)^2+(0.11^2/0.897)^2)/3-((0.34^2/0.872+0.98^2/0.974+0.11^2/0.897)/3)^2
fa<-factor(R, m=2); fa#用上面的函数计算，结果是一样的

#案例5.3 恒生金融指数数据
read.csv("hsdata.csv",header=T)->C
X<-C[1:22,2:10]
R=cor(X)
eigen(R)
fa<-factor(R, m=2); fa
E<- R-fa$loadings %*% t(fa$loadings)-diag(fa$var[,2]);sum(E^2)#求误差平方和Q(m)
vm1<-varimax(fa$loadings, normalize = F); vm1#方差最大的正交旋转
asin(-0.5305712)*180/pi#求旋转角度

#案例5.4 恒生金融指数数据因子得分
D<-as.matrix(X)
D=scale(D)
m1 <- cbind(D)
factanal(m1, factors = 2)
d=factanal(m1, factors = 2)
x=d$loadings
x=x[,] 
x
factanal(~D, factors =2,scores="Bartlett")$scores #求因子得分
fa=factanal(~D, factors =2,scores="Bartlett")$scores 
plot(fa[,1:2],type="n");text(fa[,1],fa[,2])#作图
#求残差矩阵
A=loadings(fa)
R-A%*%t(A)-diag(c(0.25,0.9,0.59,0.18,0.08,0.1,0.09,0.09,0.07))

#案例5.5
read.csv("s50.csv",header=T)->DA;R<-cor(DA[2:10])
fa<-factanal(~., factors=2, data=DA[2:10], scores="regression",rotation="none");fa#不作旋转
fa$scores[,1:2];plot(fa$scores[,1:2],type="n");text(fa$scores[,1],fa$scores[,2])

#案例5.6
read.csv("sz100.csv",header=T)->DA;R<-cor(DA[2:9])
fa<-factanal(~., factors=2, data=DA[2:9], scores="regression",rotation="none");fa#不作旋转
fa<-factanal(~., factors=2, data=DA[2:9], scores="regression");fa#旋转
fa$scores[,1:2];plot(fa$scores[,1:2],type="n");text(fa$scores[,1],fa$scores[,2])
fa<-factanal(~., factors=1, data=DA[2:6], scores="regression",rotation="none");fa#不作旋转
fa$scores