hbmreg2.gss

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

/*
************************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*	HBMREG2.GSS
*
*	HBMREG1*.GSS
*		Has a matrix X_{ij} for each (ij) observation and a vector beta_i.
*	HBMREG2*.gss
*		has a vector x_{ij} for each (ij) obsevations and a matrix B_i.
*
*	HB Multiple Regession Model
*		Y_i 	= X_i*B_i + U_i
*
*		Y_i is a n_i by nyvar matrix.
*		X_i is a n_i by rankx design matrix
*		B_i is rankx by nyvar matrix of regression coefficients.
*
*		Define:
*			yv_i = vec(Y_i')
*			beta_i = vec(B_i')
*			epsilon_i = vec(U_i') is N(0,I(n_i).*.Sigma)
*		Then
*			yv_i = (X_i.*.I(m))*beta_i + epsilon_i
*		Put
*			beta_i	= Theta'Z_i + delta_i
*
*			delta_i	is N(0,Lambda)
*
*
*	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		= "results1.dat";	@ 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		= 1000;		@ number of iterations to save for analysis 				@
skip		= 1;		@ Save every skip iterations								@
nblow		= 1000;		@ Initial transition iterations 							@
nmcmc		= nblow + skip*smcmc;	@ total number of iterations					@

/*
********************************************************************
*	Get data
********************************************************************
*/
load nyvar = nyvar;			@ Number of y variables @
load yrows = yrows;			@ Number of observations per subject @
nsub	= rows(yrows);		@ Number of subjects 				@
ntot	= sumc(yrows);
@ 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);
xynames		= setvars(inxy);			@ Get the variable names that accompnay X, Y data @
nxy			= cols(xydata);
xdim		= nxy - nyvar;
rankx		= xdim + 1;
ynames		= xynames[xdim+1:nxy];	
xnames		= "CS"|xynames[1:xdim];



@ Last nyvar rows  of xydata is Y. @
xdata		= xydata[.,1:xdim];			@ Get independent variables 	@
ydata		= xydata[.,xdim+1:nxy];		@ Get dependent variable		@

open f1		= ^inz;
zdata		= readr(f1,rowsf(f1));	
ci			= close(f1);
znames		= setvars(inz);
znames		= "CS"|znames;


@ Create pointer based on yrows to access xdata and ydata @
b	= cumsumc(yrows);	@ cumsumc is the cumulative sum of yrows @
a	= 1|(1+b[1:nsub-1]);
iptxy	= a~b;
@ To use iptxy to get the design matrix and data for subject i: @
@ x_i = xdata[iptxy[i,1]:iptxy[i,2],.] @
@ y_i = ydata[iptxy[i,1]:iptxy[i,2]] @


zdim	= cols(zdata);
@ add intercepts to xdata and zdata @
xdata 	= ones(ntot,1)~xdata;
zdata	= ones(nsub,1)~zdata;
nparb	= rankx*nyvar;				@ Number of regression coefficients @
rankz	= cols(zdata);
thdim	= nparb*rankz;				@ dimension of vec(theta) @

@ Get names for beta @
bnames = zeros(nparb,1);
ic		= 0;
for i0 (1,rankx,1); i = i0;
	for fj (1,nyvar,1); j = fj;
		ic = ic + 1;
		bnames[ic] = xnames[i] $+ "(" $+ ynames[j] $+ ")";
	endfor;
endfor;

