/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GProbit2.GSS
*-->GProbit1.GSS assumes common design matrix for all subjects
*-->GProbit2.GSS allows for different design matrices.
*		Generats data for HB PROBIT Regression Model
*		Select one of mvar+1 alternatives where the last alternative is
*		the base brand.  
*
*		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}, 0} for l \= mvar+1	
*		Select mvar+1 if max(Y) < 0.
*
*		beta_i is rankx vector
*		epsilon_{ij} is N(0,Sigma)
*       Identification:
*           sigma[mvar,mvar] = 1
*
*
*		X_{ij} is mvar x rankx
*			X will have brand intercepts for the first mvar-1 brands, 
*			and coefficients for price and advertising.
*			To identify the model, we fix the intercept for the last brand to zero.
*
*
*		beta_i	= Theta'Z_i + delta_i
*		delta_i	is N(0,Lambda)
*			Z1 is ln(income) and z2 = family size.
*
*****************************************************************
*/
new;
flagplot = 0;							@ 1 -> do a bunch of plots						@
nsub	= 100;							@ Number of subjects 							@
mvar	=  3;							@ Y_{ij} is mvar vector. Eg: mvar+1 brands 		@
										@ Choice mvar+1 is the base brand				@
yrows	= 10 + ceil(10*rndu(nsub,1));	@ Gives the number of observations per subject 	@

ntot	= sumc(yrows);					@ total number of observations        			@

rankx	= mvar + 2;			@ Brand 1, Brand 2, Brand 3, Price, Advert	@
							@ Intercept for Brand 4 = 0					@

rankz	=  3;				@ # rows of Z_i														@
							@ intercept, ln(income), and family size							@

@ Define some variable names for Gauss file @
@ Use string arrays							@
a			= seqa(1,1,mvar);
brands		= 0 $+ "Brand " $+ ftocv(a,1,0);
brands2		= brands|"Brand 4";




@ ftocv converts a numeric variable to alpha and $+ is character addition 		@
@ so brands is Brand 1, Brand 2, ..., Brand mvar								@

xyname		= brands|"Price"|"Advert"|"Choice";  @ | stacks matrices on top of each other @
zname		= {"Constant", "lnIncome",  "HH Size" };
@ To print a character string use: print $ zname; @


@ define two pointers to access xdata matrix @
@ xdata = { x_{11}, ... x_{1n_1}, x_{21}, ..., x_{2,n_2}, ..., x_{nsub,1} ... x_{nsub,n_{nsub}} } @
@ xij = xdata[lxy[i,j]:uxy[i,j],.] @
lxy		= zeros(nsub,maxc(yrows));		@ gives lower subscript @
uxy		= lxy;							@ gives upper subscript @
s1		= 0;
for i0 (1,nsub,1); i = i0;
	for fj (1,yrows[i],1); j = fj;
		uxy[i,j]	= mvar*(s1+j);
	endfor;
	s1 = s1 + yrows[i];
endfor;
lxy 	= uxy - mvar + 1;
lxy	= (0-lxy).*(lxy .< 0) + lxy;			@ zero-out the negative entries. @

@ Generate error variance Sigma @

sigmat	= { .2  .1  -.1,
            .1  .3  -.05,
            -.1   -.05  1
};
sigt12	= chol(sigmat);


@ Generate error variance Lambda @

@			b1		b2		b3			Price 	Advert		@
lambdat	= {
			1		.3		-.1			0	0	,		@ b1 	    @
			.3		.8		-.05		0	0	,		@ b2 	    @
			-.1		-.05	.5			0	0	,		@ b3 	    @
			0		0		0			.2	.1	,		@ price		@
			0		0		0			.1	.5			@ advert	@
		};
lambdat	= 0.01*lambdat;
lbd12	= chol(lambdat);


@ generate Z variables  @
@Z1 is ln(income) @
m		= ln(40000);
s		= (ln(120000) - m)/1.96;
z1		= m + s*rndn(nsub,1);
@ Mean center z1@
z1		= z1 - meanc(z1);

@Z2 is family size @
z2		= floor(rndnab(nsub,1,3,4,1,10));  
@rndnab is my truncated normal(rows, cols, mean, std, a, b)@
zdata	= z1~(z2-meanc(z2));
zdata	= ones(nsub,1)~zdata;




@ Generate theta @



@		B1 - B4		B2-B4	B3-B4	Price	Advertising					@

thetat	= {
			1		.5		 0		-2		1,		@ Intercept			@
			1		.5		 -.3	1		-.5,	@ Ln(Income)		@
			-.3		-.2	   	0		-.3		.5		@ Family Size		@
		};

@ Generate partworths beta @
betat	= zdata*thetat	+ rndn(nsub,rankx)*lbd12;


