%graphit:
% script to handle graphing from estimates 

clear all; 
close all;
printgraphs = 0; 
format bank; 

load base1;  % loads workspace from base case estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);


% 1. graph numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxip = px; 
pyip = py; % saved for joint graph later. 

% 2. Graph returns and pr(ipo/value) 

figure; 
% select which dates to show in months */

datesel = [1 4 8 20 40]'; 
py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));
px = lvgrid; 
plot(100*px,py(:,datesel)./max(max(py(:,datesel))'),'-k'); 
hold on; 
plot(100*px,findpip(lvgrid,a,b),'--k');
text(90,0.7,'3 mo.');
text(110,0.55,'1 yr.'); 
text(160,0.45,'2 yr.'); 
text(240,0.4,'5, 10 yr.'); 
text(300,0.7,'Pr(IPO,acq.|V)'); 
axis([ -400 500 -inf inf]);
%set(gca,'ytick',[])
xlabel('Log returns (%)');
ylabel('Scale for Pr(IPO,acq.|V)'); 
%set(gca,'YTick',0);
%title('Distribution of log returns at ipo');
if printgraphs;
   print -depsc2 m3.eps;
end;

meanrets =  sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel));
sdrets =  (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanrets.^2);


% graph and characterize returns to investors as a function of horizon,
% with and without complete loss if bankrupt

% call sim with pi = 0 (so we see actual values) and with big option on, so
% we see values at bankruptcy etc. 
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,0,c,d,0,logrf,logspret,lvgrid,stockind,dok,1);

         
         % it would make more sense to do continuous compounding, but
         % previous numbers used quarterly discrete compounding and
         % multiply by 400, so same here. 
