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]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 }