Implementing Deep Smoothing for Implied Volatility Surfaces (IVS)

This project is a python-based implementation of the methodologies presented in the paper Deep Smoothing of the Implied Volatility Surface by Ackerer et al (2020) ¹. A volatility surface is a representation of implied volatility across different strike prices and maturities, crucial for pricing and hedging options accurately. Ensuring the surface is arbitrage-free is essential to prevent inconsistencies that could lead to riskless profit opportunities, which would undermine the model’s reliability in real markets.

Table of Contents

• Introduction
• Notation and Initial Values
• Bates Model
  • Characteristic Function
  • Fast Fourier Transform
• Loss Function
  • Calendar and Butterfly
  • Large Moneyness
  • Around the Money
  • Total Loss
• ATM Total Variance
• Surface Stochastic Volatility Inspired (SSVI)
• Neural Network
• Convergence
• Results
• Final Remarks
• Appendix: Arb-Free Volatility Surfaces
• Citation
• References

Introduction

Deep smoothing focuses on applying deep learning methods to generate smooth, arbitrage-free implied volatility surfaces. For someone unfamiliar with quantitative finance, this problem can be summarized as follows: Imagine you are given a set of points \((k, \tau, iv)\) representing market data. These points form part of a 3D surface that we aim to construct with two key objectives:

• We need to fit the market data.
• We must ensure that the surface meets certain curvature properties. In quantitative finance, this means creating an arbitrage-free surface. Mathematically, this translates to minimizing a specified loss function across the entire surface.

Deep smoothing’s main idea is to use feedforward neural networks as a corrective tool to modify the prior models considered for volatility surfaces (See ² to learn more about prior models). By letting

we enrich the space of parameters used for fitting the volatility surface. Here \(\omega\) is the implied variance and \(\theta_{\text{prior}}\) and \(\theta_{\text{nn}}\) are two disjoint set of parameters. The figure below displays a neural network with 4 hidden layers, each containing 40 units. The weights of this network are represented by \(\theta_{\text{nn}}\).

To ensure our volatility surface is free of arbitrage, we use the ideas from Arbitrage-free SVI Volatility Surfaces (2013) ² where they argue that if the following are satisfied then the call price surface is free of Calendar & Butterfly arbitrage resp.

In this blog, we’ll walk through the main steps of the deep smoothing framework. You can find the code for this project here.

Notation and Initial Values

European call options with the following table of notation and values are used:

Table: Parameters and Definitions

Bates model

Bates dynamic is displayed below.

\(dW_{1t}dW_{2t} = \rho dt\) and \(N_t\) is a compound Poisson process with intensity \(\lambda\) and independent jumps \(J\) with

Characteristic Function

The characteristic function of the log-strike in the Bates model is given by:

We will use this characteristic function to price options in Bates model. To this end, we utilize FFTs.

Fast Fourier Transform

The option price calculation for the Bates model can be efficiently computed using an analytical formula that leverages FFT. Following Carr and Madan ³, we apply a smoothing technique to be able to compute the FFT integral. Recall that the option value is given by

where

Since the integrand is singular at the required evaluation point \(u = 0\), FFT cannot be used to evaluate call price \(C_\tau(k)\). To offset this issue, we consider the modified call price \(c_\tau(k) := \exp(\alpha k) C_\tau(k)\) for \(\alpha > 0\). Denote the Fourier transform of \(c_\tau(k)\) by

It can be shown that

We set up the FFT calculation as follows:

• Log strike levels range from \(-b\) to \(b\) where \(b = \frac{Ndk}{2}\)

• \(\Psi_\tau(u)\) is computed at the following \(v\) values: \(v_j = (j-1)du \text{ for } j = 1, \cdots, N\)

• Option prices are computed at the following \(k\) values:

• To apply FFT, we need to set \(dk \cdot du = \frac{2\pi}{N}\)

• Simpson weights are used: \(3 + (-1)^j - \delta_{j-1} \text{ for } j = 1, \cdots, N\)

Having this setup ready, call prices are obtained as follows:

