%Program: dpfigure63.m (MATLAB) for D-P figure 6.3 page 192  * 
% from "Investment Under Uncertainty", Dixit & Pindyck, 2nd  *
% edition, Princeton University Press.                       *
% Specify inputs and find P* and create graphs               *
% For Educational purposes only                              * 
% calculate P*. fzero is fussy about the range see below     *
% Written: 7/24/2002  Matlab Version: 6.0.0.88 (R12) on PCWIN*
%*********** SPECIFY INPUTS inelegant ************************
r=0.04;
delta = 0.04;
sigma = 0.2;
I = 100;
C = 10;
Beta1=0.5 - (r-delta)/sigma^2 + ...
    sqrt(((r-delta)/sigma^2 - 0.5)^2+2*r/(sigma^2))
Beta2=0.5 - (r-delta)/sigma^2 - ...
    sqrt(((r-delta)/sigma^2 - 0.5)^2+2*r/(sigma^2))
B2 = ( C^(1-Beta2)/(Beta1 - Beta2) )* ((Beta1/r) - ...
    (Beta1 - 1)/delta )
K1 = (C^(1-Beta1)/(Beta1 - Beta2))* (Beta2/r - ...
    (Beta2 -1)/delta)
%******************************************************************
% build the string function to be used by Matlab's fzero.         *
% Must be a string function. the range in whhich fzero is to look *
% is 20 to 40. The signs at the endpoints of the range must differ*
% Pstar = 23.798 for this example                                 *
%******************************************************************
firstcon = (Beta1 - Beta2)*B2;
secondcon = (Beta1 - 1)/delta;
thirdcon = Beta1*(C/r + I);
nufunct = strcat(num2str(firstcon),'*','x','^',num2str(Beta2), ...
    '+',num2str(secondcon),'*','x','-',num2str(thirdcon) );
Pstar = fzero(nufunct,[20 40] )
%*******************************************************************
% calculated  P* = 23.798 above now use it to calculate A1. 7/11   * 
% see DP p.191 -- will sub P* into equation (15) and solve for A1  *
% See dixit pindyck p.191 equation (14) and (15).                  *
%*******************************************************************
rtside = B2*Pstar^Beta2 + Pstar/delta - C/r - I;%right hand side eq. 14
A1 = rtside/Pstar^Beta1
%******************************************
% equation to value the the option at P*  * 
%******************************************
FP = A1*Pstar^Beta1;
%******************************************************************
%create table to store OutputPrice, ProjectNPV and, OptionVAlue   *
% resultis will be a  matrix ( build as needed ).                 *
% Column 1 has OutputPrice                                        *
% Column 2 has ProjectNPV                                         *
% Column 3 has OptionValue                                        *
% Create resultis matrix and initialize                           *                                                                                                                                                                                                                      
%******************************************************************
resultis = [ 0 0 0 ]; % a 1 by 3
%******************************************************************
%calculate FP, value of option to invest at intervals of $4 for P *
% calculate Project NPV, as V(P) - I based on P>C else P<=C       *
%********************* V(P) - I ***********************************
for P = 0:4:28;
    if P > C;
    VP = B2*P^Beta2 + P/delta - C/r - 100; %For P>C 
    FP = A1*P^Beta1;
else
    VP = K1*P^Beta1 - 100;
    FP = A1*P^Beta1;
end
if P >= Pstar; % once P* is reached the option cannot be greater than 
               % the underlying asset. in this case the project value 
    FP = VP ;
end
interim = [ P VP FP ];  
resultis = [ resultis ; interim ];
end
resultis(1,:) = [ ]; % deletes first row leaving just calculated values
% Print the results in the Matlab workspce on screen
fprintf(1,'%10s', '   Price','        Project VP-I','      F(P)')
fprintf(1,'\n')
fprintf(1,'%11.3f%15.3f%15.3f\n',resultis')
%**************************************************************
% Reproduce and create a plot for Figure 6.3 using P=0,4, etc *       
% values are stored in matrix 'resultis' . Create 2 plots     *                      
%**************************************************************
plot(resultis(:,1), resultis(:,2),'-+') 
projectlabel = char(num2str(resultis(:,2),'%8.3f'));
text(resultis(1:5,1), resultis(1:5,2),projectlabel(1:5,1:end), ...
    'Verticalalignment','top')
%******************************************************************
% note that col 1 has the 'x' value, column 2 has the ProjectValue*
% project value is a solid line. Col 3 has F(P) value.            *
%******************************************************************
hold
% the original plot is held and a second plot is added for OptionValue
plot(resultis(:,1), resultis(:,3),':squarer' )
FPlabel = char(num2str(resultis(:,3),'%8.3f'));
text(resultis(2:end,1), resultis(2:end,3),FPlabel(2:end,1:end), ...
    'Verticalalignment','bottom')
title('Dixit Pindyck Figure 6.3 pg 192','Fontweight', 'Bold')
xlabel('Price','Fontweight','Bold')
ylabel('F(P)with squares,"V(P)-I" solid with +','Fontweight','Bold')
text(0,0,'F(P)','fontsize',14,'verticalalignment','bottom', ...
    'Fontweight','Bold')
%***********************************************************
% Should get these Calculated Values:  Beta1=2 Beta2=-1    *
%           Price        Project VP-I      F(P)            *
%             0.000       -100.000          0.000          *
%             4.000        -86.667          7.909          *
%             8.000        -46.667         31.638          *
%            12.000         19.444         71.185          *
%            16.000        102.083        126.551          *
%            20.000        191.667        197.736          *
%            24.000        284.722        284.722          *
%            28.000        379.762        379.762          *
%***********************************************************