MATLAB Examples


Compute normalized probability density function from characteristic function. Part of the CFH Toolbox.


[F X] = CF2PDF(CF)

Given a discounted characteristic function CF, returns normalized probability density function F and corresponding sampling points X.

Input Arguments

The characteristic function CF should expect the real argument u and return the corresponding discounted characteristic function.

AUX is a structure containing optional parameters for the FFT:

  • aux.N is the number of points for FRFT evaluation, default 8192
  • aux.uMax denotes the upper limit of integration of the characteristic function, default 200
  • aux.dx describes discretization of the log strike range, default value 3/N
  • aux.x0 is the log of spot underlying, default zero
  • aux.x is a vector of pdf sampling points. If this field is used, the values in dx, x0 and N are replaced with corresponding values obtained from x.


Example 1: Black Scholes

In the Black Scholes model, the dynamics of the logarithmic spot process are:

$dy=(r_f - \frac{1}{2}\sigma^2)dt + \sigma dW$

resulting in the characteristic function

$\phi(u)=E(e^{uiX_T})=\exp(-r_f\tau +iux_t+ iu\tau(r_f-\frac{1}{2}\sigma^2)-\frac{1}{2}\tau u \sigma^2)$

where $\tau=T-t$.

The characteristic function of the Black Scholes model is also included in cflib, using the argument type='BS'.

rf          = 0.05;
tau         = 1;
sigma       = 0.25;
S0          = 1;
x0          = log(S0);
cfB= @(u) exp(-rf*tau + u*i*x0 + i*u*(rf-1/2*sigma^2)*tau - 1/2*u.^2*sigma^2*tau);
[fB xB] = cf2pdf(cfB);

Example 2: Heston's Stochastic volatility model with Merton style jumps

In this framework, the process dynamics are:

$dy = (r_f - \frac{1}{2}v_t-\lambda m)dt + \sqrt{v_t}dW_1 + JdZ$

$dv = \kappa(\theta-v_t)dt + \sigma\sqrt{v_t}dW_2$

$E(dW1dW2)=\rho dt$

and $dZ$ is a Poisson jump process with constant intensity $\lambda$.

The drift adjustment is $m = E(e^J-1)$, where the jump distribution is normal $N(\mu_J,\sigma_J)$

The corresponding characteristic function is part of the cflib toolbox:

par.rf      = 0.05;
par.q       = 0;
par.kappa   = 0.85;
par.theta   = 0.25^2;
par.sigma   = 0.10;
par.rho     = -0.8;
par.lambda  = 0.1;
par.muJ     = -0.20;
par.sigmaJ  = 0.10;
par.x0      = 0;
par.v0      = 0.25^2;
tau         = 1;
cfH         = @(u) cflib(u,tau,par,'HestonJump');
[fH xH]     = cf2pdf(cfH);
legend('Heston','Black Scholes');

Example 3: Probability density of relative asset price

In this example, we will discuss the probability density function of the relative performance of asset $X_T$ over asset $Y_T$, whose log dynamics are assumed to be:

$dx = (r_f -\frac{1}{2}\sigma_x^2 -\lambda m_x)dt + \sigma_x dW_1 + J_x dZ$

$dy = (r_f -\frac{1}{2}\sigma_y^2 \lambda m_y)dt + \sigma_y dW_2 + J_y dZ$

and constant jump intensity $\lambda$, whereas jumps in $x$ and $y$ are bivariate normal:

$[J_x , J_y]^T~\sim N(\mu,\Sigma)$

The jump transform is

$\theta(c)=\int_{R^n}exp(cJ)df(J)=\exp(c^T\mu_J + \frac{1}{2}c^T\Sigma_Jc)$

resulting in the drift adjustments

$m_x = \theta([1,0]^T)-1, m_y=\theta([0,1]^T)-1$

We are interested in the pdf, and hence the characteristic function, of of $Z_T=X_T/Y_T$. Noting that $Z_T=\exp(z_T)=\exp(x_T-y_T)$, we can rewrite the characteristic function of $Z_T$ as:

$\phi(u) = E(\exp(iuz_T))=E(\exp(iu(x_T-y_T)))$

This feat can be easily implemented via cfaffine. Let us assume correlated jumps of opposing sizes:

X0          = 100;
Y0          = 110;
rf          = 0.05;
sigmaX      = 0.20;
sigmaY      = 0.25;
muJ         = [-0.25 ; 0.15];
SigmaJ      = [0.2*0.2 0.2*0.1*0.7 ; 0.2*0.1*0.7 0.1*0.1];
lambda      = 0.15;
jump        = @(c) exp(muJ'*c + 1/2*diag(SigmaJ)'*c.^2 + SigmaJ(1,2)*c(1,:).*c(2,:));

Translating these into the AJD coefficients:

x0          = log([X0 ; Y0]);
z0          = x0(1)-x0(2);
m(1,:)      = jump([1;0])-1;
m(2,:)      = jump([0;1])-1;
K0          = rf-1/2*[sigmaX^2;sigmaY^2]-lambda*(m);
H0          = [sigmaX^2 0 ; 0 sigmaY^2];
L0          = lambda;
R0          = rf;

We will now define the characteristic function via cfaffine, which allows for an array of u values. Foreknowing that the function cf2pdf will evaluate the supplied characteristic functions at the column vector u=[0;...;uMax], we tell cfaffine that the input will be a NU x NX matrix by setting the ND argument to 2:

cf          = @(u) cfaffine(u*[1 -1],x0,tau,K0,[],H0,[],R0,[],L0,[],jump,2);

We will also compare the result to the case where no jumps are present:

K0Diffusive = rf-1/2*[sigmaX^2;sigmaY^2];
cfDiffuse   = @(u) cfaffine(u*[1 -1],x0,tau,K0Diffusive,[],H0,[],R0,[],[],[],[],2);

The initial value for cf2pdf is log(Z0), thus

aux.x0      = z0;
[f x]       = cf2pdf(cf,aux);
[fD xD]     = cf2pdf(cfDiffuse,aux);
title('Probability distribution of relative asset price, with and without jumps');
legend('pdf with jumps','pdf without jumps');
xlabel('log price ratio');
ylabel('normalized pdf');
xlim(1.5*[-1 1]);