%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bayesian Multinomial Probit model with Identity COV matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set random number generator, start stop watch, open output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch

% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
load c:\klaus\aaec6564\mlab\worksp\MNP_I_data;
M=3;
n=length(y);

% create an Xdiff matrix with "losers" subtracted from the winner
% Naturally, this is individual-specific
% We need also adjust the VCOV of the errors accordingly

k=size(X2,2);
Xcell=cell(n,1);
Ecell=cell(n,1);
D=-eye(M);
        
for i=1:n
    j=y(i); %index of winning option
    if j==1
        D1=D;
        D1(:,1)=1;
        D1(1,:)=[];
        int1=X1(i,:)-X2(i,:);
        int2=X1(i,:)-X3(i,:);
        Xcell{i}=[int1;int2];
        Ecell{i}=D1*D1';
    elseif j==2
        D2=D;
        D2(:,2)=1;
        D2(2,:)=[];
        int1=X2(i,:)-X1(i,:);
        int2=X2(i,:)-X3(i,:);
        Xcell{i}=[int1;int2];
        Ecell{i}=D2*D2';
    else
        D3=D;
        D3(:,3)=ones(M,1);
        D3(3,:)=[];
        int1=X3(i,:)-X1(i,:);
        int2=X3(i,:)-X2(i,:);
        Xcell{i}=[int1;int2];
        Ecell{i}=D3*D3';
    end
end
X=Xcell;
E=Ecell;
%note: y is an n by 1 vector indexing the chosen alternative for each
%person; X is M-1*n by k

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=1000; % burn-in
r2=1000; % keepers
R=r1+r2;
rg=30; %iterations for Gibbs-within-Gibbs

% elements for beta
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta

sig=1; %scalar for E - needed on GS

[betamat]=gs_MNP_I(X,y,k,M,n,r1,r2,mu0,V0,rg,E,sig);
'GS done'


[fid]=fopen('c:\klaus\aaec6564\mlab\logs\mod12_MNP_I.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% put all draws together & run diagnostics
allmat=[betamat];
kdiag=klausdiagnostics_greater0(allmat);

fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'\n');

% beta stuff
out=kdiag(1:k,:)';
fprintf(fid,'Output table for betas \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


save c:\klaus\aaec6564\mlab\worksp\MNP_I betamat;

finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);


st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;










