This blog post explores the process of calibrating the Heston model using market prices of European call and put options, building on the foundational Black-Scholes-Merton (BSM) framework. We start with a brief overview of the BSM model, explaining how it calculates option prices based on fixed parameters such as volatility, risk-free rate, and the underlying asset price, and providing its core equations. While the BSM model is widely used, it assumes constant volatility, which limits its accuracy in capturing the complexities of real markets.

To address this, we introduce the Heston model, an extension of the BSM framework that allows volatility to be stochastic rather than constant. This feature makes the Heston model more flexible and realistic, particularly for assets that exhibit volatility clustering or mean reversion—common characteristics in financial markets. By incorporating stochastic volatility, the Heston model provides a more dynamic framework for option pricing.

The calibration problem involves estimating the parameters of the Heston model so that it aligns closely with observed market prices. This process is essential for accurate option pricing in a stochastic volatility environment. In this post, we walk through a calibration approach using the Levenberg-Marquardt algorithm, detailing the steps required to match model prices to market data. Furthermore, we discuss how the calibration could be enhanced by incorporating a Kalman filter, which allows for real-time parameter updates as new data becomes available, improving the model’s adaptability to changing market conditions.

All code used in this blog post is open source and available in the stochastic-rs library, which provides a comprehensive Rust implementation of the Heston model along with other tools for simulating and analyzing stochastic processes. This library offers a solid foundation for those looking to explore advanced quantitative finance models and serves as a valuable resource for further experimentation and customization.

Black-Scholes-Merton (BSM) Model

The Black-Scholes-Merton (BSM) model is a foundational framework for pricing European options, assuming constant volatility and risk-free interest rates. Given an option's underlying price $S$ , strike price $K$ , risk-free interest rate $r$ , volatility $\sigma$ , and time to maturity $T$ , the BSM model calculates the theoretical prices of European call and put options.

BSM Call and Put Pricing Formulas

The BSM formulas for call and put options are as follows:

$C_{\text{BSM}} = S \cdot N(d_1) - K e^{-r T} \cdot N(d_2)$

$P_{\text{BSM}} = K e^{-r T} \cdot N(-d_2) - S \cdot N(-d_1)$

where

$d_1 = \frac{\ln(S / K) + (r + \frac{\sigma^2}{2}) T}{\sigma \sqrt{T}}$

$d_2 = d_1 - \sigma \sqrt{T}$

In these equations:

$S$ : Current price of the underlying asset
$K$ : Strike price of the option
$r$ : Risk-free interest rate
$\sigma$ : Volatility of the underlying asset’s returns
$T$ : Time to maturity
$N(\cdot)$ : Cumulative distribution function of the standard normal distribution

Heston Model for Stochastic Volatility

The Heston model expands on the BSM model by allowing the volatility of the underlying asset to be stochastic rather than constant. This model assumes that the variance of the underlying asset price follows a square-root process, leading to more realistic dynamics in options pricing, especially for assets exhibiting volatility clustering or mean reversion.

Heston Model Dynamics

The Heston model is described by the following stochastic differential equations:

$dS_t = \mu S_t \, dt + \sqrt{v_t} S_t \, dW_t^S$

$dv_t = \kappa (\theta - v_t) \, dt + \xi \sqrt{v_t} \, dW_t^v$

where:

$S_t$ : Price of the underlying asset at time $t$
$v_t$ : Stochastic variance of the asset at time $t$
$\mu$ : Drift of the asset's price
$\kappa$ : Rate of mean reversion of the variance
$\theta$ : Long-term average variance
$\xi$ : Volatility of the variance (volatility of volatility)
$W_t^S$ and $W_t^v$ : Two Wiener processes with correlation $\rho$

The main challenge in calibrating the Heston model is to determine the optimal parameters $(\kappa, \theta, \xi, \rho, v_0)$ such that the model-derived option prices closely match observed market prices.