yans = zeros(54,1); 
for tmax = 1:54; 
    ylo = 1;
    yhi = 2; 
    while (yhi-ylo)>0.00001
        y = ylo+(0:0.1:1)'*(yhi-ylo);
        pv = 0*y; 
        for t=1:tmax;
            pv = pv + sum((exp(lvgrid)'.*(pseeipt(t,:)+pnosipt(t,:)+pseebkt(t,:)+pnosbkt(t,:)))')./(y.^t); 
        end; 
        pv = pv +  sum((exp(lvgrid)'.*pnoseet(t,:))')./(y.^tmax);
        [miny,indx] = min(abs(pv-1));
        if (indx > 1)&(indx < size(y,1)); 
            ylo = y(indx-1);
            yhi = y(indx+1); 
        else; 
            disp('graphit error in yield calculation'); 
            stop; 
        end;
    end; 
    yans(tmax) = y(indx); 
end; 
    
% now, same if bankrupt are zero payoff 
    
yans2 = 0*yans; 
for tmax = 1:54; 
    ylo = 1;
    yhi = 2; 
    while (yhi-ylo)>0.00001
        y = ylo+(0:0.1:1)'*(yhi-ylo);
        pv = 0*y; 
        for t=1:tmax;
            pv = pv + sum((exp(lvgrid)'.*(pseeipt(t,:)+pnosipt(t,:)))')./(y.^t); 
        end; 
        pv = pv +  sum((exp(lvgrid)'.*pnoseet(t,:))')./(y.^tmax);
        [miny,indx] = min(abs(pv-1));
        if (indx > 1)&(indx < size(y,1)); 
            ylo = y(indx-1);
            yhi = y(indx+1); 
        else; 
            disp('graphit error in yield calculation'); 
            stop; 
        end;
    end; 
    yans2(tmax) = y(indx); 
end; 


figure; 
plot((1/4:1/4:54/4)',[400*(yans-1) 400*(yans2-1)]);
title('expected returns with and without payoff at bankruptcy'); 
disp('horizon, expected returns with and without payoff at bankruptcy'); 
horiz = (1:54)'/4; 
tbl = [horiz 400*(yans-1) 400*(yans2-1)]; 
disp(tbl([1 2 3 4 8 12 16 20 24 28 32 36 40],:)); 

% simulation with alpha set to zero and other parameters at same value. 

mlogrf = mean(logrf);
mum = mean(logspret);
sigm = std(logspret);
gama0 = -log( exp(delt*(mum-mlogrf)+1/2*delt^2*sigm^2+1/2*sig^2)*...
                     (1+(exp(delt*sigm^2)-1)*(exp(-(mum-mlogrf)-1/2*sigm^2)-1)/(exp(sigm^2)-1)));

[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gama0,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% save numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
pxipa0only = px*h; 
pyipa0only = [0 0 100 ; 100*py]; 


% -----------------------------------------------
% same thing, but for round to round sample

load base2;  % loads workspace from base case estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% 1. graph numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
pxall = px*h; 
pyall = [0 0 100 ; 100*py]; 

% *******GRAPHS of histories for model/data 

% NEW: each of IPO and all, model and data together. 

load histories; %  px1 py1 px1ip py1ip;  from datachar

figure;
indxs = (3:3:size(pxip,1))';
plot(pxip,pyip,'-k',px1ip,100*py1ip,'--k'); 
%/*Title("Outcome history");*/
xlabel('Years since investment');
ylabel('Percentage');
axis([0 8 0 100]); 
text(3,75,'IPO, acquired');
text(0.5,85,'Still private'); 
text(1,10,'Out of business'); 

% EDIT AND SAVE AS f1.eps

figure; 
plot(pxall,pyall,'-k',px1,100*py1,'--k'); 
%/*Title("Outcome history");*/
xlabel('Years since investment');
ylabel('Percentage');
axis([0 8 0 100]); 
text(3,75,'IPO, acquired, or new round');
text(0.5,85,'Still private'); 
text(1,10,'Out of business'); 

% edit and save as m2.eps. 

% OLD: joint graph of both models. 

% figure;
% indxs = (3:3:size(pxip,1))';
% plot(pxall,pyall,'--k',pxip,pyip,'-k',pxip(indxs),pyip(indxs,:),'vk'); 
% %/*Title("Outcome history");*/
% xlabel('Years since investment');
% ylabel('Percentage');
% axis([0 8 0 100]); 
% text(3,75,'IPO, acquired, or new round');
% text(0.5,85,'Still private'); 
% text(1,10,'Out of business'); 
% m2 is name of graph

% 2. Graph returns and pr(ipo/value) 

figure; 
% select which dates to show in months */

datesel = [1 4 8 20 40]'; 
py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));
px = lvgrid; 
plot(100*px,py(:,datesel)./max(max(py(:,datesel))'),'-k'); 
hold on; 
plot(100*px,findpip(lvgrid,a,b),'--k');
text(85,0.9,'3 mo.');
text(110,0.55,'1 yr.'); 
text(160,0.45,'2 yr.'); 
text(220,0.3,'5, 10 yr.'); 
text(220,0.6,'Pr(New fin.|V)'); 
axis([ -400 500 -inf inf]);
%set(gca,'ytick',[])
xlabel('Log returns (%)');
ylabel('Scale for Pr(new fin.|V)'); 
%set(gca,'YTick',0);
%title('Distribution of log returns at ipo');
if printgraphs;
   print -deps2 m3rnd.eps;
end;

disp('NOTE: add arrows to graphs and resave'); 

meanretsall =  sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel));
sdretsall =  (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanretsall.^2);


% *************************** alpha = zero cases -- IPO sample

load a0ip;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0; % graph actual returns 
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% 1. graph numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
pxipa0 = px*h; 
pyipa0 = [0 0 100 ; 100*py]; 

% 2. Graph returns and pr(ipo/value) 

figure; 
% select which dates to show in months */
datesel = [1 4 8 20 40]'; 
py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));
px = lvgrid; 
plot(100*px,py(:,datesel)./max(max(py(:,datesel))'),'-k'); 
hold on; 
plot(100*px,findpip(lvgrid,a,b),'--k');
text(90,0.7,'3 mo.');
text(110,0.55,'1 yr.'); 
text(160,0.45,'2 yr.'); 
text(240,0.4,'5, 10 yr.'); 
text(300,0.7,'Pr(IPO,acq.|V)'); 
axis([ -400 500 -inf inf]);
%set(gca,'ytick',[])
xlabel('Log returns (%)');
ylabel('Scale for Pr(IPO,acq.|V)'); 
%set(gca,'YTick',0);
%title('Distribution of log returns at ipo');

meanrets =  [meanrets; sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdrets =  [sdrets; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanrets(end,:).^2)];

% same for ER=15

load ER15ip;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxiper = px; 
pyiper = py; % saved for joint graph later. 


 %2. Graph returns and pr(ipo/value) 

figure; 
% select which dates to show in months */

datesel = [1 4 8 20 40]'; 
py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));
px = lvgrid; 
plot(100*px,py(:,datesel)./max(max(py(:,datesel))'),'-k'); 
hold on; 
plot(100*px,findpip(lvgrid,a,b),'--k');
text(90,0.7,'3 mo.');
text(110,0.55,'1 yr.'); 
text(160,0.45,'2 yr.'); 
text(240,0.4,'5, 10 yr.'); 
text(300,0.7,'Pr(IPO,acq.|V)'); 
axis([ -400 500 -inf inf]);
%set(gca,'ytick',[])
xlabel('Log returns (%)');
ylabel('Scale for Pr(IPO,acq.|V)'); 
%set(gca,'YTick',0);
%title('Distribution of log returns at ipo');
 

meanrets = [meanrets;  sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdrets =  [sdrets; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanrets(end,:).^2)];

% now the no pi cases 

load a0ipnopi;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% 1. graph numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
pxipa0nopi = px*h; 
pyipa0nopi = [0 0 100 ; 100*py]; 

% 2. Graph returns and pr(ipo/value) 

py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));
px = lvgrid; 

meanrets =  [meanrets; sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdrets =  [sdrets; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanrets(end,:).^2)];

% same for ER=15

load ER15ipnopi;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxipernopi = px; 
pyipernopi = py; % saved for joint graph later. 


% 2. Graph returns and pr(ipo/value) 

py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));
px = lvgrid; 

