function [IP, beta] = TrQUmax(X, XDistr, IPPrecisionTolerance)
% Copyright (c) 2001-2009 by Ales Cerny
%usage [IP,beta] = AAQUmax(X,XDistr,gama,IPPrecisionTolerance)
%computes IP and optimal portfolio for multiple risky assets with excess 
%return X using normalized truncated quadratic utility. See Section 3.8.1.

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

%*******************%
%   initialization  %
%*******************%
beta = zeros(size(X,1),1);
u = 1;
du = beta;
ddu = eye(size(X,1));
IPPrecision=2*IPPrecisionTolerance;

%*******************%
%   the main loop   %
%*******************%

% repeat iterations until the desired precision in certainty equivalent is
% attained.
while (and(abs(IPPrecision) >= IPPrecisionTolerance, u >= 10^-6))

    beta  = beta - ddu\du;
    wealth = 1 - X*beta;                      
    wealth = max([wealth ; zeros(1,length(X))]); 
    u = (wealth.^2)*XDistr';
    du = -(X.*(wealth)) * XDistr';
    ddu = (X.*(wealth > 0)) * ((ones(size(X,1),1)*XDistr).*X)';
    if (ddu == 0)
      IPPrecision = 0;
    else
    % precision in certainty equivalent wealth 
    IPPrecision = 1/4 * du * (ddu\du) / (2-u);
    end 
end

if (u < 10^-6)
        IP = 1;
else
        IP = 1 - sqrt(u);
end