*		var		= updated variance
*		vari		= updated inverse of sigma
*
*	Calling Statement:
{parmat, var, vari} = getmulreg(yd,xd,xdtxd,parmat,var,vari,v0i,v0iu0,f0n,g0i);
****************************************************************
*/
PROC (3) = getmulreg(yd,xd,xdtxd,parmat,var,vari,v0i,v0iu0,f0n,g0i);
local vb12, ubn, par, pdim, resid, gni, gn, rp, cp;

	rp		= rows(parmat);
	cp		= cols(parmat);
	pdim	= rp*cp;

	/*
	***********************************************************
	* Generate parmat from N_{rp x cp}(M,v)
	* par = vecr(parmat)
	* par is N(u,V) whee u = vec(M');
	*	V = (Xd'Xd.*.Var^{-1} + V_0^{-1})^{-1}
	*	u = V*( (Xd'.*.Var^{-1})*vec(Yd') + V_0^{-1}u_0 )
	*************************************************************
	*/
 
	vb12	= chol(xdtxd.*.vari + v0i);
	ubn		= ( (xd').*.vari )*vecr(yd) + v0iu0;
	par		= cholsol(ubn + vb12'rndn(pdim,1), vb12);
	parmat	= reshape(par,rp,cp);

	/*
	*********************************************************
	* Generate Var
	* 	Var^{-1} is Wishart df = f0n, scale matrix = gn
	*********************************************************
	*/
	resid 			= yd - xd*parmat;
	gni				= g0i + resid'resid;
	gn  			= invpd(gni);
	{vari, var} 	= wishart(cp,f0n,gn);

retp(parmat,var,vari);
endp;


/*
*******************************************************************
* CNDCOV
* Generates some quantities used in computing conditional normal distribtions.
* INPUT:
*	sigma	= Cov(Y)
*
* OUPUT:
*	smigni	= used in E(Y_i|Y_{not i}).
*	signi	= STD(Y_i |Y_{not i}).
*
* smigni is a mvar x mvar-1 matrix and 
* smigni[i,.] = sigma_{i, not i}*sigma_{not i, not i)^{-1}
* signi is a mvar vector and
* signi[i] = 
* sqrt(sigma_{ii} - sigma_{i,not i}*sigma_{not i, not i}^{-1} sigma_{not i,i})
******************************************************************
*/
PROC (2) = cndcov(sigma);
local smigni, signi, sqmvar, ei, einot, s2noti, sinoti, s2notii, i0, i;

sqmvar	= seqa(1,1,mvar);
smigni	= zeros(mvar,mvar-1);		@ Sigma_{ij}*Sigma_{jj}^{-1}	@
@ smuigj is used in computing mean of i given not i					@
signi	= zeros(mvar,1);			@ sigma_{ii} - sigma_{ij}*Sigma_{jj}^{-1}*sigma_{ji} @
@ signi is STD of Y[i] given not i @
for i0 (1,mvar,1); i = i0;
	ei		= sqmvar .== i;			@ ei[i] = 1, and ei[j] = 0 for j /= i 				@
	einot	= 1 - ei;				@ einot[i] = 0, and einot[j] = 1 for j /= i			@
	@ selif(x,e) selects rows of x with e[j] = 1 @
	s2noti	= selif( selif(sigma,einot)',einot)';	
	@ s2noti = sigma[j,k] with j \= i and k \= i @
	sinoti	= selif( selif(sigma',ei)',einot);
	@ sinoti is the ith column of sigma with the ith row removed @

	s2notii	= invpd(s2noti);
	smigni[i,.]	= sinoti's2notii;
	signi[i]	= sqrt(sigma[i,i] - smigni[i,.]*sinoti);
endfor;
retp(smigni, signi);
endp;		

/*
********************************************************************
* RNDNIGTJ.GSS
*	Generate a vector normal where the ith elememt is greater than
*	the rest.
*	INPUT
*		Y		= current value of  Y form MCMC Iterations
*		IPICK	= which component is greater than the rest, and zero
*		mu		= mean vector
*		smigni	= used in computing mean of y[i] given not y[i].
*		signi	= std of y[i] given not y[i].
		nygen	= number of generations before returning random utility.
*	OUTPUT
*		Y		= mvar vector where Y[i] > max( Y[j]  for j not equal to i )
*********************************************************************************
*/

PROC rndnigtj(y,ipick,mu, smigni, signi, nygen);
local mvar, sqmvar,ei,einot,munoti,ynoti, 
mui, sii, ybot,ytop,yi, i, fori, j, fj ;

mvar 	= rows(mu);					@ Dimension of problem					@
sqmvar	= seqa(1,1,mvar);

for fj (1, nygen, 1); j = fj;		@ Do muliple loops of generating Y					@
for fori (1,mvar,1); i = fori;
	ei		= sqmvar .== i;			@ ei[i] = 1, and ei[j] = 0 for j /= i 				@
	einot	= 1 - ei;				@ einot[i] = 0, and einot[j] = 1 for j /= i			@

	@ selif(x,e) selects rows of x with e[j] = 1 @
	ynoti	= selif(y,einot);		@ mvar - 1 vector that does not have y[i]			@
	munoti	= selif(mu,einot);		@ munoti is mvar-1 vector that does not have mu[i] 	@

	mui		= mu[i] + smigni[i,.]*(ynoti - munoti);
	sii		= signi[i];		@ STD of y[i] given rest @
	if i == ipick;			@ i is the selected choice, yi > max(yj,0)@
		ybot		= maxc(ynoti|0);
		y[i]		= rndtnb(mui,sii,ybot);	
		@ generates a random normal, trucnated below (in plrand.src) @
		@ rndtnb(mu,std,bot): mu = vector or means, std = vector of STD, bot = lower limit @
	else;					@ i was not selected 	@
		if ipick <= mvar;
			ytop		= maxc(y[ipick]);				@ Maximum element @
		else;		@ Base brand mvar+1 selectd, so y < 0 @
			ytop		= 0;
		endif;
		y[i]		= rndtna(mui,sii,ytop);
		@ generates a random normal, truncated above (in plrand.src) @
		@ rndtna(mu,std,top): mu = vector of mean, std = vector of STD, top = upper limit @
	endif;		
endfor;		@ End loop over generating one Y vector 	@
endfor;		@ End loop over generating nygen Y vectors 	@	
retp(y);
endp;




/*
*****************************************************************************************************
*  OUTPUTANAL
*	Does analysis of output and save some results 
****************************************************************************************************
*/
PROC (0) = outputanal;
format 10,5;
local  bout, sout, ebeta, sbeta, cb, rmse, fmts1,fmts2, fmtn1, fmtn2, a,b, flag, i0, i;



let fmtn1[1,3]	= "*.*lf" 10 5;			@ Format for printing numeric variable 		@
let fmtn2[1,3]	= "*.*lf" 10  0;		@ Format for numeric variable, no decimal	@

let fmts1[1,3]	= "-*.*s" 10 9;			@ Format for alpha, left justify		 	@
let fmts2[1,3]  = "*.*s" 10 9;			@ Format for alpha, right justify		  	@
format 10,5;							@ Default pring format						@


output file = ^outfile reset;	@ outfile is the file handle for the output file 	@
								@ Route printed output to the defined by outfile 	@
print "Results from PROBIT1.GSS";
print "Hierarchical Bayes Multivariate Probit";
print;
print "Select one of mvar+1 alternatives.";
print "Alternative mvar+1 is the base alternative.";
print "Observe choice k if Y_{ij}[k] >= max(Y_{ij},0).";
print "Choose base if all the y's < 0.";
print;
print "Latent Utility Model:";
print "Y_{ij} = X_{ij}*beta_i + epsilon_{ij}";
print "Y_{ij} is a mvar vector";
print "epsilon_{ij} is N(0,Sigma).";
print "Identification: sigma[mvar,mvar] = 1";
print;
print "beta_i = Theta'z_i + delta_i";
print "delta_i is N(0, Lambda)";
print;
print "Ouput file: " getpath(outfile);		@ File assigned to file handle outfile @
datestr(date);								@ Print the current data				@
print;
print;
print "-----------------------------------------------------------------------------------";
print;
print "MCMC Analysis";
print;
print "Total number of MCMC iterations           	 = " nmcmc;
print "Number of iterations used in the analysis 	 = " smcmc;
print "Number in transition period               	 = " nblow;
print "Number of iterations between saved iterations = " skip-1;
print;
print "-----------------------------------------------------------------------------------";
print "Number of subjects	                = " nsub;
print "Number of observations per subject = " nchoice;
print;
print "Total number of Y observations     = " ntot;
print "Dimension of Y_{ij}                = " mvar;
print "Number of alternatives             = " mvar+1;
print "Number of dependent variables X    = " rankx " (including intercept)";
print "Number of dependent variables Z    = " rankz " (including intercept)";
print;
print;
print "First level equation for Latent Utilities";
print "Y_{ij} = X_{ij}*beta_i + epsilon_{ij}";
print;
print "		Summary Statistics for X";
call sumstats(xnames,xdata,fmts1,fmts2,fmtn1);
print;
print "Second level equation:";
print "beta_i = Theta'z_i + delta_i";
print;
print "		Summary Statistis for Z:";
call sumstats(znames,zdata,fmts1,fmts2,fmtn1);
print;


if flagtrue == 1;		@ Print some fit for utilities @
	print "-----------------------------------------------------------------------------------";
	print "Fit between true and estimated utilities for each dimension.";
	sout	= {"Variable" "Multi-R" "R-Sqr" "ErrorSTD"};
	call outitle(sout,fmts1,fmts2);
	bout	= ynames~multir~rsquare~stderr;
	call outmat(bout,fmts1,fmtn1);
endif;
print "-----------------------------------------------------------------------------------";
print;
print "Statistics for Individual-Level Regression Coefficients";
if flagtrue == 1;
	ebeta	= meanc(betat);
	sbeta	= stdc(betat);
	print "True Beta";
	sout	= {"Variable" "Mean" "STD"};
	call outitle(sout,fmts1,fmts2);
	bout = xnames~ebeta~sbeta;
	call outmat(bout,fmts1,fmtn1);
endif;

print "HB Estimates of Beta";
sout	= {"Variable" "PostMean" "PostSTD" };
call outitle(sout,fmts1,fmts2);
ebeta	= meanc(betam);
sbeta	= sqrt( meanc( (betas^2)) + stdc( betam)^2);
bout = xnames~ebeta~sbeta;
call outmat(bout,fmts1,fmtn1);
print;
print "-----------------------------------------------------------------------------------";

if flagtrue == 1;
	print "Comparison of True Beta to Individual Level Estimates";
	for i0 (1,rankx,1); i = i0;
		print "Variable " $ xnames[i];
		cb 		= corrx( betat[.,i]~betam[.,i] );
		rmse 	= betat[.,i] - betam[.,i];
		rmse	= rmse'rmse;
		rmse	= sqrt(rmse/nsub);
		print "Correlation between true and HB  = " cb[1,2];
		print "RMSE between true and HB         = " rmse;
		print;
	endfor;
endif;
print "-----------------------------------------------------------------------------------";
print "Estimation of the error covariance Sigam";
sout	= "   "~(ynames');

if flagtrue == 1;
	print "True Sigma";
	call outitle(sout,fmts1,fmts2);
	bout	= ynames~sigmat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif; 
print;
print "Posterior Mean of Sigma";
call outitle(sout,fmts1,fmts2);
bout	= ynames~sigmam;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior STD  of Sigma";
call outitle(sout,fmts1,fmts2);
bout	= ynames~sigmas;
call outmat(bout,fmts1,fmtn1);
print;
print "-----------------------------------------------------------------------------------";
print;
print "HB Estimates of Theta";
sout	= "  "~(xnames');
if flagtrue == 1;
	print "True Theta";
	call outitle(sout,fmts1,fmts2);
	bout = znames~thetat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif;
print "Posterior Mean of Theta";
print outitle(sout,fmts1,fmts2);
bout = znames~thetam;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior STD of Theta";
call outitle(sout,fmts1,fmts2);
bout	= znames~thetas;
call outmat(bout,fmts1,fmtn1);
print;
print "-----------------------------------------------------------------------------------";
print;
sout = "  "~(xnames');
print "HB Estimate of Lambda";
if flagtrue == 1;
	print "True Lambda";
	call outitle(sout,fmts1,fmts2);
	bout = xnames~lambdat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif;
print "Posterior Mean of Lambda";
call outitle(sout,fmts1,fmts2);
bout = xnames~lambdam;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior STD of Lambda";
call outitle(sout,fmts1,fmts2);
bout = xnames~lambdas;
call outmat(bout,fmts1,fmtn1);
print;
print "===================================================================================";



output off;
closeall;
endp;


/*
****************************************************************************
*  TR(X)
*  Compute the trace of X. Functional call (FN). 
***************************************************************************
*/
fn tr(x) = sumc(diag(x));



/*
*****************************************************************************************
* OUTITLE
*	Prints header for columns of numbers.
*	INPUT
*		a 	= character row vector of column names
*		fmts1	= format for first column
*		fmts2	= format for second column
*	OUTPUT
*		None
******************************************************************************************
*/
PROC (0) = outitle(a,fmt1,fmt2);
local mask, fmt, flag, ncols;
ncols	= cols(a);
mask	= zeros(1,ncols);
fmt		= fmt1|(ones(ncols-1,1).*.fmt2);
flag	= printfm(a,mask,fmt);
print;
endp;
/*
***************************************************************************************
* OUTMAT
*	Outputs a matrix:
*	(Character Vector)~(Numeric matrix);
*	The entries in the numeric matrix have the same format
*	INPUT
*		bout		= matrix to be printed
*		fmts		= format for string
*		fmtn		= format for numeric matrix
*	OUTPUT
*		None
******************************************************************************************
*/
PROC (0) = outmat(bout,fmts,fmtn);
local fmt,mask,flag,ncols, nrows;
ncols		= cols(bout);
nrows		= rows(bout);
fmt			= fmts|(ones(ncols-1,1).*.fmtn);
mask		= zeros(nrows,1)~ones(nrows,ncols-1);
flag		= printfm(bout,mask,fmt);
print;
endp;

/*
*****************************************************************************************
* SUMSTATS
*	Prints summary statistics for a data matrix
*	INPUT
*		names		= charater vector of names
*		data		= data matrix
*		fmts1		= format for string
*		fmts2		= format for string
*		fmtn		= format for numbers
*	OUTPUT
*		None
********************************************************************************************
*/
PROC (0) = sumstats(names,data,fmts1,fmts2,fmtn);
local a, bout;
a	= {"Variable" "Mean" "STD" "MIN" "MAX"};
call outitle(a,fmts1,fmts2);
bout	= names~meanc(data)~stdc(data)~minc(data)~maxc(data);
call outmat(bout,fmts1,fmtn);
endp;