Forum Navigation: Select WHO'S ON: 07:25 AM Wilmott / BookMap Competition General Forum Technical Forum Economics Forum Numerical Methods Forum Trading Forum The Quantitative Finance FAQs Proje... Student Forum Book And Research Paper Forum Programming and Software Forum The Quantitative Finance Code Libra... Careers Forum Events Board Brainteaser Forum Off Topic Forum and Website Bugs and Suggesti...

 new topic
 search
 bookshop
 jobs board
 news
 help
 file share
 audio visual
 articles
 blogs
 wiki
 forums
 home
 login

 FORUMS > Numerical Methods Forum < refresh >
 Topic Title: Levenberg-Marquardt in matlab Created On Fri Jun 13, 08 03:54 PM Topic View: Branch View Threaded (All Messages) Threaded (Single Messages) Linear

stt106
Member

Posts: 53
Joined: Nov 2007

Fri Jun 13, 08 03:54 PM

I have got a matlab code for L-M algorithm but i can't apply it for the model i want to minimise.

The model i want to minimise is EWMA

sigma(n)= lambda * sigma(n-1) + (1- lambda)* R(n)^2

where everything is known except lambda.
sigma is volatility and R is the log-return of stock price.
I think the code is correct and well written. If anyone can give me a hand on how to apply the code to the above model, it will be highly appreciated.

function [xf, S, cnt] = LMFsolve(varargin)
% Solve a Set of Overdetermined Nonlinear Equations in Least-Squares Sense.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A solution is obtained by a Fletcher version of the Levenberg-Maquardt
% algoritm for minimization of a sum of squares of equation residuals.
%
% [Xf, Ssq, CNT] = LMFsolve(FUN,Xo,Options)
% FUN     is a function handle or a function M-file name that evaluates
%         m-vector of equation residuals,
% Xo      is n-vector of initial guesses of solution,
% Options is an optional set of Name/Value pairs of control parameters
%         of the algorithm. It may be also preset by calling:
%         Options = LMFsolve('default'), or by a set of Name/Value pairs:
%         Options = LMFsolve('Name',Value, ... ), or updating the Options
%                   set by calling
%         Options = LMFsolve(Options,'Name',Value, ...).
%
%    Name   Values {default}         Description
% 'Display'     integer     Display iteration information
%                            {0}  no display
%                             k   display initial and every k-th iteration;
% 'Jacobian'    handle      Jacobian matrix function handle; {@finjac}
% 'FunTol'      {1e-7}      norm(FUN(x),1) stopping tolerance;
% 'XTol'        {1e-7}      norm(x-xold,1) stopping tolerance;
% 'MaxIter'     {100}       Maximum number of iterations;
% 'ScaleD'                  Scale control:
%               value        D = eye(m)*value;
%               vector       D = diag(vector);
%                {[]}        D(k,k) = JJ(k,k) for JJ(k,k)>0, or
%                                   = 1 otherwise,
%                                     where JJ = J.'*J
% Not defined fields of the Options structure are filled by default values.
%
% Output Arguments:
% Xf        final solution approximation
% Ssq       sum of squares of residuals
% Cnt       >0          count of iterations
%           -MaxIter,   did not converge in MaxIter iterations

% Example:
% The general Rosenbrock's function has the form
%    f(x) = 100(x(1)-x(2))^2 + (1-x(1))^2
% Optimum solution gives f(x)=0 for x(1)=x(2)=1. Function f(x) can be
% expressed in the form

%    f(x) = f1(x)^2 + f2(x)^2,

% where f1(x) = 10(x(1)-x(2)) and f2(x) = 1-x(1).

% Values of the functions f1(x) and f2(x) can be used as residuals.
% LMFsolve finds the solution of this problem in 5 iterations.

%The more complicated problem sounds:
% Find the least squares solution of the Rosenbrock valey inside a circle
% of the unit diameter centered at the origin. It is necessary to build
% third function, which is zero inside the circle and increasing outside it.
% This property has, say, the next function:
%    f3(x) = sqrt(x(1)^2 + x(2)^2) - r, where r is a radius of the circle.

