function [val1, val2] = HARAmax(X,XDistr,gama,IPTolerance,outputflag)
% Copyright (c) 2001-2009 by Ales Cerny
% usage [IP,alpha1] = HARAmax(X,XDistr,gama,IPPrecisionTolerance,outputflag)
% computes IP and optimal portfolio for multiple risky assets with excess 
% return X using normalized HARA utility (local risk aversion = 1). When
% outputflag is not specified or it is not equal to 1 then no output is 
% printed. When IPTolerance is not specified the default precision is 10^-6

% Copyright (c) 2001-2009 by Ales Cerny

%************************************************************************%
% HARAmax.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 HARAmax.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.

alpha = zeros(size(X,1),1);
iteration  = 0;
if (nargin<5);
    outputflag=0;
    if nargin<4
        IPTolerance = 10^-6;
    end
end;
IPPrecision= 2*IPTolerance;

% repeat iterations until the desired precision in certainty equivalent is
% attained.
while (abs(IPPrecision) >= IPTolerance )
    
    wealth = 1 + alpha'*X;
    u = (wealth.^(1-gama))*XDistr';
    du = (1-gama)*X*((wealth.^(-gama)).*XDistr)';
    ddu = gama*(gama-1)*X*(X.*(ones(size(X,1),1)*((wealth.^(-gama-1)).*XDistr)))';
    CE=u^(1/(1-gama)); 
    IPPrecision = -1/2/(1-gama)*du'*(ddu\du)/u*gama*CE/(1+gama*(CE-1));
    if (outputflag == 1)
        disp(sprintf('iteration                     %12.0g               ', iteration));
        disp(sprintf('investment potential          %12.3f               ', gama*(CE-1)));
        disp(sprintf('estimated precision in 1+IP   %12.1e               ', IPPrecision));
        disp(sprintf('%s',['alpha1                               ' num2str(gama*alpha','%12.3f')]));
        disp(sprintf('%s',['estimated precision in alpha1     ' num2str(gama*(-ddu\du)','%12.1e')]));
        disp(' ');
        disp(' ');
    end

    iteration=iteration+1;
    alpha = alpha - ddu\du;
end
val1 = gama*(CE - 1);
val2 = gama*alpha;