@ Compute some sufficient statistics @
@ Get point to access stacked matrices of xtx  @
b	= seqa(rankx,rankx,nsub);
a	= 1|(1+b[1:nsub-1]);
iptxt	= a~b;
xtx	= zeros(rankx*nsub,rankx);
@ Pooled estimated of beta @
bpool 	= invpd(xdata'xdata)*xdata'ydata;
bhat	= zeros(nsub,nparb);		@ MLE of beta_i @
sse	= zeros(nyvar,nyvar);
for i0 (1,nsub,1);	i = i0;
	xi	= xdata[iptxy[i,1]:iptxy[i,2],.];
	yi	= ydata[iptxy[i,1]:iptxy[i,2],.];
	xitxi	= xi'xi;
	if rank(xitxi) >= rankx;		@ Individual level mle's exist @
		xitxii	= invpd(xitxi);
		xityi	= xi'yi;
		xtx[iptxt[i,1]:iptxt[i,2],.] = xitxi;
		bhati = xitxii*xityi;
	else;			@ Use pooled estimate for a common beta @
		bhati = bpool;
	endif;
	bhat[i,.]	= vecr(bhati)';
	resid	= yi - xi*bhati;
	sse	= sse + resid'resid;
endfor;
s2hat	= sse/ntot;		@ MLE of sigma @

ztz	= zdata'zdata;



/*
********************************************************************
*	Initialize Priors 
********************************************************************
*/

@ Prior for Sigma is IW(sf0,sg0)  @
sf0	= nyvar+4; sfn = sf0+ntot;
sg0	= eye(nyvar);
sg0i = invpd(sg0);


@ 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*vec(u0);          @ used in updating theta @

@ Lambda^{-1} is W_nparb(f0,g0 ) @
@ f0 = prior df, g0 = prior scale matrix @
f0 		= nparb+2;  f0n = f0 + nsub;
g0i 	= eye(nparb);   @ g0^{-1} @

/*
*******************************************************************
*	Initialize MCMC
******************************************************************
*/

beta 		= zeros(nsub,nparb);
sigma		= eye(nyvar);
sigmai	= invpd(sigma);

theta		= zeros(rankz,nparb);
lambda	= eye(nparb);
lambdai	= invpd(lambda);

@ Define data structures for saving iterates & computing posterior means & std @
betam	= zeros(nsub,nparb);	@ posterior mean of beta 	@
betas	= zeros(nsub,nparb);	@ posterior std of beta 	@
c		= nyvar*(nyvar+1)/2;
sigmag	= zeros(smcmc,c);
thetag	= zeros(smcmc,thdim);
c	= nparb*(nparb+1)/2;
lambdag = zeros(smcmc,c);	@ save unique elements of lambda @


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

for i1 (1,smcmc,1);	imcmc = i1;		@ Save smcmc iterations 			@
	for i2 (1,skip,1); jmcmc = i2;		@ Save every skip iterations 	@
		call gethbmreg2;
	endfor;
	sigmag[imcmc,.]	= vech(sigma)';
	thetag[imcmc,.]	= vecr(theta)';
	betam			= betam + beta;
	betas			= betas + beta^2;
	lambdag[imcmc,.]	= vech(lambda)'; @ vech gets unique elements of symmetric matrix @
endfor;


/*
******************************************************************
*	Compute Posterior Means and STD
******************************************************************
*/
@ Compute Posterior Means @
betam			= betam/smcmc;
sigmam		= xpnd(meanc(sigmag));
thetam		= meanc(thetag);
thetam		= reshape(thetam,rankz,nparb);
lambdam		= xpnd(meanc(lambdag));	@ xpnd is opposite of vech @

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

@ Predict yi @
yhat	= zeros(ntot,nyvar);
br		= zeros(nyvar,1);		@ multiple correlation for each y @	
estd	= br;					@ error std			@
sse		= zeros(nyvar,nyvar);
for i0 (1, nsub, 1); i = i0;
	xi	= xdata[iptxy[i,1]:iptxy[i,2],.];
	yi	= ydata[iptxy[i,1]:iptxy[i,2],.];
	betai	= betam[i,.]';
	bmati	= reshape(betai,rankx,nyvar);
	yhati	= xi*bmati;
	yhat[iptxy[i,1]:iptxy[i,2],.] = yhati;
	resid		= yi - yhati;
	sse		= sse + resid'resid;
endfor;
estd = sqrt(  diag( sse/ntot)  );

for fj (1, nyvar,1); j = fj;
	cr = corrx(ydata[.,j]~yhat[.,j]);
	br[j]	= cr[1,2];
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);
title("MLE & HB for " $+ xnames[2]$+ " versus " $+ xnames[1]);
_plctrl = -1;
xy(bhat[.,1]~betam[.,1],bhat[.,2]~betam[.,2]);
graphset;			@ Return to default settings @


end;

/*
****************************************************************
* GETHBMREG2
*	Does one iteration of the HB regression model.
*	INPUT
*		Global Variables
*	OUTPUT
*		Global Variables
****************************************************************
*/
PROC (0) = gethbmreg2;
local vibn, vibn12, ebin, yhat, sse, sgn, sgni,
 resid, gni, gn,  w, i0, i, limub, yvi, bmati, bi;

	limub = zdata*theta*lambdai;		@ Used in posterior mean of beta @
	/*
	***********************************************************
	* Generate beta
	* beta_i is N(mbin, vbn)
	* vbn 	= ( X_i'X_i.*.Sigma^{-1}  + Lambda^{-1} }^{-1}
	* mbin	= vbn*( (X_i'.*.Sigma^{-1})Y_i + Lambda^{-1}*Theta*Z_i)
	**********************************************************
	*/
	sse	= zeros(nyvar,nyvar);
	for i0 (1,nsub,1); i = i0;
		xi 		= xdata[iptxy[i,1]:iptxy[i,2],.];
		yi 		= ydata[iptxy[i,1]:iptxy[i,2],.];
		yvi		= vecr(yi);
		xitxi 	= xtx[iptxt[i,1]:iptxt[i,2],.];
		vibn	= xitxi.*.sigmai + lambdai;
		vibn12	= chol(vibn);
		ebin	= ((xi').*.sigmai)*yvi + limub[i,.]';
		bi 	= cholsol(ebin + vibn12'rndn(nparb,1), vibn12);
		beta[i,.]	= bi';
		bmati	= reshape(bi,rankx,nyvar);
		yhati	= xi*bmati;
		resid	= yi - yhati;
		sse	= sse + resid'resid;
	endfor;

	/*
	***************************************************************
	* Generate sigma
	* 
	* sn = s0 + sum_{i=1}^{nsub} (Y_i - X*beta_i)'(Y_i - X*beta_i)
	****************************************************************
	*/
	sgni	= sg0i + sse;
	sgn	= invpd(sgni);

	{sigmai, sigma} 	= wishart(nyvar,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