# 6 time series
# 22 obs per day, 991 days = 21802 observations

y = read.table("https://hedibert.org/wp-content/uploads/2024/10/bovespa.txt",header=TRUE)
y = y[,c(1,2,3,5,6)]
p = ncol(y)
n = nrow(y)
names = names(y)

names = c("vale","itau","bbrasil","cyrela","redecard")
range = range(y)
par(mfrow=c(2,3))
for (i in 1:p)
  ts.plot(y[,i],main=names[i],ylim=range,ylab="Return")


# PCA = S = u%*%d%*%v
# In my notes W=u and lambda=d
S       = cov(y)
pca0    = svd(S)
k       = 2
W0      = pca0$u[,1:k]
lambda0 = pca0$d[1:k]
x0      = t(t(W)%*%t(y))

# components x are orthogonal with variance lambda
round(t(x0)%*%x0/n,3)
round(diag(lambda0),3)

# optimal linear reconstruction
y1     = t(W0%*%t(x0))
par(mfrow=c(2,3))
for (i in 1:p){
  plot(y[,i],y1[,i],xlab=names[i],ylab="Linear reconstruction")
  abline(0,1,col=2,lwd=2)
}

# Probabilistic PCA
Lambda = diag(pca0$d[1:k])
sig2   = mean(pca0$d[(k+1):p])
W      = W0%*%sqrt(Lambda-sig2*diag(1,k))
M      = t(W)%*%W + sig2*diag(1,k)
iM     = solve(M)
Ex1    = t(iM%*%t(W)%*%t(y))
y2     = t(W%*%t(Ex1))                          # Not orthogonal reconstruction
y3     = t(W%*%solve(t(W)%*%W)%*%t(M)%*%t(Ex1)) # Orthogonal reconstruction

par(mfrow=c(3,5))
for (i in 1:p){
  plot(y[,i],y1[,i],xlab=names[i],ylab="Linear")
  abline(0,1,col=2,lwd=2)
}
for (i in 1:p){
  plot(y[,i],y2[,i],xlab=names[i],ylab="Not orthogonal")
  abline(0,1,col=2,lwd=2)
}
for (i in 1:p){
  plot(y[,i],y3[,i],xlab=names[i],ylab="Orthogonal")
  abline(0,1,col=2,lwd=2)
}

# Factor analysis
k  = 2
fa = factanal(y,factors=k,scores="regression")

fa

beta  = fa$loadings
Sig   = diag(fa$uniqueness)
iSig  = solve(Sig)
Omega = beta%*%t(beta)+Sig
Vf    = solve(diag(1,k)+t(beta)%*%iSig%*%beta)
Ef    = t(Vf%*%t(beta)%*%iSig%*%t(y))
round(Ef-fa$scores,4)

par(mfrow=c(2,2))
plot(x0[,1],Ex1[,1])
plot(x0[,2],Ex1[,2])
plot(Ef[,1],-Ex1[,1])
plot(Ef[,2],Ex1[,2])