%% Workflow for Plotting an Option Sensitivity Surface Using the Merton76 Model
% Use optSensByMertonFFT to calibrate the FFT strike grid for sensitivities, compute option sensitivities, and plot option sensitivity surfaces.
%% Define Option Variables and Merton76 Model Parameters

AssetPrice = 80;
Rate = 0.03;
DividendYield = 0.02;
OptSpec = 'call';

Sigma = 0.16;
MeanJ = 0.02;
JumpVol = 0.08;
JumpFreq = 2;

%
%% Compute the Option Prices for the Entire FFT (or FRFT) Strike Grid, Without Specifying "Strike"
% Compute option sensitivities and also output the corresponding strikes. If the Strike input is empty ( [] ), option sensitivities will be computed on the entire FFT (or FRFT) strike grid. The FFT (or FRFT) strike grid is determined as exp(log-strike grid), where each column of the log-strike grid has NumFFT points with LogStrikeStep spacing that are roughly centered around each element of log(AssetPrice). The default value for NumFFT is 2^12. In addition to the sensitivities in the first output, the optional last output contains the corresponding strikes.

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, 6);
Strike = []; % Strike is not specified

[Delta, Kout] = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta");

% Show the lowest and highest strike values on the FFT strike grid
format
[Kout(1) Kout(end)]

% Show a subset of the strikes and corresponding option sensitivities
Range = (2046:2052);
[Kout(Range) Delta(Range)]

%
%% Change the Number of FFT (or FRFT) Points and Compare with optSensByMertonNI
% Try a different number of FFT (or FRFT) points, and compare the results with numerical integration. Unlike optSensByMertonFFT, which uses FFT (or FRFT) techniques for fast computation across the whole range of strikes, the optSensByMertonNI function uses direct numerical integration and it is typically slower, especially for multiple strikes. However, the values computed by optSensByMertonNI can serve as a benchmark for adjusting the settings for optSensByMertonFFT.

% Try a smaller number of FFT points 
% (e.g. for faster performance or smaller memory footprint)
NumFFT = 2^10; % Smaller than the default value of 2^12
Strike = []; % Strike is not specified (will use the entire FFT strike grid)
[Delta, Kout] = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT);

% Compare with numerical integration method
Range = (510:516);
Strike = Kout(Range);
DeltaFFT = Delta(Range);
DeltaNI = optSensByMertonNI(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta");
Error = abs(DeltaFFT-DeltaNI);
table(Strike, DeltaFFT, DeltaNI, Error)

%
%% Make Further Adjustments to FFT (or FRFT)
% If the values in the output DeltaFFT are significantly different from those in DeltaNI, try making adjustments to optSensByMertonFFT settings, such as CharacteristicFcnStep, LogStrikeStep, NumFFT, DampingFactor, and so on. Note that if (LogStrikeStep * CharacteristicFcnStep) is 2*pi/ NumFFT, FFT is used. Otherwise, FRFT is used.

Strike = []; % Strike is not specified (will use the entire FFT or FRFT strike grid)
[Delta, Kout] = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001);

% Compare with numerical integration method
Strike = Kout(Range);
DeltaFFT = Delta(Range);
DeltaNI = optSensByMertonNI(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta");
Error = abs(DeltaFFT-DeltaNI);
table(Strike, DeltaFFT, DeltaNI, Error)

% Save the final FFT (or FRFT) strike grid for future reference. For
% example, it provides information about the range of Strike inputs for
% which the FFT (or FRFT) operation is valid.
FFTStrikeGrid = Kout;
MinStrike = FFTStrikeGrid(1) % Strike cannot be less than MinStrike
MaxStrike = FFTStrikeGrid(end) % Strike cannot be greater than MaxStrike

%
%% Compute the Option Sensitivity for a Single Strike
% Once the desired FFT (or FRFT) settings are determined, use the Strike input to specify the strikes rather than providing an empty array. If the specified strikes do not match a value on the FFT (or FRFT) strike grid, the outputs are interpolated on the specified strikes.

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, 6);
Strike = 80;

