ghbmreg1.gss

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

/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GHBMREG1.GSS
*		Generats data for HB Multivariate Regression Model
*		Y_{ij} 	= X_{ij}*beta_i + epsilon_{ij} 
*			for i = 1, ..., I and j = 1, ..., n_i
*		Y_{ij} is mvar vector
*		beta_i is rankx vector
*		epsilon_{ij} is N(0,Sigma)
*
*			Think of i as household, j as purchase occasion,
*			mvar as # of brands under consideration, and
*			Y as "utility" for the brand.
*
*		X_{ij} is mvar x rankx
*			X will have brand intercepts, 
*			and coefficients for price and advertising.
*
*
*		beta_i	= Theta'Z_i + delta_i
*		B		= Z*Thetwa + D
*		delta_i	is N(0,Lambda)
*			Z1 is ln(income) and z2 = family size.
*****************************************************************
*/
flagplot = 0;							@ 1 -> do a bunch of plots						@
nsub	= 100;							@ Number of subjects 							@
mvar	=  4;							@ Y_{ij} is mvar vector. Eg: mvar brands  		@
@yrows	= 10 + ceil(10*rndu(nsub,1));@	   @ Gives the number of observations per subject 	@
yrows   = 2*ones(nsub,1);
ntot	= sumc(yrows);					@ total number of observations        			@



rankx	=  mvar + 2;  		@ # columns of X_{ij}: think of mvar brands, price and advertising 	@
rankz	=  3;				@ # rows of Z_i														@
							@ Think of 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);		
@ ftocv converts a numeric variable to alpha and $+ is character addition 		@
@ so brands is Brand 1, Brand 2, ..., Brand mvar								@

xynames		= brands|"Price"|"Advert"|"Utility";  @ | stacks matrices on top of each other @
znames		= "Constant"|"lnIncome"|"HH Size";
@ To print a character string use: print $ znames; @


@ 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 @
/***********************
sig12	= {
			1 	.1	0	1,
			0	2	0	2,
			0	0	3	0,
			0	0	0	4
		};

sigmat	= sig12'sig12;
*****************************/

sigmat = {
40  0    7   0,
0   .1   0   0,
7    0   10   -3,
0    0   -3    5
};

sig12 = chol(sigmat);
@ error variance Lambda @
/************************
lbd12	= {
			1	-2	3	0	.5	-.6,
			0	2	1	0	-1	-2,
			0	0	3	0	2	.1,
			0	0	0	4	0	0,
			0	0	0	0	5	2,
			0	0	0	0	0	6
		};
lbd12 = lbd12/2;
lambdat = lbd12'lbd12;
*************************/

lambdat = { 
2   1  0  -1   0   0,
1   3  0  -2   0   0,
0   0  .5   0   0   0,
-1 -2  0   9   0   0,
0   0  0   0  16  -4,
0   0  0   0  -4   8
};

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      = rndnab(nsub,1,3,4,1,10);
z2		= floor(z2);  @rndnab is my truncated normal(rows, cols, mean, std, a, b)@

zdata	= z1~z2;
zdata	= ones(nsub,1)~zdata;




@ Generate theta @

@ High income brands like brands 1 and 2   					@
@ Small families perfer brands 1 and 2    					@
@ Low income & large familes prefer brand 4 (generic brand) @
@ Price partworth is negative, increases with income, and decreases with family size. @
@ Advertising effectiveness decreases with income, and increases with family size.   @

@		x		B1		B2		B3		B4		Price	Adverting 	@ @ Z		@
thetat	= {		-15		-5		10		20		-5		10,			@ Intercept @
				2		1		-1		-2		2		-3,		    @ Ln(Income) @
				-1		-.5		 2		3		-1		2			@ Family Size @
			};


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



@ generate X & Y data @
ydata 	= zeros(ntot*mvar,1);
xdata	= zeros(ntot*mvar,rankx);
ipick	= zeros(mvar,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);						@ Brands 									@
		@ Do price. Brand 1 & 2 are premimum brands with sharpe promotions					@
		@ Brand 3 is a want-to-be, and brand 4 is generic.									@
        @     	   Regular Price            Price Promotion			@
		p1		= 5.50 + .1*rndn(1,1) - (rndu(1,1) < .25)*.5;
		p2		= 5.25 + .1*rndn(1,1) - (rndu(1,1) < .4)*.25;
		p3		= 5.10 + .2*rndn(1,1) - (rndu(1,1) < .2)*.10;
		p4		= 5.00 + .1*rndn(1,1);
		price	= p1|p2|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|a2|a3|a4;
		xij		= xij~advert;
		@ Generate utility for the four brands @
		yij	= xij*(betat[i,.]') + sig12'rndn(mvar,1);
		bchoice	= maxindc(yij);					@ bchoice is the brand with the maximum utility @
		ipick[bchoice]	= ipick[bchoice] + 1;
		@ Save it in the data matrix 			@
		ydata[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;
@ Create a Gauss file.  f1 is the file handle.  @
create f1 = xydata with ^xynames, 0, 8;
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 ^znames, 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 xdata 		= xdata;

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)~ones(mvar,1);	@ 0 for alpha, and 1 for numeric     		@
fmt1	= fmtsb|fmt1a;					@ Format for columns of output				@
bout	= brands~(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 (1,rankx,1); j = fj;
	atitle = "Partworth for " $+ xynames[j] $+ " versus " $+ znames[2];
	title(atitle);
	xy(zdata[.,2], betat[.,j]);
	atitle = "Partworth for " $+ xynames[j] $+ " versus " $+ znames[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); 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(ydata, bd);
	atitle = "Utility for " $+ xynames[j] $+ " versus " $+ xynames[mvar+1];
	title(atitle);
	xy(xj[.,mvar+1], yj);
	atitle = "Utility for " $+ xynames[j] $+ " versus " $+ xynames[mvar+2];
	title(atitle);
	xy(xj[.,mvar+2], yj);
	wait;
endfor;
graphset;
endif;

end;