
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% POSTESTIMATION SCRIPT for fish stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% this is for data space 2 (both activities)

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;
load c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_fishing_Mprobs;


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

fprintf(fid,'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \n\n');
fprintf(fid,'data space = D2 \n\n');
fprintf(fid,'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \n\n');

pM=relcount;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare data
%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;

%Common parameters for all models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=10000; % burn-in
r2=10000; % keepers
R = r1+r2;

xp=[1 1.55 11]'; %settings for predictions

ypwmat=[];

for i=1:M    
    if pM(i)>0
        'Model ='
        i

        f=find(Mmat(i,:)==0);
        int=Z;
        int(:,f)=[];
        X=[Xm int];
        k=size(X,2);
      
        % generic OLS
        bols=inv(X'*X)*X'*y;
        res=y-X*bols;
        s2=(res'*res)/(n-k);
        betadraw=bols;

        % elements for beta
        mu0=zeros(k,1); % prior mean for all parameters
        V0=100*eye(k);

        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]=...
        gs_normal_independent(X,y,n,k,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw);    
       
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % GENERATE PREDICTIONS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for r=1:r2
            br=betamat(:,r);
            br=br(1:k1); %we're only using the baseline variables
            % (we're basically setting the interaction terms to zero
            sig2r=sig2mat(r);
            ypmat(r)=normrnd(xp'*br,sqrt(sig2r)); %careful - normrnd requires 
            %the standard deviation as input, not the variance
        end
        ypemat=exp(ypmat);

        fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
        fprintf(fid,'\n');
        kdiag=klausdiagnostics(ypmat);
        out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
        fprintf(fid,'Output table for log-predictions Model\t%6.0f\n', i);
        fprintf(fid,'mean\t\tstd\t\tnse\n\n');
        fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
        fprintf(fid,'\n');

        kdiag=klausdiagnostics(ypemat);
        out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
        fprintf(fid,'Output table for predictions in $$s Model\t%6.0f\n', i);
        fprintf(fid,'mean\t\tstd\t\tnse\n\n');
        fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
        fprintf(fid,'\n');   
        
        %Take pM(i)*r2 draws from that PPD
        if pM(i)>0
            int=randsample(ypmat,round(pM(i)*r2));
            ypwmat=[ypwmat int];
        end
        
        % save results
        filename = strcat('c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_fishing_MA',int2str(i));
        save(filename, 'betamat', 'sig2mat','ypmat');
        clear betamat sig2mat ypmat; %make some room
%             
        
    end % end model-specific loop
end % end main loop
ypewmat = exp(ypwmat);

% capture statistics for model-weighted predictions
kdiag=klausdiagnostics(ypwmat);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Output table for log-predictions, BMA\n\n');
fprintf(fid,'mean\t\tstd\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

kdiag=klausdiagnostics(ypewmat);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Output table for predictions in $$s, BMA\n\n');
fprintf(fid,'mean\t\tstd\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


save c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_fishing_MA ypwmat ypewmat;

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;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create plot for baseline model, full model, and BMA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%First, define a color for each plot
% Each color is expressed in numrical triplets
c1=randsample([0 0.25 .5 1],M,true)';
c2=randsample(c1,M);
c3=randsample(c1,M);
ColMap=[c1 c2 c3];


figure(1);
subplot(1,2,1)

for i=1:M
     filename = strcat('c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_fishing_MA',int2str(i));
        load(filename, 'ypmat');
     [f1,x1]=ksdensity(ypmat,'kernel','epanechnikov','npoints',300);
     plot(x1,f1, 'Color',ColMap(i,:));
     set(gca,'Xlim',[-2 10]);
     set(gca,'Ylim',[0 0.45]);
     hold on;
end
[f1,x1]=ksdensity(ypwmat,'kernel','epanechnikov','npoints',300);
plot(x1,f1,'--r');
legend('BMA');
title('PPD for individual models and BMA in log form');
hold off;

subplot(1,2,2)
for i=1:M
     filename = strcat('c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS_fishing_MA',int2str(i));
        load(filename, 'ypmat');
     ypemat=exp(ypmat);
        [f1,x1]=ksdensity(ypemat,'kernel','epanechnikov','npoints',300);
     plot(x1,f1, 'Color',ColMap(i,:));
     set(gca,'Xlim',[0 1000]);
     set(gca,'Ylim',[0 0.01]);
     hold on;
end
[f1,x1]=ksdensity(ypewmat,'kernel','epanechnikov','npoints',300);
plot(x1,f1,'--r');
title('PPD for individual models and BMA in $$s');
legend('BMA');

  