Currently, the correlation is considered constant over time, but recent research has highlighted that correlation can vary over time and may be modeled by an Ornstein-Uhlenbeck (OU) or Fractional Ornstein-Uhlenbeck (FOU) process. However, this approach is beyond the scope of this post, and there may not always be a clear solution for it.

BS Price under the Heston model

In the Heston model, the price of a European call option can be calculated using the formula:

$C(t, S_t, v_t, K, T) = S_t P_1 - K e^{-r(T-t)} P_2$

where $P_1$ and $P_2$ are probabilities defined as follows for $j = 1, 2$ :

$P_j(x_t, v_t; x_T, \ln K) = \frac{1}{2} + \frac{1}{\pi} \int_0^{\infty} \text{Re} \left( \frac{e^{-i\phi \ln K} f_j(\phi; t, x, v)}{i\phi} \right) d\phi$

The function $f_j(\phi; v_t, x_t)$ is given by:

$f_j(\phi; v_t, x_t) = \exp \left( C_j(\tau, \phi) + D_j(\tau, \phi) v_t + i \phi x_t \right)$

where $C_j(\tau, \phi)$ and $D_j(\tau, \phi)$ are defined by the following equations:

$C_j(\tau, \phi) = (r - q) i \phi \tau + \frac{a}{\sigma^2} \left( (b_j - \rho \sigma i \phi + d_j) \tau - 2 \ln \frac{1 - g_j e^{d_j \tau}}{1 - g_j} \right)$

$D_j(\tau, \phi) = \frac{b_j - \rho \sigma i \phi + d_j}{\sigma^2} \left( \frac{1 - e^{d_j \tau}}{1 - g_j e^{d_j \tau}} \right)$

with parameters defined as follows:

$g_j = \frac{b_j - \rho \sigma i \phi + d_j}{b_j - \rho \sigma i \phi - d_j}$

$d_j = \sqrt{(b_j - \rho \sigma i \phi)^2 - \sigma^2 (2 i u_j \phi - \phi^2)}$

and

$u_1 = \frac{1}{2}, \quad u_2 = -\frac{1}{2}, \quad a = \kappa \theta, \quad b_1 = \kappa + \lambda - \rho \sigma, \quad b_2 = \kappa + \lambda$

To compute the option price, you need to:

Calculate $P_1$ and $P_2$ by integrating over the characteristic functions $f_j(\phi; v_t, x_t)$ for $j = \{1,2\}$ .
Substitute $P_1$ and $P_2$ into the price formula.

This approach leverages the Fourier transform to compute the probabilities $P_1$ and $P_2$ , which encapsulate the complex dynamics of the Heston model, including the stochastic volatility component.

Implementation in Rust

In this section, we’ll take a closer look at the Rust code implementing the Heston model. Here’s the full code with explanations on the important parts. The double_exponential::integrate method is used here for numerical integration, which is particularly effective for handling oscillatory integrals, making it well-suited for complex functions like the ones in the Heston model.

Heston Model Code in Rust

#[derive(ImplNew, Clone)]
pub struct HestonPricer {
    /// Stock price
    pub s: f64,
    /// Initial volatility
    pub v0: f64,
    /// Strike price
    pub k: f64,
    /// Risk-free rate
    pub r: f64,
    /// Dividend yield
    pub q: Option<f64>,
    /// Correlation between the stock price and its volatility
    pub rho: f64,
    /// Mean reversion rate
    pub kappa: f64,
    /// Long-run average volatility
    pub theta: f64,
    /// Volatility of volatility
    pub sigma: f64,
    /// Market price of volatility risk
    pub lambda: Option<f64>,
    /// Time to maturity
    pub tau: Option<f64>,
    /// Evaluation date
    pub eval: Option<chrono::NaiveDate>,
    /// Expiration date
    pub expiry: Option<chrono::NaiveDate>,
}