The Python snippet below computes option prices using the Fast Fourier Transform (FFT). We use np.vectorize multiple times to vectorize the computations across various maturities and the corrective parameter, \(\alpha\).

def cf_bates(tau,u):


    u2, sigma2, alpha2, ju = u**2,sigma**2, alpha**2, 1j*u

    def log_phi_j():
        term1 = -0.5*u2*alpha2+ju*(np.log(1+beta)-0.5*alpha2)
        term2 = np.exp(term1)-1
        return tau*lam*term2

    a = kappa-ju*rho*sigma
    gamma = np.sqrt((sigma2)*(u2+ju)+a**2)
    b = 0.5*gamma*tau
    c = kappa*theta/(sigma2)

    coshb, sinhb = np.cosh(b), np.sinh(b)
    
    term1 = c*tau*a+ju*(tau*(rate-lam*beta)+np.log(spot))
    term2 = coshb+(a/gamma)*sinhb
    term3 = (u2+ju)*v0/(gamma*(coshb/sinhb)+a)
    res = log_phi_j()+term1-term3-2*c*np.log(term2)
    
    return np.exp(res)

N = 2**14  # number of knots for FFT
B = 1000  # integration limit

alpha_corrs = np.array([[[0.4,0.5,0.7,1]]])
n_a = alpha_corrs.shape[2]

du = B / N
u = np.array(np.arange(N)* du).reshape(1,-1,1) 

w_fft = 3 + (-1) ** (np.arange(N) + 1)   # Simpson weights
w_fft[0] = 1 # [1, 4, 2, ..., 4, 2, 4]
dk = 2*np.pi/(du*N) # dk * du = 2pi/N
upper_b = N * dk / 2
lower_b = -upper_b
kus = lower_b + dk * np.arange(N)

taus_0_np = np.array(taus_0).reshape(-1,1,1) #n_t*1*1
w_fft = w_fft.reshape(1,-1,1) # 1*N*1
kus = kus.reshape(1,-1,1) # 1*N*1

cust_interp1d = lambda x,y: interp1d(x, y, kind='linear')
fn_vec = np.vectorize(cust_interp1d, signature='(n),(n)->()')

term1 = np.exp(-rate*taus_0_np)/(alpha_corrs**2+alpha_corrs-u**2+1j*(2*alpha_corrs+1)*u)
term2 = cf_bates(taus_0_np,u-(alpha_corrs+1)*1j) 
mdfd_cf = term1*term2

integrand = np.exp(1j * upper_b * u)*mdfd_cf*w_fft* du / (3*np.pi)
vectorized_fft = np.vectorize(lambda i, j: fft(integrand[i,:, j]), signature='(),()->(n_k)')
tau_indices, corr_indices = np.indices((n_t, n_a))
integrand_fft = vectorized_fft(tau_indices, corr_indices).transpose(0, 2, 1)

Ck_u =np.exp(-alpha_corrs*kus)*np.real(integrand_fft )

f1 = Ck_u.transpose(1, 0, 2).reshape(N, -1)
fn = fn_vec([kus[0,:,0] for _ in range(n_t*n_a)],[f1[:,i] for i in range(n_t*n_a)])

spline = fn.reshape(n_t, n_a) 
vectorized_spline = np.vectorize(lambda i, j, k: spline[i, j](ks[k]), signature='(),(),()->()')

tau_indices, corr_indices, ks_indices  = np.indices((n_t, n_a, n_k))
prices = vectorized_spline(tau_indices, corr_indices, ks_indices)
prices_ = prices.copy()
prices = np.median(prices,axis=1)

After computing prices, we use Quantlib library to obtain the volatility surface. We are now ready to explain the loss functions through which the arb-free conditions are enforced!

Loss Function

We define four different loss functions and construct the total loss function as a linear combination of these four, with coefficients being the penalty parameters. The first loss function is the prediction error, which is defined as below:

Here, \(\mathcal{I}_0\) is the set of log-moneyness, implied volatility, and maturities for each observed market option. For future use, denote

The other three loss functions are auxiliary, and consequently, we need to introduce auxiliary maturities and log-moneyness values. This ensures desired features for the constructed implied variance \(\omega\) across unobserved market data. Define

