##########################################
#
# Skewed regression model
#
##########################################

# Data
y = scan("https://web.ics.purdue.edu/~jltobias/wages.raw")
n = length(y)
hist(y,breaks=seq(0,max(y),length=30))

# Prior hyperparameters
iB0 = diag(1/10,2)
e0 = 3
f0 = 2

# MCMC set-up
delta  = 0
beta   = 0
sigma2 = 1
M0     = 2000
M      = 3000
niter = M0+M
draws = matrix(0,niter,3)

# MCMC algorithm
for (iter in 1:niter){
	
	 # Sampling the z's
  # z = rep(0,n)
  # for (i in 1:n){
    # m = delta*(y[i]-beta)/(sigma2+delta^2)
    # v = sigma2/(sigma2+delta^2)
    # z[i] = rtruncnorm(1,a=0,b=Inf,mean=m,sd=sqrt(v))
  # }
  m = delta*(y-beta)/(sigma2+delta^2)
  v = sigma2/(sigma2+delta^2)
  z = rtruncnorm(n,a=0,b=Inf,mean=m,sd=sqrt(v))

  # Sampling beta and delta
  X = cbind(1,z)
  B = solve(t(X)%*%X/sigma2+iB0)
  b = B%*%t(X)%*%y/sigma2
  draw = b + t(chol(B))%*%rnorm(2)
  beta = draw[1]
  delta = draw[2]

  # Sampling sigma2
  e1 = e0 + n/2
  f1 = f0 + sum((y-beta-delta*z)^2)/2
  sigma2 = 1/rgamma(1,e1,f1)

  # Storing the results
  draws[iter,] = c(beta,delta,sigma2)
}
draws = draws[(M0+1):niter,]

# Posterior summaries
par(mfrow=c(3,3))
ts.plot(draws[,1],xlab="iteration",ylab="",main=expression(beta))
ts.plot(draws[,2],xlab="iteration",ylab="",main=expression(beta))
ts.plot(draws[,3],xlab="iteration",ylab="",main=expression(beta))
acf(draws[,1],main="")
acf(draws[,2],main="")
acf(draws[,3],main="")
hist(draws[,1],prob=TRUE,main="",xlab="")
hist(draws[,2],prob=TRUE,main="",xlab="")
hist(draws[,3],prob=TRUE,main="",xlab="")