impl Pricer for HestonPricer {
    /// Calculate the price of a European call and put option using the Heston model
    fn calculate_call_put(&self) -> (f64, f64) {
        let tau = self.tau().unwrap_or(1.0);

        let call = self.s * (-self.q.unwrap_or(0.0) * tau).exp() * self.p(1, tau)
            - self.k * (-self.r * tau).exp() * self.p(2, tau);
        let put = call + self.k * (-self.r * tau).exp() - self.s * (-self.q.unwrap_or(0.0) * tau).exp();

        (call, put)
    }
}

impl Time for HestonPricer {
    fn tau(&self) -> Option<f64> {
        self.tau
    }

    fn eval(&self) -> chrono::NaiveDate {
        self.eval.unwrap()
    }

    fn expiration(&self) -> chrono::NaiveDate {
        self.expiry.unwrap()
    }
}

impl HestonPricer {
    pub(self) fn u(&self, j: u8) -> f64 {
        match j {
            1 => 0.5,
            2 => -0.5,
            _ => panic!("Invalid j"),
        }
    }

    pub(self) fn b(&self, j: u8) -> f64 {
        match j {
            1 => self.kappa + self.lambda.unwrap_or(1.0) - self.rho * self.sigma,
            2 => self.kappa + self.lambda.unwrap_or(1.0),
            _ => panic!("Invalid j"),
        }
    }

    pub(self) fn d(&self, j: u8, phi: f64) -> Complex64 {
        ((self.b(j) - self.rho * self.sigma * phi * Complex64::i()).powi(2)
            - self.sigma.powi(2) * (2.0 * Complex64::i() * self.u(j) * phi - phi.powi(2)))
        .sqrt()
    }

    pub(self) fn g(&self, j: u8, phi: f64) -> Complex64 {
        (self.b(j) - self.rho * self.sigma * Complex64::i() * phi + self.d(j, phi))
            / (self.b(j) - self.rho * self.sigma * Complex64::i() * phi - self.d(j, phi))
    }

    pub(self) fn C(&self, j: u8, phi: f64, tau: f64) -> Complex64 {
        (self.r - self.q.unwrap_or(0.0)) * Complex64::i() * phi * tau
            + (self.kappa * self.theta / self.sigma.powi(2))
                * ((self.b(j) - self.rho * self.sigma * Complex64::i() * phi + self.d(j, phi)) * tau
                    - 2.0
                        * ((1.0 - self.g(j, phi) * (self.d(j, phi) * tau).exp()) / (1.0 - self.g(j, phi))).ln())
    }

    pub(self) fn D(&self, j: u8, phi: f64, tau: f64) -> Complex64 {
        ((self.b(j) - self.rho * self.sigma * Complex64::i() * phi + self.d(j, phi))
            / self.sigma.powi(2))
            * ((1.0 - (self.d(j, phi) * tau).exp())
                / (1.0 - self.g(j, phi) * (self.d(j, phi) * tau).exp()))
    }

    pub(self) fn f(&self, j: u8, phi: f64, tau: f64) -> Complex64 {
        (self.C(j, phi, tau) + self.D(j, phi, tau) * self.v0 + Complex64::i() * phi * self.s.ln()).exp()
    }

    pub(self) fn re(&self, j: u8, tau: f64) -> impl Fn(f64) -> f64 {
        let self_ = self.clone();
        move |phi: f64| -> f64 {
            (self_.f(j, phi, tau) * (-Complex64::i() * phi * self_.k.ln()).exp() / (Complex64::i() * phi))
                .re
        }
    }

    pub(self) fn p(&self, j: u8, tau: f64) -> f64 {
        0.5 + FRAC_1_PI * double_exponential::integrate(self.re(j, tau), 0.00001, 50.0, 10e-6).integral
    }
}

Why Use `double_exponential::integrate`

The double_exponential::integrate function is particularly useful for the Heston model because it efficiently handles oscillatory integrals, which are typical in option pricing models with complex-valued characteristic functions. This integration method offers high accuracy even with oscillating functions, making it well-suited for calculating the real part of the characteristic functions $f_j(\phi)$ over a range of values.

