/*
******************************************************************
*   (C) Copyright 1999, Peter Lenk. All Rights Reserved.
*  GHBMREG2.GSS
*		Generats data for HB Regression Model
*		Y_i 	= X_i*B_i + U_i
*
*		Y_i is a n_i by nyvar matrix.
*		X_i is a n_i by rankx design matrix
*		B_i is rankx by nyvar matrix of regression coefficients.
*
*		Define:
*			yv_i = vec(Y_i')
*			beta_i = vec(B_i')
*			epsilon_i = vec(U_i') is N(0,I(n_i).*.Sigma)
*		Then
*			yv_i = (X_i.*.I(m))*beta_i + epsilon_i
*		Put
*			beta_i	= Theta'Z_i + delta_i
*			delta_i	is N(0,Lambda)
*
*****************************************************************
*/
new;
nobs	= 150;							@ Number of subjects 					@
yrows	= 10 + 	ceil(10*rndu(nobs,1));	@ Number of observations per subject		@
ntot	= sumc(yrows);
nyvar	= 3;
rankx	= 2;							@ Rank of design matrix (includes intercept) @
nparb	= nyvar*rankx;
sig12	= {
			1   -4    0,
			0   2     0,
			0	0	   3
};

sigmat	= sig12'sig12;							@ True error STD					@
lbd12	= {	1	-2	-3	3	0	0,
			0	2	-1	1	0	0,
			0	0	3	1	0	0,
			0	0	0	4	0	0,
			0	0	0	0	5	-3,
			0	0	0	0	0	6
			};
lbd12 	= lbd12/2;
lambdat	= lbd12'lbd12;

thetat	= {
			2	-1	-3	4	1	-3,
			-1	0	-2	3	1	-2,
			3	2	1	0	2	-1
		};
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);
ynames	= 0 $+ "Y" $+ ftocv(seqa(1,1,nyvar),1,0);
a	= seqa(1,1,zdim);
znames	= 0 $+ "Z " $+ ftocv(a,1,0);
xynames	= xnames|ynames;

xdata	= rndn(ntot,xdim);
xmat	= ones(ntot,1)~xdata;
zdata	= rndn(nobs,zdim);
zmat	= ones(nobs,1)~zdata;
betat	= zmat*thetat + rndn(nobs,nparb)*lbd12;
ydata	= zeros(ntot,nyvar);
@ Generate and store y data @
for i0 (1,nobs,1); i = i0;
	xi	= xmat[iptxy[i,1]:iptxy[i,2],.];
	betai	= betat[i,.]';
	bmati	= reshape(betai,rankx,nyvar);
	yi		= xi*bmati + rndn(yrows[i],nyvar)*sig12;
	ydata[iptxy[i,1]:iptxy[i,2],.] = yi;
endfor;
xydata	= xdata~ydata;
@ Save the data @
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;
save nyvar 		= nyvar;