Code covered by the BSD License

# Improving MATLABĀ® performance when solving financial optimization problems

### Jorge Paloschi (view profile)

Jorge Paloschi,PHD and Sri Krishnamurthy,CFA May 2011, http://www.wilmott.com/magazine.cfm

OptimizationWithSymbolicToolboxDemo()
```%% OptimizationWithSymbolicToolboxDemo
% In this example we will demonstrate solving an optimization problem using
% both numeric derivatives and analytical derivatives. We will also
% leverage parallel computing if the Parallel Toolbox is available.
% The Case study involves computing volatility estimates using the implied
% volatility method approach. Refer to our paper in the Wilmott magazine
% for more details.

% Authors:
% Jorge Paloschi (jpalosch@mathworks.com)
% Sri Krishnamurthy(skrishna@mathworks.com)

function OptimizationWithSymbolicToolboxDemo()
%% Problem parameters
% We have implemented two cases:
% 1. Creates Random Data sets
% 2. Uses Google Options Data as of March 7th 2011

[S,K,T,CP,r] = getProblemParameters(2);% parameter indicates case study
n = numel(S);
x0 = ones(n,1); %initial guess

% Pass constant parameters to numeric evaluation function
numObj(x0,S,K,T,CP,r); %#ok<NASGU>

%% Using parallel programming constructs.
% This sets it up also in other processors if parallel run
% To open MATLAB pools use
% matlabpool open <pool size> before running this script
parfor ii=1:matlabpool('size')
numObj(x0,S,K,T,CP,r); %#ok<NASGU>
end

%% Optimize using fd approximations
fprintf('\nFirst optimize using numeric derivatives\n');
optimize( false );

%% Using the Symbolic Toolbox, obtain the Symbolic gradient and hessian
tic
toc

%% Optimize using analytical information
fprintf('\nNow optimize using analytical derivatives\n');
optimize( true );

%--------------------------------------------------------------------------------

%% Helper functions
%% function: optimize
% Input : Flag: true OR false ;
% true indicates optimize using analytical derivatives
% false indicates optimize using numeric derivatives
function x = optimize(useSymb)
% Optimization method parameters
opts = optimset('Algorithm','interior-point');
opts.MaxFunEvals = 100000;
opts.Display = 'off';
opts.UseParallel = 'always';
switch useSymb
case true %Using symbolic data
objFunc = @symObj;
opts.Hessian = 'user-supplied';
opts.HessFcn = @hess;
Title = 'Using symbolic data';
case false %Using numerical data
objFunc = @numObj;
opts.Hessian = 'bfgs';
Title = 'Using numerical data';
end

%% Solve optimization problem using fmincon
tic %Start measuring time
lb = zeros(n,1);
ub = 10*ones(n,1);
[x,f,exitFlag,output] = fmincon(objFunc,x0,[],[],[],[],lb,ub,[],opts);
toc %Stop measuring time
%% Display results
fprintf('Exit flag=%d ||x||=%g obj=%g, %d func eval\n',...
exitFlag,norm(x),f,output.funcCount);
% Display surface fit graph
createSurfaceFit(K, T, x, Title);

end

function [o,g] = symObj(x)
o = numObj(x);
end

function H = hess(x,varargin)
H = symHess(x);
end

%% Define symbolic objects
reset(symengine);
n = numel(S);

%Define symbolic residuals
x = evalin(symengine,['n:=',num2str(n),';[x[j] \$ j = 1..n]']);
d1 = (log(S./K) + (r+x.'.^2/2).*T)./(x.'.*sqrt(T));
d2 = d1 - x.'.*sqrt(T);
v = S*0.5.*(1+erf(d1/sqrt(2))) - ...
K.*exp(-r*T) * 0.5.*(1+erf(d2/sqrt(2)));
f = v - CP;
obj = f.'*f;

%Transform symbolics into MATLAB function handles
symHess = matlabFunction(hess);
end
end

%% function numObj : Define numerical residual evaluation
% This is needed to make the code more efficient when run in parallel
% Rather than passing them as parameters or inheriting from workspace
% constant parameters are kept in static workspace.
% The function also returns the handle to the objective function to be
% evaluated
function [obj] = numObj(x , varargin)

persistent S;
persistent K;
persistent T;
persistent CP;
persistent r;
if nargin>1
S = varargin{1};
K = varargin{2};
T = varargin{3};
CP = varargin{4};
r = varargin{5};
end

n = numel(x);
f = zeros(1,n);
for jj=1:n
d1 = (log(S(jj)/K(jj)) + (r+x(jj)^2/2)*T(jj))/(x(jj)*sqrt(T(jj)));
d2 = d1 - x(jj)*sqrt(T(jj));
val = S(jj)*normcdf(d1) - K(jj)*exp(-r*T(jj)) * normcdf(d2);
f(jj) = (val-CP(jj));
end
obj = f*f';
end

%% function getProblemParameters
% Input:Flag 1 or 2
% Outputs
% S = Current Share price
% K = Strike price
% T = Time to Maturity
% CP = Call price
% r = Risk Free Rate

function [S,K,T,CP,r] = getProblemParameters( icase )
r = 0.06; % risk free IR
switch icase
% Ensure reproducibility of results
stream = RandStream('mt19937ar','Seed',1);
RandStream.setDefaultStream(stream);
S = [85:3:115,85:3:115,85:3:115]; % Share price
K = [90:2:110,90:2:110,90:2:110]; % Strike price
T = [0.5*ones(1,11),0.75*ones(1,11),1*ones(1,11)]; % Maturities
n = numel(S); %size
volatility = 0.2*rand(1,n);
CP = zeros(n,1);
for ii=1:n
CP(ii) = blsprice(S(ii), K(ii), r, T(ii), volatility(ii)); % Call prices using MATLAB
end
case 2
%% Sample Call option data of Google from Yahoo as of March 7 2011
T = googData(:,1);
K = googData(:,2);
CP = googData(:,3);
n = numel(K);
S= repmat(591.66, 1,n); % Closing price as of March 7 2011
otherwise
error('Not implemented');
end
S = S(:);
T = T(:);
K = K(:);
CP = CP(:);
end

function [fitresult, gof] = createSurfaceFit(K, T, x, Title)
%CREATESURFACEFIT(K,T,X)
%  Fit surface to data.
%
%  Data for 'untitled fit 1' fit:
%      X Input : K
%      Y Input : T
%      Z Output: x
%  Output:
%      fitresult : an sfit object representing the fit.
%      gof : structure with goodness-of fit info.
%

%  Auto-generated by MATLAB on 08-Mar-2011 09:37:39

%% Fit: 'untitled fit 1'.
[xInput, yInput, zOutput] = prepareSurfaceData( K(:), T(:), x(:) );

% Set up fittype and options.
ft = 'cubicinterp';
opts = fitoptions( ft );
opts.Normalize = 'on';

% Fit model to data.
[fitresult, gof] = fit( [xInput, yInput], zOutput, ft, opts );

% Plot fit with data.
figure( 'Name', [Title,':Fitted Volatility'] );
h = plot( fitresult, [xInput, yInput], zOutput );
legend( h, 'Fitted Vol', 'Vol vs. K, T', 'Location', 'NorthEast' );
% Label axes
xlabel( 'Strike' );
ylabel( 'Maturity' );
zlabel( 'Vol' );
grid on
view( -28, 35 );
end```