glogit1.gss

gauss 离散型计量估计源代码,直接下载下来就可以使用

/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GLOGIT1.GSS
*		Generats data for HB LOGIT Regression Model
*-----------> Uses common design matrix
*
*		Select one of mvar+1 alternatives where the last alternative is
*		the base brand.  
*
*		P(C_{ij} = k|X_{ij}) = exp(X_{ij}[k]'beta_i)/(1 + sum exp(X_{ij}[l]'beta_i))

*		for 	i = 1, ..., nsub
*				j = 1, ..., nchoice
*				k = 1, ..., mvar
*		Alternative mvar+1 is the base vector.
*		
*
*		beta_i is rankx vector
*
*		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 = 1;							@ 1 -> do a bunch of plots						@
nsub	= 100;							@ Number of subjects 							@
mvar	=  3;							@ mvar+1 brands									@
										@ Choice mvar+1 is the base brand				@
nchoice = 20;							@ Number of choice sets per subject				@
a1 = seqa(1,mvar,nchoice);
a2 = seqa(mvar,mvar,nchoice);
iptx = a1~a2;							@ iptx is a pointer into the stacked design matrix @

ntot	= nchoice*nsub;					@ total number of observations        			@

rankx	= mvar + 2;			@ Brand 1, Brand 2, Brand 3, Price, Advert	@
							

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|"Base";



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 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);

@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-meanc(z1))~(z2-meanc(z2));
zdata	= ones(nsub,1)~zdata;




@ Generate theta @
@ beta_i = Theta'z_i + delta)i @

@		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;

@ Get common design matrix @
for j0 (1,nchoice,1); j = j0;
		@ 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)*.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);
		xij 	= xij~price;
		@ Do advertising 			@
		@ Brand 1 heavily advertises, followed by the other four		@
		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);
		xij		= xij~advert;
		if j == 1;
			xdata = xij;
		else;
			xdata = xdata|xij;
		endif;
endfor;


cdata 	= zeros(nsub,nchoice);				@ 0/1 choice  	@
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 @

		@ Generate utility for  brands @
		xij = xdata[iptx[j,1]:iptx[j,2],.];
		uij	= xij*(betat[i,.]');
		uij	= uij|0;				@ utility for last brand is 0 @
		pij	= exp(uij);
		pij	= pij./sumc(pij);
		zij	= rndzmn(pij');			@ zij takes values 1 to mvar+1 @
		cdata[i,j] = zij;
	

		ipick[zij]	= ipick[zij] + 1;	

			

	endfor;
endfor;


create f1 = xdata with ^xname, 0, 8;
if writer(f1,xdata) /= rows(xdata);
		errorlog "Conversion of XDATA to Gauss File did not work";
endif;
closeall f1;

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 cdata 		= cdata;
save iptx		= iptx;



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 Share (%)";
flag		= printfm(bout,mask,fmt1);	@ Formated print.							@
print;

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]);

endfor;


graphset;			@ Set plots back to default values @

endif;
end;