/*
************************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*	PROBIT1.GSS
*	HB Probit Regression Model.
*----->PROBIT1 uses common choice set matrix for all subjects.
*----->PROBIT2 allows for different design matrices for subjects.
*
*       Use McCullough and Rossi's method of post-MCMC normalizing to identify model
*
*		Select one of mvar+1 alternatives.
*
*		Y_{ij} 	= X_{j}*beta_i + epsilon_{ij} 
*			for i = 1, ..., I and j = 1, ..., n_i
*		Y_{ij} is mvar vector
*		Y_{ijk} is the utility from subject i, choice set j, and alternative k
*		for 	i = 1, ..., nsub
*				j = 1, ..., nchoice
*				k = 1, ..., mvar
*
*		Alternative mvar+1 is the base vector.
*		
*		
*		Select alternative k if:
*			Y_{ijk} > max( Y_{ijl} } for l \= k < mvar+1.
*		Select base brand if max(Y) < 0.
*
*		Observe the choices, not the utilities Y_{ij}.
*		Pick_{ij} is a mvar vector of 0/1.  
*
*
*		beta_i is rankx vector
*		epsilon_{ij} is N(0,Sigma)
*       Error variance for brand mvar is one.
*
*
*		X_{j} is mvar x rankx for choice j.
*
*		beta_i	= Theta'Z_i + delta_i
*		delta_i	is N(0,Lambda)
*			Z1 is ln(income) and z2 = family size.
*	PRIORS
*		Sigma  is Inverted Wishart(sf0,sg0)
*		Theta is maxtrix normal(u0,v0).
*		That is, vec(Theta) is N(vec(u0),v0).
*		vec(theta) stacks the columns of theta.
*		Lambda is Inverted Wishart(f0, g0) 
***********************************************************************
*/
new;
outfile		= "results1.dat";	@ Specify output file for saving results 					@
								@ outfile is a string variable that contains a file name 	@
inx	 		= "xdata";			@ Name of Gauss file with choice set design matrix			@

inz	 		= "zdata";			@ Name of Gauss file with Z data							@
flagtrue	= 1;				@ 1 -> knows true parameters from simulation				@

/*
********************************************************************
*	Initialize parameters for MCMC 
********************************************************************
*/
smcmc		= 100;		@ number of iterations to save for analysis 				@
skip		= 1;			@ Save every skip iterations							@
nblow		= 100;		@ Initial transition iterations 							@
nmcmc		= nblow + skip*smcmc;	@ total number of iterations					@
nygen		= 1;		@ Do nygen generations of Y for each MCMC iteration.		@

/*
********************************************************************
*	Get data
********************************************************************
*/

@ Get dimensions and pointers @
load iptx = iptx;
mvar = iptx[1,2];						@ # choices = mvar + 1 @
load cdata = cdata;



@ Input Gauss files @
open f1 	= ^inx;					@ Get Gauss file for X data 				@
xdata		= readr(f1,rowsf(f1));		@ "p" for picks								@	
ci			= close(f1);



nsub		= rows(cdata);				@ Number of subjects @
nchoice		= cols(cdata);				@ Number of choices @
ntot			= nsub*nchoice;
xnames		= setvars(inx);			      @ Get the variable names that accompnay X, Y data @
ynames		= xnames[1:mvar];			@ Use names of interceps for names of components of Y @
ynames2		= ynames|" Base ";
rankx		= cols(xdata);

open f1		= ^inz;
zdata		= readr(f1,rowsf(f1));		@ First column of zdata is a vector of ones 	@
ci			= close(f1);
znames		= setvars(inz);
rankz		= cols(zdata);					@ # of Z variables (includes intercept) @
thdim		= rankx*rankz;					@ dimension of vec(theta) @

@ Compute some sufficient statistics @
ztz		= zdata'zdata;

/*
********************************************************************
*	Initialize Priors 
********************************************************************
*/


