logit1.gss

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

/*
************************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*	LOGIT1.GSS
*------>Common choice sets for each subject.
*		Select one of mvar+1 alternatives where the last alternative is
*		the base brand.  
*
*		P(C_{ij} = k|X_{j}) = exp(X_{j}[k]'beta_i)/(1 + sum exp(X_{j}[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_{j} is mvar x rankx
*			X will have brand intercepts for the first mvar brands. 
*
*
*		beta_i	= Theta'Z_i + delta_i
*		delta_i	is N(0,Lambda)
*
*	Priors:
*		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 	@
inx 		= "xdata";			@ 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		= 1000;		@ number of iterations to save for analysis 				@
skip			= 1;			@ Save every skip iterations								@
nblow		= 500;		@ Initial transition iterations 							@
nmcmc		= nblow + skip*smcmc;	@ total number of iterations					@
metstd		= 0.1;					@ STD for random walk metropolis				@

/*
********************************************************************
*	Get data
********************************************************************
*/
load cdata 	= cdata;					@ Gives choices @
load iptx	= iptx;						@ Pointer into x matrix					@
mvar		= iptx[1,2];				@ mvar + 1 brands						@
nchoice		= cols(cdata);				@ number of choice sets					@
nsub		= rows(cdata);				@ number of subjects					@
ntot		= nsub*nchoice;				@ total number of choices				@
@ Input Gauss files @
open f1 	= ^inx;					@ Get Gauss file for X, Y data 				@
										@ Opens Gauss file & assigns file handle f1 @
xdata		= 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);
xnames		= setvars(inx);			@ Get the variable names that accompnay X, Y data @
ynames		= xnames[1:mvar];			@ Use names of interceps for names of components of Y @
rankx		= cols(xdata);


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 @
	


@ 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
******************************************************************
*/

beta 	= zeros(nsub,rankx);
{llbeta, umean}	= loglike(beta);			@ Compute log-likelihood at beta 	@
									@ llbeta is used in metropolis		@
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 		@
thetag	= zeros(smcmc,thdim);
c		= rankx*(rankx+1)/2;
lambdag = zeros(smcmc,c);		@ save iterations for lambda 	@
c2		= nchoice*(mvar+1);
umeanm	= zeros(nsub, c2);	@ matrix of estimated expected utilities @



/*
********************************************************************
*	Do MCMC
********************************************************************
*/
mixp	= 0;
@ Do the initial transition period @
for i1 (1,nblow,1);	imcmc = i1;
	call getlogit;
endfor;


for i1 (1,smcmc,1);	imcmc = i1;		@ Save smcmc iterations 			@
	for i2 (1,skip,1); jmcmc = i2;		@ Save every skip iterations 	@
		call getlogit;
	endfor;
	@ Save the random deviates @
	thetag[imcmc,.]	= vecr(theta)';
	betam			= betam 	+ beta;
	betas			= betas 	+ beta^2;
	lambdag[imcmc,.]= vech(lambda)';
	umeanm		= umeanm + umean;
endfor;


/*
******************************************************************
*	Compute Posterior Means and STD
******************************************************************
*/
mixp		= mixp/(nmcmc*nsub);
umeanm	= umeanm/smcmc;
betam		= betam/smcmc;
thetam		= reshape(meanc(thetag),rankz,rankx);
lambdam		= xpnd(meanc(lambdag));


betas		= sqrt( abs(betas 	- smcmc*betam^2)/smcmc);
thetas		= reshape(stdc(thetag),rankz,rankx);
lambdas		= xpnd(stdc(lambdag));


	

/*
****************************************************************
*	Do some output
****************************************************************
*/
call outputanal;


@ Plot saved iterations against iterations number @
t 	= seqa(nblow+skip,skip,smcmc);		@ saved iteration number @
title("Theta versus Iteration");
xy(t,thetag);
title("Lambda versus Iteration");
xy(t,lambdag);
graphset;

@ Get names for exporting estimated expected utilities @
c =seqa(1,1,nchoice).*.ones(mvar+1,1);  	@ Choices @
a = ones(nchoice,1).*.seqa(1,1,mvar+1);	@ Alternatives @
unames = 0 $+ "C" $+ ftocv(c,3,0) $+ "A" $+ ftocv(a,3,0);

@ Output Estimated Expected Utility = x*beta to EXCEL file: Alternatives nested in choice Sets @
ok = export(umeanm, "utility.xls", unames);
@ Output Posterior Mean of Beta @
ok = export(betam, "betaMean.xls", xnames);
@ Output Posterior STD of Beta @
ok = export(betas, "betaSTD.xls", xnames);
end;

/*
****************************************************************
* GETLOGIT
*	Does one iteration of the HB regression model.
*	INPUT
*		Global Variables
*	OUTPUT
*		Global Variables
****************************************************************
*/
PROC (0) = getlogit;
local 
	thz,bnew,llbnew,i0,i,mui,testp,vibn12,ebin,resid,gni,gn,gn12,w, u;

	thz		= zdata*theta;	
	/*
	***********************************************************
	* Generate beta
	* Use symmetric random walk metropolis
	**********************************************************
	*/
	bnew	= beta + metstd*rndn(nsub,rankx);
	{llbnew, umean}  = loglike(bnew);					@ Get log likelihood at new beta @
	u	= ln(rndu(nsub,1));
	@ Do metropolis step  @
	for i0 (1,nsub, 1); i = i0;
		mui  = thz[i,.]';
		testp = llbnew[i] - llbeta[i] 								@log Likelihood @
			- ( bnew[i,.]' -  mui )'lambdai*( bnew[i,.]' - mui )/2	@ "Priors"		@	
			+ ( beta[i,.]' -  mui )'lambdai*( beta[i,.]' - mui )/2;

		if u[i] < testp;  @ Accept Candidate @
			mixp		= mixp + 1;
			beta[i,.] 	= bnew[i,.];
			llbeta[i]  	= llbnew[i];
		endif;			@ else reject candidate & stay pat @
	endfor;
   
	/*
	***********************************************************
	* 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);
****************************************************************
*/
PROC (3) = getmulreg(yd,xd,xdtxd,parmat,var,vari,v0i,v0iu0,f0n,g0i);
local vb12, ubn, par, pdim, resid, gni, gn, rp, cp;

	rp		= rows(parmat);
	cp		= cols(parmat);
	pdim	= rp*cp;

	/*
	***********************************************************
	* Generate parmat from N_{rp x cp}(M,v)
	* par = vecr(parmat)
	* par is N(u,V) whee u = vec(M');
	*	V = (Xd'Xd.*.Var^{-1} + V_0^{-1})^{-1}
	*	u = V*( (Xd'.*.Var^{-1})*vec(Yd') + V_0^{-1}u_0 )
	*************************************************************
	*/