meanrets = [meanrets;  sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdrets =  [sdrets; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanrets(end,:).^2)];



% ------------------------- alpha = zero cases -- ALL sample

% 1. alpha zero 
load a0all;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxalla0 = px; 
pyalla0 = py; % saved for joint graph later. 

% Graph returns and pr(ipo/value) 

py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));

meanretsall =  [meanretsall; sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdretsall =  [sdretsall; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanretsall(end,:).^2)];

% 2.  ER=15

load ER15all;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxaller = px; 
pyaller = py; % saved for joint graph later. 

% returns and pr(ipo/value) 

py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));

meanretsall =  [meanretsall; sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdretsall =  [sdretsall; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanretsall(end,:).^2)];

% 3. alpha zero nopi 

load a0allnopi;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

% numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxalla0nopi = px; 
pyalla0nopi = py; % saved for joint graph later. 

%  returns and pr(ipo/value) 

py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));

meanretsall =  [meanretsall; sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdretsall =  [sdretsall; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanretsall(end,:).^2)];

% same for ER=15 no pi

load ER15allnopi;   % loads workspace from alhpa zero  estimation 
clear xlabel; % remove this line after next mlest3 run; won't include xlabel as a variable. 
pim = 0;
[pnoseet,pseeipt,pnosipt,pseebkt,pnosbkt] = ...
             sim3(gam,delt,sig,k,a,b,pim,c,d,0,logrf,logspret,lvgrid,stockind,dok,0);

%  numbers in each state as a function of time */

