yx <-read.table("c:/jorge/business/bayes/trans3_1_4g1.txt")
y<-as.matrix(yx[,1])
x<-as.matrix(yx[,2:16])
nobs<-nrow(x) 								#(1)number of firms
nvar<-ncol(x) 								#(2)number of explanatory variables
sig2<-0.05 								#i(3)prior value for the variance
rstar<-0.85	 							#(4)prior for the median efficiency level
v<-matrix(data=0.05, nrow=nobs, ncol=1)  					#(5)a vector of prior values for the inefficiency parameter
ones<-matrix(data=1, nrow=nobs, ncol=1) 					#(6)a vector containing 1
samp<-matrix(data=NA, nrow=10000, ncol=nobs+nvar+2) 				#(7)a big matrix where each iteration outcome will be dumped

#(1)STARTING BIG SAMPLING LOOP

for(i in 1:11000)
{

	ystar<- y - v 							#(8)define ystar as the initial y minus the prior inefficiency level (0.05)
	betahat<-as.matrix( solve(t(x) %*% x) %*% t(x) %*% ystar)  		#(9)estimate of initial regression coefficients using ystar
	covb<-sig2*solve(t(x)%*%x) 						#(10)estimate of covariance matrix using prior variance (0.05)
	beta <- as.matrix(mvrnorm(n=1, betahat, covb, tol=1e-6, empirical=FALSE)) 	#(11)sampling regression coefficients from a multivariate 
									#  normal(betahat, covb)

	#(2)SAMPLE THE PRECISION ESTIMATE FROM AN INVERSE CHI-SQUARE
	# CONDITIONAL ON BETA AND 
	#ESTIMATE SIG2 AS THE INVERSE OF THE PRECISION

	e<-ystar-x%*%beta 						#(12)estimating the vector of errors (residuals)
	sse<-t(e)%*%e 							#(13)estimating the sum of square errors (residuals)
	h<-rgamma(1,shape=(nobs-2)/2,scale=sse/2) 				#(14)sampling the precision estimate (h) from an inverse chi-square 
									# using the gamma distribution
	sig2<-1/h 							#(15)sigma square is defined as 1/precision

	#(3) SAMPLE LAMBDA, THE MEAN OF THE INFFICIENCY
	#ERROR FROM ITS CONDITIONAL

	df1<- nobs+1 							#(16)the shape parameter of lambda inverse 
	df2<-1/(t(v)%*%ones-log(rstar)) 					#(17)the scale parameter of lambda inverse
	lambdainv<-rgamma(1,shape=df1,scale=df2) 				#(18)sampling lambda inverse from an inverse chi-square 
									#using the gamma distribution lambdainv~Gamma(df1,df2) and the 
									#mean is defined as df1*df2
	lambda<-1/lambdainv	   					#(19)obtaining lambda, the mean of the inefficiency error

	#(4)FIND THE CONDITIONAL MEAN OF THE TRUNCATED
	#NORMAL, I.E., SAMPLE V FROM ITS CONDITIONAL
	#SAMPLE FROM THE TRUNCATED NORMAL

	meanv<-y-x%*%beta-(sig2/lambda)*ones 				#(20)generating a vector of inefficiency errors
	varv<-sig2*ones							#(21)creating a vector containing the variance estimate 
	b<-meanv/sqrt(varv)						#(22) endpoint of half-line interval 

	p<-matrix(data=NA, ncol=1, nrow=nobs)				#(23) The for... loop generates a random number
	for(j in 1:nobs)							#p from a truncated to the left normal distribution
	{								#this was taken from the function x = nmlt_rnd(b)
		c1<- b[j]							#written for MATLAB by James P. LeSage, Dept of Economics
		t4<-0.450							#University of Toledo
		if(c1>t4)							#2801 W. Bancroft St,
		{							#Toledo, OH 43606
			z<- -log(runif(1)/c1)				# jpl@jpl.econ.utoledo.edu
			while(runif(1)>exp(-0.5*z*z))
			{
				z<- -log(runif(1))/c1
			}
			p[j]<-c1+z
		}
		else
		{
			p[j]<- runif(1)
			while(p[j]<c1)
			{
				p[j]<- runif(1)
			}
		}
	}
	v<-meanv+p*sqrt(varv)						#(24)once p has been drawn, it can be used to estimate
									#inefficiency error
	effic<- exp(-v)							#(25) estimating the level of efficiency for firm i

	#DROP INITIAL 1000 POINTS AND PUT REST IN SAMP
	if(i>1000)
	{
		samp[i-1000, ]<-c( t(effic),t(beta),sig2,lambda)			#(26)dumping results in matrix samp for further analysis
	}								#after dropping 1000 burn-in iterations
}