are1.gss

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

smat	= sigma2/(1-rho^2)*cor;
smatin	= invpd(smat);
smati12n = chol(smatin);
detsign = det(smat);
detsiga = sigma2^(nobs)/(1-rho^2);

smati12 = sqrt(1-rho^2)*cori12/sigma;
*********/


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



vibn	= xdata'smati*xdata + vib0;
vibn12	= chol(vibn);						@ vibn12'vibn12 = vibn							@
ebn		= xdata'smati*ydata + vieb0;
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) 
************************************************************************************************************
*/

resid	= ydata - xdata*beta;
sn		= s0 + (1-rho^2)*resid'cori*resid;

sigma2	= sn/(2*rndgam(1,1,rn/2));

sigma	= sqrt(sigma2);			@ Error STD					@





/*
*****************************************************
* Generate rho 
*	Use Damien's method to handle the determinant factor.
*   Truncated normal on -1 to 1
*
* (See program TEST2.GSS to test the following identities.)
* Smat = sigma2*cor/(1-rho^2)  where
* cor  is Topelitz with (i,j) element rho^{|i-j|}.
* Det|Smat|^{-1/2} = sqrt((1 - rho^2))/(sigma2^{n/2})

*****************************************************
*/

/*
****************************************************
* To include the determinate:
* V < sqrt(1 - rho^2)
* - sqrt( 1 - V^2 ) < rho < sqrt( 1 + V^2)
****************************************************
*/
v2		= (1 - rho^2)*rndu(1,1);
rhotop	= sqrt(1 - v2);
rhobot 	= - rhotop;

@ Keep rho between -1 and 1 @
rhobot 	= maxc(rhobot|-1);
rhotop	= minc(rhotop|1);


arho	= resid[2:nobs-1]'resid[2:nobs-1];
brho	= resid[2:nobs]'resid[1:nobs-1];
rv		= sigma2/arho;

rm		= brho/arho;

rho		= rndnab(1,1,rm,rv,rhobot,rhotop);		@ truncated normal from -1 to 1 @




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 ARE1.GSS";
print "Bayesian analysis of linear regression model using MCMC.";
print "AR Error Terms.";
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 "epsilon_t = rho*epsilon_{t-1} + z_t";
print "Var(Z_t) = sigma^2";
print;
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 "Error Correlation";
print "Posterior mean of rho = " rhom;
print "Posterior STD  of rho = " rhos;
print;
print "Regression Coefficients";
sout	= {"Variable" "PostMean" "PostSTD"};
call outitle(sout,fmts1,fmts2);

bout	= xnames~betam~betas;
call outmat(bout,fmts1,fmtn1);
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;