hbmreg2.gss

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

*		
*	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 )
	*************************************************************
	*/
 
	vb12	= chol(xdtxd.*.vari + v0i);
	ubn		= ( (xd').*.vari )*vecr(yd) + v0iu0;
	par		= cholsol(ubn + vb12'rndn(pdim,1), vb12);
	parmat	= reshape(par,rp,cp);

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




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

if flagtrue == 1;		@ Did a simulation @
	load betat = betat;
	load sigmat = sigmat;
	load thetat = thetat;
	load lambdat = lambdat;
endif;

@ Define formats for fancy printing @
@ Used to print a matrix of alpha & numeric variables @

format 10,5;							@ Default pring format						@
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		   	@
output file = ^outfile reset;			@ outfile is the file handle for the output file @
			@ Route printed output to the defined by outfile 	@
print "Results from HBMReg2.GSS";
print "Hierarchical Bayesian Multivariate Regression Model using MCMC.";
print "Different design matrices for each subject";
print "Y_i = X_i*B_i + U_i";
print "yv_i = (X_i.*.I(nyvar))*beta_i + epsilon_i";
print "		yv_i	= vec(Y_i') = vecr(Y_i)";
print "     beta_i = vec(B_i') = vecr(B_i)";
print "		epsilon_i = vecr(U_i)";
print;
print "beta_i = Theta'z_i + delta_i";
print "epsilon_i is N(0,I.*.Sigma)";
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;
print "Number of subjects	                 = " nsub;
print "Mean # of observations per subject  = " meanc(yrows);
print "STD  # of observations per subject  = " stdc(yrows);
print "MIN  # of observations per subject  = " minc(yrows);
print "MAX  # of observations per subject  = " maxc(yrows);
print "Total number of observations        = " ntot;
print "Number of Y variables               = " nyvar;
print "Number of dependent variables X     = " xdim " (excluding intercept)";
print "Number of dependent variables Z     = " zdim " (excluding intercept)";
print;
print;
print "Variables in first level equation: Y_i = X_i*beta_i + epsilon_i";
print;
print "       Summary Statistics for Y";
call sumstats(ynames,ydata,fmts1,fmts2,fmtn1);
print;
print "       Summary Statistics for X";
call sumstats(xnames,xdata,fmts1,fmts2,fmtn1);
print;
print "Variables in second level equation: beta_i = Theta*z_i + delta_i";
print;
print "      Summary Statistics for Z";
print sumstats(znames,zdata,fmts1,fmts2,fmtn1);
print;
print "-----------------------------------------------------------------------------------";
print;
print "Statistics of Fit Measures for each Y";
sout = {"Variable" "Mult-R" "R-Sq" "ErrSTD"};
print outitle(sout,fmts1,fmts2);
bout	= ynames~br~(br^2)~estd;
print outmat(bout,fmts1,fmtn1);
print;
print;
print "-----------------------------------------------------------------------------------";
print "Estimation of the error covariance Sigma";
sout = "Variable"~(ynames');
if flagtrue == 1;
	print "True Sigma";
	call outitle(sout,fmts1,fmts2); 
	bout	= ynames~sigmat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif; 
print "MLE of Sigma";
call outitle(sout,fmts1,fmts2);
bout = ynames~s2hat;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior Mean of Sigma ";
call outitle(sout,fmts1,fmts2);
bout = ynames~sigmam;
call outmat(bout,fmts1,fmtn1);
print "Posterior STD  of Sigma";
call outitle(sout,fmts1,fmts2);
bout = ynames~sigmas;
call outmat(bout,fmts1,fmtn1);
print;

print "-----------------------------------------------------------------------------------";
print;
print "Statistics for Individual-Level Regression Coefficients";
sout = "Variable"~(ynames');
if flagtrue == 1;
	ebeta	= reshape(meanc(betat),rankx,nyvar);
	sbeta	= reshape(stdc(betat),rankx,nyvar);
	print "Mean of True Beta";

	call outitle(sout,fmts1,fmts2);
	bout = xnames~ebeta;
	call outmat(bout,fmts1,fmtn1);
	print;
	print "STD of True Beta";
	call outitle(sout,fmts1,fmts2);
	bout = xnames~sbeta;
	call outmat(bout,fmts1,fmtn1);
	print;
endif;

@ MLE @
	ebeta	= reshape(meanc(bhat),rankx,nyvar);
	sbeta	= reshape(stdc(bhat),rankx,nyvar);
	print "Mean of MLE for Beta";

	call outitle(sout,fmts1,fmts2);
	bout = xnames~ebeta;
	call outmat(bout,fmts1,fmtn1);
	print;
	print "STD of MLE for Beta";
	call outitle(sout,fmts1,fmts2);
	bout = xnames~sbeta;
	call outmat(bout,fmts1,fmtn1);
	print;

@ HB @
	ebeta	= reshape(meanc(betam),rankx,nyvar);
	sbeta	= reshape(stdc(betam),rankx,nyvar);
	print "Mean of HBE for Beta";

	call outitle(sout,fmts1,fmts2);
	bout = xnames~ebeta;
	call outmat(bout,fmts1,fmtn1);
	print;
	print "STD of HBE for Beta";
	call outitle(sout,fmts1,fmts2);
	bout = xnames~sbeta;
	call outmat(bout,fmts1,fmtn1);
	print;


print;
print "-----------------------------------------------------------------------------------";
print;
if flagtrue == 1;
	print "Comparison of True Beta to Individual Level Estimates";
	ic 	= 0;
	for i0 (1,rankx,1); i = i0;
		for j0 (1,nyvar,1); j = j0;
		ic = ic + 1;
		sout = "Beta for X variable " $+ xnames[i] $+ " and Y variable " $+ ynames[j];
		print $ sout;
		cb 		= corrx( betat[.,ic]~betam[.,ic] );
		rmse 	= betat[.,ic] - betam[.,ic];
		rmse	= rmse'rmse;
		rmse	= sqrt(rmse/nsub);
		print "Correlation between true and HB  = " cb[1,2];
		print "RMSE between true and HB         = " rmse;
		print;
		cb 		= corrx( betat[.,ic]~bhat[.,ic] );
		rmse 	= betat[.,ic] - bhat[.,ic];
		rmse	= rmse'rmse;
		rmse	= sqrt(rmse/nsub);
		print "Correlation between true and MLE = " cb[1,2];
		print "RMSE between true and MLE        = " rmse;
		print;
		endfor;
	endfor;
endif;
print "-----------------------------------------------------------------------------------";
print;
print "HB Estimates of Theta";
sout	= "  "~(bnames');
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 "Posterior Mean/Posterior STD";
call outitle(sout,fmts1,fmts2);
bout	= znames~(thetam./thetas);
call outmat(bout,fmts1,fmtn1);
print;
print;
print "-----------------------------------------------------------------------------------";
print;
sout = "  "~(bnames');
print "HB Estimate of Lambda";
if flagtrue == 1;
	print "True Lambda";
	call outitle(sout,fmts1,fmts2);
	bout = bnames~lambdat;
	call outmat(bout,fmts1,fmtn1);
	print;
endif;
print "Posterior Mean of Lambda";
call outitle(sout,fmts1,fmts2);
bout = bnames~lambdam;
call outmat(bout,fmts1,fmtn1);
print;
print "Posterior STD of Lambda";
call outitle(sout,fmts1,fmts2);
bout = bnames~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;