%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic Probit model using Fair's (1978) affairs data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;  % start stop watch

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


% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
% load c:\klaus\AAEC6564\mlab\worksp\Fair_data.txt;
% data=Fair_data;
% clear Fair_data;
% save c:\klaus\AAEC6564\mlab\worksp\Fair_data data;
% 
% 'ok'
% break
load c:\klaus\AAEC6564\mlab\worksp\Fair_data;


% Contents of Data (601 observations)

%1       affair    1=yes, 0=no
%2       male      1=male, 0=female
%3       years married
%4       kids      1=yes, 0=no (from this marriage)
%5       religious 1=yes, 0=no 
%6       education (years of completed schooling)
%7       happy     1=yes, 0=no (marriage = happier than the average
                                %marriage)
                                
y=data(:,1);
n=length(y);
X=[ones(n,1) data(:,2:end)];

f0=find(y==0);
f1=find(y==1);

l0=length(f0);
l1=length(f1);
                                 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=15000; % burn-ins
r2=10000; % keepers
R=r1+r2;
k=size(X,2);

% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(n-k);

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

load c:\klaus\AAEC6564\mlab\worksp\mod5_probit_Fair betamat;

% 'ok'
% return

[logpmarg,lp,llf,lpost]...
    =gs_probit_chib(X,y,k,n,r1,r2,mu0,V0,f0,f1,betamat);

'GS done'


% Open log file 
[fid]=fopen('c:\Klaus\AAEC6564\mlab\logs\mod5_probit_Fair_chib.txt','w');
if fid==-1;
   warning('File could not be opened');
else;
   disp('File opened successfully');
end;

fprintf(fid,'logpmarg=\t%6.4f \n',logpmarg);
fprintf(fid,'lp=\t%6.4f \n',lp);
fprintf(fid,'llf=\t%6.4f \n',llf);
fprintf(fid,'lpost=\t%6.4f \n',lpost);


save c:\klaus\AAEC6564\mlab\worksp\mod5_probit_Fair_chib logpmarg lp llf lpost;

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