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 ([e-mail])
% Sri Krishnamurthy([e-mail])
% Copyright 2011 MathWorks, Inc.

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
[symGrad,symHess] = defineSymbolics();
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.GradObj = 'on';
                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);
        g = symGrad(x);
    end

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


%% Define symbolic objects
    function [symGrad,symHess] = defineSymbolics()
        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;
        
        %Define objective gradient and hessian
        grad = jacobian(obj,x);
        hess = jacobian(grad,x);
        
        %Transform symbolics into MATLAB function handles
        symGrad = matlabFunction(grad);
        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
    case 1 %Made up data
        % 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
        googData=xlsread('GoogleOrig.xlsx');
        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.
%
%  See also FIT, SFIT.

%  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