% Its implementation using anonymous functions has the form
%    R  = @(x) sqrt(x'*x)-.5;    %   A distance from the radius r=0.5
%    ros= @(x) [10*(x(2)-x(1)^2); 1-x(1); (R(x)>0)*R(x)*1000];
%    [x,ssq,cnt]=LMFsolve(ros,[-1.2,1],'Display',1,'MaxIter',50)
% Solution: x = [0.4556; 0.2059],  |x| = 0.5000
% sum of squares: ssq = 0.2966,
% number of iterations: cnt = 18.
%
% Note:
% Users with older MATLAB versions, which have no anonymous functions
% implemented, have to call LMFsolve with named function for residuals.
% For above example it is
%
%   [x,ssq,cnt]=LMFsolve('rosen',[-1.2,1]);
%
% where the function rosen.m is of the form
%
%   function r = rosen(x)
%%  Rosenbrock valey with a constraint
%   R = sqrt(x(1)^2+x(2)^2)-.5;
%%  Residuals:
%   r = [ 10*(x(2)-x(1))  %   first part
%         1-x(1)            %   second part
%         (R>0)*R*1000.     %   penalty
%       ];

% Reference:
% Fletcher, R., (1971): A Modified Marquardt Subroutine for Nonlinear Least
% Squares. Rpt. AERE-R 6799, Harwell

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   OPTIONS
%%%%%%%
%               Default Options
if nargin==1 && strcmpi('default',varargin(1))
xf.Display  = 0;         %   no print of iterations
xf.Jacobian = @finjac;   %   finite difference Jacobian approximation
xf.MaxIter  = 100;       %   maximum number of iterations allowed
xf.ScaleD   = [];        %   automatic scaling by D = diag(diag(J'*J))
xf.FunTol   = 1e-7;      %   tolerace for final function value
xf.XTol     = 1e-4;      %   tolerance on difference of x-solutions
return

%               Updating Options
elseif isstruct(varargin{1}) % Options=LMFsolve(Options,'Name','Value',...)
if ~isfield(varargin{1},'Jacobian')
error('Options Structure not Correct for LMFsolve.')
end
xf=varargin{1};          %   Options
for i=2:2:nargin-1
name=varargin{i};     %   option to be updated
if ~ischar(name)
error('Parameter Names Must be Strings.')
end
name=lower(name(isletter(name)));
value=varargin{i+1};  %   value of the option
if strncmp(name,'d',1), xf.Display  = value;
elseif strncmp(name,'f',1), xf.FunTol   = value(1);
elseif strncmp(name,'x',1), xf.XTol     = value(1);
elseif strncmp(name,'j',1), xf.Jacobian = value;
elseif strncmp(name,'m',1), xf.MaxIter  = value(1);
elseif strncmp(name,'s',1), xf.ScaleD   = value;
else   disp(['Unknown Parameter Name --> ' name])
end
end
return

%               Pairs of Options
elseif ischar(varargin{1})  % check for Options=LMFSOLVE('Name',Value,...)
Pnames=char('display','funtol','xtol','jacobian','maxiter','scaled');
if strncmpi(varargin{1},Pnames,length(varargin{1}))
xf=LMFsolve('default');  % get default values
xf=LMFsolve(xf,varargin{:});
return
end
end

%   LMFSOLVE(FUN,Xo,Options)
%%%%%%%%%%%%%%%%%%%%%%%%

FUN=varargin{1};            %   function handle
if ~(isvarname(FUN) || isa(FUN,'function_handle'))
error('FUN Must be a Function Handle or M-file Name.')
end

xc=varargin{2};             %   Xo

if nargin>2                 %   OPTIONS
if isstruct(varargin{3})
options=varargin{3};
else
if ~exist('options','var')
options = LMFsolve('default');
end
for i=3:2:size(varargin,2)-1
options=LMFsolve(options, varargin{i},varargin{i+1});
end
end
else
if ~exist('options','var')
options = LMFsolve('default');
end
end

x=xc(:);
lx=length(x);

r=feval(FUN,x);             % Residuals at starting point
%~~~~~~~~~~~~~~
S=r'*r;
epsx=options.XTol(:);
epsf=options.FunTol(:);
if length(epsx)<lx, epsx=epsx*ones(lx,1); end
J=options.Jacobian(FUN,r,x,epsx);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
A=J.'*J;                    % System matrix
v=J.'*r;

D = options.ScaleD;
if isempty(D)
D=diag(diag(A));        % automatic scaling
for i=1:lx
if D(i,i)==0, D(i,i)=1; end
end
else
if numel(D)>1
D=diag(sqrt(abs(D(1:lx)))); % vector of individual scaling
else
D=sqrt(abs(D))*eye(lx);     % scalar of unique scaling
end
end

Rlo=0.25; Rhi=0.75;
l=1;      lc=.75;     is=0;
cnt=0;
ipr=options.Display;
printit(-1);                %   Table header
d=options.XTol;             %   vector for the first cycle
maxit = options.MaxIter;    %   maximum permitted number of iterations

while cnt<maxit && ...      %   MAIN ITERATION CYCLE
any(abs(d)>=epsx) && ...    %%%%%%%%%%%%%%%%%%%%
any(abs(r)>=epsf)
d=(A+l*D)\v;            % negative solution increment
xd=x-d;
rd=feval(FUN,xd);
%   ~~~~~~~~~~~~~~~~~
Sd=rd.'*rd;
dS=d.'*(2*v-A*d);       % predicted reduction
R=(S-Sd)/dS;

if R>Rhi
l=l/2;
if l<lc, l=0; end
elseif R<Rlo
nu=(Sd-S)/(d.'*v)+2;
if nu<2
nu=2;
elseif nu>10
nu=10;
end
if l==0
lc=1/max(abs(diag(inv(A))));
l=lc;
nu=nu/2;
end
l=nu*l;
end
cnt=cnt+1;
if ipr>0 && (rem(cnt,ipr)==0 || cnt==1)
printit(cnt,[S,l,R,x(:).'])
printit(0,[lc,d(:).'])
end
if Sd>S && is==0
is=1;
St=S; xt=x; rt=r; Jt=J; At=A; vt=v;
end
if Sd<S || is>0
S=Sd; x=xd; r=rd;
J=options.Jacobian(FUN,r,x,epsx);
%       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
A=J.'*J;   v=J.'*r;
else
S=St; x=xt; r=rt; J=Jt; A=At; v=vt;
end
if Sd<S, is=0; end
end
xf = x;                         %   finat solution
if cnt==maxit, cnt=-cnt; end    %   maxit reached
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   PRINTIT     Printing of intermediate results
%   %%%%%%%
function printit(cnt,mx)
%        ~~~~~~~~~~~~~~~
%  cnt     =  -1  print out the header
%           0  print out second row of results
%          >0  print out first row of results

if ipr>0
if cnt<0                 %   table header
disp('')
disp(char('*'*ones(1,100)))
disp(['  cnt   SUM(r^2)        l            R' blanks(21) 'x(i) ...'])
disp([blanks(24) 'lc' blanks(32) 'dx(i) ...'])
disp(char('*'*ones(1,100)))
disp('')
else                     %   iteration output
if cnt>0 || rem(cnt,ipr)==0
f='%12.4e ';
form=[f f f f '\n' blanks(49)];
if cnt>0
fprintf(['%4.0f ' f f f ' x = '],[cnt,mx(1:3)])
fprintf(form,mx(4:length(mx)))
else
fprintf([blanks(18) f blanks(13) 'dx = '], mx(1))
fprintf(form,mx(2:length(mx)))
end
fprintf('\n')
end
end
end
end

%   FINJAC       numerical approximation to Jacobi matrix
%   %%%%%%
function J = finjac(FUN,r,x,epsx)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lx=length(x);
J=zeros(length(r),lx);
for k=1:lx
dx=.25*epsx(k);
xd=x;
xd(k)=xd(k)+dx;
rd=feval(FUN,xd);
%   ~~~~~~~~~~~~~~~~
J(:,k)=((rd-r)/dx);
end
end
end

 Reply Quote Top Bottom

GypCasino
Junior Member

Posts: 13
Joined: Jun 2008

Sat Jun 14, 08 12:35 AM

Before you try that routine, I would suggest using one of the built-in nonlinear optimizers in MATLAB.

Basically just write a .m file that accepts lambda as an argument and returns your objective function. And pass the name of that function and an initial guess for lambda to the optimizer. If that is working for you then you could start testing that coded L-M.
 Reply Quote Top Bottom

AvH
Junior Member

Posts: 20
Joined: Mar 2008

Sun Jun 15, 08 08:46 PM

The built-in function lsqnonlin (available in the optimization toolbox) can also use the LM-algorithm.

See http://www.mathworks.com/access/helpdesk/help/toolbox/optim/index.html?/access/helpdesk/help/toolbox/optim/ug/lsqnonlin.html
 Reply Quote Top Bottom

 FORUMS > Numerical Methods Forum < refresh >

 Forum Navigation: Select WHO'S ON: 07:25 AM Wilmott / BookMap Competition General Forum Technical Forum Economics Forum Numerical Methods Forum Trading Forum The Quantitative Finance FAQs Proje... Student Forum Book And Research Paper Forum Programming and Software Forum The Quantitative Finance Code Libra... Careers Forum Events Board Brainteaser Forum Off Topic Forum and Website Bugs and Suggesti...

© All material, including contents and design, copyright Wilmott Electronic Media Limited - FuseTalk 4.01 © 1999-2014 FuseTalk Inc. Terms & Conditions