function [CE, alpha] = CRRAmax(X, XDistr, gama, CEqTol, outputflag)
% Copyright (c) 2001-2009 by Ales Cerny
% CRRAmax(X, XDistr, gama, CEqTol, outputflag) returns the
% certainty equivalent (CE) computed to the precision CEqTol and the 
% optimal proportion of risky investment in safe wealth for an agent with 
% CRRA utility and the coefficient of local relative risk aversion = gama. 
% Procedure assumes single risky investment with excess returns X 
% distributed with PDF XDistr. Unless outputflag=1 no output is printed.

%************************************************************************%
% CRRAmax.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 CRRAmax.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  @
% @*******************@
% local alpha,iteration,CEqPrecision,wealth,u,du,ddu,CEq;
alpha = 0; iteration=0;
CEqPrec=2*CEqTol;
if (nargin==4);
    outputflag = 0;
end;

% @*******************@
% @   the main loop   @
% @*******************@
% repeat iterations until the desired precision in certainty equivalent is
% attained.

while ( abs(CEqPrec) >= CEqTol)
    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';
    % precision in certainty equivalent wealth @
    CE = u^(1/(1-gama)); 
    CEqPrec = -1/2/(1-gama)*du^2/ddu/u;
    if (outputflag == 1)
        disp(sprintf('iteration                               %12.0f            ', iteration));
        disp(sprintf('Certainty equivalent/Safe wealth        %12.4f            ', CE));
        disp(sprintf('estimated precision in CE               %12.1e            ', CEqPrec));
        disp(sprintf('alpha                                   %12.4f            ', alpha));
        disp(sprintf('estimated precision in alpha            %12.1e            ', -du/ddu));
        disp(' ');
    end
    iteration=iteration+1;
    alpha = alpha - du/ddu;
end;