%Program: dpfig72.m  Run on 10/25/02 correctly
%uses the optimization toolbox as licensed by UTA
% Written 7/25/2002   Matlab Version:6.0.0.88 (R12) on PCWIN
% Optimization Toolbox Version: Version 2.1   (R12)
%****************************************************************
%* creates Chapter figures 7.2, page 226 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 was tested for each individual sigma value     *
%* and subsequent system of nonlinear simultaneous equations.   *
%* When I found a guess for 'vars0' that could generate         *
%* multiple solutions I then created this integrated software.  *
%* Interim output is shown.                                                             *
%* Note that future Matlab releases could cause code errors.    *
%****************************************************************
format short
%*************************************************************
% Part I: Specify parameters and calculate Beta1 and Beta2   * 
% for al values of sigma ( the volatility )                  *
%*************************************************************
interplot = zeros(1,4);
sigmastart=0.20   
%******************************************************
% guess OK for 0.2,0.25,0.3,0.35,0.4, not for 0.15    *
% see below for "vars0". Therefore sigma begins at 0.2*
% incremnented by 0.05 up to sigma=0.4.               *
% UNPREDICTABLE RESULTS CAN BE GENERATED due to fsolve*
% AND THE NEED FOR A STARTING GUESS                   *
%******************************************************
for kk =1:5
sigma = sigmastart + (kk-1)*0.05;
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))
%pause(10)
%***********************************************************
% 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       *
%***********************************************************
%vars0is your guess
vars0 = [10 1.2  6  0.55 ];
%OPTIMSET is recommended for setting options.
opts=optimset('fsolve');
opts=optimset(opts,'Maxiter',900,'Tolx',1e-4,'tolfun',1e-4);
%**********************************************************
% next instruction calls 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".*
% **********************************************************
interplot = [interplot; sigma vars(4) vars(2) C ]
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')
hold
plot (interplot(:,1), interplot(:,4),':black' )
philabel = char(num2str(interplot(:,3),'%7.3f'))
plot (interplot(:,1), interplot (:,3),'-diamond' )
% philabel and plowlabel are strings hence column is 1:end 
text(interplot(:,1),interplot(:,3),philabel(:,1:end), ...
    'Verticalalignment','bottom')
text(interplot(2,1),interplot(2,4),'C = 0.8', ...
    'Verticalalignment','bottom', ...
    'Fontweight', 'Bold')  
xlabel ('Sigma')
ylabel ('Price( dollars per pound)')
title('Dixit & Pindyck, Fig 7.2, p226 Copper Entry & Exit Thresholds', ...
'Fontweight', 'Bold')
% *******************************************
% * Values calculated:                      *
% * at sigma = 0.2 ph = 1.336 plow =0.545   *
% * at sigma = 0.25 ph = 1.441 plow = 0.514 *
% * at sigma = 0.3  ph = 1.548 plow = 0.488 *
% * at sigma = 0.35 ph = 1.657 plow = 0.465 *
% * at sigma = 0.4  ph = 1.769 plow = 0.445 *
% *******************************************