This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

optByMertonFFT

Option price by Merton76 model using FFT and FRFT

Syntax

[Price,StrikeOut] = optByMertonFFT(Rate,AssetPrice,Settle,Maturity,OptSpec,Strike,Sigma,MeanJ,JumpVol,JumpFreq)
[Price,StrikeOut] = optByMertonFFT(___,Name,Value)

Description

example

[Price,StrikeOut] = optByMertonFFT(Rate,AssetPrice,Settle,Maturity,OptSpec,Strike,Sigma,MeanJ,JumpVol,JumpFreq) computes vanilla European option price by Merton76 model, using Carr-Madan FFT and Chourdakis FRFT methods.

example

[Price,StrikeOut] = optByMertonFFT(___,Name,Value) adds optional name-value pair arguments.

Examples

collapse all

Use optByMertonFFT to calibrate the FFT strike grid, compute option prices, and plot an option price surface.

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 prices and also output the corresponding strikes. If the Strike input is empty ( [] ), option prices 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 prices 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 (will use the entire FFT strike grid)

% Compute option prices for the entire FFT strike grid
[Call, Kout] = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield);

% Show the lowest and highest strike values on the FFT strike grid
format
MinStrike = Kout(1) % Lowest possible strike in the current FFT strike grid
MinStrike = 2.9205e-135
MaxStrike = Kout(end) % Highest possible strike in the current FFT strike grid
MaxStrike = 1.8798e+138
% Show a subset of the strikes and corresponding option prices
Range = (2046:2052);
[Kout(Range) Call(Range)]
ans = 7×2

   50.4929   29.4645
   58.8640   21.2601
   68.6231   12.2218
   80.0000    4.5600
   93.2631    0.9579
  108.7251    0.1236
  126.7505    0.0113

Change the Number of FFT (or FRFT) Points and Compare with optByMertonNI

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

% Try a smaller number of FFT (or FRFT) 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)
[Call, Kout] = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT);

% Compare with numerical integration method
Range = (510:516);
Strike = Kout(Range);
CallFFT = Call(Range);
CallNI = optByMertonNI(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield);
Error = abs(CallFFT-CallNI);
table(Strike, CallFFT, CallNI, Error)
ans=7×4 table
    Strike     CallFFT       CallNI         Error  
    ______    _________    ___________    _________

    12.696       66.328         66.696      0.36786
    23.449       55.922         56.103      0.18071
    43.312       36.481         36.536     0.055233
        80       4.7387           4.56      0.17867
    147.76     0.046602      0.0008089     0.045793
    272.93    0.0092842    -7.0709e-08    0.0092842
    504.11    0.0024041    -2.4515e-07    0.0024044

Make Further Adjustments to FFT (or FRFT)

If the values in the output CallFFT are significantly different from those in CallNI, try making adjustments to optByMertonFFT 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)
[Call, Kout] = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001);

% Compare with numerical integration method
Strike = Kout(Range);
CallFFT = Call(Range);
CallNI = optByMertonNI(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield);
Error = abs(CallFFT-CallNI);
table(Strike, CallFFT, CallNI, Error)
ans=7×4 table
    Strike    CallFFT    CallNI      Error   
    ______    _______    ______    __________

    79.76      4.674      4.674    4.9664e-10
    79.84     4.6358     4.6358    4.9651e-10
    79.92     4.5978     4.5978    4.9642e-10
       80       4.56       4.56    4.9641e-10
    80.08     4.5224     4.5224    4.9642e-10
    80.16      4.485      4.485     4.965e-10
    80.24     4.4478     4.4478     4.966e-10

% 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
MinStrike = 47.9437
MaxStrike = FFTStrikeGrid(end) % Strike cannot be greater than MaxStrike
MaxStrike = 133.3566

Compute Option Price 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; 

Call = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001)
Call = 4.5600

Compute the Option Prices 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)';

Call = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001)
Call = 5×1

    6.7411
    5.5762
    4.5600
    3.6891
    2.9551

Compute the Option Prices 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