where, for example, \(k_{\max} = k_{\max}(\mathcal{I}_0)\). Note that we consider more monyness around the money.

Calendar and Butterfly

Calendar and Butterfly loss functions are then defined as

Large Moneyness

For large moneyness behavior, we set

where \(\mathcal{I}_{asym} = \left\{ (k, \tau) : k \in \{ 6k_{\min}, 4k_{\min}, 4k_{\max}, 6k_{\max} \}, \tau \in \mathcal{T}_{Aux} \right\}\)

Around the Money

We prefer that the prior component gives the best possible fit for ATM, while the neural network corrects the prior’s limitations for OTM. To this end, the following loss function is considered:

Here \(\mathcal{I}_{\text{atm}} = \{(0, \tau) : \tau \in \mathcal{T}_{Aux}\}\).

Total Loss

The total loss function is constructed as follows:

Regularization parameters \(\lambda_{but}, \lambda_{cal}, \lambda_{atm}\) are tunable. In our experiments, we consider:

Table: Loss Function Regularization Parameters

Python snippet below shows how loss function is computed given a pair of prior and neural network model i.e., \(\theta_{\text{prior}}\) and \(\theta_{nn}\), respectively.

def loss_butcal(model_prior, model_nn, ks, ts):
    
    kts = torch.stack((ks, ts), dim=1)
    output = model_prior(ks, ts)*model_nn(kts).squeeze(1)
    
    dw = torch.autograd.grad(outputs=output, inputs=[ks,ts], grad_outputs= torch.ones_like(output), create_graph=True, retain_graph=True)
    dwdk_,dwdtau_ = dw
    d2wdk_ = torch.autograd.grad(outputs=dwdk_, inputs=ks, grad_outputs=torch.ones_like(dwdk_), create_graph=True)[0]

    loss_but = relu(-((1-ks*dwdk_/(2*output))**2-(dwdk_/4)*(1/output+0.25)+d2wdk_/2)).mean()
    loss_cal = relu(-dwdtau_).mean()

    return loss_but, loss_cal 

def loss_large_m(model_prior, model_nn, ks, ts):

    kts = torch.stack((ks,ts), dim=1)
    output = model_prior(ks, ts)*model_nn(kts).squeeze(1)
    dw= torch.autograd.grad(outputs=output, inputs=ks, grad_outputs=torch.ones_like(output), create_graph=True)
    dwdk = dw[0]
    grad_output_dwdk = torch.ones_like(dwdk)
    d2wdk = torch.autograd.grad(outputs=dwdk, inputs=ks, grad_outputs=grad_output_dwdk, create_graph=True)[0]
    loss_large_m = torch.abs(d2wdk).mean()

    return loss_large_m

def loss_atm(model_nn, ks, ts):
    
    kts = torch.stack((ks,ts), dim=1)
    output = model_nn(kts).squeeze(1)
    loss_atm = torch.sqrt(torch.pow(1-output,2).sum()/output.shape[0])
    
    return loss_atm

def loss_0(model_prior, model_nn):

    output = model_prior(df_torch_gpu[:,2], df_torch_gpu[:,3])*(model_nn(df_torch_gpu[:,2:4]).squeeze(1))
    sigma_theta = torch.sqrt(output/df_torch_gpu[:,3])
    temp = sigma_theta-df_torch_gpu[:,4]
    rmse = torch.sqrt(torch.pow(temp,2).mean())
    mape = (torch.abs(temp)/df_torch_gpu[:,4]).mean()

    return rmse, mape 

