%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SSVS model for fishing data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;  % start stop watch

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

load c:\klaus\AAEC6564\mlab\worksp\fish_define_dataspace D2 Z2;

% Open log file 
[fid]=fopen('c:\Klaus\AAEC6564\mlab\logs\mod8_SSVS_fishing.txt','w');
if fid==-1;
   warning('File could not be opened');
   return;
else;
   disp('File opened successfully');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data Set D2: cold or warm, RIVER fishing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=D2;
Z=Z2;
n=size(data,1);

y=data(:,9);

% define and standardize matrix with main effects
Xm=[ones(n,1) data(:,8) data(:,10)]; 
% constant, 
% log(catch)
% log(income)

%contents of Z
% dummy for "warm water"
% interaction of dummy with log(catch)
% interaction of dummy with log(income)

X=[Xm Z];
k1=size(Xm,2);
k2=size(Z,2);
k=k1+k2;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=25000; % burn-in
r2=10000; % keepers
R = r1+r2;

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

% elements for beta
p=0.5; %ssvs prior probability that a slope coeff comes from a distr. with
% a large variance
t2=0.01;
c2=10/t2;

mu0=zeros(k,1); % prior mean for all parameters
Vx=ones(k1,1)*(c2*t2); % prior variances for elements of X

%Create starting Variance matrix for beta
gamvec=binornd(1,p,k2,1); %draw a 0/1 value for each interaction coefficient
int=gamvec.*c2*t2 + (1-gamvec).*t2; %set starting variance for each coefficient
dvec=[Vx;int]; %combine prior variances for main effects and interactions
V0=diag(dvec); % k by k
%This prior variance will be updated throughout the sampling process

betadraw=bols; %use OLS values as starting draws for Gibbs Sampler

% elements for sig2
sig2draw=s2; % use OLS variance as starting draw for Gibbs Sampler
v0=1/2;
tau0=1/2; % diffuse prior shape and scale

[betamat,sig2mat,gammat]=...
    gs_SSVS_fishing(X,y,n,k1,k2,k,r1,r2,mu0,Vx,V0,t2,c2,p,betadraw,v0,tau0,sig2draw);

'GS done'

% Run diagnostics on beta and sigma
allmat=[betamat;sig2mat];
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');

%sig2 stuff
out=kdiag(k+1,:)';
fprintf(fid,'Output table for sig2 \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');

%gamma stuff
out=[mean(gammat,2)]; %inclusion probabilities for each interaction term
fprintf(fid,'Inclusion probabilities \n\n');
fprintf(fid,'%6.3f\n',out);
fprintf(fid,'\n');

save c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_fishing betamat sig2mat gammat X y Z;

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;


