data = read.table("http://hedibert.org/wp-content/uploads/2020/03/sp500-djia-nasdaq.txt",header=TRUE)

par(mfrow=c(2,3))
plot(data[,2],xlab="day",ylab="price",type="l",main="sp500")
plot(data[,3],xlab="day",ylab="price",type="l",main="djia")
plot(data[,4],xlab="day",ylab="price",type="l",main="nasdaq")

plot(diff(log(data[,2])),xlab="day",ylab="log-return",type="l",main="sp500")
plot(diff(log(data[,3])),xlab="day",ylab="log-return",type="l",main="djia")
plot(diff(log(data[,4])),xlab="day",ylab="log-return",type="l",main="nasdaq")


y = as.matrix(data[,2:4])
n = nrow(y)
y1 = y - matrix(apply(y,2,mean),n,3,byrow=TRUE)
y1 = y1%*%diag(sqrt(1/diag(var(y1))))

par(mfrow=c(1,2))
ts.plot(y,col=1:3,ylab="Price")
ts.plot(y1,col=1:3,ylab="Standardized price")
legend("topleft",legend=c("sp500","djia","nasdaq"),col=1:3,lty=1)


r = as.matrix(100*(data[2:n,2:4]/data[1:(n-1),2:4]-1))
n1 = nrow(r)

par(mfrow=c(2,3))
ts.plot(r[,1],ylab="Returns",main="sp500",ylim=range(r))
ts.plot(r[,2],ylab="Returns",main="djia",ylim=range(r))
ts.plot(r[,3],ylab="Returns",main="nasdaq",ylim=range(r))
plot((n1-100):n1,r[(n1-100):n1,1],ylab="Returns",main="sp500",ylim=range(r),type="l")
plot((n1-100):n1,r[(n1-100):n1,2],ylab="Returns",main="djia",ylim=range(r),type="l")
plot((n1-100):n1,r[(n1-100):n1,3],ylab="Returns",main="nasdaq",ylim=range(r),type="l")


r1 = r - matrix(apply(r,2,mean),n1,3,byrow=TRUE)
r1 = r1%*%diag(sqrt(1/diag(var(r1))))

par(mfrow=c(1,3))
ts.plot(r1[,1],ylab="Standardized log-returns",main="sp500",ylim=range(r1))
ts.plot(r1[,2],ylab="Standardized log-returns",main="djia",ylim=range(r1))
ts.plot(r1[,3],ylab="Standardized log-returns",main="nasdaq",ylim=range(r1))


# Factor analysis
fac = factanal(r1,factors=1,scores="regression")

par(mfrow=c(1,1))
ts.plot(fac$scores,ylab="Common factor")

cor(fac$scores,r1)

vd = round(cbind(fac$load^2,fac$uniq),3)
colnames(vd) = c("common factor","unique factor")
rownames(vd) = c("sp500","djia","nasdaq")
vd

install.packages("factorstochvol")
library(factorstochvol)

fsv = fsvsample(r1)

par(mfrow=c(1,3))
hist(fsv$facload[1,1,])
hist(fsv$facload[2,1,])
hist(fsv$facload[3,1,])

par(mfrow=c(1,1))
ts.plot(sqrt(fsv$runningstore$cov[,,1]),ylab="Standard deviations",col=1:3)
legend("topleft",legend=c("sp500","djia","nasdaq"),col=1:3,lty=1)

par(mfrow=c(1,1))
ts.plot(fsv$runningstore$cor[,,1],col=1:3,ylab="Correlations")
legend("bottomleft",legend=c("sp500,djia","sp500,nasdaq","djia,nasdaq"),col=1:3,lty=1)
abline(h=1,lty=2)

par(mfrow=c(1,2))
ts.plot(sqrt(fsv$runningstore$cov[(n1-200):n1,,1]),ylab="Standard deviations",col=1:3)
legend("topleft",legend=c("sp500","djia","nasdaq"),col=1:3,lty=1)
ts.plot(fsv$runningstore$cor[(n1-200):n1,,1],col=1:3,ylab="Correlations")
legend("bottomright",legend=c("sp500,djia","sp500,nasdaq","djia,nasdaq"),col=1:3,lty=1)
abline(h=1,lty=2)

par(mfrow=c(1,1))
plot(sqrt(fsv$runningstore$cov[(n1-200):n1,1,1]),fsv$runningstore$cor[(n1-200):n1,1,1],
xlab="Standard deviations (SP500)",ylab="Correlations",pch=16)
points(sqrt(fsv$runningstore$cov[(n1-200):n1,1,1]),fsv$runningstore$cor[(n1-200):n1,2,1],col=2,pch=16)
legend("bottomright",legend=c("sp500,djia","sp500,nasdaq"),col=1:2,pch=16)