By using this method, we can accurately compute the probabilities $P_1$ and $P_2$ necessary for pricing the options, capturing the full dynamics of the Heston model with complex numbers and avoiding numerical instability.

Calibration Objective

The goal of the calibration process is to minimize the difference between the market prices and model prices of options, typically using a least-squares optimization method. This enables us to fit the Heston model parameters to real market data, achieving more accurate option pricing that accounts for stochastic volatility.

Certainly! Below is the expanded calibration section that includes both the partial derivatives required for constructing the Jacobian matrix and the formulas for the initial parameter estimates using the Normal Maximum Likelihood Estimation (NMLE) method.

Calibration Using the Levenberg-Marquardt Algorithm

For calibrating the Heston model, we use the Levenberg-Marquardt algorithm to optimize the model parameters, $\theta = (\kappa, \sigma, \theta, \rho, v_0)$ , to best fit observed market option prices. This process requires an accurate initial parameter guess, which we estimate using the Normal Maximum Likelihood Estimation (NMLE) method, and a Jacobian matrix with closed-form partial derivatives of the model with respect to each parameter.

Step 1: Initial Parameter Guess Using NMLE

When the initial parameters are unknown or need refinement, we can use the NMLE method to estimate them. The initial parameter formulas are derived from Parameter Estimates of Heston Stochastic Volatility Model with MLE and Consistent EKF Algorithm by Ximei Wang, Xingkang He, Ying Bao, and Yanlong Zhao.

The NMLE formulas are as follows:

$\hat{\kappa} = \frac{2}{\delta} \left(1 + \frac{\hat{P} \delta}{2} \frac{1}{n} \sum_{k=1}^n \frac{1}{V_{k-1}} - \frac{1}{n} \sum_{k=1}^n \sqrt{\frac{V_k}{V_{k-1}}}\right)$ ,

$\hat{\sigma} = \sqrt{\frac{4}{\delta n} \sum_{k=1}^n \left( \sqrt{V_k} - \sqrt{V_{k-1}} - \frac{\delta}{2 \sqrt{V_{k-1}}}(\hat{P} - \hat{\kappa} V_{k-1}) \right)^2}$ ,

$\hat{\theta} = \hat{P} + \frac{\hat{\sigma}^2}{2 \hat{\kappa}}$ .

where:

$\delta$ : Time step (e.g., daily, monthly),
$V_k$ : Observed variance at time step $k$ ,
$\hat{P}$ : Long-term mean of variance.

These initial estimates improve the convergence speed and accuracy of the calibration by providing a reasonable starting point for the optimization algorithm.

Step 2: Jacobian Matrix Calculation with Closed-Form Derivatives

Once we have an initial parameter guess, we proceed with calculating the Jacobian matrix, which is essential for the Levenberg-Marquardt optimization. The Jacobian contains the partial derivatives of the residuals with respect to each parameter, which are calculated as follows:

The Jacobian matrix $\nabla C(\theta; K, \tau)$ for the Heston model parameters is given by:

$\nabla C(\theta; K, \tau) = \frac{e^{-r\tau}}{\pi} \left[ \int_0^{\infty} \text{Re} \left( \frac{e^{-iu \log \frac{K}{S_0}}}{iu} \nabla \phi(\theta; u - i, \tau) \right) du \right] - K \int_0^{\infty} \text{Re} \left( \frac{e^{-iu \log \frac{K}{S_0}}}{iu} \nabla \phi(\theta; u, \tau) \right) du$

where $\nabla \phi(\theta; u, \tau) = \phi(\theta; u, \tau) h(u)$ , and $h(u) = [h_1(u), h_2(u), h_3(u), h_4(u), h_5(u)]^T$ represents the derivatives of the characteristic function with respect to each parameter.

The components of $h(u)$ are as follows:

$h_1(u) = -\frac{A}{v_0}$ ,

$h_2(u) = \frac{2\kappa}{\sigma^2} D - \frac{\kappa \rho i u}{\sigma}$ ,

