linreg1.gss

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

@ shape parameters alpha and scale parameter beta = 1				@
@ mean = alpha and variance = alpha								@
sigma	= sqrt(sigma2);			@ Error STD						@
endp;


/*
*****************************************************************************************************
*  OUTPUTANAL
*	Does analysis of output and save some results 
****************************************************************************************************
*/
PROC (0) = outputanal;

local  bout, sout, fmtn1, fmtn2, fmts1, fmts2;
format 10,5;							@ Default 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 LINREG1A.GAS";
print "Bayesian analysis of linear regression model using MCMC.";
print "Ouput file: " getpath(outfile);			@ File assigned to file handle outfile @
datestr(date);									@ Print the current data				@
print;
print "Model";
print "Y = X*beta + epsilon";
print "Number of observations          = " nobs;
print "Number of independent variables = " xdim "(excluding the intercept).";
print "Summary Statistics";
call sumstats(xynames,xydata,fmts1,fmts2,fmtn1);	@ Print summary statistics @
print;
print;
print "-----------------------------------------------------------------------------------";
print "MLE Analysis";
print;
print "Number of observations        = " nobs;
print "Number of dependent variables = " xdim;
print;
print "R-Squared     = " rsquare;
print "Multiple R    = " sqrt(rsquare);
print "MLE Error STD = " sighat;
print;
print "Estimated Regression Coefficents";
sout	= {"Variable"  "MLE"  "StdError"};
call outitle(sout,fmts1,fmts2);
bout 	= xnames~bhat~bstd;
call outmat(bout,fmts1,fmtn1);

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 observations        = " nobs;
print "Number of dependent variables = " xdim " (excluding the intercept)";
@ Compute posterior means and std @
print;
print "Bayes R-Square   =   " cy2;
print "Bayes Multiple R =   " cy[1,2];
print;

print "Error Standard Deviation";
print "Posterior mean of sigma = " sigm;
print "Posterior STD  of sigma = " sigstd;
print;
print "Regression Coefficients";
sout	= {"Variable" "PostMean" "PostSTD"};
call outitle(sout,fmts1,fmts2);

bout	= xnames~bm~bs;
call outmat(bout,fmts1,fmtn1);
output off;
closeall;
endp;

/*
*****************************************************************
* PDFSIGMA
*	Computes posterior density if sigma over a grid.
*	INPUT
*		Uses global variables
*	OUTPUT
*		sgrid	= grid for the density
*		fs		= posterior pdf on sgrid
****************************************************************
*/
PROC (2) = pdfsigma;
local smax, smin, ms, sdelta, sgrid, scon, fs, s2grid, lnsgrid, i0, i,
sse, sn;
@ compute posterior density of sigma over a grid @
smax 	= maxc(sigmag);
		@ maxc(x) returns maximum of x 			@
smin	= minc(sigmag);
		@ minc(x) returns minimum of x				@
@ Get grid for plotting		@
ms		= 100;
sdelta	= (smax-smin)/ms;
sgrid	= seqa(smin, sdelta, ms+1);
scon	= ln(2) - lnfact(1+rn/2);
fs		= zeros(ms+1,1);
s2grid	= sgrid.*sgrid;
lnsgrid	= ln(sgrid);
for i0 (1,smcmc,1); i = i0;
	sse	= ydata - xdata*(betag[i,.]');
	sse = sse'sse;
	sn	= s0 + sse;
	fs  = fs + exp(scon + rn*ln(sn/2)/2 + (0.5 - rn)*lnsgrid - sn/(2*s2grid));
endfor;
fs = fs/smcmc;
fs	= fs/intsim(fs,sdelta);		@ intsim is simpson's integration @
retp(sgrid,fs);
endp;

/*
********************************************************************************
* PDFBETA
*	Computes marginal posterior density of beta over a grid.
*	These are marginal, univariate densities.
*	INPUT
*		Uses globals.
*	OUTPUT
*		bgrid	= each row corresponds to grid for beta_j
*		fb		= each row corresponds to marginal density for beta_j
********************************************************************************
*/
PROC (2) = pdfbeta;
local mb, bgrid, bmax, bmin, bdelta, fj, j, fb, bcon, sigma, sigma2, 
vibn, vbn, ebn, vbnd;
@ Plot posterior density for beta @
mb 		= 100;
bgrid 	= zeros(mb+1,rankx);
bmax 	= bm + 5*bs;
bmin	= bm - 5*bs;
bdelta	= (bmax - bmin)/mb;
for fj (1,rankx,1) ; j = fj;
	bgrid[.,j] = seqa(bmin[j],bdelta[j], mb+1);
endfor;
	
fb		= zeros(mb+1,rankx);
bcon	= -ln(2*pi)/2;
for fj (1,smcmc, 1); j = fj;
	sigma	= sigmag[j];
	sigma2	= sigma^2;
	vibn	= (xtx/sigma2 + vib0);
	vbn		= invpd(vibn);					@ Posterior variance given beta 	@
	ebn		= vbn*(xty/sigma2 + vieb0);		@ Posterior mean given beta 		@
	vbnd	= diag(vbn);					@ diagonal elements 				@
	fb		= fb + exp(bcon - 0.5*ln(vbnd)' - ((bgrid - ebn')^2)./(2*vbnd') );
endfor;
fb = fb/smcmc;
fb = fb./(intsim(fb,bdelta)');

retp(bgrid, fb);
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;