def total_loss(model_prior, model_nn):

    ###########  Butterfly & Calendar Loss function ###############

    k_but_cal   = torch.tensor(ext(0,I_butcal), requires_grad=True).to(device)
    t_but_cal   = torch.tensor(ext(1,I_butcal), requires_grad=True).to(device) 
    l_but, l_cal = loss_butcal(model_prior, model_nn,k_but_cal, t_but_cal)
    
    ###########  Large moneyness loss function ################

    k_large_m  = torch.tensor(ext(0,I_large_m), requires_grad=True).to(device)
    t_large_m  = torch.tensor(ext(1,I_large_m), requires_grad=True).to(device)
    l_large_m = loss_large_m(model_prior, model_nn, k_large_m, t_large_m)
    
    ###########  ATM loss function ###############

    k_atm = torch.tensor(ext(0,I_atm)).to(device)
    t_atm = torch.tensor(ext(1,I_atm)).to(device)
    l_atm = loss_atm(model_nn, k_atm, t_atm)
    
    ###########  L0 loss function ###############

    rmse, mape  = loss_0(model_prior, model_nn)

    ###########  Total loss ###############

    return rmse, mape, l_cal,l_but,l_large_m,l_atm

ATM Total Variance

The prior model is expected to provide a first-order approximation of the volatility surface, making it essential for the prior to accurately reproduce ATM (at-the-money) values. Since the ATM term structure is observable through market data for \(\tau \in \mathcal{T}_0\), we construct the ATM variance for each \(\tau \in \mathcal{T}_{aux}\) as described below.

Below are the steps to construct ATM variance:

• Step A: Ensure that \(\omega_{atm}(\tau_2) \geq \omega_{atm}(\tau_1)\) for every \(\tau_2 \geq \tau_1\). To achieve this, we solve the following optimization problem using quadratic programming (QP) with CVX:

• Step B: Interpolate results from Step A using scipy with a smoothing spline of degree 3

• Step C: Fit a degree-five polynomial to the result from Step B. This approach allows access to both \(\omega_{atm}(\tau)\) and its derivative \(\frac{\partial}{\partial \tau} \omega_{atm}(\tau)\) for each \(\tau \in \mathcal{T}_{aux}\), which is useful during the training process.

Python snippet below shows how ATM total variance is constructed.

import cvxpy as cv
from torch.autograd import Function
from scipy.optimize import curve_fit
from scipy.interpolate import UnivariateSpline

spot_variances = ... ## a numpy array containing atm variance 
z = cv.Variable(len(spot_variances))
objective = cv.Minimize(cv.sum_squares(z - spot_variances))
constraints = [z[i] >= z[i-1]+1e-10 for i in range(1,len(spot_variances))]
prob = cv.Problem(objective, constraints)
prob.solve(verbose=False)
spot_variances_mdfd= z.value

spline_curve = UnivariateSpline(spot_taus, spot_variances_mdfd, k=5)
spline_curve.set_smoothing_factor(0.5)

def poly(x, *coeffs):
    return sum(c * x**i for i, c in enumerate(coeffs))

def poly_derivative(x, *coeffs):
    return sum(i * c * x**(i-1) for i, c in enumerate(coeffs) if i > 0)

degree = 5
popt, _ = curve_fit(poly, taus_aux,  spline_curve(taus_aux), p0=np.ones(degree+1))
popt_tensor = torch.tensor(popt).to(device)

# The at-the-money (ATM) variance and its derivatives are used as part of the neural network training process.
# To facilitate this, we define the following custom spline class, SPLFunction, which inherits from PyTorch's Function class.
class SPLFunction(Function):
    @staticmethod
    def forward(ctx, tau):
        ctx.save_for_backward(tau)
        return poly(tau,*popt_tensor)

    @staticmethod
    def backward(ctx, grad_output):
        tau, = ctx.saved_tensors
        grad_input = poly_derivative(tau,*popt_tensor)*grad_output
        return grad_input
    
def apply_spline(tau): return SPLFunction.apply(tau)

Surface Stochastic Volatility Inspired (SSVI)

The following version of the SSVI prior model (with power-law parameterization) has been used for the entire volatility surface.

Here \(\omega_{\text{atm}}(\tau)\) is the ATM term structure, and \(\rho\), \(\gamma\), and \(\eta\) are tunable parameters. In other words, we have:

The Python snippet below shows how Wssvi is constructed. The function apply_spline, used here, is defined in the ATM Variance section above.

