%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PPP value for skewness in earnings data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

tic;  % start stop watch

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

% load data posterior results for female earnings regression
load c:\klaus\AAEC6564\mlab\worksp\mod2application;
R=size(betamat,2); %number of draws from the GS
n=length(y);
k=size(X,2);

sk_y=zeros(R,1); 
sk_ystar=zeros(R,1);

for r=1:R % start main loop
    br=betamat(:,r);
    sig2r=sig2mat(r);
    %compute sk for actual data:
    sk_y(r)=sqrt(n)*sum((y-X*br).^3)./(sum((y-X*br).^2)^(3/2));
    % draw simulated data:
    ystar=normrnd(X*br,sqrt(sig2r)); 
    %this generates N by 1 draws of yi, each corresponding to its respective xi
    %compute sk for simulated data
    sk_ystar(r)=sqrt(n)*sum((ystar-X*br).^3)./(sum((ystar-X*br).^2)^(3/2));
end

% Compute PPP
%%%%%%%%%%%%%%%%%%%%%%%%
f=find(sk_ystar<mean(sk_y));
int=length(f)/R; %proportion of sky_star's that are smaller than E(sk_y)
PPP=min(int,1-int)

fprintf(fid,'total number of predictive draws =\t%8.0f \n',R);
fprintf(fid,'\n');
fprintf(fid,'PPP =\t%10.8f \n',PPP);
fprintf(fid,'\n');
    
% Compute and plot difference 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d=sk_ystar-sk_y;
[f1,x1]=ksdensity(d,'kernel','epanechnikov','npoints',300);

figure(1);
subplot(1,2,1)
plot(x1,f1,'-k');
title('sk(y*) - sk(y)');
set(gca,'Xlim',[-2 2]);
set(gca,'Ylim',[0 2.5]);

% plot sk(y*) and relative location of sk_y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[f1,x1]=ksdensity(sk_ystar,'kernel','epanechnikov','npoints',300);
xm=repmat(mean(sk_y),1,1000);
ym=(0:1:999);

subplot(1,2,2)
plot(x1,f1,'-k',xm,ym,'--k');
title('sk(y*) and E(sk_y)');
legend('sk(y*)','E(sk_y)');
set(gca,'Xlim',[-2 2]);
set(gca,'Ylim',[0 2.5]);

text(-1,2.2,['E[sk_y] = ' num2str(mean(sk_y))]);
    
    
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;


    