hbmreg1.gss

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

/*
************************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*	HBMREG1.GSS
*	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*Theta + Delta
*		delta_i	is N(0,Lambda)
*			Z1 is ln(income) and z2 = family size.
*	PRIORS
*		Sigma is Inverted Wishart(sf0,sg0)
*		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		= "results2.txt";	@ Specify output file for saving results 					@
								@ outfile is a string variable that contains a file name 	@
inxy 		= "xydata";			@ Name of Gauss file with X,Y data							@
inz	 		= "zdata";			@ Name of Gauss file with Z data							@
flagtrue	= 1;				@ 1 -> knows true parameters from simulation				@

/*
********************************************************************
*	Initialize parameters for MCMC 
********************************************************************
*/
smcmc		= 10000;		@ number of iterations to save for analysis 				@
skip		= 1;		@ Save every skip iterations								@
nblow		= 10000;		@ Initial transition iterations 							@
nmcmc		= nblow + skip*smcmc;	@ total number of iterations					@

/*
********************************************************************
*	Get data
********************************************************************
*/
@ Get dimensions and pointers @
load yrows	= yrows;	@ Number of observations per subject 	@
load lxy	= lxy;		@ xij = xdata[lxy[i,j]:uxy[i,j],.]		@
load uxy	= uxy;
mvar	= uxy[1,2] - uxy[1,1];			@ Y_{ij} is a mvar vector. 				@

@ Input Gauss files @
open f1 	= ^inxy;					@ Get Gauss file for X, Y data 				@
										@ Opens Gauss file & assigns file handle f1 @
xydata		= readr(f1,rowsf(f1));		
			@ 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);
rankxy		= cols(xydata);
rankx		= rankxy - 1;				@ # of X variables (includes intercept) @
xynames		= setvars(inxy);			@ Get the variable names that accompnay X, Y data @
xnames		= xynames[1:rankx];
ynames		= xynames[1:mvar];			@ Use names of interceps for names of components of Y @

@ Last row of xydata is Y. @

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) @
nsub	= rows(yrows);					@ number of subjects		@
ntot	= sumc(yrows);

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

@ Prior for sigma is IW(sf0, gs0) @
sf0 	= mvar+2; sfn = sf0 + ntot;
sg0i	= eye(mvar);



@ 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);
sigma	= eye(mvar);
sigmai	= invpd(sigma);
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 		@
c		= mvar*(mvar+1)/2;
sigmag	= zeros(smcmc,c);		@ save iterations for sigma 	@
thetag	= zeros(smcmc,thdim);
c		= rankx*(rankx+1)/2;
lambdag = zeros(smcmc,c);		@ save iterations for lambda 	@



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

for i1 (1,smcmc,1);	imcmc = i1;		@ Save smcmc iterations 			@
	for i2 (1,skip,1); jmcmc = i2;		@ Save every skip iterations 	@
		call gethbmreg;
	endfor;
	sigmag[imcmc,.]	= vech(sigma)';	@ vech({1 2 3, 4 5 6, 7 8 9}) = {1, 4 5, 7 8 9} @
									@ xpnd is the inverse operator of vech			@
	thetag[imcmc,.]	= vecr(theta)';
	betam			= betam 	+ beta;
	betas			= betas 	+ beta^2;
	lambdag[imcmc,.]= vech(lambda)';
endfor;


/*
******************************************************************
*	Compute Posterior Means and STD
******************************************************************
*/

betam		= betam/smcmc;
thetam		= reshape(meanc(thetag),rankz,rankx);
sigmam		= xpnd(meanc(sigmag));		@ xpnd reconstructs symmetric matrix @
lambdam		= xpnd(meanc(lambdag));

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

@ Compute predictive value of Y_{ij} @
yhat		= zeros(rows(xydata),1);
for i0 (1, nsub, 1); i = i0;
	for fj (1, yrows[i],1); j = fj;
		xij		= xydata[lxy[i,j]:uxy[i,j],1:rankx];
		yij		= xydata[lxy[i,j]:uxy[i,j],rankxy];
		yhat[lxy[i,j]:uxy[i,j]] = xij*(betam[i,.]');
	endfor;
endfor;

@ Pick out each dimension of Y_{ij} and compute fit statistics @
multir	= zeros(mvar,1);
rsquare	= zeros(mvar,1);
stderr	= zeros(mvar,1);
for fm (1,mvar,1); m = fm;
	b 		= zeros(mvar,1); b[m] = 1;
	bn 		= ones(ntot,1).*.b;
	ym		= selif(xydata[.,rankxy],bn);
	yhatm	= selif(yhat, bn);
	cm		= corrx(ym~yhatm);
	multir[m]	= cm[1,2];
	rsquare[m]	= cm[1,2]^2;
	resid		= ym - yhatm;
	stderr[m]	= sqrt(resid'resid/ntot);
endfor;
	

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


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

end;

/*
****************************************************************
* GETHBMREG
*	Does one iteration of the HB regression model.
*	INPUT
*		Global Variables
*	OUTPUT
*		Global Variables
****************************************************************
*/
PROC (0) = gethbmreg;
local betai, vibn, vibn12, ebin, yhat, sse, sn, resid, gni, gn, gn12, w, sum1, sum2, 
i0, i, fj, j, xij, yij, sgni, sgn, sgn12,mz;

	/*
	***********************************************************
	* Generate beta
	* beta_i is N(mbin, vbn)
	* vbn 	= ( sum_{j=1}^{n_i} X_{ij}'Sigma^{-1} X_{ij} + Lambda^{-1} }^{-1}
	* mbin	= vbn*( sum_{j=1}^{n_i} X_{ij}'Sigma^{-1}Y_{ij} + Lambda^{-1}*Theta*Z_i)
	**********************************************************
	*/

	mz	= zdata*theta*lambdai;			
	for i0 (1, nsub,1); i = i0;
		sum1	= 0;
		sum2	= 0;
		for fj (1,yrows[i],1); j = fj;
			xij		= xydata[lxy[i,j]:uxy[i,j],1:rankx];
			yij		= xydata[lxy[i,j]:uxy[i,j],rankxy];
			sum1	= sum1 + xij'sigmai*xij;
			sum2	= sum2 + xij'sigmai*yij;
		endfor;
		vibn		= sum1 + lambdai;
		vibn12		= chol(vibn);
		ebin		= sum2 + mz[i,.]';
		betai		= cholsol(ebin + vibn12'rndn(rankx,1), vibn12);
		beta[i,.]	= betai';
	endfor;

	/*
	***************************************************************
	* Generate sigma2
	****************************************************************
	*/
	@ Compute SSE		@
	sse		= 0;
	for i0 (1, nsub,1); i = i0;
		for fj (1,yrows[i],1); j = fj;
			xij		= xydata[lxy[i,j]:uxy[i,j],1:rankx];
			yij		= xydata[lxy[i,j]:uxy[i,j],rankxy];
			resid	= yij - xij*(beta[i,.]');
			sse		= sse + resid*resid';
		endfor;
	endfor;
	sgni			= sg0i + sse;
	sgn			= invpd(sgni);

	{sigmai, sigma} 	= wishart(mvar,sfn,sgn);

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