linreg1.gss

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

/*
*********************************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
* LINREG*.GSS
*	Does basic, linear regression model:
*
*	Y = X*beta + epsilon
*	epsilon is N(0,sigma^2I)
*	beta is N(b0,B0)
*	sigma^2 is IG(r0/2,s0/2)
*
*	library pgraph, plbam;
***********************************************************************************
*/
new;
inxy 		= "xydata";			@ Name of Gauss file with X,Y data				@
outfile	= "results1.dat";	@ Specify output file for saving results 			@
						@ outfile is a string variable that contains a file name @
/*
********************************************************************
*	Initialize parameters for MCMC 
********************************************************************
*/
smcmc		= 100;		@ number of iterations to save for analysis 				@
skip		= 1;		@ Save every skip iterations								@
nblow		= 1;		@ Initial transition iterations 							@
nmcmc		= nblow + skip*smcmc;	@ total number of iterations					@

/*
********************************************************************
*	Get data
********************************************************************
*/
@ 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 @
ynames		= xynames[rows(xynames)];	
xnames		= "Constant"|xynames[1:rows(xynames)-1];

xdim		= cols(xydata)-1;			@ number of x variables 				@
										@ cols(x) = # of columns of x			@
rankx		= xdim + 1;					@ number of beta parameters				@
nobs		= rows(xydata);				@ number of observations				@
xdata		= ones(nobs,1)~xydata[.,1:xdim];	
										@ design matrix, includes an intercept 	@
										@ ~ sticks two matrices side by side	@
										@ ones(i,j) = i x j matrix of ones		@
										@ x[i,j] is the i,j element of x		@
										@ x[.,j] is the column j of  x			@
										@ x[.,j1:j2] are columns j1 to j2 of x	@
										@ x[i,.] is row i of x					@
										@ x[i1:i2,.] are rows i1 to i2 of x		@
ydata		= xydata[.,cols(xydata)];

@ Sufficient statistics used in MCMC @
xtx			= xdata'xdata;
xty			= xdata'ydata;

@ Get MLE @
{bhat, bstd, sighat, rsquare} = regmle(xdata,ydata);

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

@ Beta is N(eb0, vb0) @
eb0		= zeros(rankx,1);		@ prior mean							@
vb0		= 100*eye(rankx);		@ prior variance						@	
							@eye(m) = m x m identity matrix 			@
							@ Generally, use big variance			@
vib0		= invpd(vb0);			@ invpd(x) = inverse of positive definite x 	@
vieb0	= vib0*eb0;			@ used in full conditional of beta			@

@ Sigma2 is Inverted Gamma(r0/2, s0/2) 			@
@ E(1/sigma2) = r0/s0 and Var(1/sigma2) = 2*r0/s0^2 	@
@ Usually pick r0 and s0 small and positive 		@
r0 	= 2; s0 = 2;
rn	= r0 + nobs;			@ Posterior shape parameters @	

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

@ Initial parameters for MCMC 	@
@ Could initialize to MLE.		@
beta 	= zeros(rankx,1);		@ regression coefficients 	@
sigma	= 1;					@ error std				@
sigma2	= sigma*sigma;		@ error variance			@

@ Define matrics for saving MCMC iterations @
betag	= zeros(smcmc, rankx);
sigmag	= zeros(smcmc,1);

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

for i1 (1,smcmc,1);	imcmc = i1;		@ Save smcmc iterations 			@
	for i2 (1,skip,1); jmcmc = i2;		@ Save every skip iterations 	@
		call getbetasigma;
	endfor;
	betag[imcmc,.] 	= beta';		@ Save regression coefficients 		@
	sigmag[imcmc,.] 	= sigma;		@ Save error std					@
endfor;


/*
**************************************************************************************
*  Analyze Results
**************************************************************************************
*/
@ Compute posterior means and std from MCMC iterations @
sigm	= meanc(sigmag);				
sigstd	= stdc(sigmag);
bm		= meanc(betag);
bs		= stdc(betag);	
yhat	= xdata*bm;
cy		= corrx(ydata~yhat);			@ Correlation between Y and Yhat @
cy2		= cy[1,2]*cy[1,2];
call outputanal;
@ Plot saved iterations against iterations number @
t 	= seqa(nblow+skip,skip,smcmc);		@ saved iteration number @
title("Error STD versus Iteration");
xy(t,sigmag);
title("Coefficients versus Iteration");
xy(t,betag);


