%Program: dpfig74.m
%uses the optimization toolbox as licensed by UTA
% Written 7/26/2002   Matlab Version:6.0.0.88 (R12) on PCWIN
% Optimization Toolbox Version: Version 2.1   (R12)
%****************************************************************
%* creates Chapter figures 7.4, page 227 Dixit & Pindyck        *
%*"Investment Under Uncertainty", 2nd edition, Princeton Press  * 
%* % "fsolve" uses a starting guess. A solution MAY NOT BE FOUND*
%* based on the starting guess. Thus creating a single m-file   *
%* that creates multiple solutions can pose problems.           *
%* This software replace each 'guess' with the previous         *
%* set of calculations. The initial guess was set to values     *
%* that were close to the E= -5 prices. A1 & B2 were guesses    * 
%* Note that future Matlab releases could cause code errors.    *
%*                                                              *
%* Bottom line, finding values at C=0 are difficult             *                                                 
%****************************************************************
format short
%*************************************************************
% Part I: Specify parameters and calculate Beta1 and Beta2   * 
% for al values of sigma ( the volatility )                  *
%*************************************************************
interplot = zeros(1,3); %Initialize
%******************************************************
% UNPREDICTABLE RESULTS CAN BE GENERATED due to fsolve*                                                                                       
% AND THE NEED FOR A STARTING GUESS                   *
%******************************************************
sigma = 0.2;
r=0.04;
delta = 0.04;
II = 20;
C = 0.8;
Q=10;
%E=2.0; E varies -- see loop below
r=0.04; delt = 0.04;
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))
%***********************************************************
% Part II: Create parameters to solve the non linear system*
% and call the function. Must start with a guess in "vars0"*                                    
% the function "fff" returns the vars vector with:         *
% vars(1,1)= A1  vars(1,2)= ph vars(1,3)= B2 vars(1,4)=pl  *
% where ph is the high price and pl is the low price       *
%***********************************************************
%guessA1 = vars(1,1) 
%guessPH = vars(1,2)
%guessB2 = vars(1,3)
%guessPL = vars(1,4)
guessA1 = 2%36.05
guessPH = 1%0.4%0.33; changed from 0.4 for C=0
guessB2 = 0.4%0.057
guessPL = 0.5%0.02% 0.05; changed from 0.02 for C=0
vars0 = [ guessA1 guessPH  guessB2  guessPL]; % initial guess
for E = -5:5:15
%C = kk*0.1 %vars0is your guess% c=.2,.3, 0.4 OK
%OPTIMSET is recommended for setting options.
% seems that moving from le-4 to le-6 "loosened" tolerance
% for 'good' roots of the equations
opts=optimset('fsolve');
opts=optimset(opts,'Maxiter',1100,'Tolx',1e-6,'tolfun',1e-6);
%was 900 e-4 then e-6
%**********************************************************
%* Next instruction calls function fff.m                  *
%**********************************************************
vars = fsolve(@fff,vars0,opts,Q,Beta1,Beta2,delt,C,r,II,E)
%***********************************************************
% Part III: Populate interplot by row; each row has the    *
% variables necessary to generate the figure including "C".*
% **********************************************************
guessA1 = vars(1,1); % creates a new 'guess' for the 
guessPH = vars(1,2); % next value of "C"
guessB2 = vars(1,3); % Once a good initial guess for C=0.1 
guessPL = vars(1,4); % this 'replacement' procedure worked
vars0 = [ guessA1 guessPH  guessB2  guessPL];
interplot = [interplot; E vars(4) vars(2)  ];
clear functions 
clear A1 B2 ph pl
end
interplot(1,:) = [ ]% delete first row leaving just calculated values
%********************************************************************
% Print the results in the Matlab workspce on screen
fprintf(1,'%11s', 'E ','       Price-Low ','  Price-High')
fprintf(1,'\n')
fprintf(1,'%11.3f%15.3f%15.3f\n',interplot')% transpose interplot
%********************************************************************
plowlabel = char(num2str(interplot(:,2),'%7.3f'));
plot (interplot(:,1), interplot(:,2), ':squarer' );
% remember philabel and plowlabel are strings hence column is 1:end 
text(interplot(:,1),interplot(:,2),plowlabel(:,1:end), ...
    'Verticalalignment','bottom','Fontweight', 'Bold','Fontsize',9);
hold
philabel = char(num2str(interplot(:,3),'%7.3f'));
plot (interplot(:,1), interplot (:,3),'-diamond','Linewidth',1.2 );
% philabel and plowlabel are strings hence column is 1:end 
text(interplot(:,1),interplot(:,3),philabel(:,1:end), ...
    'Verticalalignment','bottom','Fontweight', 'Bold','Fontsize',9);
xlabel ('E(millions of Dollars)','Fontweight', 'Bold');
ylabel ('Price( dollars per pound)','Fontweight', 'Bold');
title('Dixit&Pindyck,Fig 7.4 p227 Thresholds as f(Abandonment Cost) sigma=0.2', ...
'Fontweight', 'Bold');
% ********************************
%  Values calculated:            *
%     E         P-Low     P-High *
%  -5.0000    0.5869    1.2754   *
%        0    0.5563    1.3202   *
%   5.0000    0.5296    1.3582   *
%  10.0000    0.5057    1.3914   *
%  15.0000    0.4837    1.4208   *
%*********************************