## ----setup, include=FALSE------------------------------------------------

## ----install-------------------------------------------------------------
list.of.packages=c("ggfortify","psych")
if(length(which(!list.of.packages %in% installed.packages()))){
  install.packages(list.of.packages[!list.of.packages %in% installed.packages()])
}
setwd("C:/Users/dujuan/Documents/workshop_MV/workshop_MV/workshopcode")
## ----Example 6.1---------------------------------------------------------

   R=matrix(c(1.00,0.02,0.96,0.42,0.01,
              0.02,1.00,0.13,0.71,0.85,
              0.96,0.13,1.00,0.50,0.11,
              0.42,0.71,0.50,1.00,0.79,
              0.01,0.85,0.11,0.79,1.00),nrow=5)
   # Customized Codes
   R.eigen=eigen(R);
   R.eigenvalue=R.eigen$values;
   R.eigenvector=R.eigen$vectors;
   m=2;
   L=R.eigenvector[,1:m]%*%diag(sqrt(R.eigenvalue[1:m]))    
   h2=rowSums(L^2)
   psi=diag(R)-h2
   cum.prop=cumsum(R.eigenvalue)/sum(diag(R))
   result1=cbind(L,h2,psi)
   dimnames(result1)=list(NULL,c("F1","F2","Community","Specific Variance"))
   result1
   
   
## ----Example 6.2---------------------------------------------------------
   stock=read.table("T8-4.DAT")
   
   # Principal Component Method: Eigen Decomposition
   R=cor(stock)
   R.eigen=eigen(R);
   R.eigenvalue=R.eigen$values;
   R.eigenvector=R.eigen$vectors;
   m=2;
   L=R.eigenvector[,1:m]%*%diag(sqrt(R.eigenvalue[1:m]))    
   h2=rowSums(L^2)
   psi=diag(R)-h2
   cum.prop=cumsum(R.eigenvalue)/sum(diag(R))
   result2=cbind(L,h2,psi)
   nam=c("JP Morgan", "Citibank", "Wells Fargo", "Royal Dutch Shell","ExxonMobil")
   dimnames(result2)=list(paste(nam,sep=""),c("F1","F2","Community","Specific Variance"))
   result2
   
   
   # Principal Component Method:Using prcomp function 
   
   R.prcomp=prcomp(stock,scale=T)
   screeplot(R.prcomp,type="l")
   library(ggfortify)
   autoplot(R.prcomp)
   m=2;
   L=R.prcomp$rotation[,1:m]%*%diag(R.prcomp$sdev[1:m])
   h2=rowSums(L^2)
   psi=diag(R)-h2
   result3=cbind(L,h2,psi)
   dimnames(result3)=list(NULL,c("F1","F2","Community","Specific Variance"))
   result3
   
   
   # MLE method: factanal is a function in package stats
   
   result=factanal(~V1+V2+V3+V4+V5, factors=2,data=stock,rotation="none")
   
   psi=result$uniquenesses
   L=result$loadings
   h2=rowSums(L^2)
   result4=cbind(L,h2,psi)
   dimnames(result4)=list(paste(nam,sep=""),c("F1","F2","Community","Specific Variance"))
   result4
   
   # MLE method+Rotation: factanal is a function in package stats
   
   result=factanal(~V1+V2+V3+V4+V5, factors=2,data=stock,rotation="varimax")
   
   Psi=result$uniquenesses
   L=result$loadings
   h2=rowSums(L^2)
   result5=cbind(L,h2,psi)
   dimnames(result5)=list(paste(nam,sep=""),c("F1","F2","Community","Specific Variance"))
   result5
   

## ----Example 6.3---------------------------------------------------------
   R=matrix(c(1, 0.439, 0.410, 0.288, 0.329, 0.248,
              0.439, 1, 0.351, 0.354, 0.320, 0.329,              
              0.410, 0.351, 1, 0.164, 0.190, 0.181,
              0.288, 0.354, 0.164, 1, 0.595, 0.470,
              0.329, 0.320, 0.190, 0.595, 1, 0.464,
              0.248, 0.329, 0.181, 0.470, 0.464, 1),nrow=6);

   # No rotation
   result=factanal(factors=2,covmat=R,rotation="none")
   Psi=result$uniquenesses
   L=result$loadings
   h2=rowSums(L^2)
   result6=cbind(L,h2,Psi)
   dimnames(result6)=list(NULL,c("F1","F2","Community","Specific Variance"))
   result6
   
   # Rotation with varimax criterion
   result=factanal(factors=2,covmat=R,rotation="varimax")
   Psi=result$uniquenesses
   L=result$loadings
   h2=rowSums(L^2)
   result7=cbind(L,h2,Psi)
   dimnames(result7)=list(NULL,c("F1","F2","Community","Specific Variance"))
   result7
   
   # Rotation with certain degree