/*
*************************************************************************
*	GARE1.GSS
*	Generates data for linear regression + AR errors:
*
*	Y = x*beta + epsilon
*	epsilon_t = rho*epsilon_{t-1} + zeta_t
*	zeta_t is N(0,sigma^2)
*
*	Issue the command
*	library pgraph,plbam;
*	before running.
*************************************************************************
*/
new;
nobs	= 100;		@ number of observations 			@
betat	= {
5, 3
};
sigmat	= 3;
rhot	= 0.7;

sigma2	= sigmat*sigmat;
rankx = rows(betat);				@ rows(x) = number of rows for x 			@
xdim	= rankx - 1;

xnames	= 0 $+ "X" $+  ftocv(seqa(1,1,xdim),1,0);


ynames	= "Y";
xynames = xnames|ynames;

/*
*************************************************************************
* Generate Toplitz matrix for error covariance
**************************************************************************
*/
cor = zeros(nobs,nobs);
for fori (1,nobs,1); i = fori;
	for forj (1,nobs,1); j = forj;
		cor[i,j]	= rhot^(abs(i-j));
	endfor;
endfor;
/*
****************************************************************
*  Cor^{-1} is tri-diagonal
*  cori = (1/1-rhot^2)*C. 
*  C is tri-digaonal with
*  (1, 1-3*rhot^2, ..., 1-3*rhot^2, 1) on main diagonal.
*  -rhot on minor diagonals.
****************************************************************
*/

subc	= -rhot*eye(nobs-1)~zeros(nobs-1,1);
subc	= zeros(1,nobs)|subc;
cori	= (1+rhot^2)*eye(nobs) + subc + subc';
cori[1,1] = 1;
cori[nobs,nobs] = 1;
cori	= cori/(1-rhot^2);
/*
*********************************************************************
* cori12 
********************************************************************
*/
cori12 = eye(nobs)/sqrt(1-rhot^2);
cori12[nobs,nobs] = 1;
subc	= -rhot*eye(nobs-1)/sqrt(1-rhot^2);
subc	= zeros(nobs-1,1)~subc;
subc	= subc|zeros(1,nobs);
cori12	= cori12 + subc;
detcorn	= det(cor);
@ det of cor = (1-rhot^2)^(nobs-1) @

smat	= sigma2/(1-rhot^2)*cor;

@ Det of smat = sigma2^(nobs)/(1-rhot^2) @

smati	= (1-rhot^2)*cori/sigma2;
smati12 = sqrt(1-rhot^2)*cori12/sigmat;

smat12	= chol(smat);




xdata	= rndn(nobs,xdim);			@ rndn -> N(0,1) random numbers				@
xmat	= ones(nobs,1)~xdata;		@ design matrix								@
									@ ~  pastes two matrics side-by-side		@

@ Generate depenent observations @
errors	= smat12'rndn(nobs,1);
ydata	= xmat*betat + errors;

xydata	= xdata~ydata;


save sigmat	= sigmat;
save rhot	= rhot;
save betat	= betat;

/*
**************************************************************************
* Create a Gauss file.  f1 is the file handle.  
* The Gauss file will be called "XYDATA."       
* The column will be named by the strings in the character array xyname.  
* ^xyname means use the names in the character string.	
* 0, 8 gives double precision real numbers.			
**************************************************************************
*/
create f1 = xydata with ^xynames, 0, 8;
/*
**************************************************************************
* Next read data into the Gauss file by using the writer command.	
* f1 is the file handle defined in previous command.			
* xydata is the data matrix that we just created.			
* writer returns the number of rows read to f1.			
* If it is not rows(xydata), something bad happended.		
**************************************************************************
*/
if writer(f1,xydata) /= rows(xydata);
		errorlog "Conversion of XYDATA to Gauss File did not work";
endif;
closeall f1;


@ Plot Y versus X_i @
_plctrl 	= -1;					@ Plot symbols and no lines 				@

for fj (1,xdim,1); j = fj;
	title(ynames $+ " versus " $+ xnames[j]);
	xy(xdata[.,j],ydata);			@ Plot y versus x_i 						@
endfor;

graphset;							@ Return to default graphs					@