#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter4: Section4.2.4
#
# Gibbs sampling for seemingly unrelated regression model
#
# Author : Tomohiro Ando
#
################################################################################


# Generates samples from the multivariate normal distribution

prnorm <- function(mu,A,n,p) {
   U <- svd(A)$u
   V <- svd(A)$v
   D <- diag(sqrt(svd(A)$d))
   B <- U %*% D %*% t(V)
   w <- matrix(0,n,p)
   for (i in 1:n){w[i,] <- mu + B%*%rnorm(p)}
   return(w)
 }

# library

library(bayesm)
library(MCMCpack)

# Data generation

n  <- 100
p  <- 2

X1 <- matrix(runif(n*p,-4,4),ncol=p)
X2 <- matrix(runif(n*p,-4,4),ncol=p)
B1 <- c(3,-2); B2 <- c(2,1)
SS <- 2*matrix(c(0.1,-0.05,-0.05,0.2),2,2)
er <- prnorm(rep(0,len=2),SS,n,2)
y1 <- X1%*%B1+er[,1]
y2 <- X2%*%B2+er[,2]



# Find the generalized least squared estimate

X <- matrix(0,ncol=p*2,nrow=n*2)
X[1:n,1:p] <- X1
X[(n+1):(2*n),(p+1):(2*p)] <- X2
y <- c(y1,y2)

B <- solve(t(X)%*%X)%*%t(X)%*%y
Beta1 <- B[1:p]
Beta2 <- B[-(1:p)]
R <- matrix(0,2,2)
R[1,1] <- sum( (y1-X1%*%Beta1)^2 )
R[1,2] <- R[2,1] <- sum( t(y1-X1%*%Beta1)%*%(y2-X2%*%Beta2) )
R[2,2] <- sum( (y2-X2%*%Beta2)^2 )
S <- R/n
O <- kronecker(solve(S),diag(1,n))

for(i in 1:100){
Beta1 <- B[1:p]
Beta2 <- B[-(1:p)]
B <- solve(t(X)%*%O%*%X)%*%t(X)%*%O%*%y
R <- matrix(0,2,2)
R[1,1] <- sum( (y1-X1%*%Beta1)^2 )
R[1,2] <- R[2,1] <- sum( t(y1-X1%*%Beta1)%*%(y2-X2%*%Beta2) )
R[2,2] <- sum( (y2-X2%*%Beta2)^2 )
S <- R/n
O <- kronecker(solve(S),diag(1,n))
}






# Gibbs sampling

MHn <- 6000 # Numb. of iterations

AB <- matrix(0,MHn,p*2)
AS <- matrix(0,MHn,4)

for(j in 1:MHn){

Bold <- B
Sold <- S
Rold <- R

R[1,1] <- sum( (y1-X1%*%Beta1)^2 ); R[2,2] <- sum( (y2-X2%*%Beta2)^2 )
R[1,2] <- R[2,1] <- sum( t(y1-X1%*%Beta1)%*%(y2-X2%*%Beta2) )
O <- kronecker(solve(S),diag(1,n))
B <- solve(t(X)%*%O%*%X)%*%t(X)%*%O%*%y
Bnew <- prnorm(B,solve(t(X)%*%O%*%X),1,4)
Beta1 <- Bnew[1:p]; Beta2 <- Bnew[-(1:p)]

AB[j,] <- Bnew

R[1,1] <- sum( (y1-X1%*%Beta1)^2 ); R[2,2] <- sum( (y2-X2%*%Beta2)^2 )
R[1,2] <- R[2,1] <- sum( t(y1-X1%*%Beta1)%*%(y2-X2%*%Beta2) )
S <- R/n

Snew <- rwishart(n,solve(S)/n)$IW

AS[j,] <- as.vector(Snew)

}



# Try also the same procedures below with respect to the remaining parameters

# Geweke's test

geweke.diag(mcmc(AB[,1]))$z

# Plot Sample pass of beta11

plot(AB[,1],xlab="",ylab="",type="l")

# Plot autocorrelation function of beta11

autocorr.plot(mcmc(AB[,1]),lag.max=50)

# Plot posterior density of beta11

densplot(mcmc(AB[,1]))