@ generate X & Y data @
ydata 	= zeros(ntot*mvar,1);				@ 0/1 choice  	@
ydatat	= ydata;							@ true utilities of Brand j - Brand mvar+1 @
xdata	= zeros(ntot*mvar,rankx);
ipick	= zeros(mvar+1,1);		@ keep track of the number of choices @
for i0 (1,nsub,1); i = i0;
	for fj (1,yrows[i],1); j = fj;
		@ Do subject i, purchase j @
		@ xij = brands, price, advertising @
		xij	= eye(mvar);						@ Brand Intercepts @
		@ Generate Prices											@
        @     	   Regular Price            Price Promotion			@
		p1		= 5.5 + .1*rndn(1,1) - (rndu(1,1) < .4)*.3;
		p2		= 5.5 + .1*rndn(1,1) - (rndu(1,1) < .4)*.3;
		p3		= 5.2 + .1*rndn(1,1) - (rndu(1,1) < .2)*.2;
		p4		= 5.0 + .05*rndn(1,1);
		price	= p1-p4|p2-p4|p3-p4;
		xij 	= xij~price;
		@ Do advertising 			@
		@ Brand 1 heavily advertises, followed by the other four		@
		a1		= rndu(1,1) < .4;					
		a2		= rndu(1,1) < .4;
		a3		= rndu(1,1) < .2;
		a4		= rndu(1,1) < .1;
		advert	= a1-a4|a2-a4|a3-a4;
		xij		= xij~advert;

		@ Generate utility for the four brands @
		bi			= betat[i,.]';
		yij			= xij*bi + sigt12'rndn(mvar,1);
		choice 		= zeros(mvar,1);	
		ib			= maxindc(yij|0);
		if ib <= mvar;
			choice[ib] 	= 1;
		endif;
		ipick[ib]	= ipick[ib] + 1;	

			
		@ Save it in the data matrix 			@
		ydata[lxy[i,j]:uxy[i,j],.] 		= choice;
		ydatat[lxy[i,j]:uxy[i,j],.] 	= yij;
		xdata[lxy[i,j]:uxy[i,j],.] 		= xij;
		@ rows of ydata2 are purchase occassions and columns are brands @
		@ ydata2 will be used for computing summary statiscs for utilities @
		if i == 1 and j == 1;
			ydata2 = yij';
		else;
			ydata2 = ydata2|(yij');
		endif;
	endfor;
endfor;
xydata			= xdata~ydata;

/*
*****************************************************************************
* Output to a Gauss file.
* Gauss files are faster to read, and they assign variable names to 
* the columns.  Also, Gauss has a number of special commands that operate
* on Gauss files, such as file merges, variable recoding, and summary statistics.
*
* First, create a Gauss file.  f1 is the file handle.  
* The Gauss file will be called "XPDATA."       
* The column will be named by the strings in the character array xyname.  	
* ^xyname means use the names in the character string.						
* 0, 8 gives double precision real numbers.								
********************************************************************************
*/

create f1 = xpdata with ^xyname, 0, 8;

/*
********************************************************************************
* Next read data into the Gauss file by using the writer command.			
* f1 is the file handle defined in previous command.						
* xydata is the data matrix that we just created.							
* writer returns the number of rows read to f1.								
* If it is not rows(xydata), something bad happended.						
********************************************************************************
*/

if writer(f1,xydata) /= rows(xydata);
		errorlog "Conversion of XYDATA to Gauss File did not work";
endif;
closeall f1;
@ We do not need f1 anymore, so close it up.								@
@ Do the same for zdata														@
create f1 = zdata with ^zname, 0, 8;
if writer(f1,zdata) /= rows(zdata);
	errorlog "Conversion of ZDATA to Gauss File did not work";
endif;
closeall f1;

save yrows		= yrows;
save lxy		= lxy;
save uxy		= uxy;

save ydatat		= ydatat;
save sigmat		= sigmat;
save betat		= betat;
save thetat		= thetat;
save lambdat	= lambdat;

@ Define formats for fancy printing @
@ Used to print a matrix of alpha & numeric variables @
let fmt1a[1,3]	= "*.*lf" 10 5;			@ Format for printing numeric variable 		@
let fmtsb[1,3]	= "*.*s" 8 8;			@ Format for printing character variable 	@
mask	= zeros(mvar+1,1)~ones(mvar+1,1);	@ 0 for alpha, and 1 for numeric     		@
fmt1	= fmtsb|fmt1a;					@ Format for columns of output				@
bout	= brands2~(ipick/ntot*100);
print "    Brand Market Shares";
print " Brand    Market Share (%)";
flag		= printfm(bout,mask,fmt1);	@ Formated print.							@
print;
mask	= zeros(mvar,1)~ones(mvar,4);
fmt2	= fmtsb|fmt1a|fmt1a|fmt1a|fmt1a;
bout	= brands~meanc(ydata2)~stdc(ydata2)~minc(ydata2)~maxc(ydata2);
print "    Summary Statistics for Brand by Purchase Occasion Utilities";
print " Brand     Mean      STD       MIN       MAX";
flag	= printfm(bout,mask,fmt2);

if flagplot == 1;
_plctrl = -1;			@ use symbols only in plots @
@ plot beta versus ln(income) and family size @
for fj (mvar,rankx,1); j = fj;
	atitle = "Partworth for " $+ xyname[j] $+ " versus " $+ zname[2];
	title(atitle);
	xy(zdata[.,2], betat[.,j]);
	atitle = "Partworth for " $+ xyname[j] $+ " versus " $+ zname[3];
	title(atitle);
	xy(zdata[.,3], betat[.,j]);
	wait;							@ Hit any key to continue @
endfor;

@ plot y versus price for the mvar brands @
for fj (1,mvar-1,1); j = fj;
	@ create a vector to select brands from data matrices @
	b		= zeros(mvar,1); b[j] = 1; 
	bd	 	= ones(ntot,1).*.b;  @ .*. is Kronecker product. bd is a vector of 0 & 1 @
	xj		= selif(xdata, bd);	 @ selif extracts the rows where bd = 1 @
	yj		= selif(ydatat, bd);
	atitle = "Utility for " $+ brands[j] $+ " versus " $+ xyname[mvar];
	title(atitle);
	xy(xj[.,mvar], yj);
	atitle = "Utility for " $+ brands[j] $+ " versus " $+ xyname[mvar+1];
	title(atitle);
	xy(xj[.,mvar+1], yj);
	wait;
endfor;

graphset;			@ Set plots back to default values @

endif;
end;