$h_3(u) = -\frac{\partial A}{\partial \rho} + \frac{2 \kappa \theta}{\sigma^2 d} \left( \frac{\partial d}{\partial \rho} - \frac{d A_2}{\partial \rho} \right) - \kappa \theta \rho i u$ ,

$h_4(u) = \frac{1}{\sigma i u} \frac{\partial A}{\partial \rho} + \frac{2 \theta}{\sigma^2} D + \frac{2 \kappa \theta}{\sigma^2 d} \left( \frac{\partial d}{\partial \kappa} - \frac{d A_2}{\partial \kappa} \right)$ ,

$h_5(u) = -\frac{\partial A}{\partial \sigma} - \frac{4 \kappa \theta}{\sigma^3} D + \frac{2 \kappa \theta}{\sigma^2 d} \left( \frac{\partial d}{\partial \sigma} - \frac{d A_2}{\partial \sigma} \right) + \frac{\kappa \theta \rho i u}{\sigma^2}$ .

These partial derivatives allow us to construct the Jacobian matrix efficiently, which is crucial for the optimization algorithm to converge quickly and accurately.

Rust Code for Calibration

Here's the Rust code that performs the calibration, incorporating both the initial parameter estimation with NMLE and the calculation of the Jacobian matrix using the closed-form derivatives:

/// Heston model parameters
#[derive(Clone, Debug)]
pub struct HestonParams {
    pub v0: f64,
    pub theta: f64,
    pub rho: f64,
    pub kappa: f64,
    pub sigma: f64,
}

impl From<HestonParams> for DVector<f64> {
    fn from(params: HestonParams) -> Self {
        DVector::from_vec(vec![
            params.v0,
            params.theta,
            params.rho,
            params.kappa,
            params.sigma,
        ])
    }
}

impl From<DVector<f64>> for HestonParams {
    fn from(params: DVector<f64>) -> Self {
        HestonParams {
            v0: params[0],
            theta: params[1],
            rho: params[2],
            kappa: params[3],
            sigma: params[4],
        }
    }
}

/// A calibrator for the Heston model.
#[derive(ImplNew, Clone)]
pub struct HestonCalibrator {
    pub params: HestonParams,             // Parameters to calibrate
    pub c_market: DVector<f64>,           // Market option prices
    pub s: DVector<f64>,                  // Asset price vector
    pub k: DVector<f64>,                  // Strike price vector
    pub tau: f64,                         // Time to maturity
    pub r: f64,                           // Risk-free rate
    pub q: Option<f64>,                   // Dividend yield
    pub option_type: OptionType,          // Option type (Call or Put)
    derivates: RefCell<Vec<Vec<f64>>>,    // Derivative matrix
}

impl HestonCalibrator {
    pub fn calibrate(&self) {
        println!("Initial guess: {:?}", self.params);

        let (result, ..) = LevenbergMarquardt::new().minimize(self.clone());

        // Print market prices and model prices after calibration
        println!("Market prices: {:?}", self.c_market);
        let residuals = result.residuals().unwrap();
        println!("Model prices: {:?}", self.c_market.clone() + residuals);
        println!("Calibration report: {:?}", result.params);
    }

    /// Set initial parameters based on Normal Maximum Likelihood Estimation (NMLE)
    pub fn set_initial_params(&mut self, s: Array1<f64>, v: Array1<f64>, r: f64) {
        self.params = nmle_heston(s, v, r);
    }
}

