new; /* erases everything in memory and closes files and windows 1 */

/*      One-sided inequality constraints and central credible inteverals.
**
**	This is an annotated Gauss batch program for Bayesian 
**      estimation of a normal linear regression,  based on part of
**      Example 24.1.2 in Mittelhammer et al. (2000).
**
**	The purpose of this program is to demonstrate the computation of 
**      central (equal-tail) credible intervals for regression parameters. The
**      likelihood function is based on a normal distribution, and the prior 
**      distribution is in the form of inequality contraints on the beta 
**      parameters of the model. In contrast to BayesIC.e, this program uses
**      ONE-SIDED inequality constraints. The standard weakly informative 
**      prior on sigma2 is used.
**
**	For present purposes, the Bayes analysis is based on outcomes 
**      of a linear regression model of the form,
**
**				y = b1 + b2*x2 + b3*x3 + e
**
**      The Y variable is specified to be distributed as 
**      N(x*Beta, sigma2*eye(n)), i = 1,..., n.  
**
**	The observations for this example are the textile data provided
**      by Theil (Principles of Econometrics, John Wiley and Sons, 1971,
**      p. 102).  The data relate to textile consumption in the 
**      Netherlands. The variables consist of indices measuring
**      the per-capita consumption of textiles, the per-capita level 
**      of real disposable income, and the real price of textiles. All 
**      of the variables in the model are converted, following
**      Theil's example, to natural logarithms.  Thus the model takes the
**      form of a constant elasticity demand function.  The data are
**      found loaded internally to the code of this program.
**
**      In this example it is reasonable to assume that the price
**      elasticity was negative and the income elasticity was positive 
**      but neither was very large in absolute value. 
*/

/*	Set some graphics and program parameters  */
gausset; /* Resets global control variables declared in gauss.dec */ 
/* The next 3 lines set up an output file */
format /rd 11,5; /* right justified numbers of form xxxxxx.xxxxx */
outwidth 110; /* 110 columns in output */
output file = BayesICCI.out reset; /* Overwrite previous version of file */

/* Normal model,  Y ~ Normal(x*Beta, sigma2*eye(n))  */

/* Prior information on Beta in this example is in the form of inequality
** constraints */
lowinc = 0;    /* Set lower bound on income elasticity of demand */
upinc  = 1.5;   /* Set upper bound on income elasticity of demand */
lowprice = -1.5; /* Set lower bound on price elasticity of demand */
upprice  = 0;   /* Set upper bound on price elasticity of demand */

/*	Load part of Theil's data set and transform the data to natural 
**      logarithms. The data represent indices of per-capita textile 
**      consumption, real per-capita disposable income, and the real 
**      price of textiles. We use only the first five observations to 
**      illustrate a case in which the data are not sufficiently 
**      informative to yield good estimates of all parameters of interest
**      unless supplemented by theoretically grounded prior information.
*/

data = { 99.2   96.7  101.0,
         99.0   98.1  100.1,
 	100.0  100.0  100.0,
 	111.6  104.9   90.6,
        122.2  104.9   86.5};
/* 	117.6  109.5   89.7,
** 	121.1  110.8   90.6,
**      136.0  112.3   82.8,
** 	154.2  109.3   70.1,
** 	153.6  105.3   65.4,
** 	158.5  101.7   61.3,
** 	140.6   95.4   62.5,
** 	136.2   96.4   63.6,
** 	168.0   97.6   52.6,
** 	154.3  102.4   59.7,
**	149.0  101.6   59.5,
**	165.5  103.8   61.3 };
*/

/*	Generate the other model components  */
n = rows(data);  /* number of observations */
k = 3;   /* number of regressors, including constant term */
v = n - k;  /* degrees of freedom  */
y = ln(data[.,1]);  /* take natural logs of all entries in first column
**  of data matrix */
x = ones(n,1)~ln(data[.,2 3]);  /* Create nx3 design matrix by concatenation */

