logit2.gss

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

	/*
	*********************************************************
	* Generate Var
	* 	Var^{-1} is Wishart df = f0n, scale matrix = gn
	*********************************************************
	*/
	resid 			= yd - xd*parmat;
	gni				= g0i + resid'resid;
	gn  			= invpd(gni);
	{vari, var} 	= wishart(cp,f0n,gn);

retp(parmat,var,vari);
endp;


/*
*********************************************************************************
* LOGLIKE 
*	Computes multinomial-logit log-likelihood at parameter beta
*	INPUT:
*		beta
*	OUTPUT:
*		Log likelihood
*
********************************************************************************
*/
PROC (1) = loglike(beta);
local  llbeta,i0,i,fj,j,xij,cij,zij,muij,emuij;
llbeta    = zeros(nsub,1);				@ ln(Likelihood x prior) of beta for each subject @
for i0 (1, nsub, 1);	i = i0;					@ loop over subjects   @
	for fj (1,yrows[i],1);	j = fj;				@ loop over selections @
		xij		= xpdata[lxy[i,j]:uxy[i,j],1:rankx];	@ Get design matrix 	@
		cij		= xpdata[lxy[i,j]:uxy[i,j],rankxp];		@ Get vector of choices @
		if maxc(cij) == 0;		@ Picked base brand @
			zij	= mvar+1;
		else;					@ Picked one of brands 1 to mvar @
			zij = maxindc(cij);
		endif;
		muij 		= xij*(beta[i,.]');		@ logits   @
		muij 		= muij|0;
		emuij 		= exp(muij);
		emuij		= sumc(emuij);
		llbeta[i]	= llbeta[i] + muij[zij] -  ln(emuij);
	endfor;			
endfor;
retp(llbeta);
endp;


/*
*****************************************************************************************************
*  OUTPUTANAL
*	Does analysis of output and save some results 
****************************************************************************************************
*/
PROC (0) = outputanal;
format 10,5;
local  bout, sout, ebeta, sbeta, cb, rmse, fmts1,fmts2, fmtn1, fmtn2, a,b, flag, i0, i,
betat, thetat, lambdat;

@ Get true parameters if simulation @
if flagtrue == 1;
	load betat 		= betat;
	load thetat		= thetat;
	load lambdat	= lambdat;
endif;


let fmtn1[1,3]	= "*.*lf" 10 5;			@ Format for printing numeric variable 		@
let fmtn2[1,3]	= "*.*lf" 10  0;		@ Format for numeric variable, no decimal	@

let fmts1[1,3]	= "-*.*s" 10 9;			@ Format for alpha, left justify		 	@
let fmts2[1,3]  = "*.*s" 10 9;			@ Format for alpha, right justify		  	@
format 10,5;							@ Default pring format						@


output file = ^outfile reset;	@ outfile is the file handle for the output file 	@
								@ Route printed output to the defined by outfile 	@
print "Results from LOGIT1.GSS";
print "Hierarchical Bayes Multivariate LOGIT";
print;
print "Select one of " mvar+1 " alternatives.";
print;
print "beta_i = Theta'z_i + delta_i";
print "delta_i is N(0, Lambda)";
print;
print "Ouput file: " getpath(outfile);		@ File assigned to file handle outfile @
datestr(date);								@ Print the current data				@
print;
print;
print "-----------------------------------------------------------------------------------";
print;
print "MCMC Analysis";
print;
print "Total number of MCMC iterations           	   = " nmcmc;
print "Number of iterations used in the analysis 	   = " smcmc;
print "Number in transition period               	   = " nblow;
print "Number of iterations between saved iterations = " skip-1;
print "Proportion of Metroplis Steps Accepted        = " mixp;
print;
print "-----------------------------------------------------------------------------------";
print "Number of subjects	              = " nsub;
print "Number of observations per subject:";
print "     Average  = " meanc(yrows);
print "     STD      = " stdc(yrows);
print "     MIN      = " minc(yrows);
print "     MAX      = " maxc(yrows);
print;
print "Total number of Choices            = " ntot;
print "Number of Alternatives             = " mvar+1;
print "Number of dependent variables X    = " rankx " (including intercept)";
print "Number of dependent variables Z    = " rankz " (including intercept)";
print;



print "Dependent variables are " $xpnames[rankxp];
print "		Summary Statistics for Choices (excluding Base)";
call sumstats(ynames,reshape(xpdata[.,rankxp],ntot,mvar),fmts1,fmts2,fmtn1);
print;
print "First level equation for Logits";
print "logit_{ij} = X_{ij}*beta_i + epsilon_{ij}";
print;
print "		Summary Statistics for X";
call sumstats(xnames[mvar+1:rankx],xpdata[.,mvar+1:rankx],fmts1,fmts2,fmtn1);
print;
print "Second level equation:";
print "beta_i = Theta*z_i + delta_i";
print;
print "		Summary Statistis for Z:";
call sumstats(znames[2:rankz],zdata[.,2:rankz],fmts1,fmts2,fmtn1);
print;



print "-----------------------------------------------------------------------------------";
print;
print "Statistics for Individual-Level Regression Coefficients";
if flagtrue == 1;
	ebeta	= meanc(betat);
	sbeta	= stdc(betat);
	print "True Beta";
	sout	= {"Variable" "Mean" "STD"};
	call outitle(sout,fmts1,fmts2);
	bout = xnames~ebeta~sbeta;
	call outmat(bout,fmts1,fmtn1);
endif;

print "HB Estimates of Beta";
sout	= {"Variable" "PostMean" "PostSTD" };
call outitle(sout,fmts1,fmts2);
ebeta	= meanc(betam);
sbeta	= sqrt( meanc( (betas^2)) + stdc( betam)^2);
bout = xnames~ebeta~sbeta;
call outmat(bout,fmts1,fmtn1);
print;
print "-----------------------------------------------------------------------------------";

if flagtrue == 1;
	print "Comparison of True Beta to Individual Level Estimates";
	for i0 (1,rankx,1); i = i0;
		print "Variable " $ xnames[i];
		cb 		= corrx( betat[.,i]~betam[.,i] );
		rmse 	= betat[.,i] - betam[.,i];
		rmse	= rmse'rmse;
		rmse	= sqrt(rmse/nsub);
		print "Correlation between true and HB  = " cb[1,2];
		print "RMSE between true and HB         = " rmse;
		print;
	endfor;
endif;
print "-----------------------------------------------------------------------------------";
print;
print "HB Estimates of Theta";
sout	= "  "~(xnames');
if flagtrue == 1;
	print "True Theta";
	call outitle(sout,fmts1,fmts2);
	bout = znames~thetat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif;
print "Posterior Mean of Theta";
print outitle(sout,fmts1,fmts2);
bout = znames~thetam;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior STD of Theta";
call outitle(sout,fmts1,fmts2);
bout	= znames~thetas;
call outmat(bout,fmts1,fmtn1);
print;
print "-----------------------------------------------------------------------------------";
print;
sout = "  "~(xnames');
print "HB Estimate of Lambda";
if flagtrue == 1;
	print "True Lambda";
	call outitle(sout,fmts1,fmts2);
	bout = xnames~lambdat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif;
print "Posterior Mean of Lambda";
call outitle(sout,fmts1,fmts2);
bout = xnames~lambdam;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior STD of Lambda";
call outitle(sout,fmts1,fmts2);
bout = xnames~lambdas;
call outmat(bout,fmts1,fmtn1);
print;
print "===================================================================================";



output off;
closeall;
endp;




/*
*****************************************************************************************
* OUTITLE
*	Prints header for columns of numbers.
*	INPUT
*		a 	= character row vector of column names
*		fmts1	= format for first column
*		fmts2	= format for second column
*	OUTPUT
*		None
******************************************************************************************
*/
PROC (0) = outitle(a,fmt1,fmt2);
local mask, fmt, flag, ncols;
ncols	= cols(a);
mask	= zeros(1,ncols);
fmt		= fmt1|(ones(ncols-1,1).*.fmt2);
flag	= printfm(a,mask,fmt);
print;
endp;
/*
***************************************************************************************
* OUTMAT
*	Outputs a matrix:
*	(Character Vector)~(Numeric matrix);
*	The entries in the numeric matrix have the same format
*	INPUT
*		bout		= matrix to be printed
*		fmts		= format for string
*		fmtn		= format for numeric matrix
*	OUTPUT
*		None
******************************************************************************************
*/
PROC (0) = outmat(bout,fmts,fmtn);
local fmt,mask,flag,ncols, nrows;
ncols		= cols(bout);
nrows		= rows(bout);
fmt			= fmts|(ones(ncols-1,1).*.fmtn);
mask		= zeros(nrows,1)~ones(nrows,ncols-1);
flag		= printfm(bout,mask,fmt);
print;
endp;

/*
*****************************************************************************************
* SUMSTATS
*	Prints summary statistics for a data matrix
*	INPUT
*		names		= charater vector of names
*		data		= data matrix
*		fmts1		= format for string
*		fmts2		= format for string
*		fmtn		= format for numbers
*	OUTPUT
*		None
********************************************************************************************
*/
PROC (0) = sumstats(names,data,fmts1,fmts2,fmtn);
local a, bout;
a	= {"Variable" "Mean" "STD" "MIN" "MAX"};
call outitle(a,fmts1,fmts2);
bout	= names~meanc(data)~stdc(data)~minc(data)~maxc(data);
call outmat(bout,fmts1,fmtn);
endp;