impl<'a> LeastSquaresProblem<f64, Dyn, Dyn> for HestonCalibrator {
    type JacobianStorage = Owned<f64, Dyn, Dyn>;
    type ParameterStorage = Owned<f64, Dyn>;
    type ResidualStorage = Owned<f64, Dyn>;

    fn set_params(&mut self, params: &DVector<f64>) {
        self.params = HestonParams::from(params.clone());
    }

    fn params(&self) -> DVector<f64> {
        self.params.clone().into()
    }

    fn residuals(&self) -> Option<DVector<f64>> {
        let mut c_model = DVector::zeros(self.c_market.len());
        let mut derivates = Vec::new();

        for (idx, _) in self.c_market.iter().enumerate() {
            let pricer = HestonPricer::new(
                self.s[idx],
                self.params.v0,
                self.k[idx],
                self.r,
                self.q,
                self.params.rho,
                self.params.kappa,
                self.params.theta,
                self.params.sigma,
                None,
                Some(self.tau),
                None,
                None,
            );
            let (

call, put) = pricer.calculate_call_put();

            match self.option_type {
                OptionType::Call => c_model[idx] = call,
                OptionType::Put => c_model[idx] = put,
            }

            derivates.push(pricer.derivatives());
        }

        let _ = std::mem::replace(&mut *self.derivates.borrow_mut(), derivates);
        Some(c_model - self.c_market.clone())
    }

    fn jacobian(&self) -> Option<DMatrix<f64>> {
        let derivates = self.derivates.borrow();
        let derivates = derivates.iter().flatten().cloned().collect::<Vec<f64>>();

        let jacobian = DMatrix::from_vec(derivates.len() / 5, 5, derivates);
        Some(jacobian)
    }
}

/// Maximum likelihood estimation for Heston model parameters.
///
/// Based on the methodology from:
/// [Parameter estimates of Heston stochastic volatility model with MLE and consistent EKF algorithm](http://scis.scichina.com/en/2018/042202.pdf)
///
/// # Arguments
/// * `s` - Stock prices as an `Array1<f64>` (vector).
/// * `v` - Volatility values as an `Array1<f64>` (vector).
/// * `r` - Risk-free rate as a `f64`.
///
/// # Returns
/// * `HestonParams` - Struct containing estimated parameters for the Heston model.
pub fn nmle_heston(s: Array1<f64>, v: Array1<f64>, r: f64) -> HestonParams {
    let n = v.len();
    let delta = 1.0 / n as f64;
    let mut sum = [0.0; 4];

    // Accumulate required summations for MLE estimation
    for i in 1..n {
        // Sum of (V_i / V_{i-1})
        sum[0] += v[i] / v[i - 1];

        // Sum of V_i
        sum[1] += v[i];

        // Sum of V_{i-1}
        sum[2] += v[i - 1];

        // Sum of 1 / V_{i-1}
        sum[3] += 1.0 / v[i - 1];
    }

    // Calculate intermediate beta estimates
    let beta_hat1 = ((n as f64).powi(-2) * sum[1] * sum[3] - (n as f64).powi(-1) * sum[0])
        / ((n as f64).powi(-2) * sum[2] * sum[3] - 1.0);
    let beta_hat2 = ((n as f64).powi(-1) * sum[0] - beta_hat1) 
        / ((1.0 - beta_hat1) * (n as f64).powi(-1) * sum[3]);

    // Calculate beta_hat3 for further parameter estimation
    let sum_beta_hat3 = {
        let mut sum = 0.0;
        for i in 1..n {
            sum += (v[i] - v[i - 1] * beta_hat1 - beta_hat2 * (1.0 - beta_hat1).powi(2)) 
                * (1.0 / v[i - 1]);
        }
        sum
    };
    let beta_hat3 = (n as f64).powi(-1) * sum_beta_hat3;

    // Compute Heston model parameters
    let kappa_hat = -(1.0 / delta) * beta_hat1.ln();
    let theta_hat = beta_hat2;
    let sigma_hat = (2.0 * kappa_hat * beta_hat3) / (1.0 - beta_hat1.powi(2)).sqrt();

    // Calculate the correlation (rho) estimate
    let mut sum_dw1dw2 = 0.0;
    for i in 1..n {
        let dw1_i = (s[i].ln() - s[i - 1].ln() - (r - 0.5 * v[i - 1]) * delta) / v[i - 1].sqrt();
        let dw2_i = (v[i] - v[i - 1] - kappa_hat * (theta_hat - v[i - 1]) * delta)
            / (sigma_hat * v[i - 1].sqrt());
        
        sum_dw1dw2 += dw1_i * dw2_i;
    }
    let rho_hat = sum_dw1dw2 / (n as f64 * delta);

    // Return estimated parameters as HestonParams struct
    HestonParams {
        v0: v[0],         // Initial volatility
        theta: theta_hat, // Long-run average volatility
        rho: rho_hat,     // Correlation between stock price and variance
        kappa: kappa_hat, // Mean reversion rate
        sigma: sigma_hat, // Volatility of volatility
    }
}

