 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data for treatment effect model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% data generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=1000; %sample size
m=2; % number of correlated equations
Etrue=[1 0.75;0.75 1.8]; %error varcov, with a "healthy" covariance term
mu=zeros(m,1); % error mean

beta1true=[1 2]'; %coefficients for treatment equation
beta2true=[6 0.5 1]';%coefficients for outcome equation

X1=[ones(n,1) -0.5+1.3*randn(n,1)];
e=mvnrnd(mu,Etrue,n); % will be n by 2
y1=X1*beta1true+e(:,1);
f=find(y1>0);
fi=length(f)/n % stop at about 50% treated

treat=zeros(n,1);
treat(f)=1;
X2=[ones(n,1) 4+2*randn(n,1) treat];
y2=X2*beta2true+e(:,2);

%also, all of y2>0 - convenient if we want to take logs down the line.

save c:\klaus\AAEC6564\mlab\worksp\treat_data;