/*	Compute the LS estimates  */
{bhat, ehat, yhat} = olsqr2(y, x); /* Compute OLS estimates using QR 
** decompostion */
s2 = ehat'ehat/(n - k);  /* Calculate LS estimate of sigma^2 */
Covbhat = s2*invpd(x'x);  /* LS estimate of cov. matrix via inverse of 
** symmetric pos. def. matrix */

/*	Calculate precision h for use later in the normal - MultT case  */
h = x'x/s2;  

/* Create procedure to draw from a multivariate t distribution */
proc rndmultT(u, h, v, n);
	local t, Zsqr, w, k;
	
	k = rows(u);                  /* count the rows in u */
        /* multiply the transpose of the Cholesky factor of the covariance matrix by 
        ** a matrix of standard normal draws */
	w = chol(invpd(h))'rndn(k,n); 
	Zsqr = 2*rndgam(1,n,v/2);  /* twice a 1xn vector of draws from gamma 
        ** distribution with parameter v/2 */
	t = u + w.*sqrt(v./Zsqr); /* 1xk vector of variables with t-distribution */
	
	retp(t');
endp;

/* Choose a sample size (m) for the Monte Carlo simulation. The sampling process
** will continue until this number of Monte Carlo outcomes have been obtained from the
** Truncated multivariate T-distribution. Choose a multiple of 100. */
m = 17000;
/* Set seed of random number generator for replicability */
rndseed 1;
/*	Initialize variables for the simulation  */
Bmeans = zeros(3,1);
TotKept = 0;
numkept = 0;
bkept = zeros(m,3);
TotReps = 0;
do while numkept < m;
	TotReps = TotReps + 1;

	Btry = rndMultT(bhat, h, v, 1);

	/*	Identify truncated outcomes by zeros  */
	select = (Btry[.,2] .>= lowinc).*(Btry[.,3] .<= upprice);
	select = prodc(select');
	if select == 1;
/*		locate(5,1);  */
		numkept = numkept + 1;
		bkept[numkept,.] = btry;
		probineq = numkept/TotReps;
		 if numkept/int(numkept/5000) == 5000;
			format /mat /rd 5,0;
			print /flush "Monte Carlo Outcomes =" numkept " of a target of " m;
			print;
			format /mat /rd 7,4;
			print /flush "Estimated Posterior Probability of Inequality Constraints =" probineq;
		endif;
	else;
		continue;
	endif;
endo;
/*	Now, calculate the credible regions as well as classical confidence
**	intervals and bounds and print them to compare  */
MLest = bhat;
MLstd = sqrt(diag(covbhat));
tcrit = cdftci(0.05, 14);

format /m1 /rd 8,4;
print "Comparison of ML Confidence Intervals and Confidence Bounds, and Bayes"; 
print "    Credible Regions Under Beta Parameter Inequality Restrictions";
print;
print "     Maximum Likelihood (unconstrained):";
print " Parameter       90% Confidence Intervals";
print;
print "Constant         [" bhat[1]-tcrit*MLstd[1] " , " bhat[1]+tcrit*MLstd[1] "]";
print "Income           [" bhat[2]-tcrit*MLstd[2] " , " bhat[2]+tcrit*MLstd[2] "]";
print "Price            [" bhat[3]-tcrit*MLstd[3] " , " bhat[3]+tcrit*MLstd[3] "]";
print;
print;
format /mat /rd 5,0;
print "       BAYES(MC Repetitions =" m "):";
print "Parameter     90% Credible Region (equal tail probability)";
print;
format /m1 /rd 8,4;
print "Constant         [" quantile(bkept[.,1],.05) " , " quantile(bkept[.,1],.95) "]";
print "Income           [" quantile(bkept[.,2],.05) " , " quantile(bkept[.,2],.95) "]";
print "Price            [" quantile(bkept[.,3],.05) " , " quantile(bkept[.,3],.95) "]";
print;
print "Estimated Posterior Probability of Inequality Constraints =" numkept/TotReps;
print;	
	
/* ndpclex;  clears math coprocessor exception flags (under DOS only) */
/* system;   quit Gauss and return to the operating system */
end; /* Terminate program, revert to command mode, close open files, turn screen on */