Code covered by the BSD License  

Highlights from
Modern Pricing Method using Transforms

image thumbnail

Modern Pricing Method using Transforms

by

 

25 Jul 2012 (Updated )

COS, CONV, Lewis Option Pricing Methods including Bermudan and American Options.

FFTCONV_B_Fast(n, L, alpha, cp, model, S0, t, r, q, ...
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
% 
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license 
%
% This code is being provided solely for information and general 
% illustrative purposes. The authors will not be responsible for the 
% consequences of reliance upon using the code or for numbers produced 
% from using the code. 



function price = FFTCONV_B_Fast(n, L, alpha, cp, model, S0, t, r, q, ...
    strike, Nex, varargin)
% n           -> points used for the grid construction, 2^n
% L           -> width of the interval approx the density
% alpha       -> dampening factor
% cp          -> call cp=1; put cp=-1;
% model       -> string for finding model from CF lib
% t           -> maturity
% r           -> riskfree rate
% q           -> dividend
% strike      -> strike of the option
% parameters  -> model parameters
% Nex         -> number of exercise dates

Ngrid = 2^n;                                    % num grid points
dt = t / Nex;                                   % time step
rdt = r * dt;                                   % riskfree times dt

Delta_y = L/Ngrid;                              % discrete y
Delta_x = Delta_y;                              % discrete x
Delta_u = 2 * pi / L;                           % discrete u

adj_y = log(strike/S0) ...
    - (ceil(log(strike/S0)/Delta_y) * Delta_y); % adjustment y-grid

% Construct the Grids
Grid_i = (0:Ngrid-1)';                          % index nums for grids
Grid_m = (-1).^(0:Ngrid-1)';                    % (-1)^Grid_i

y = adj_y + (Grid_i .* Delta_y) ...
    - (.5*Ngrid * Delta_y);                      % Adjusted y-grid
x = y;                                          % x-grid from y-grid

u = zeros(Ngrid, 1);                            % the base u-grid 
u = u + (Grid_i .* Delta_u) ...
    - (.5*Ngrid * Delta_u);                      % Set u-grid

w = ones(Ngrid,1); w(1) = 0.5; w(Ngrid) = 0.5;  % coefficients

V = max(cp .* (S0*exp(y) - strike), 0);         % option value at t_Nex
v = V .* exp(alpha .* y);                       % dampened option value

for m = Nex-1:-1:1                              % backward induction
    FT_Vec = ifft( (Grid_m) .* w .* v );% direct FFT for v
    FT_Vec_tbt = exp( 1i .* Grid_i .* (y(1) - x(1)) .* Delta_u ) ...
    .* exp(feval(@CF, model,-(u - (1i*alpha)), dt,r,q,varargin{:})) ...
    .* FT_Vec;                                 % vector to be transf.
    C = abs(exp(-rdt-(alpha .* x) + (1i .* u .* (y(1) - x(1))) ) ...
        .* (Grid_m) .* fft(FT_Vec_tbt));% continuation value
    h = max(cp .* (S0*exp(x) - strike), 0);     % payoff at t_m
    V = max(C, h);                              % option value at t_m    
  
    l = find(V==h,1,'last');                    % find exercise boundary
    if cp == 1                                  % case of a call           
        l = l-1;
    end
    xstar = ( x(l+1)*(C(l) - h(l)) - x(l)*(C(l+1) - h(l+1)) ) ...
        / ( (C(l) - h(l)) - (C(l+1) - h(l+1)) );% x^star
    adj_x = xstar-(ceil(xstar/Delta_x)*Delta_x);% adjust x-grid
    x = adj_x + (Grid_i .* Delta_x) ...
        - (.5*Ngrid * Delta_x);                  % new x-grid
    
    % calculate with new x-grid
    FT_Vec_tbt = exp( 1i .* Grid_i .* (y(1) - x(1)) .* Delta_u ) ...
    .* exp(feval(@CF, model,-(u - (1i*alpha)), dt,r,q,varargin{:})) ...
    .* FT_Vec;                                 % vector to be transf.

    % use new grid
    C = abs(exp( -rdt - (alpha .* x) + (1i .* u .* (y(1) - x(1))) ) ...
        .* (Grid_m) .* fft(FT_Vec_tbt)); % continuation
    V = max(C,max(cp .* (S0*exp(x) - strike), 0));% option value
    
    y = x;                                      % y-grid = x-grid
    
    v = V .* exp(alpha .* y);                   % dampened option value   
end

% Last Step
x = (Grid_i .* Delta_x)-(.5*Ngrid * Delta_x);   % final x-grid

% final transforms
FT_Vec = ifft( (Grid_m) .* w .* v );
FT_Vec_tbt = exp( 1i .* Grid_i .* (y(1) - x(1)) .* Delta_u ) ...
    .* exp(feval(@CF, model,-(u - (1i*alpha)), dt,r,q,varargin{:})) ...
    .* FT_Vec;                                 % vector to be transf.

C = abs(exp( -rdt - (alpha .* x) + (1i .* u .* (y(1) - x(1))) ) ...
    .* (Grid_m) .* fft(FT_Vec_tbt));    % final value

price = C(.5*Ngrid + 1, 1);                      % return price
end

Contact us