Tests

Here’s a table that compares the model prices generated by the calibrated parameters against the actual market prices. Each row corresponds to a different initial guess for $v_0$ and shows how the calibrated Heston model parameters evolve and how the resulting model prices approximate the market prices.

Initial $v_0$	Theta	Rho	Kappa	Sigma	Market Prices	Model Prices	Difference (Market - Model)
0.0	6.47e-5	-0.00198	0.00657	0.000509	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[33.00, 28.05, 22.12, 15.57, 9.31, 4.01, 0.35, -1.47, -1.65, -0.54, 1.02, 2.09]	[-2.25, -2.17, -1.12, 0.93, 2.57, 3.68, 4.09, 3.57, 2.43, 0.79, -0.92, -1.99]
0.0011	0.1029	-0.0617	0.1097	-0.0188	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.21, 26.47, 21.81, 17.31, 13.24, 9.67, 6.78, 4.58, 3.00, 1.89, 1.25, 0.92]	[-0.46, -0.59, -0.81, -0.81, -1.36, -1.98, -2.34, -2.48, -2.22, -1.64, -1.15, -0.82]
0.0022	0.0909	-0.0570	0.0991	-0.0185	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.21, 26.46, 21.76, 17.20, 13.08, 9.48, 6.57, 4.37, 2.83, 1.77, 1.18, 0.88]	[-0.46, -0.58, -0.76, -0.70, -1.20, -1.79, -2.13, -2.27, -2.05, -1.52, -1.08, -0.78]
0.0033	0.0793	-0.0528	0.0889	-0.0184	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.22, 26.45, 21.70, 17.09, 12.91, 9.26, 6.33, 4.15, 2.63, 1.63, 1.11, 0.85]	[-0.47, -0.57, -0.70, -0.59, -1.03, -1.57, -1.89, -2.05, -1.85, -1.38, -1.01, -0.75]
0.0044	0.0683	-0.0491	0.0792	-0.0186	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.24, 26.44, 21.65, 16.97, 12.72, 9.02, 6.07, 3.89, 2.42, 1.49, 1.03, 0.83]	[-0.49, -0.56, -0.65, -0.47, -0.84, -1.33, -1.63, -1.79, -1.64, -1.24, -0.93, -0.73]
0.0056	0.0577	-0.0459	0.0698	-0.0190	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.26, 26.45, 21.60, 16.84, 12.52, 8.76, 5.78, 3.62, 2.19, 1.34, 0.96, 0.81]	[-0.51, -0.57, -0.60, -0.34, -0.64, -1.07, -1.34, -1.52, -1.41, -1.09, -0.86, -0.71]
0.0067	0.0476	-0.0431	0.0608	-0.0196	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.31, 26.47, 21.56, 16.72, 12.30, 8.47, 5.46, 3.31, 1.95, 1.18, 0.90, 0.81]	[-0.56, -0.59, -0.56, -0.22, -0.42, -0.78, -1.02, -1.21, -1.17, -0.93, -0.80, -0.71]
0.0078	0.0379	-0.0406	0.0521	-0.0204	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.36, 26.50, 21.53, 16.58, 12.07, 8.15, 5.10, 2.98, 1.68, 1.01, 0.83, 0.82]	[-0.61, -0.62, -0.53, -0.08, -0.19, -0.46, -0.66, -0.88, -0.90, -0.76, -0.73, -0.72]
0.0089	0.0285	-0.0386	0.0437	-0.0214	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.44, 26.55, 21.51, 16.45, 11.82, 7.81, 4.72, 2.61, 1.39, 0.84, 0.78, 0.85]	[-0.69, -0.67, -0.51, 0.05, 0.06, -0.12, -0.28, -0.51, -0.61, -0.59, -0.68, -0.75]
0.0100	0.0269	-0.0274	0.0369	-0.0113	[30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10]	[31.48, 26.59, 21.50, 16.39, 11.70, 7.64, 4.53, 2.44, 1.24, 0.76, 0.76, 0.86]	[-0.73, -0.71, -0.50, 0.11, 0.18, 0.05, -0.09, -0.34, -0.46, -0.51, -0.66, -0.76]

