/*
************************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*	PROBIT2.GSS
*	HB Probit Regression Model.
*-->PROBIT1.GSS has common design matrix for all subjects
*-->PROBIT2.GSS allows different design matrices.
*
*   Uses McCulloch & Rossi's method for handling identification.
*
*		Select one of mvar+1 alternatives.
*
*		Y_{ij} 	= X_{ij}*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, ..., yrows[i]
*				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)
*       Sigma is full
*       Divide by last element after MCMC to identify coefficients.
*
*
*		X_{ij} is mvar x rankx
*			X will have brand intercepts, 
*			and coefficients for price and advertising relative to the base brand.
*
*
*		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 	@
inxy 		= "xpdata";			@ Name of Gauss file with X, Choice data					@
inz	 		= "zdata";			@ Name of Gauss file with Z data							@
flagtrue	= 1;				@ 1 -> knows true parameters from simulation				@

/*
********************************************************************
*	Initialize parameters for MCMC 
********************************************************************
*/
smcmc		= 200;		@ 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 yrows	= yrows;	@ Number of observations per subject 	@
load lxy	= lxy;		@ xij = xdata[lxy[i,j]:uxy[i,j],.]		@
load uxy	= uxy;
	
nsub	= rows(yrows);					@ number of subjects 					@
mvar	= uxy[1,2] - uxy[1,1];			@ Y_{ij} is a mvar vector. 				@
ntot	= sumc(yrows);

@ Input Gauss files @
open f1 	= ^inxy;					@ Get Gauss file for X, Y data 				@
										@ Opens Gauss file & assigns file handle f1 @
xpdata		= readr(f1,rowsf(f1));		@ "p" for picks								@	
			@ readr reads in Gauss file with file handle f1.			@
			@ rowsf(f1) returns the number of rows in the Gauss file. 	@
			@ readr reads rowsf(f1) rows, which is the entir dataset.	@
ci			= close(f1);
xpnames		= setvars(inxy);			@ Get the variable names that accompnay X, Y data @
ynames		= xpnames[1:mvar];			@ Use names of interceps for names of components of Y @
ynames2		= ynames|" Base ";
rankxp		= cols(xpdata);
rankx		= rankxp - 1;					@ # of X variables (includes intercept) @
xnames		= xpnames[1:rankx];
@ Last row of xpdata is the choice vector. @

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	= xpdata[.,rankxp];		@ latent y variables		@
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(mvar*ntot,1);	@ 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 parameter, divide by sqrt(sigma[mvar,mvar]) or sigma[mvar,mvar] @
	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		= 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;
		b 		= zeros(mvar,1); b[m] = 1;
		bn 		= ones(ntot,1).*.b;
		ym		= selif(ydatat,bn);
		yhatm	= selif(ydatam, bn);
		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;

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,yrows[i],1); j = fj;
			xij		= xpdata[lxy[i,j]:uxy[i,j],1:rankx];
			cij		= xpdata[lxy[i,j]:uxy[i,j],rankxp];		@ Choice vector		@
			yij		= ydata[lxy[i,j]:uxy[i,j]];
			muij	= xij*(beta[i,.]');						@ Mean for y_{ij} 	@

			if maxc(cij) == 0;		@ selected base brand mvar + 1@
				ic = mvar+1;
			else;				@ selected one of brand 1 to mvar @
				ic = maxindc(cij);
			endif;
			yij	= rndnigtj(yij, ic, muij, smigni, signi, nygen);
			ydata[lxy[i,j]:uxy[i,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,yrows[i],1); j = fj;
			xij		= xpdata[lxy[i,j]:uxy[i,j],1:rankx];
			yij		= ydata[lxy[i,j]:uxy[i,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,yrows[i],1); j = fj;
			xij		= xpdata[lxy[i,j]:uxy[i,j],1:rankx];
			yij		= ydata[lxy[i,j]:uxy[i,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
*		var		= updated variance
*		vari		= updated inverse of sigma
*
*	Calling Statement:
{parmat, var, vari} = getmulreg(yd,xd,xdtxd,parmat,var,vari,v0i,v0iu0,f0n,g0i);
****************************************************************