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

/*      One-sided inequality constraints.
**
**	This is an annotated Gauss batch program for Bayesian 
**      estimation of a normal linear regression,  based on part of
**      Example 23.3.1 in Mittelhammer et al. (2000).
**
**	The purpose of this program is to demonstrate the use of Monte
**      Carlo methods for calculating properties of the posterior 
**      distribution of the parameters of a linear model. 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 = BayesIC2.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  = 2;   Set upper bound on income elasticity of demand */
/* lowprice = -2; 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;  

/* Set the number of Monte Carlo trials to a postivie interger multiple of 100 */
m =  100000; 

/*	Set print format for matrices.
**  	/m1 means print 1 carriage return/line feed pair before each row of a
**      matrix. /rd 16,4 means numbers should be right-justified with field 
**      width 16 and 4.0 digits right of the decimal place.
*/
format /m1 /rd 16,4;

/* 	We now generate m independent outcomes from the appropriate multivariate
**	T distribution that characterizes the UNTRUNCATED posterior distribution 
**      of the model parameters.  We then truncate this posterior according to the
**      inequality-based prior information on the parameters.  The truncation is 
**      accomplished by discarding all of the outcomes of the multivariate 
**      T-distribution that do not satisfy the inequality constraints.
**      We also form the ratio of the number of retained outcomes divided by the total
**      number of MC outcomes, which is an estimate of the posterior probability that 
**      the inequality constraints are valid based on the data information.
**
**	As a numerical measure of the accuracy with which the posterior means are
**      calculated, we also report the value of the sample standard deviation of the
**      sample mean for each of the three posterior means (intercept, income parameter,
**      price parameter) calculated.
**
**	Initialize some variables used in the simulation loop
*/
Bmeans = zeros(3,1);
Bsqrs = zeros(3,1);
TotKept = 0;

/* Create procedure to draw from a multivariate t distribution, using
** the algorithm given by Mittelhammer, pp. 684-85.  The first argument
** of the procedure (u) is the location (mean) vector.  The second
** argument (h) is the precision matrix.  The third argument (v) is the
** degrees of freedom. The fourth argument (n) is the dimension of the vector
** returned (the number of draws).
*/
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.  Each column of w will have a 
        ** multivariate normal distribution with mean 0 and precision h. */
	w = chol(invpd(h))'rndn(k,n); 
	Zsqr = 2*rndgam(1,n,v/2);  /* A 1xn vector of draws from gamma 
        ** distribution with scale parameter 2 and shape parameter v/2.
        ** These are equivalent to draws from a Chi-squared
        ** distribution with v degrees of freedom. */
	t = u + w.*sqrt(v./Zsqr); /* Use Gauss's element-by-element operators to 
        ** create a kxn matrix.  Each column is drawn from  multivariate
        ** t-distribution with location u, precision h, and degrees of freedom v */
	retp(t');  /* Return a nxk matrix */
endp;

/*	Start the simulation loop  */
rndseed 1;  /* set seed of random number generator for replicability */
numloop = m/100;  
/* set format for matrix to right-justified fields of width 6 and no decimal */
format /mat /rd 6,0; 
/* 	Begin a loop with counter variable j with initial value 1, final value numloop, 
** 	and increments of size 1.
*/
for j (1, numloop, 1);
/*	locate(4,1);
**	print /flush  "	Trial" j*100 "of" m;  */ 
/* Draw from a multivariate t distribution with mean bhat, precision h, degrees of 
** freedom v, and dimension 100)
*/


	Btry = rndmultT(bhat, h, v, 100); 
	/* Identify truncated outcomes by zeros. select will be a 100x1 vector of 
        ** zeros and ones.  */
	select = (Btry[.,2] .>= lowinc).*(Btry[.,3] .<= upprice);
	select = prodc(select'); /* Now select is a 1x100 vector of zeros and ones. */
	numkept = sumc(select);

	/*	Accumulate information on the number of draws kept and the
	**	moments of the accepted draws  */
	TotKept = TotKept + numkept;
/* 	Select a row from matrix btry if the corresponding element of select is a 1. 
**  	Then sum the elements in each column of the matrix formed from the selected
**      rows.  Also sum the corresponding squares.  
*/
	bmeans = bmeans + sumc(selif(btry,select));
	bsqrs = bsqrs + sumc(selif(btry.*btry,select));
endfor;
cls;

/*	Calculate the posterior results and the accuracy measures  */
bmeans = bmeans/TotKept;  /* Mean of accepted draws of betas */
bsqrs = bsqrs/TotKept;    /* Mean of squares of accepted draws */
fracacc = TotKept/m;     /* Fraction of draws accepted */

PostVars = bsqrs-bmeans^2;  /*posterior var(beta) = post E(beta^2)-(post E(beta))^2 */
StdMean = sqrt(PostVars/TotKept); /* Standard deviation of posterior mean */

MLest = bhat;
MLstd = sqrt(diag(covbhat));

lowbound = bmeans - 2*StdMean;
upbound = bmeans + 2*StdMean;

/*	Report the results */
print;
format /mat /rd 7,4;
print "Comparison of ML and Bayes Estimates Under Beta Parameter Inequality Restrictions";
print;
print "Maximum Likelihood (unconstrained):";
print "Parameter       Estimate         Std Error";
print;
print "Constant         " MLest[1,1] "          " MLstd[1,1];
print "Income           " MLest[2,1] "          " MLstd[2,1];
print "Price            " MLest[3,1] "          " MLstd[3,1];
print;
format /mat /rd 6,0;
print "BAYES:  (Attempted MC Repetitions =" m ")";
print "Parameter          Estimate    Posterior StdDev      StdDev of Mean";
print;
format /mat /rd 10,4;
print "Constant         " Bmeans[1,1] "         " sqrt(PostVars[1,1]) "         " Stdmean[1,1];
print "Income           " Bmeans[2,1] "         " sqrt(PostVars[2,1]) "         " Stdmean[2,1];
print "Price            " Bmeans[3,1] "         " sqrt(PostVars[3,1]) "         " Stdmean[3,1];
print;
format /mat /rd 7,4;
print "Fraction of draws accepted = " fracacc;


/* 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 */