@ Prior for theta is N(u0,v0) @

u0 		= zeros(thdim,1);
v0   	= 100*eye(thdim);		@ thdim = rankx*rankz    @
v0i  	= invpd(v0);    		@ used in updating theta @
v0iu0 	= v0i*u0;          		@ used in updating theta @

@ Prior for sigma is IW(sf0, gs0) @
sf0 	= mvar+2; sfn = sf0 + ntot;
sg0i	= eye(mvar);



@ Lambda^{-1} is W_rankx(f0,g0 ) @
@ f0 = prior df, g0 = prior scale matrix @
f0 		= rankx+2;  f0n = f0 + nsub;
g0i 	= eye(rankx);   @ g0^{-1} @

/*
*******************************************************************
*	Initialize MCMC
******************************************************************
*/
ydata	= zeros(ntot,mvar);

beta 	= zeros(nsub,rankx);
sigma	= eye(mvar);			
sigmai	= invpd(sigma);

theta	= zeros(rankz,rankx);
lambda	= eye(rankx);
lambdai	= invpd(lambda);

@ Define data structures for saving iterates & computing posterior means & std @
betam	= zeros(nsub,rankx);	@ posterior mean of beta 		@
betas	= zeros(nsub,rankx);	@ posterior std of beta 		@
c		= mvar*(mvar+1)/2;
sigmag	= zeros(smcmc,c);		@ save iterations for sigma 	@
thetag	= zeros(smcmc,thdim);
c		= rankx*(rankx+1)/2;
lambdag = zeros(smcmc,c);		@ save iterations for lambda 	@
ydatam	= zeros(ntot,mvar);		@ posterior mean utilities		@
ydatas	= ydatam;				@ posterior std utilities		@



/*
********************************************************************
*	Do MCMC
********************************************************************
*/
etime = hsec;
@ Do the initial transition period @
for i1 (1,nblow,1);	imcmc = i1;
	call getprobit;
	if imcmc == 100*floor(imcmc/100);
		dtime = (hsec -etime)/(60*100);
		print "TP Iteration = " imcmc " D.time = " dtime;
		etime = hsec;
	endif;
endfor;

etime = hsec;
for i1 (1,smcmc,1);	imcmc = i1;		@ Save smcmc iterations 			@
	for i2 (1,skip,1); jmcmc = i2;		@ Save every skip iterations 	@
		call getprobit;
	endfor;
	if imcmc == (100/skip)*floor(skip*imcmc/100);
		dtime 	= (hsec - etime)/(60*100);
		tit		= nblow + skip*imcmc;
		print "Iteration = " tit " D.time = " dtime;
		etime 	= hsec;
	endif;
	@ When saving iterates, divide by sigma[mvar,mvar] or its square root to identify model @
	sqrtsiglast    = sqrt(sigma[mvar,mvar]);
	sigmag[imcmc,.]	= vech(sigma/sigma[mvar,mvar])';	
	                                @ vech({1 2 3, 4 5 6, 7 8 9}) = {1, 4 5, 7 8 9} @
									@ xpnd is the inverse operator of vech			@
	thetag[imcmc,.]	= vecr(theta/sqrtsiglast)';
	betam			= betam 	+ beta/sqrtsiglast;
	betas			= betas 	+ (beta/sqrtsiglast)^2;
	lambdag[imcmc,.]= vech(lambda/sigma[mvar,mvar])';
	ydatam			= ydatam + ydata/sqrtsiglast;
	ydatas			= ydatas + (ydata/sqrtsiglast)^2;
endfor;


/*
******************************************************************
*	Compute Posterior Means and STD
******************************************************************
*/
ydatam	= ydatam/smcmc;
betam		= betam/smcmc;
thetam		= reshape(meanc(thetag),rankz,rankx);
sigmam		= xpnd(meanc(sigmag));		@ xpnd reconstructs symmetric matrix @
lambdam		= xpnd(meanc(lambdag));

