# Factor Analysis for Stock data via PCA
#
e<-eigen(var(normdf(stocks)))

# 1 faktor modell
t(t(e$vectors)*sqrt(e$values))[,1]

# LL'
(t(t(e$vectors)*sqrt(e$values))) %*% (t(e$vectors)*sqrt(e$values))
cor(stocks)

# psi_i
1-t(t(e$vectors)*sqrt(e$values))[,1]**2

# 2 faktor modell
t(t(e$vectors)*sqrt(e$values))[,1:2]

# psi_i
1 - apply(t(t(e$vectors)*sqrt(e$values))[,1:2]**2, 1, sum)

# Residual Matrix
round(
cor(stocks) - 
(t(t(e$vectors[,1:2])*sqrt(e$values[1:2]))) %*% (t(e$vectors[,1:2])*sqrt(e$values[1:2])) -
diag(1 - apply(t(t(e$vectors)*sqrt(e$values))[,1:2]**2, 1, sum))
, 3)


# Factor Analysis for Stock data via MLE
#
sf <- factanal(stocks, 2, rotation="none")

# psi_i
1-apply(sf$loadings[,1:2]**2, 1, sum)

# Residual Matrix
round(
cor(stocks) - 
(sf$loadings[,1:2]) %*% t(sf$loadings[,1:2]) -
diag(1 - apply(sf$loadings[,1:2]**2, 1, sum))
, 3)

# Factor Analysis for Stock data via MLE with varimax Rotation
#
sfr <- factanal(stocks, 2)

#plots
plot(sfr$loadings[,1:2], xlim=c(0,1), ylim=c(-0.5,1), pch="R")
points(sf$loadings[,1:2], pch="O")
abline(h=0)

#test
det((sf$loadings[,1:2]) %*% t(sf$loadings[,1:2]) + diag(1-apply(sf$loadings[,1:2]**2, 1, sum)))/
det(cor(stocks))