class Wssvi(nn.Module):
    def __init__(self):
        super().__init__()
        self.initialize_weights()

    def initialize_weights(self):
        
        self.gamma = nn.Parameter(torch.tensor(0.5))
        self.eta   = nn.Parameter(torch.tensor(1.0))
        self.rho   = nn.Parameter(torch.tensor(0.0))

    def forward(self, k, tau):

        t1 = apply_spline(tau)
        phi = self.eta*torch.pow(t1,-self.gamma**2)*torch.pow(1+t1,self.gamma**2-1)
        t2 = 1+self.rho*phi*k+torch.sqrt((phi*k+self.rho)**2+1-self.rho**2)
        return t1*t2/2 

Neural Network

For the feed-forward neural network, we use 4 layers, each with 40 units. The model is trained using the Adam optimizer in conjunction with the optim.lr_scheduler.ReduceLROnPlateau scheduler. To escape local minima, we reinitialize weights when necessary. See the Convergence section below for more details. The Python snippet below provides the code for constructing the neural network class.

class NN(nn.Module):
    def __init__(self, input_size, hidden_sizes, output_size):
        super().__init__()
        layers = []
        in_size = input_size
        for hidden_size in hidden_sizes:

            layers.append(nn.Linear(in_size, hidden_size).to(torch.float64))
            layers.append(nn.Softplus())
            in_size = hidden_size
        layers.append(nn.Linear(in_size, output_size).to(torch.float64))
        layers.append(oneplustanh())
        self.layers = nn.ModuleList(layers)
        self.network = nn.Sequential(*layers)
        self.initialize_weights()

    def initialize_weights(self):
        for layer in self.layers:
            if isinstance(layer, nn.Linear):
                n_in = layer.in_features
                n_out = layer.out_features
                std = (1 / (n_in + n_out))**0.5
                nn.init.normal_(layer.weight, mean=0, std=std)
                nn.init.normal_(layer.bias, mean=0, std=std)
    
    def forward(self, x):
        # return torch.tensor([[1]]).to(device)
        return self.network(x)

def reinitialize_model(model,weight_init,scheduler_class, scheduler_args, lr=learning_rate):

    if weight_init:
        model.initialize_weights()

    optimizer = optim.Adam(model.parameters(), lr=lr)
    scheduler = scheduler_class(optimizer, **scheduler_args)
    return optimizer, scheduler

Convergence

For learning rate scheduling, a slightly different approach is taken compared to Ackerer et al. The following table summerizes the convergence techniques used for training:

Results

We have two kinds of training: training a reference model and training models using the reference model for initialization. Training the reference model takes significantly more time, whereas training other models based on it takes much less time. The following table shows the parameter values of the reference model as well as an example of those of a subsequent model trained based on the reference model. It is emphasized that training the subsequent model is much faster.

Table: Parameter Values

Numerical results show that the enhanced model, incorporating a neural network with the loss function \(\omega(k,\tau; \theta)\) with SSVI as prior, fits the Bates model data perfectly and produces an arbitrage-free volatility surface. The figure below shows how well the model fits the training data.

The following also is an example of the training trajectory where a feedforward neural network with 4 hidden layers, with 40 units in each layer.

Final Remarks

Several numerical experiments are still to be conducted using the project code. These experiments include applying deep smoothing to real market data rather than synthetic data, performing backtesting, examining W-shaped smile patterns, investigating different sizes for neural network layers, and comparing results with a stand-alone prior (which can be easily implemented by setting the neural network’s forward pass output to 1). We encourage interested readers to explore these avenues further.

Appendix: Arb-Free Volatility Surfaces

The figure below illustrates 29 trained implied variance surfaces obtained using the deep smoothing algorithm. Different values for the Bates model parameters are used in each case. Displayed parameters are \(\alpha, \beta, \kappa, v_0, \theta, \rho, \sigma, \lambda\) respectively. .

Bibtex

@misc{baghal2024implementing,
  title={Implementing Deep Smoothing for Implied Volatility Surfaces (IVS)},
  author={Baghal, Sina},
  journal={https://sinabaghal.github.io/deepsmoothing/},
  url={https://sinabaghal.github.io/deepsmoothing/},
  year={2024}
}

References