ydatas		= sqrt( abs(ydatas 	- smcmc*ydatam^2)/smcmc);
betas		= sqrt( abs(betas 	- smcmc*betam^2)/smcmc);
thetas		= sqrt( reshape(stdc(thetag),rankz,rankx) );
sigmas		= xpnd(stdc(sigmag));
lambdas		= xpnd(stdc(lambdag));

if flagtrue == 1;			@ Did a simulation, so we have the true utilities. 	@
	@ Get true parameters if simulation @

	load ydatat		= ydatat;
	load betat 		= betat;
	load sigmat		= sigmat;
	load thetat		= thetat;
	load lambdat	= lambdat;



	@ Pick out each dimension of Y_{ij} and compute fit statistics @
	multir	= zeros(mvar,1);
	rsquare	= zeros(mvar,1);
	stderr	= zeros(mvar,1);
	for fm (1,mvar,1); m = fm;
		ym		= ydatat[.,m];
		yhatm	= ydatam[.,m];
		cm		= corrx(ym~yhatm);
		multir[m]	= cm[1,2];
		rsquare[m]	= cm[1,2]^2;
		resid		= ym - yhatm;
		stderr[m]	= sqrt(resid'resid/ntot);
	endfor;
endif;
	

/*
****************************************************************
*	Do some output
****************************************************************
*/
call outputanal;


@ Plot saved iterations against iterations number @
t 	= seqa(nblow+skip,skip,smcmc);		@ saved iteration number @
title("Latent Error Cov vs Iteration");
xy(t,sigmag);
title("Theta vs Iteration");
xy(t,thetag);
title("Lambda vs Iteration");
xy(t,lambdag);
graphset;

@ Get ydatam ready for output @
@ Add column of 0's for base brand @
utilitym = ydatam~zeros(rows(ydatam),1);
utilitym = reshape(utilitym, nsub, nchoice*(mvar+1));

@ Get names for exporting estimated expected utilities @
c =seqa(1,1,nchoice).*.ones(mvar+1,1);  	@ Choices @
a = ones(nchoice,1).*.seqa(1,1,mvar+1);	@ Alternatives @
unames = 0 $+ "C" $+ ftocv(c,3,0) $+ "A" $+ ftocv(a,3,0);

@ Output Estimated Expected Utility = x*beta to EXCEL file: Alternatives nested in choice Sets @
ok = export(utilitym, "utility.xls", unames);
@ Output Posterior Mean of Beta @
ok = export(betam, "betaMean.xls", xnames);
@ Output Posterior STD of Beta @
ok = export(betas, "betaSTD.xls", xnames);

end;

/*
****************************************************************
* GETPROBIT
*	Does one iteration of the HB regression model.
*	INPUT
*		Global Variables
*	OUTPUT
*		Global Variables
****************************************************************
*/
PROC (0) = getprobit;
local zbl, bi, vibn, vibn12, ebin, yhat, sse, sn, resid, gni, gn, gn12, w, sum1, sum2, 
i0, i, fj, j, xij, yij, sgni, sgn, sgn12, muij, cij, ic,
v, sig11, sig11i, smigni, signi;


	/*
	*******************************************************************
	* Compute quantities used in conditional normal distribution.
	* Need to run cndcov(sigma) before generating the random utilities.
	********************************************************************
	*/
	{smigni, signi} = cndcov(sigma);
	/*
	*******************************************************************
	* smigni is a mvar x (mvar-1) matrix and 
	* used in the conditional mean of Y_{i} given Y_{not i}:  
	* smigni[i,.] = sigma_{i, not i}*sigma_{not i, not i)^{-1}
	* signi is a mvar matrix and
	* signi[i] = STD(Y_{i}| Y_{not i}) 
	*          = sqrt(sigma_{ii} - sigma_{i,not i}*sigma_{not i, not i}^{-1} sigma_{not i,i})
	******************************************************************
	*/



	/*
	********************************************************************
	* Generate Y_{ij}, the utility.
	*
	* If alternative k (k = 1, .., mvar) was selected, then
	* Y_{ij} is N(X_{ij}*beta_i, Sigma) and Y_{ij}[k] >= max(Y_{ij})
	**********************************************************************
	*/
	@ Do multiple loops of generating the Utilities for each MCMC Iteration @

	for i0 (1, nsub, 1); i = i0;
		for fj (1,nchoice,1); j = fj;
			xij		= xdata[iptx[j,1]:iptx[j,2],.];
			ic		= cdata[i,j];				@ Index of selected brand		@
			yij		= ydata[(i-1)*nchoice+j,.]';
			muij		= xij*(beta[i,.]');						@ Mean for y_{ij} 	@



			yij	= rndnigtj(yij, ic, muij, smigni, signi, nygen);
			ydata[(i-1)*nchoice+j,.] = yij';					@ store the utility				@
		endfor;
	endfor;

	/*
	***********************************************************
	* Generate beta
	* beta_i is N(mbin, vbn)
	* vbn 	= ( sum_{j=1}^{n_i} X_{ij}'Sigma^{-1} X_{ij} + Lambda^{-1} }^{-1}
	* mbin	= vbn*( sum_{j=1}^{n_i} X_{ij}'Sigma^{-1}Y_{ij} + Lambda^{-1}*Theta*Z_i)
	**********************************************************
	*/
	zbl		= zdata*theta*lambdai;
	for i0 (1, nsub,1); i = i0;
		sum1	= 0;
		sum2	= 0;
		for fj (1,nchoice,1); j = fj;
			xij		= xdata[iptx[j,1]:iptx[j,2],.];
			yij		= ydata[(i-1)*nchoice+j,.]';
			sum1	= sum1 + xij'sigmai*xij;
			sum2	= sum2 + xij'sigmai*yij;
		endfor;
		vibn		= sum1 + lambdai;
		vibn12		= chol(vibn);
		ebin		= sum2 + zbl[i,.]';
		bi			= cholsol(ebin + vibn12'rndn(rankx,1), vibn12);
		beta[i,.] 	= bi';

	endfor;

	/*
	***************************************************************
	* Generate sigma
	****************************************************************
	*/
	@ Compute SSE		@
	sse	= zeros(mvar,mvar);
	for i0 (1, nsub,1); i = i0;
		for fj (1,nchoice,1); j = fj;
			xij		= xdata[iptx[j,1]:iptx[j,2],.];
			yij		= ydata[(i-1)*nchoice+j,.]';
			resid	= yij - xij*(beta[i,.]');
			sse		= sse + resid*resid';
		endfor;
	endfor;
	sgni				= sg0i + sse;
	sgn					= invpd(sgni);
	{sigmai, sigma}     = wishart(mvar,sfn, sgn);
   
	/*
	***********************************************************
	* Generate Theta and Lambda from multivariate model:
	*	B = Z*Theta + N(0,Lambda)
	************************************************************
	*/
	{theta, lambda, lambdai} = 
	getmulreg(beta,zdata,ztz,theta,lambda,lambdai,v0i,v0iu0,f0n,g0i);

endp;


/*
****************************************************************
* GETMULREG
*	Generate multivariate regression parameters.
*		Yd		= Xd*parmat + epsilon
*	
*	INPUT
*		yd		= dependent variables
*		xd		= independet variables
*		xdtxd		= xd'xd
*		
*		parmat	= current value of coefficient matrix
*		var		= current value of covariance matrix
*		vari		= its inverse
*		v0i		= prior precisions for bmat
*		v0iu0		= prior precision*prior mean for bmat
*		f0n		= posterior df for sigma
*		g0i		= prior scaling matrix inverse for sigma

*		
*	OUTPUT
*		parmat	= updated rankx x mvar coefficient matrix