Call = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001) % Five values in vector output
Call = 5×1

    8.5589
    8.9439
    9.2316
    9.4653
    9.6565

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.

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

    8.5589    9.9675   11.1343   12.1492   13.0464
    7.4844    8.9439   10.1481   11.1939   12.1181
    6.5125    8.0023    9.2316   10.2999   11.2449
    5.6401    7.1402    8.3827    9.4653   10.4249
    4.8630    6.3545    7.5990    8.6881    9.6565

Kout = 5×5

    76    76    76    76    76
    78    78    78    78    78
    80    80    80    80    80
    82    82    82    82    82
    84    84    84    84    84

Compute the Option Prices 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

Call = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001, ...
    'ExpandOutput', true) % (5 x 6) matrix output
Call = 5×6

    6.7411    8.5589    9.9675   11.1343   12.1492   13.0464
    5.5762    7.4844    8.9439   10.1481   11.1939   12.1181
    4.5600    6.5125    8.0023    9.2316   10.2999   11.2449
    3.6891    5.6401    7.1402    8.3827    9.4653   10.4249
    2.9551    4.8630    6.3545    7.5990    8.6881    9.6565

Compute the Option Prices 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

Call = optByMertonFFT(Rate, ManyAssetPrices, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001, ...
    'ExpandOutput', true) % (5 x 4) matrix output
Call = 5×4

    3.4187    5.6579    8.5589   12.0417
    2.8538    4.8401    7.4844   10.7343
    2.3718    4.1205    6.5125    9.5230
    1.9635    3.4922    5.6401    8.4090
    1.6198    2.9476    4.8630    7.3921

Plot an Option Price Surface

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 surface as a NStrikes-by-NMaturities matrix.

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;

Call = optByMertonFFT(Rate, AssetPrice, Settle, Maturity, OptSpec, Strike, ...
    Sigma, MeanJ, JumpVol, JumpFreq, 'DividendYield', DividendYield, 'NumFFT', NumFFT, ...
    'CharacteristicFcnStep', 0.065, 'LogStrikeStep', 0.001, 'ExpandOutput', true);

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

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

Input Arguments

collapse all

Continuously compounded risk-free interest rate, specified as a scalar decimal value.

Data Types: double

Current underlying asset price, specified as numeric value using a scalar or a NINST-by-1 or NColumns-by-1 vector.

For more information on the proper dimensions for AssetPrice, see the name-value pair argument ExpandOutput.

Data Types: double

Option settlement date, specified as a NINST-by-1 or NColumns-by-1 vector using serial date numbers, date character vectors, datetime arrays, or string arrays. The Settle date must be before the Maturity date.

For more information on the proper dimensions for Settle, see the name-value pair argument ExpandOutput.

Data Types: double | char | datetime | string

Option maturity date, specified as a NINST-by-1 or NColumns-by-1 vector using serial date numbers, date character vectors, datetime arrays, or string arrays.

For more information on the proper dimensions for Maturity, see the name-value pair argument ExpandOutput.

Data Types: double | char | datetime | string

Definition of the option, specified as a NINST-by-1 or NColumns-by-1 vector using a cell array of character vectors or string arrays with values 'call' or 'put'.

For more information on the proper dimensions for OptSpec, see the name-value pair argument ExpandOutput.

Data Types: cell | string

Option strike price value, specified as a NINST-by-1, NRows-by-1, NRows-by-NColumns vector of strike prices.

If this input is an empty array ([]), option prices are computed on the entire FFT (or FRFT) strike grid, which is determined as exp(log-strike grid). Each column of the log-strike grid has'NumFFT' points with 'LogStrikeStep' spacing that are roughly centered around each element of log(AssetPrice).

For more information on the proper dimensions for Strike, see the name-value pair argument ExpandOutput.

Data Types: double

Volatility of the underling asset, specified as a scalar numeric value.

Data Types: double

Mean of the random percentage jump size (J), specified as a scalar decimal value where log(1+J) is normally distributed with mean (log(1+MeanJ)-0.5*JumpVol^2) and the standard deviation JumpVol.

Data Types: double