px = (0:rows(pnoseet))';
py = [cumsum(sum(pseeipt')'+pnosipt) cumsum(pseebkt+pnosbkt) pnoseet];
px = px*h; 
py = [0 0 100 ; 100*py]; 
pxallernopi = px; 
pyallernopi = py; % saved for joint graph later. 

%  Graph returns and pr(ipo/value) 
py = pseeipt'./(ones(size(pseeipt',1),1)*sum(pseeipt'));

meanretsall =  [meanretsall; sum((lvgrid*ones(1,size(datesel,1))).*py(:,datesel))];
sdretsall =  [sdretsall; (sum(((lvgrid.^2)*ones(1,size(datesel,1))).*py(:,datesel))-meanretsall(end,:).^2)];



 % ************************ Finally, graph of all these histories 

 
figure;
indxs = (4:4:size(pxip,1))';
mindxs = (11:12:size(px1ip,1))';

plot( px1ip,100*py1ip(:,2),'--k',...   % data
      px1ip,100*py1ip(:,1),'-k',...
      px1ip(mindxs),100*py1ip(mindxs,1),'sk',... 
      px1ip(mindxs),100*py1ip(mindxs,2),'sk',...
     pxipa0,pyipa0(:,1),'--k',...  % alpha zero
     pxipa0,pyipa0(:,2),'-k',...
     pxipa0(indxs),pyipa0(indxs,1),'vk',...  % alpha zero
     pxipa0(indxs),pyipa0(indxs,2),'vk',...
     pxiper,pyiper(:,1),'--k',... % ER=15
     pxiper,pyiper(:,2),'-k',... 
     pxiper(indxs),pyiper(indxs,1),'ok',... 
     pxiper(indxs),pyiper(indxs,2),'ok',... 
     pxipa0only,pyipa0only(:,1),'--k',... % alpha zero no other parameters
     pxipa0only,pyipa0only(:,2),'-k'); 

xlabel('Years since investment');
ylabel('Percentage');
axis([0 8 0 65]); 

% Edit and save as ma0hist

%text(3,75,'IPO, acquired, or new round');
%text(0.5,85,'Still private'); 
%text(1,10,'Out of business'); 

% old figure of all simulations. Replaced by new figure of model and data
% for each of ipo, round

figure;
indxs = (4:4:size(pxip,1))';
plot(pxalla0,pyalla0(:,1:2),'-k',pxalla0nopi,pyalla0nopi(:,1:2),'-k',pxall(indxs),pyall(indxs,1:2),'vk',pxall,pyall(:,1:2),'-k'); 
% plot(pxalla0,pyalla0(:,1:2),'-k',pxaller,pyaller(:,1:2),'--k',pxalla0nopi,pyalla0nopi(:,1:2),'-k',pxallernopi,pyallernopi(:,1:2),'--k',...
%                      pxall(indxs),pyall(indxs,1:2),'vk',pxall,pyall(:,1:2),'-k'); 

                
xlabel('Years since investment');
ylabel('Percentage');
axis([0 8 0 100]); 
text(3,75,'IPO, acquired, or new round');
%text(0.5,85,'Still private'); 
text(1,10,'Out of business'); 
title('effects of measurement error'); 

  head = [];
 fmt = ' %8.0f '; 
 row = { 'horizon'; 'base IP, mean'; '$\alpha=0$ IP, mean'; 'ER=15\% IP, mean';'$\alpha=0, \pi=0$ IP, mean'; 'ER=15\%, $\pi=0$, IP, mean';...
                    'base IP, s.d.' ;'$\alpha=0$ IP, s.d.'; 'ER=15\% IP, s.d.';'$\alpha=0, \pi=0$ IP, s.d.'; 'ER=15\%, $\pi=0$, IP, s.d.';...
                    'base rnd, mean'; '$\alpha=0$ rnd, mean'; 'ER=15\% rnd, mean';'$\alpha=0, \pi=0$ rnd, mean'; 'ER=15\%, $\pi=0$, rnd, mean';...
                    'base rnd, s.d.' ;'$\alpha=0$ rnd, s.d.'; 'ER=15\% rnd, s.d.';'$\alpha=0, \pi=0$ rnd, s.d.'; 'ER=15\%, $\pi=0$, rnd, s.d.'};
 texprint([datesel'/4;100*meanrets;100*sdrets; 100*meanretsall; 100*sdretsall], fmt, head, row);