%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multinomial Probit Data Generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% data generation with E=I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assume "0" utility for first option

n=1000; % number of observations
M=3; %number of original options

bet=[1 -1]';
X1=[0.1+2*randn(n,1) 0.5+1.5*randn(n,1)];
X2=[0.1+2*randn(n,1) 0.5+1.5*randn(n,1)];
X3=[0.1+2*randn(n,1) 0.5+1.5*randn(n,1)];

% VCOV for error differences
Etrue=eye(M);
mu=zeros(M,1);
e=mvnrnd(mu,Etrue,n); % will be n by M
e1=e(:,1);
e2=e(:,2);
e3=e(:,3);

U1=X1*bet+e1;
U2=X2*bet+e2;
U3=X3*bet+e3;

U=reshape([U1';U2';U3'],M*n,1);
Ucell=mat2cell(U,M*ones(n,1),1);

y=zeros(n,1);
for i=1:n
    Ui=Ucell{i};
    [j jj]=max(Ui);
    if jj==1
        y(i)=1;
    elseif jj==2
        y(i)=2;
    else
        y(i)=3;
    end
end

l1=length(find(y==1));
l2=length(find(y==2));
l3=length(find(y==3));

[l1 l2 l3]
%stop at 323, 350, 327

save c:\klaus\aaec6564\mlab\worksp\mod12_MNP_I_data y X1 X2 X3;
