#
# Daily Closing Prices of Major European Stock Indices, 1991–1998
# Contains the daily closing prices of major European stock indices: 
#
# Germany DAX (Ibis)
# Switzerland SMI
# France CAC
# UK FTSE. 

library("astsa")
library("stochvol")
library("factorstochvol")

data("EuStockMarkets")
price = EuStockMarkets
ret   = scale(apply(log(price),2,diff))
n     = nrow(ret)
names = c("Germany DAX","Swiss SMI","France CAC","UK FTSE")
names1 = c("DAX","SMI","CAC","FTSE")

par(mfrow=c(1,1))
ts.plot(price,col=1:4,ylab="Daily Closing Prices")
legend("topleft",legend=names,col=1:4,bty="n",lty=1)

par(mfrow=c(1,1))
ts.plot(ret,col=1:4,ylab="Log returns (standardized)")

cor(ret)

# Fitting univariate stochastic volatility models
qstd = array(0,c(4,n,3))
for (i in 1:4){
  fit = svtlsample(ret[,i])
  qstd[i,,] = t(apply(exp(fit$latent[[1]]/2),2,quantile,c(0.05,0.5,0.95)))
}

par(mfrow=c(2,2))
for (i in 1:4){
  ts.plot(qstd[i,,],col=c(3,2,3),main=names[i],ylim=c(0,max(qstd[i,,])))
  lines(0.1*abs(ret[,i]),type="h")
}

par(mfrow=c(1,1))
ts.plot(t(qstd[,,2]),col=1:4)
legend("topleft",legend=names,col=1:4,bty="n",lty=1)




# Standardized returns 
# --------------------
ret1 = ret/t(qstd[,,2])

par(mfrow=c(2,2))
for (i in 1:4)
  ts.plot(ret1[,i],main=names[i],ylab="")

cor(ret1)


# Principal components analysis
# -----------------------------
pca = princomp(ret1)

summary(pca)

pairs(pca$scores)

par(mfrow=c(2,2))
for (i in 1:4)
  plot(pca$scores[,1],ret1[,i],xlab="1st principal component",ylab=names[i])


# Factor analysis
# ---------------
fa = factanal(ret1,factors=1,scores="regression")

summary(fa)

ts.plot(fa$scores[,1],ylab="Common factor")

par(mfrow=c(1,1))
plot(fa$scores[,1],pca$scores[,1],xlab="1st principal component",ylab="Common factor")

cbind(fa$load,fa$uniq)


# Factor stochastic volatility
fsv = fsvsample(ret)

par(mfrow=c(1,1))
ts.plot(exp(fsv$runningstore$logvar[,5,1]/2),ylim=c(0,3),col=2,ylab="Standard deviation")
title("Common factor")
lines(0.2*abs(fsv$runningstore$fac[,,1]),type="h")

par(mfrow=c(2,2))
for (i in 1:4){
  hist(fsv$facload[i,1,],main=names[i],prob=TRUE,xlab="")
  abline(v=fa$load[i],lwd=2)
}

par(mfrow=c(1,1))
ts.plot(fsv$runningstore$sd[,,1],col=1:5)
legend("top",legend=names,col=1:4,bty="n",lty=1)

par(mfrow=c(3,2))
ts.plot(fsv$runningstore$cor[,1,1],ylim=c(0,1),ylab="")
abline(h=cor(ret)[2,1],col=2)
title(paste("corr(",names1[1],",",names1[2],")",sep=""))
ts.plot(fsv$runningstore$cor[,2,1],ylim=c(0,1),ylab="")
abline(h=cor(ret)[3,1],col=2)
title(paste("corr(",names1[1],",",names1[3],")",sep=""))
ts.plot(fsv$runningstore$cor[,3,1],ylim=c(0,1),ylab="")
abline(h=cor(ret)[4,1],col=2)
title(paste("corr(",names1[1],",",names1[4],")",sep=""))
ts.plot(fsv$runningstore$cor[,4,1],ylim=c(0,1),ylab="")
abline(h=cor(ret)[3,2],col=2)
title(paste("corr(",names1[3],",",names1[2],")",sep=""))
ts.plot(fsv$runningstore$cor[,5,1],ylim=c(0,1),ylab="")
abline(h=cor(ret)[4,2],col=2)
title(paste("corr(",names1[4],",",names1[2],")",sep=""))
ts.plot(fsv$runningstore$cor[,6,1],ylim=c(0,1),ylab="")
abline(h=cor(ret)[4,3],col=2)
title(paste("corr(",names1[4],",",names1[3],")",sep=""))