function [val1, val2] = HARAmax1(X,XDistr,gama,CEqTolerance)
%Copyright (c) 2001-2009 by Ales Cerny
%usage [IP,alpha1] = HARAmax1(X,XDistr,gama,CertaintyEqTolerance)
%computes IP and optimal portfolio for a single risky asset with excess
%return X using normalized HARA utility (local risk aversion = 1)

%************************************************************************%
% HARAmax1.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 HARAmax1.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 = 0;
CEqPrecision=2*CEqTolerance;
%*******************%
%   the main loop   %
%*******************%
while abs(CEqPrecision) >= CEqTolerance;
    wealth = 1 + X*alpha;
    u = (wealth.^(1-gama))*XDistr';
    du = (1-gama)*(X.*(wealth.^(-gama)))*XDistr';
    ddu = gama*(gama-1)*((X.^2).*(wealth.^(-gama-1)))*XDistr';
    CE = u.^(1/(1-gama));
    CEqPrecision = -1/2/(1-gama)*du^2/u/ddu*gama*CE/(1+gama*(CE-1));
    alpha = alpha - du/ddu;        
end;
val1=gama*(CE-1);
val2=gama*alpha;