/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GHBREG2.GSS
*		Generats data for HB Regression Model
*		Y_i 	= X_i*beta_i + epsilon_i
*		beta_i	= Theta'Z_i + delta_i
*		B		= Z*Theta + D
*		epsilon_i is N(0,sigma2*I)
*		delta_i	is N(0,Lambda)
*--->	Different Design Matrices
*****************************************************************
*/
new;
nobs	= 100;							@ Number of subjects 					@
yrows	= 5 + floor(10*rndu(nobs,1));	@ Number of observations per subject		@
ntot	= sumc(yrows);
sigmat	= 5;							@ True error STD					@

lbd12	= { 
			5	.5	-1	.1,
			0	4	-2	0,
			0	0	3	2,
			0	0	0	1
		};
lbd12	= lbd12/2;
lambdat	= lbd12'lbd12;

thetat	= {
			2	-1	-3	4,
			-1	0	2	-3,
			3	2	1	0
		};
rankx	= rows(lambdat);
rankz	= rows(thetat);
@ Get pointer into stacked xy matrices @
b	= cumsumc(yrows);
a	= 1|(1+b[1:nobs-1]);
iptxy	= a~b;

xdim	= rankx - 1;
zdim	= rankz - 1;
a	= seqa(1,1,xdim);
xnames	= 0 $+ "X " $+ ftocv(a,1,0);
a	= seqa(1,1,zdim);
znames	= 0 $+ "Z " $+ ftocv(a,1,0);
xynames	= xnames|"Y ";

xdata	= rndn(ntot,xdim);
xmat	= ones(ntot,1)~xdata;
zdata	= rndn(nobs,zdim);
zmat	= ones(nobs,1)~zdata;
betat	= zmat*thetat	+ rndn(nobs,rankx)*lbd12;
ydata	= zeros(ntot,1);
@ Generate and store y data @
for i0 (1,nobs,1); i = i0;
	xi	= xmat[iptxy[i,1]:iptxy[i,2],.];
	yi	= xi*(betat[i,.]') + sigmat*rndn(yrows[i],1);
	ydata[iptxy[i,1]:iptxy[i,2],.] = yi;
endfor;
xydata	= xdata~ydata;
/*
**************************************************************************
* Create & Read a Gauss file.  f1 is the file handle.  			
**************************************************************************
*/
create f1 = xydata with ^xynames, 0, 8;
if writer(f1,xydata) /= rows(xydata);
		errorlog "Conversion of XYDATA to Gauss File did not work";
endif;
closeall f1;
create f1 = zdata with ^znames, 0, 8;
if writer(f1,zdata) /= rows(zdata);
	errorlog "Conversion of ZDATA to GAUSS File did not work";
endif;
closeall f1;
save yrows		= yrows;
save sigmat		= sigmat;
save betat		= betat;
save thetat		= thetat;
save lambdat	= lambdat;