@ Compute the posterior density of sigma over a grid @
{sgrid,fs} = pdfsigma;

@ Compute the marginal posterior density of beta_j over a grid @
{bgrid, fb} = pdfbeta;
title("Posterior Density of the Error STD");
xy(sgrid,fs);


for fj (1,rankx,1); j = fj;
	title("Posterior Density of Coefficient for " $+ xnames[j]);
	xy(bgrid[.,j],fb[.,j]);
endfor;
	
graphset;			@ Return to default graphs @
end;

/*
*****************************************************************************************
*   REGMLE
*	Compute MLE for simple regression
*	INPUT:
*		XDATA	= Design matrix
*		YDATA	= Dependent Variable
*	OUTPUT:
*		BHAT	= MLE for regression coefficients
*		BSTD	= STD Error of beta
*		Sighat	= Error STD
*		Rsqure	= R-Squared
***********************************************************************************************
*/
PROC (4) = regmle(x,y);
local xtx, xtxi, xty, b, yhat, resid, sse, s, sst, r2, bstd;
	xtx		= x'x;
	xtxi		= invpd(xtx);				@ Inverse of xtx			@
	xty		= x'y;
	b		= xtxi*xty;				@ MLE of beta 			@
	yhat		= x*b;					@ Predicted y values 		@
	resid	= y - yhat;				@ Residuals				@
	sse		= resid'resid;				@Sums of Squares Error 	@
	s		= sqrt(sse/nobs);			@ Error STD				@
	sst		= y - meanc(y);
	sst		= sst'sst;					@ SS Total				@
	r2		= 1 - sse/sst;				@ R-squared				@
	bstd		= s*sqrt(diag(xtxi));
retp(b,bstd, s,r2);						@ Return to main program	@
endp;

/*
*******************************************************************************************
*  GETBETASIGMA
*	Generate one MCMC random deviate of beta and sigma
*	No input or output.
*	Data & values are passed through global variables
*******************************************************************************************
*/
PROC (0) = getbetasigma;
local  ebn, vibn, vibn12, vbn, vbn12, sse, sn;
/*
*******************************************************************************************************
*  Generate beta:
*	beta is N(ebn, vbn);
*	vbn = (X'X/sigma2 + Vb0^{-1} )^{-1}		is posterior variance 
*	ebn = vbn*(X'Y/sigma2 + Vb0^{-1}*eb0)	is posteriro mean
*********************************************************************************************************
*/
ebn		= xty/sigma2 + vieb0;		@ Posterior variance * ebn = posterior mean of beta  	@
vibn		= xtx/sigma2 + vib0;		@ Inverse of posterior variance of beta			 	@
@vibn12	= chol(vibn);@				@ vibn12'vibn12 = vibn							@
vbn		= invpd(vibn);
vbn12	= chol(vbn);
ebn		= vbn*ebn;
beta		= ebn + vbn12'rndn(rankx,1);
@beta		= cholsol(ebn+vibn12'rndn(rankx,1), vibn12);@
		@ cholsol is efficient method of solving linear equations if you have chol decompostion 	@
		@ Suppose C'C = A															@
		@ Need to find x to solve A*x = b												@
		@ x = cholsol(b, C) does it 													@

/*
***********************************************************************************************************
*  Generate sigma2
*	sigma2 is IG(rn/2, sn/2)
*	rn 	= r0 + nobs
*	sn 	= s0 + (y - X*beta)'(y - X*beta) 
************************************************************************************************************
*/
sse		= ydata - xdata*beta;
sse		= sse'sse;
sn		= s0 + sse;
sigma2	= sn/(2*rndgam(1,1,rn/2));
@ rndgam(i,j,alpha) generates a matrix of gamma random variables with 	@