Standard deviation of log(1+J) where J is the random percentage jump size, specified as a scalar decimal value.

Data Types: double

Annual frequency of Poisson jump process, specified as a scalar numeric value.

Data Types: double

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: [Price,StrikeOut] = optByMertonFFT(Rate,AssetPrice,Settle,Maturity,OptSpec,Strike,Sigma,MeanJ,JumpVol,JumpFreq,'Basis',7)

Day-count of the instrument, specified as the comma-separated pair consisting of 'Basis' and a scalar using a supported value:

  • 0 = actual/actual

  • 1 = 30/360 (SIA)

  • 2 = actual/360

  • 3 = actual/365

  • 4 = 30/360 (PSA)

  • 5 = 30/360 (ISDA)

  • 6 = 30/360 (European)

  • 7 = actual/365 (Japanese)

  • 8 = actual/actual (ICMA)

  • 9 = actual/360 (ICMA)

  • 10 = actual/365 (ICMA)

  • 11 = 30/360E (ICMA)

  • 12 = actual/365 (ISDA)

  • 13 = BUS/252

For more information, see basis.

Data Types: double

Continuously compounded underlying asset yield, specified as the comma-separated pair consisting of 'DividendYield' and a scalar numeric value.

Data Types: double

Number of grid points in the characteristic function variable and in each column of the log-strike grid, specified as the comma-separated pair consisting of 'NumFFT' and a scalar numeric value.

Data Types: double

Characteristic function variable grid spacing, specified as the comma-separated pair consisting of 'CharacteristicFcnStep' and a scalar numeric value.

Data Types: double

Log-strike grid spacing, specified as the comma-separated pair consisting of 'LogStrikeStep' and a scalar numeric value.

Note

If (LogStrikeStep*CharacteristicFcnStep) is 2*pi/NumFFT, FFT is used. Otherwise, FRFT is used.

Data Types: double

Damping factor for Carr-Madan formulation, specified as the comma-separated pair consisting of 'DampingFactor' and a scalar numeric value.

Data Types: double

Type of quadrature, specified as the comma-separated pair consisting of 'Quadrature' and a single character vector or string array with a value of 'simpson' or 'trapezoidal'.

Data Types: char | string

Flag to expand the outputs, specified as the comma-separated pair consisting of 'ExpandOutput' and a logical:

  • true — If true, the outputs are NRows-by- NColumns matrices. NRows is the number of strikes for each column and it is determined by the Strike input. For example, Strike can be a NRows-by-1 vector, or a NRows-by-NColumns matrix. If Strike is empty, NRows is equal to NumFFT. NColumns is determined by the sizes of AssetPrice, Settle, Maturity, and OptSpec, which must all be either scalar or NColumns-by-1 vectors.

  • false — If false, the outputs are NINST-by-1 vectors. Also, the inputs Strike, AssetPrice, Settle, Maturity, and OptSpec must all be either scalar or NINST-by-1 vectors.

Data Types: logical

Output Arguments

collapse all

Option prices, returned as a NINST-by-1, or NRows-by-NColumns, depending on ExpandOutput.

Strikes corresponding to Price, returned as a NINST-by-1, or NRows-by-NColumns, depending on ExpandOutput.

More About

collapse all

Merton Jump Diffusion Model

The Merton jump diffusion model (Merton (1976)) is a different extension of the Black-Scholes model, where sudden asset price movements (both up and down) are modeled by adding the jump diffusion parameters with the Poisson process.

The stochastic differential equation is:

dSt=(rqλpμj)Stdt+σStdWt+JStdPtprob(dPt=1)=λpdt

where

r is the continuous risk-free rate.

q is the continuous dividend yield.

Wt is the Weiner process.

J is the random percentage jump size conditional on the jump occurring, where ln(1+J) is normally distributed with mean ln(1+μJ)δ22 and the standard deviation δ, and (1+J) has a lognormal distribution:

