% Copyright (c) 2001-2009 by Ales Cerny

%************************************************************************%
% chapter3sect10.m - supplementary program to                            %
% Ales Cerny (2009) Mathematical Techniques in Finance (2nd ed.)         % 
% Princeton University Press http://press.princeton.edu/titles/9079.html %
%************************************************************************%
    
% This code is provided 'as-is', without any express or implied warranty.  
% 
% Permission is granted to anyone to use this code for any purpose, 
% subject to the following restrictions:
% 
% 1. The origin of this code must not be misrepresented; you must not
%    claim that you wrote the original code.
% 2. Modified code versions must be plainly marked as such, and must not 
%    be misrepresented as being the original code.
% 3. This notice may not be removed from any source distribution.

% NOTICE TO STUDENTS: To avoid accusations of plagiarism, if you use this 
% code or its modifications in assessed work you should prepend it with a 
% note stating: 
%   "This is the original/modified version of the code chapter3sect10.m by 
%    Ales Cerny (2009), Mathematical Techniques in Finance (2nd ed.),  
%    Princeton University Press. The original version is available from 
%    http://www.martingales.info/mtfweb2".
% A similar acknowledgement should appear prominently inside your written 
% report.

% One period optimal investment with HARA utility, incomplete markets, two
% risky assets, one year horizon, NIKKEI 225 real returns.

clear;
clc;

%*******************%
%    parameters     %
%*******************%
Risky1 = -0.35:0.05:0.90;       % row vector of rates of return from -35% to 90% in 5% increments
XDistr = [0.011                 % relative frequency of annual real returns on NIKKEI 225, monthly data 1960-2000
          0.023
          0.021
          0.041
          0.034
          0.061
          0.112
          0.095
          0.097
          0.081 
          0.091
          0.068
          0.074
          0.049
          0.047
          0.028
          0.015
          0.008
          0.013
          0.021
          0.000
          0.000 
          0.004
          0.000
          0.002
          0.004 ]'; 
strike    = 1.0;
SafeRate  = 0.02;
vol       = sqrt((log(1+Risky1).^2)*XDistr'-(log(1+Risky1)*XDistr')^2);
Risky2    = max([zeros(1,length(Risky1));((1+Risky1)-strike)]);
BSprice2  = BScall(1,strike,SafeRate,1,vol);    
gamavect  = [1.000001; 2; 5; 15; 10000; 100000000];
disp('                   *********************************************************************');
disp('                    TYPE IN THE OPTION PRICE as a FRACTION of FAIR PRICE and PRESS ENTER');
disp('                   *********************************************************************');
pricefract= input('Enter the option price as a fraction of fair price (0,1): ');
clc;
price2    = pricefract*BSprice2;
X         =(Risky2/price2-1-SafeRate);

%*********************************************************%
%        global parameters for procedure AASR             %
% (this procedure computes Arbitrage-Adjusted Sharpe R.)  %
%							                              %
%   X, RiskyDistr, PrecisionTolerance                     %
%*********************************************************%
IPPrecisionTolerance = 10^(-12);
%_output=0;

%****************************%
% main body of the programme %
%****************************%
HARAIP      = zeros(length(gamavect),1);
HARAbeta    = zeros(length(gamavect),1);
suboptHARAIP= zeros(length(gamavect),1);
[IPQm,betaQm]=TrQUmax(X,XDistr,IPPrecisionTolerance);

for ii = 1 : length(gamavect)
    gama=gamavect(ii);
    [HARAIP(ii),HARAbeta(ii)] = HARAmax(X,XDistr,gama,IPPrecisionTolerance);
    suboptHARAIP(ii) = HARAIPvalue(X, XDistr, betaQm, gama);
end



disp(' ______________________________________________________________________________________________________________________ ');
disp(' ');
disp(sprintf('                         option price/Black-Scholes value =  %5.2f', pricefract));
disp(' ');                                           
disp(' ');
disp(sprintf('  HARA parameter gamma      %8.0f    %8.0f    %8.0f    %8.0f    %8.0f    %8.0e   %8.0e', [-1 gamavect']));
disp(' ');
disp(sprintf('  HARA beta                 %8.3f    %8.3f    %8.3f    %8.3f    %8.3f    %8.3f   %8.3f', [betaQm HARAbeta']));
disp(' ');
disp(' ');
disp(sprintf('  HARA IP of TrQU beta      %8.3f    %8.3f    %8.3f    %8.3f    %8.3f    %8.3f   %8.3f', [IPQm suboptHARAIP']));
disp(' ');
disp(' ');
disp(sprintf('  maximum HARA IP           %8.3f    %8.3f    %8.3f    %8.3f    %8.3f    %8.3f   %8.3f', [IPQm HARAIP']));
disp(' ______________________________________________________________________________________________________________________ ');
disp(' ');
disp(' ');
disp(' ');
disp('                   *********************************************************************');
disp('                   PRESS ANY KEY TO SEE WEALTH DISTRIBUTION WITH LOCAL RISK AVERSION = 5');
disp('                   *********************************************************************');
pause;
clc;
disp(' ______________________________________________________________________________________________________________________ ');
disp(sprintf('              option price/Black-Scholes value =  %5.2f ', pricefract));
disp(' ');
disp(' ');
disp(' NIKKEI rate   probability               total return on safe wealth with local risk aversion 5  '); 
disp('   of return     NIKKEI >=    ');
disp(sprintf('                this level   gama = -1    %8.0f   %8.0f   %8.0f   %8.0f   %8.0e   %8.0e ', gamavect'));
disp(' ');
TailDistr = 1 - cumsum(XDistr);
TotRet=1+X'*[betaQm ; HARAbeta]'/5;
for i = 1 : length(Risky1)
    disp(sprintf('     %8.1f    %8.1f    %8.1f     %8.1f   %8.1f   %8.1f   %8.1f    %8.1f   %8.1f   %8.1f',...
        100*Risky1(i), 100*TailDistr(i), 100*TotRet(i,:)));
end
disp(' ');
disp(' ______________________________________________________________________________________________________________________ ');disp(' ');
disp(' ');
disp(' ');
disp(' ');