Analysis

Model Fit: As we see from the differences between the market and model prices, the calibration starts off with larger discrepancies when $v_0$ is initialized at very low values (e.g., 0.0). However, as $v_0$ increases, the calibrated model parameters tend to improve the fit, reducing the differences with market prices.
Parameter Adjustments: The values of theta, rho, kappa, and sigma vary with each initial guess of $v_0$ , showing how sensitive the calibration is to the initial starting values.
Best Fit: Generally, smaller values for $v_0$ (closer to 0.0011 - 0.0056) yield better fits, as evidenced by the smaller differences in market and model prices.

Conclusion

The Heston model is an industry-standard for options pricing, particularly valued for its ability to capture stochastic volatility, a crucial feature in real-world markets where asset volatility is not constant. As such, the Heston model is widely used by practitioners and researchers alike, making it a foundational tool worth mastering for those involved in financial engineering and quantitative analysis.

In this study, we demonstrated the calibration of the Heston model using the Levenberg-Marquardt algorithm, showing how different initial values for $v_0$ influence the calibration outcome. The process revealed that with appropriately chosen initial parameters, the model can be fine-tuned to closely approximate observed market prices. The calibration outcomes highlight the importance of selecting suitable initial values and illustrate the sensitivity of the model to parameter adjustments.

However, as markets evolve, the static nature of traditional calibration methods like Levenberg-Marquardt may limit their effectiveness in capturing time-varying dynamics. To address this, the research paper by Wang, He, Bao, and Zhao suggests an extension using the Kalman filter, which allows for dynamic updating of parameters as new data becomes available. By incorporating a Kalman filter, the calibration can adapt to shifts in market conditions in real-time, leading to more responsive and accurate pricing over time.

Integrating the Kalman filter into the calibration of the Heston model would not only enhance accuracy but also offer a "rolling" calibration approach, especially beneficial in high-frequency trading environments. This would enable the model to adjust continuously, better reflecting current market volatility and improving risk management and hedging strategies.

In summary, the Heston model’s flexibility and adaptability in capturing stochastic volatility make it an indispensable tool in options pricing. Enhancing its calibration with dynamic methods like the Kalman filter holds promise for the future, providing a robust framework for managing financial risk in rapidly changing markets. For those in quantitative finance, understanding and working with the Heston model is essential, as it remains one of the most reliable and widely adopted models in the industry.

References

Yiran Cui, Sebastian del Baño Rollin, and Guido Germano. Full and Fast Calibration of the Heston Stochastic Volatility Model. European Journal of Operational Research, 2017. DOI: 10.1016/j.ejor.2017.05.018
Ximei Wang, Xingkang He, Ying Bao, and Yanlong Zhao. Parameter Estimates of Heston Stochastic Volatility Model with MLE and Consistent EKF Algorithm. Science China Information Sciences, 2018. DOI: 10.1007/s11432-017-9215-8

These references provide both the analytical forms of the partial derivatives for constructing the Jacobian and the methods for estimating initial parameters via Normal Maximum Likelihood Estimation (NMLE), enhancing calibration accuracy and efficiency in the Heston model.