%Program: dpfig73.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.3, 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 C=0.1 prices. The initial B2 and A1   * 
%* values were guesses based on their behavior for C=0.2 and 0.3*
%*            A starting C=0 value was not used due to           *
%*            difficulties finding a solution                   * 
%* 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;
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 = 36.05%
guessPH = 0.4%
guessB2 = 0.057%
guessPL = 0.02% 
vars0 = [ guessA1 guessPH  guessB2  guessPL]; % initial guess
for kk = 1:10;
C = kk*0.1 %vars0is your guess% c=.2,.3, 0.4 OK
%C=0.0
%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',1300,'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; C vars(4) vars(2)  ];
clear functions 
clear A1 B2 ph pl
end
interplot(1,:) = [ ]% delete first row leaving just calculated values
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 ('C(dollars per pound)','Fontweight', 'Bold');
ylabel ('Price( dollars per pound)','Fontweight', 'Bold');
title('Dixit&Pindyck,Fig 7.3 p227 Entry & Exit Thresholds as f(Oper cost)', ...
'Fontweight', 'Bold');
% ********************************
%  Values calculated:            *
%     C         P-Low     P-High *
%   0.1000    0.0518    0.3399   *
%   0.2000    0.1151    0.4978   *
%   0.3000    0.1823    0.6468   *
%   0.4000    0.2519    0.7904   * 
%   0.5000    0.3233    0.9304   *
%   0.6000    0.3962    1.0677   *
%   0.7000    0.4702    1.2028   *
%   0.8000    0.5452    1.3361   *
%   0.9000    0.6210    1.4680   *
%   1.0000    0.6976    1.5987   *
%*********************************