1(1+J)δ2πexp{[ln(1+J)(ln(1+μJ)δ22]2δ22}

μJ is the mean of J for (μJ > -1).

δ is the standard deviation of ln(1+J) for (δ≥ 0).

ƛp is the annual frequency (intensity) of Poisson process Ptfor (ƛp ≥ 0).

σ is the volatility of the asset price for (σ > 0).

The characteristic function fMerton76j(ϕ) for j = 1 (asset prices measure) and j = 2 (risk-neutral measure) is:

fMerton76j=fBSjexp(λpτ(1+μj)mj+12[(1+μj)iϕeδ2(mjiϕ+(iϕ)22)1]λpτμjiϕ)where for j=1,2:fBS1(ϕ)=fBS2(ϕi)fBS2(i)fBS2(ϕ)=exp(iϕ[lnSt+(rqσ22)τ]ϕ2σ22τ)m1=12,m2=12

where

ϕ is the characteristic function variable.

τ is the time to maturity (τ = T- t).

i is the unit imaginary number ( i2 = -1).

Carr-Madan Formulation

The Carr and Madan (1999) formulation is a popular modified implementation of Heston (1993) framework.

Rather than computing the probabilities P1 and P2 as intermediate steps, Carr and Madan developed an alternative expression so that taking its inverse Fourier transform gives the option price itself directly.

Call(k)=eαkπ0Re[eiukψ(u)]duψ(u)=erτf2(ϕ=(u(α+1)i))α2+αu2+iu(2α+1)Put(K)=Call(K)+KerτSteqτ

where

r is the continuous risk-free rate.

q is the continuous dividend yield.

St is the asset price at time t.

τ is time to maturity (τ = T-t).

Call(K) is the call price at strike K.

Put(K) is the put price at strike K

i is a unit imaginary number (i2= -1)

ϕ is the characteristic function variable.

α is the damping factor.

u is the characteristic function variable for integration, where ϕ = (u - (α+1)i).

f2(ϕ) is the characteristic function for P2.

P2 is the probability of St > K under the risk-neutral measure for the model.

To apply FFT or FRFT to this formulation, the characteristic function variable for integration, u, is discretized into NumFFT(N) points with the step size CharacteristicFcnStepu), and the log-strike k is discretized into N points with the step size LogStrikeStepk).

The discretized characteristic function variable for integration, uj(for j = 1,2,3,…,N), has a minimum value of 0 and a maximum value of (N-1) (Δu), and it approximates the continuous integration range from 0 to infinity.

The discretized log-strike grid, kn(for n = 1, 2, 3, N) is approximately centered around ln(St), with a minimum value of

ln(St)N2Δk

and a maximum value of

ln(St)+(N21)Δk

Where the minimum allowable strike is

Stexp(N2Δk)

and the maximum allowable strike is

Stexp[(N21)Δk]

As a result of the discretization, the expression for the call option becomes

Call(kn)=Δueαknπj=1NRe[eiΔkΔu(j1)(n1)eiuj[NΔk2ln(St)]ψ(uj)]wj

where

Δu is the step size of discretized characteristic function variable for integration.

Δk is the step size of discretized log-strike.

N is the number of FFT or FRFT points.

wj is the weights for quadrature used for approximating the integral.

FFT is used to evaluate the above expression if Δk and Δu are subject to the following constraint:

ΔkΔu=(2πN)

otherwise, the functions use the FRFT method described in Chourdakis (2005).

References

[1] Bates, D. S. “Jumps and Stochastic Volatility: Exchange Rate Processes Implicit in Deutsche Mark Options.” The Review of Financial Studies. Vol 9. No. 1. 1996.

[2] Carr, P., and D.B. Madan. “Option Valuation Using the Fast Fourier Transform.” Journal of Computational Finance. Vol 2. No. 4. 1999.

[3] Cont, R. and P. Tankov. Financial Modeling with Jump Processes. Chapman & Hall/CRC Press, 2004.

[4] Chourdakis, K. “Option Pricing Using Fractional FFT.” Journal of Computational Finance. 2005.

[5] Merton, R. “Option Pricing When Underlying Stock Returns are Discontinuous.” Journal of Financial Economics. Vol 3. 1976.

Introduced in R2018a