#install.packages("MCMCpack")
library(MCMCpack)

data(swiss)

?swiss

dim(swiss)

pairs(swiss)

# MLE factor analysis
mle = factanal(~Agriculture+Examination+Education+Catholic
                     +Infant.Mortality, factors=2,data=swiss)

mle

# Bayes factor analysis                    
bayes = MCMCfactanal(~Agriculture+Examination+Education+Catholic
                     +Infant.Mortality, factors=2,
                      lambda.constraints=list(Examination=list(1,"+"),
                      Examination=list(2,"-"), Education=c(2,0),
                      Infant.Mortality=c(1,0)),
                      verbose=0,store.scores=FALSE,a0=1,b0=0.15,
                      data=swiss,burnin=5000,mcmc=50000,thin=20)
summary(bayes)
   
par(mfrow=c(3,3))
hist(bayes[,1],xlab="",prob=TRUE,main="B11")
abline(v=mle$load[1,1],col=2)
hist(bayes[,2],xlab="",prob=TRUE,main="B12")
abline(v=mle$load[1,2],col=2)
hist(bayes[,9],xlab="",prob=TRUE,main="Psi1")
abline(v=mle$uniq[1],col=2)

hist(bayes[,3],xlab="",prob=TRUE,main="B21")
abline(v=mle$load[2,1],col=2)
hist(bayes[,4],xlab="",prob=TRUE,main="B22")
abline(v=mle$load[2,2],col=2)
hist(bayes[,10],xlab="",prob=TRUE,main="Psi2")
abline(v=mle$uniq[2],col=2)

hist(bayes[,5],xlab="",prob=TRUE,main="B31")
abline(v=mle$load[3,1],col=2)
plot(0,col=0,axes=FALSE,xlab="",ylab="",main="B32")
hist(bayes[,11],xlab="",prob=TRUE,main="Psi3")
abline(v=mle$uniq[3],col=2)

hist(bayes[,6],xlab="",prob=TRUE,main="B41")
abline(v=mle$load[4,1],col=2)
hist(bayes[,7],xlab="",prob=TRUE,main="B42")
abline(v=mle$load[4,2],col=2)
hist(bayes[,12],xlab="",prob=TRUE,main="Psi4")
abline(v=mle$uniq[4],col=2)

plot(0,col=0,axes=FALSE,xlab="",ylab="",main="B51")
hist(bayes[,8],xlab="",prob=TRUE,main="B52")
abline(v=mle$load[5,2],col=2)   
hist(bayes[,13],xlab="",prob=TRUE,main="Psi5")
abline(v=mle$uniq[5],col=2) 


##########################################
# BAYESIAN FACTOR STOCHASTIC VOLATILITY  #
##########################################
#install.packages("factorstochvol")
library("factorstochvol")

data(exrates, package = "stochvol")

exrates$date  =  NULL

?exrates

dim(exrates)

# Compute the de-meaned percentage log returns:
dat  =  100 * logret(exrates, demean = TRUE)

par(mfrow=c(1,1))
corr = vech(cor(dat))
corr = corr[corr<1]
hist(corr,xlim=c(-1,1),prob=TRUE,xlab="",main="All correlations")

res  =  fsvsample(dat, factors = 1, draws = 2000, burnin = 1000, runningstore = 6)

par(mfrow=c(1,1))
boxplot(t(res$facload[,1,]),names=colnames(dat),cex.axis=0.75)


par(mfrow=c(1,1))
voltimeplot(res)

par(mfrow=c(1,1))
corimageplot(res, nrow(dat), plotCI = 'circle')