/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GPR0BIT1.GSS
*		Generats data for HB PROBIT Regression Model
*------->GPR0BIT1 uses common design matrix.
*------->GPROBIT2 allows each subject to have own design matrix.
*
*		Select one of mvar+1 alternatives where the last alternative is
*		the base brand.  
*
*		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, ..., 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)
*       sigma[1,1] = 1
*		That is, the last brand has fixed variance = 1 
*
*
*		X_{j} 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				@
nchoice = 20;				@ Number of choices per subject 	@
a1 = seqa(1,mvar,nchoice);
a2 = seqa(mvar,mvar,nchoice);
iptx = a1~a2;				@ Pointer into design matrix for each choice set @


ntot	= nsub*nchoice;					@ 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							@

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



@ Generate error variance Sigma @
sigmat	= { .2  .1  -0.1,
            .1  .3  -0.2,
           -0.1 -0.2  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)@
z2		= z2 - meanc(z2);		@ mean center z2 @
zdata	= z1~z2;
zdata	= ones(nsub,1)~zdata;




@ Generate theta @



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

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


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


@ generate X & Y data @
cdata	= zeros(nsub,nchoice);				@ Gives the choice data @

for j0 (1,nchoice,1); j = j0;
		xj	= 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)*.5;
		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;
		xj 		= xj~price;
		@ Do advertising 			@
		a1		= rndu(1,1) < .4;					
		a2		= rndu(1,1) < .3;
		a3		= rndu(1,1) < .2;
		a4		= rndu(1,1) < .1;
		advert	= a1-a4|a2-a4|a3-a4;
		xj		= xj~advert;
		if j == 1;
			xdata = xj;
		else;
			xdata = xdata|xj;
		endif;
endfor;

ipick	= zeros(mvar+1,1);		@ keep track of the number of choices @
for i0 (1,nsub,1); i = i0;
	for fj (1,nchoice,1); j = fj;
		@ Do subject i, purchase j @
		xij = xdata[iptx[j,1]:iptx[j,2],.];

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

			
		@ Save it in the data matrix 			@

		@ 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;
			ydatat = yij';
		else;
			ydatat = ydatat|(yij');
		endif;
	endfor;
endfor;


/*
*****************************************************************************
* 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 = xdata with ^xname, 0, 8;
if writer(f1,xdata) /= rows(xdata);
		errorlog "Conversion of XYDATA to Gauss File did not work";
endif;
closeall f1;


@ 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 iptx			= iptx;			@ Pointer into xdata to get choice set design matrix @
save cdata 		= cdata;
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(ydatat)~stdc(ydatat)~minc(ydatat)~maxc(ydatat);
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 " $+ xname[j] $+ " versus " $+ zname[2];
	title(atitle);
	xy(zdata[.,2], betat[.,j]);
	atitle = "Partworth for " $+ xname[j] $+ " versus " $+ zname[3];
	title(atitle);
	xy(zdata[.,3], betat[.,j]);
	wait;							@ Hit any key to continue @
endfor;



graphset;			@ Set plots back to default values @

endif;
end;