Delta = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001)

%
%% Compute the Option Sensitivities for a Vector of Strikes
% Use the Strike input to specify the strikes.

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, 6);
Strike = (76:2:84)';

Delta = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001)

%
%% Compute the Option Sensitivities for a Vector of Strikes and a Vector of Dates of the Same Lengths
% Use the Strike input to specify the strikes. Also, the Maturity input can be a vector, but it must match the length of the Strike vector if the ExpandOutput name-value pair argument is not set to "true".

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, [12 18 24 30 36]); % Five maturities
Strike = [76 78 80 82 84]'; % Five strikes

Delta = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, ...
    'LogStrikeStep', 0.001) % Five values in vector output

%
%% Expand the Outputs for a Surface
% Set the ExpandOutput name-value pair argument to "true" to expand the outputs into NStrikes-by-NMaturities matrices. In this case, they are square matrices.

[Delta, Kout] = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, ...
    'LogStrikeStep', 0.001, 'ExpandOutput', true) % (5 x 5) matrix output

%
%% Compute the Option Sensitivities for a Vector of Strikes and a Vector of Dates of Different Lengths
% When ExpandOutput is "true", NStrikes do not have to match NMaturities. That is, the output NStrikes-by-NMaturities matrix can be rectangular.

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, 12*(0.5:0.5:3)'); % Six maturities
Strike = (76:2:84)'; % Five strikes

Delta = optSensByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, ...
    'LogStrikeStep', 0.001, 'ExpandOutput', true) % (5 x 6) matrix output

%
%% Compute the Option Sensitivities for a Vector of Strikes and a Vector of Asset Prices
% When ExpandOutput is "true", the output can also be a NStrikes-by-NAssetPrices rectangular matrix by accepting a vector of asset prices.

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, 12); % Single maturity
ManyAssetPrices = [70 75 80 85]; % Four asset prices
Strike = (76:2:84)'; % Five strikes

Delta = optSensByMertonFFT(Rate, ManyAssetPrices, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'OutSpec', "delta", ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, ...
    'LogStrikeStep', 0.001, 'ExpandOutput', true) % (5 x 4) matrix output

%
%% Plot Option Sensitivity Surfaces
% Use the Strike input to specify the strikes. Increase the value for NumFFT to support a wider range of strikes. Also, the Maturity input can be a vector. Set ExpandOutput to "true" to output the surfaces as NStrikes-by-NMaturities matrices.

Settle = datenum('29-Jun-2017');
Maturity = datemnth(Settle, 12*[1/12 0.25 (0.5:0.5:3)]');
Times = yearfrac(Settle, Maturity);
Strike = (2:2:200)';

% Increase 'NumFFT' to support a wider range of strikes
NumFFT = 2^13;

[Delta, Gamma, Rho, Theta, Vega] = optSensByMertonFFT(...
    Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, ...
    'NumFFT', NumFFT, 'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001, ...
    'OutSpec', ["delta", "gamma", "rho", "theta", "vega"], ...
    'ExpandOutput', true);

[X,Y] = meshgrid(Times,Strike);

figure;
surf(X,Y,Delta);
title('Delta');
xlabel('Years to Option Expiry');
ylabel('Strike');
view(-112,34);
xlim([0 Times(end)]);

figure;
surf(X,Y,Gamma)
title('Gamma')
xlabel('Years to Option Expiry')
ylabel('Strike')
view(-112,34);
xlim([0 Times(end)]);

figure;
surf(X,Y,Rho)
title('Rho')
xlabel('Years to Option Expiry')
ylabel('Strike')
view(-112,34);
xlim([0 Times(end)]);

figure;
surf(X,Y,Theta)
title('Theta')
xlabel('Years to Option Expiry')
ylabel('Strike')
view(-112,34);
xlim([0 Times(end)]);

figure;
surf(X,Y,Vega)
title('Vega')
xlabel('Years to Option Expiry')
ylabel('Strike')
view(-112,34);
xlim([0 Times(end)]);

% Copyright 2012 The MathWorks, Inc.