function [Y X Z] = DataGeneration(model,n,m,psi,gamma,parameter)
% Generate correlated fixed effects with mean zero and covar matrix omega
%rho =-.25; mu = zeros(n,2); omega = [1 rho;rho, 1]; log_fixed_effects = mvnrnd(mu,omega); 
%alpha = exp(log_fixed_effects(:,1)-1/2); 
%gamma = exp(log_fixed_effects(:,2)-1/2);

% generate covariates

 z = randn(n,m)-1/2; u = randn(n,m)-1/2; x = z+gamma*u; 
 X{1} = x; Z{1} = z;


% generate conditional-mean function
varphi = exp(x*psi).*exp(u);

% generate outcome variable
if strcmp(model,'Poisson')==1,
   Y = random('Poisson',varphi); % poisson with mean = var = varphi
end
if strcmp(model,'Negbin2')==1,
    Y = random('nbin',parameter,parameter./(parameter+varphi)); % negbin2  with mean = varphi and var = mean(1+parameter*varphi)
end
if strcmp(model,'Lnormal')==1,
    if parameter==1, variance = 1             ; end % log-normal errors are standard normal
    if parameter==2, variance =   1./varphi   ; end % Poisson-type variance
    if parameter==3, variance = 1+1./varphi   ; end % Negbin2-type variance with theta = 1;
    if parameter==4, variance =   1./varphi.^2; end % outcomes are homoskedastic outcomes
    
    sigma2 = log(1+variance); mu = -sigma2/2; e = mu+sqrt(sigma2).*randn(n,m); epsilon = exp(e); 
    Y = varphi.*epsilon; % log-normal model  
    
%    Y = round(Y); % introduce rounding error
    
end
if strcmp(model,'Mixture')==1,
    M = random('nbin',parameter,parameter./(parameter+varphi)); Y = zeros(n,m);
    for i = 1:n,
        for j = 1:m,
            Q = chi2rnd(1,M(i,j),1); 
            Y(i,j) = sum(Q);
        end
    end
    
    
end


