/*
*********************************************************************************
*	ARE1.GSS
*	linear regression model + AR Errors
*
*	Y = X*beta + epsilon
*	epsilon_t = rho*epsilon_{t-1} + zeta_t
*	zeta_t iid N(0,sigma^2);
*
*	Priors
*
*
*	beta is N(b0,B0)
*	sigma^2 is IG(r0/2,s0/2)
*	rho is uniform on [-1,1];
*	
*	Issue the command:
*	library pgraph, plbam;
*	before running.
***********************************************************************************
*/
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		= 1000;		@ number of iterations to save for analysis 				@
skip		= 1;		@ Save every skip iterations								@
nblow		= 100;		@ 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)];	



xdim		= cols(xydata)-1;			@ number of x variables 				@
										@ cols(x) = # of columns of x			@
rankx		= xdim+1;					@ number of beta parameters				@
xnames		= "Const"|xynames[1:xdim];

nobs		= rows(xydata);				@ number of observations				@
xmat		= xydata[.,1:xdim];	
xdata		= ones(nobs,1)~xmat;
										@ 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		= 10*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) 			@
smati = 0;
@ 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 + 1;			@ Posterior shape parameters: Add one for the epsilon_0 @	

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

rho		= 0;



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


/*
********************************************************************
*	Do MCMC
********************************************************************
*/
@ Do the initial transition period @
icount = 0;
for i1 (1,nblow,1);	imcmc = i1;
	call getall;
	icount = icount + 1;

endfor;

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

endfor;


/*
**************************************************************************************
*  Analyze Results
**************************************************************************************
*/
@ Compute posterior means and std from MCMC iterations @
sigm	= meanc(sigmag);				
sigstd	= stdc(sigmag);
rhom	= meanc(rhog);
rhos	= stdc(rhog);
betam	= meanc(betag);
betas	= stdc(betag);	
yhat	= xdata*betam;
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 @






	

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;

/*
*******************************************************************************************
*  GETALL
*	Generate one MCMC random deviate of beta and sigma
*	No input or output.
*	Data & values are passed through global variables
*******************************************************************************************
*/
PROC (0) = getall;
local  xstar, resid, rstar, ebn, vibn, vibn12, es, cor, cori,
subc, 
sse, sn, cmat, ccmat, subcmat, xivi, xivi12, exim, sntau, rv, rm, 
v2, rhotop, rhobot, arho, brho;


/*
****************************************************************
* Sigma = sigma^2/(1-rho^2)*Cor 
*	where Cor is Topelitz with i,j entry rho^{|i-j|}
*  Compute Sigma^{-1}
*  Cor^{-1} is tri-diagonal
*  cori = (1/1-rho^2)*C. 
*  C is tri-digaonal with
*  (1, 1-3*rho^2, ..., 1-3*rho^2, 1) on main diagonal.
*  -rho on minor diagonals.
****************************************************************
*/

subc	= -rho*eye(nobs-1)~zeros(nobs-1,1);
subc	= zeros(1,nobs)|subc;
cori	= (1+rho^2)*eye(nobs) + subc + subc';
cori[1,1] = 1;
cori[nobs,nobs] = 1;
cori	= cori/(1-rho^2);
smati	= (1-rho^2)*cori/sigma2;
/*
*********************************************************************
* cori12 and sigam^{-1/2}
********************************************************************
*
cori12 = eye(nobs)/sqrt(1-rho^2);
cori12[nobs,nobs] = 1;
subc	= -rho*eye(nobs-1)/sqrt(1-rho^2);
subc	= zeros(nobs-1,1)~subc;
subc	= subc|zeros(1,nobs);
cori12	= cori12 + subc;
detcorn	= det(cor);
detcora = (1-rho^2)^(nobs-1);