Computational Reconstruction of the Volatility Term Structure in the General Hull–White Model

Abstract

Volatility recovery is of paramount importance in contemporary finance. Volatility levels are heavily used in risk and portfolio management. We employ the Hull–White one- and two-factor models to describe the market condition. We computationally recover the volatility term structure as a piecewise-linear function of time. For every maturity, a cost functional, defined as the squared differences between theoretical and market prices, is minimized and the respective linear part is reconstructed. On the last time steps, before each maturity, the derivative price is decomposed in order to make the minimization problem analytically solvable. The procedure works fast since only scalar values are obtained on each minimization. However, the predictor–corrector nature of the algorithm allows for the precise recovery of very complex volatility functions. An implicit scheme is used to solve the PDEs on bounded domains. The computational simulations with artificial and real data show that the proposed algorithm is stable, accurate and efficient.

1. Introduction

Since its introduction by Hull and White back in 1990, the Hull–White model [1] has become a cornerstone in the modeling of interest rates and credit risk within the financial industry. Its widespread adoption among practitioners can be attributed to its simplicity and flexibility, particularly its ability to accommodate negative interest rates—a feature once regarded as a limitation but now recognized as a significant advantage in modern markets. The model’s framework allows for the calibration of observed term structures of interest rates, as well as of market prices or implied volatilities of various derivatives, including caps, floors and European swaptions. This adaptability stems from its mean-reversion and volatility parameters, which are explicitly time-dependent. Furthermore, the Hull–White model offers analytical solutions for the pricing of zero-coupon bonds and European vanilla options. For more sophisticated and exotic derivative instruments, numerical methods need to be employed in order to obtain pricing solutions.

In addition to pricing, the accurate and robust calibration of financial models to real-world market data remains one of the most critical challenges in modern finance. As models evolve beyond the foundational Black–Scholes framework [2], the complexity of the calibration process increases significantly. Empirical and theoretical investigations consistently highlight the inadequacies of constant-volatility models in capturing the nuances of market price behavior [3]. Such models fail to account for critical phenomena, including heavy-tailed distributions, volatility clustering and the presence of volatility smiles and skews.

In the Black–Scholes framework, volatility, denoted as $\sigma$, is a key parameter. However, assuming $\sigma$ to be constant frequently leads to inaccurate predictions of option prices, as the actual volatility of the underlying asset exhibits significant temporal variation. Extensive research has further demonstrated that the volatility term structure is rarely flat, particularly for instruments with short maturities. This suggests a pronounced dependence of volatility on the time remaining to maturity, as supported by both empirical observations and theoretical analyses [4,5].

To address these limitations, we focus on modeling volatility as a time-dependent parameter alongside other dynamic features [6]. This approach provides a more comprehensive and effective framework for accurately pricing and hedging options, delivering insights that align closely with real-world market behavior.

The calibration of financial models is as critical as the construction of the models themselves. Calibration involves determining parameter values that enable the model to replicate observed market prices with the highest possible accuracy. Both precision and computational efficiency in the calibration process are essential, as practitioners rely on the calibrated parameters to value complex derivative instruments and to design high-frequency trading strategies.

Volatility, which quantifies the degree of randomness in price movements, plays a pivotal role in this context. Along with other parameters such as drift and mean-reversion speed, volatility is widely utilized in risk management and portfolio optimization. These parameters, however, are not directly observable in financial markets. Instead, they are typically inferred from the prices of traded options. As a result, accurately reconstructing these hidden parameters from market data is of significant interest to researchers and practitioners [7,8]. This reconstruction is fundamental to ensure that models are not only theoretically sound but also practically applicable in real-world financial scenarios.

The determination of implied volatility and other model parameters can be interpreted as an inverse problem for partial differential equations (PDEs). Such problems are typically ill-posed, requiring careful treatment, often through regularization techniques. Various approaches to address these challenges have been developed within this field [9].

The classical problem of calibrating the Black–Scholes model or related frameworks to vanilla option prices has been extensively studied in the literature [10]. Similar efforts have also been devoted to calibrating interest rate models. One study [11] introduces a calibration approach by treating the short rate as an unobservable variable, simultaneously estimating it along with the model parameters. This method addresses the challenges of market distortions, such as speculative influences on overnight rates, ensuring a robust calibration of the Vašíček model. The authors in [12] present a detailed calibration framework for Ueno’s model, which adapts the Vašíček process to incorporate negative shadow rates, reflecting the modern reality of negative interest rates. By employing maximum likelihood estimation for parameter determination and numerically solving a parabolic PDE to align the market price of risk with long-term rates, the paper suggests a method for calibrating shadow rate models to real-world data. Reference [13] introduces a computationally efficient approach for calibrating the Heston model with a piecewise constant in time parameters, enabling time-dependent market parametrization. By deriving semi-analytical formulae and using Gauss–Kronrod quadrature with control variates, the study calibrates the model to FX options data, demonstrating improved performance over the standard Heston model in specific scenarios.

The authors in [14] extend the Hull–White framework by introducing an SDE with state-dependent coefficients to incorporate market-implied skew and smile features. Employing the Randomized Affine Diffusion (RAnD) technique and preserving the analytic tractability of the Hull–White dynamics, the study provides a semi-analytic calibration method for European swaptions and demonstrates the impact of skew and smile on Valuation Adjustment (xVA) metrics through regression-based Monte Carlo simulations. A model order reduction (MOR) approach for high-dimensional problems of financial risks is presented in [15]. It is based on a proper orthogonal decomposition (POD) method. The authors tested the designed algorithms on the Hull–White model, which is utilized to value a future stochastic rate today. A calibration approach for one-factor and two-factor Hull–White models, using swaptions under a market-consistent framework, is presented in [16]. Reference [17] suggests a framework for calibrating the Hull–White model by employing finite difference methods to solve its associated PDEs. By validating the numerical results, obtained by operator splitting, against analytical solutions, the study ensures accuracy and demonstrates the efficacy of the approach for pricing one- and two-asset bond options under the Hull–White model.

In contrast, the body of research addressing the calibration of time-dependent parameters remains relatively sparse compared to studies focusing on constant parameters [8,18]. Building on this gap, our previous work [7] proposed a robust algorithm for reconstructing time-dependent volatility as a piecewise-linear function of time, based on option price data. This methodology offers a systematic and effective approach to capturing the dynamic nature of volatility due to a number of reasons. First of all, the option quotes, although vast in volume, are not continuous, but, rather, they are specified in the form $V=V(K,T)$, i.e., prices are associated with pairs of strike levels and maturities. Such shaped data are often called market observations. What is more, the algorithm avoids inverting an ill-posed formula or solving a complicated inverse equation. Also, on each step only a scalar is recovered rather that a function, contributing to the robustness of the approach. In that way, the main novelty of the study is achieved, resulting in a time-dependent calibration of the dynamics parameters of the general Hull–White model.

The paper is continued as follows. In the next section, the general 1D and 2D Hull–White models are defined and commented on. The numerical solution to the direct problem is given as well. In Section 3, the solutions to the inverse problems are thoroughly explained, and the algorithm is presented in detail. The decomposition, which greatly improves the algorithm’s performance, is also discussed. Section 4 is dedicated to the numerical experiment for both direct and inverse problems, and the cases for a zero-coupon bond option are considered and finalized with a discussion. The last section concludes the study.

2. Direct Problem

In this section, we will recall the general Hull–White models in the PDE framework [19]. Further, we will introduce the computational solution to the direct problem, which consists of (numerically) obtaining the ‘fair’ price of the financial derivative instrument under consideration, given the market conditions.

2.1. One-Factor General Hull–White Model

The main stochastic model, assumed to model the interest rate dynamics $r\left(t\right)$ in the Hull–White model, is

$\mathrm{d}{r}_t=(\eta \left(t\right)-\gamma \left(t\right)r_t)\mathrm{d}t+\sigma \left(t\right)\mathrm{d}{W}_t,$

where $\eta \left(t\right)$ is the deterministic drift, $\gamma \left(t\right)>0$ is the reversion speed, $\sigma \left(t\right)$ is the short rate volatility and $W_t$ is a standard Wiener process [6].

This short rate model is a Vašíček extended model.

2.2. Two-Factor General Hull–White Model

The two-factor model introduces another source of uncertainty by incorporating a stochastic reversion level:

$\begin{array}{c}\mathrm{d}{r}_t=(\theta \left(t\right)+u_t-a\left(t\right)r_t)\mathrm{d}t+{\sigma }_1\left(t\right)\mathrm{d}{W}_t^1,\hfill \\ \mathrm{d}{u}_t=-b\left(t\right)u_t\mathrm{d}t+{\sigma }_2\left(t\right)\mathrm{d}{W}_t^2,\hfill \end{array}$

where $\theta \left(t\right)+u_t$ is the reversion level, $a\left(t\right)$ is the reversion speed, $b\left(t\right)$ is the reversion speed of $u_t$ and ${\sigma }_i\left(t\right),i=1,2$ are the volatilities [6].

The Wiener processes $W_t^i$ are correlated with instantaneous correlation $\rho \left(t\right)$, i.e.,

$\mathbb{E}\left[\mathrm{d}{W}_t^1\mathrm{d}{W}_t^2\right]=\rho \left(t\right)\mathrm{d}t.$

2.3. Definitions

We consider the European interest rate derivative premium $V=V(r,\tau ;T)$ in a Hull–White setting, which satisfies

(1)${\displaystyle \frac{\partial V}{\partial \tau }}={\displaystyle \frac{1}{2}}{\sigma }^2\left(\tau \right){\displaystyle \frac{{\partial }^{2}V}{\partial r^2}}+(\eta \left(\tau \right)-\gamma \left(\tau \right)r){\displaystyle \frac{\partial V}{\partial r}}-r\left(\tau \right)V,$

$(r,\tau )\in \mathbb{R}\times [0,T)$

with an initial condition of

(2)$V(r,0)=1\mathrm{for} \mathrm{a} \mathrm{zero}-\mathrm{coupon} \mathrm{bond}$

and

(3)$V(r,0)={\left(\mathrm{ZCB}(T,S)-K\right)}^{+} \mathrm{for} \mathrm{a} \mathrm{European} \mathrm{call}$

on a zero-coupon bond, maturing at time S, option maturity $T<S$ and strike K.

When applying a numerical method to solve equations like (1), (2) or (3), the infinite computational domain must be appropriately truncated to a bounded one. Part of it consists of imposing suitable boundary conditions on the left and right boundary of the truncated domain [20].

Neumann boundary conditions are usually imposed of the form

(4)${\displaystyle \frac{\partial V}{\partial r}(r_{min},\tau )=\frac{\partial V}{\partial r}(r_{max},\tau )=0,}$

where $r_{min}\to -\infty$ and $r_{max}\to +\infty$ in a theoretical setting and are a small enough and a large enough number in computational simulations.

We use the linearity boundary conditions

(5)${\displaystyle \frac{{\partial }^{2}V}{\partial r^2}(r_{min},\tau )=\frac{{\partial }^{2}V}{\partial r^2}(r_{max},\tau )=0.}$

We assume that we know the interest rate derivative price V. Then, we find that volatility, $\sigma \left(\tau \right)$, for which the theoretical result coincides with the observed quoted price on the market. This volatility is called implied volatility, i.e., IP1:

$V^{\mathrm{obs}}=V^{\mathrm{HW}}\left(r_0,t;K,T,\eta \left(\tau \right),\gamma \left(\tau \right),{\sigma }^{\mathrm{imp}}\left(\tau \right)\right).$

Now, we assume again that we know the interest rate derivative price V. Then, we find the parameters $\sigma \left(\tau \right)$, $\eta \left(\tau \right)$ and $\gamma \left(\tau \right)$ for which the theoretical result coincides with the observed quoted price on the market. These are the implied parameter values, i.e., IP2:

$V^{\mathrm{obs}}=V^{\mathrm{HW}}\left(r_0,t;K,T,{\eta }^{\mathrm{imp}}\left(\tau \right),{\gamma }^{\mathrm{imp}}\left(\tau \right),{\sigma }^{\mathrm{imp}}\left(\tau \right)\right).$

We also consider the two-dimensional general Hull–White equation

(6)$\begin{array}{c}{\displaystyle \frac{\partial V}{\partial \tau }}={\displaystyle \frac{1}{2}}{\sigma }_1^2\left(\tau \right){\displaystyle \frac{{\partial }^{2}V}{\partial r^2}}+\rho \left(\tau \right){\sigma }_1\left(\tau \right){\sigma }_2\left(\tau \right){\displaystyle \frac{{\partial }^{2}V}{\partial r\partial u}}+{\displaystyle \frac{1}{2}}{\sigma }_2^2\left(\tau \right){\displaystyle \frac{{\partial }^{2}V}{\partial u^2}}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left(\theta \left(\tau \right)+u-a\left(t\right)r\right){\displaystyle \frac{\partial V}{\partial r}}-b\left(t\right)u{\displaystyle \frac{\partial V}{\partial r}}-rV.\end{array}$

with linear boundary conditions

(7)${\displaystyle \frac{{\partial }^{2}V}{\partial r^2}(r_{min},u,\tau )=\frac{{\partial }^{2}V}{\partial r^2}(r_{max},u,\tau )=\frac{{\partial }^{2}V}{\partial u^2}(r,u_{min},\tau )=\frac{{\partial }^{2}V}{\partial u^2}(r,u_{max},\tau )=0}$

and the initial condition

(8)$V(r,u,0)=1.$

It is possible to define similar inverse problems using (IP1) and (IP2) for models (6)–(8).

2.4. Finite Difference Scheme

Let $f:[r_{min},r_{max}]\to \mathbb{R}$ and, if $r_{min}=r_0<r_1<\dots <r_{I+1}=r_{max}$ is the spatial grid, $h_i=r_i-r_{i-1}$, $H_i=h_i+h_{i+1}$ and $(1\le i\le I)$, then the first derivative $f^{\prime }\left(r_i\right)$ could be approximated in the following ways:

(9)$f^{\prime }\left(r_i\right)\approx {\displaystyle \frac{f\left(r_i\right)-f\left(r_{i-1}\right)}{h_i}}\equiv {\mathrm{D}}^l{f}_i,$

(10)$f^{\prime }\left(r_i\right)\approx {\displaystyle \frac{f\left(r_{i+1}\right)-f\left(r_i\right)}{h_{i+1}}}\equiv {\mathrm{D}}^r{f}_i$

The second derivative $f^{″}\left(r_i\right)$ could be approximated as follows:

(11)$f^{″}\left(r_i\right)\approx {\displaystyle \frac{2}{h_i{H}_i}}f\left(r_{i-1}\right)-{\displaystyle \frac{2}{h_i{h}_{i+1}}}f\left(r_i\right)+{\displaystyle \frac{2}{h_{i+1}{H}_i}}f\left(r_{i+1}\right)\equiv {\mathrm{D}}^2{f}_i.$

Let $▵{\tau }^j={\tau }^{j+1}-{\tau }^j,(0\le j\le J)$ and $v_i^j\approx V(r_i,{\tau }^j)$. Applying (9)–(11) to (1), we have the upwind implicit scheme

${\displaystyle \frac{v_0^{j+1}-v_0^j}{▵{\tau }^j}}=({\eta }^{j+1}-{\gamma }^{j+1}{r}_0){\mathrm{D}}^l{v}_0^{j+1}-r_0{v}_0^{j+1},$

(12)$\begin{array}{c}{\displaystyle \frac{v_i^{j+1}-v_i^j}{▵{\tau }^j}}={\displaystyle \frac{1}{2}}{\left({\sigma }^{j+1}\right)}^2{\mathrm{D}}^2{v}_i^{j+1}+\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\ast \left\{\begin{array}{c}{\mathrm{D}}^l{v}_i^{j+1}\mathrm{for}\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)\le 0,\hfill \\ {\mathrm{D}}^r{v}_i^{j+1}\mathrm{for}\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)>0\hfill \end{array}\right\}-r_i{v}_i^{j+1},(1\le i\le I),\end{array}$

${\displaystyle \frac{v_{I+1}^{j+1}-v_I^j}{▵{\tau }^j}}=({\eta }^{j+1}-{\gamma }^{j+1}{r}_{I+1}){\mathrm{D}}^r{v}_{I+1}^{j+1}-r_{I+1}{v}_{I+1}^{j+1},(0\le j<J).$

If we define $v^j={\left[v_0^j,v_1^j,\dots ,v_{I+1}^j\right]}^{\top }$, we could write the three-diagonal linear system in a matrix form as

(13)$\mathbf{M}{v}^j=v^{j+1},$

where

${\mathbf{M}}^l=\left[\begin{array}{ccccc}{C}_0& B_0& & & \\ A_1^l& C_1^l& B_1^l& & \\ & \ddots & \ddots & \ddots & \\ & & A_{I-1}^l& C_{I-1}^l& B_{I-1}^l\\ & & & A_I& C_I\end{array}\right],$

${\mathbf{M}}^r=\left[\begin{array}{ccccc}{C}_0& B_0& & & \\ A_1^r& C_1^r& B_1^r& & \\ & \ddots & \ddots & \ddots & \\ & & A_{I-1}^r& C_{I-1}^r& B_{I-1}^r\\ & & & A_I& C_I\end{array}\right],$

where

$\mathbf{M}={\mathbf{M}}^l.\ast \left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\le 0\right)+{\mathbf{M}}^r.\ast \left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i>0\right),$

$A_i^l=-{\displaystyle \frac{{\left({\sigma }^{j+1}\right)}^2▵{\tau }^j}{h_i{H}_i}}+{\displaystyle \frac{\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)▵{\tau }^j}{h_i}}$

$B_i^l=-{\displaystyle \frac{{\left({\sigma }^{j+1}\right)}^2▵{\tau }^j}{h_{i+1}{H}_i}},$

$C_i^l={\displaystyle \frac{{\left({\sigma }^{j+1}\right)}^2▵{\tau }^j}{h_i{h}_{i+1}}}-{\displaystyle \frac{\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)▵{\tau }^j}{h_i}}+r_i▵{\tau }^j+1,$

and

$A_i^r=-{\displaystyle \frac{{\left({\sigma }^{j+1}\right)}^2▵{\tau }^j}{h_i{H}_i}},$

$B_i^r=-{\displaystyle \frac{{\left({\sigma }^{j+1}\right)}^2▵{\tau }^j}{h_i{H}_i}}-{\displaystyle \frac{\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)▵{\tau }^j}{h_{i+1}}},$

$C_i^r={\displaystyle \frac{{\left({\sigma }^{j+1}\right)}^2▵{\tau }^j}{h_i{h}_{i+1}}}+{\displaystyle \frac{\left({\eta }^{j+1}-{\gamma }^{j+1}{r}_i\right)▵{\tau }^j}{h_{i+1}}}+r_i▵{\tau }^j+1.$

Solving (13) allows us to make a transition to the new temporal layer, and thus to obtain the premium at the ‘initial’ time $\tau =T$.

3. Inverse Problem

In this section, we will describe in great detail the algorithm to solve inverse problems (IP1) and (IP2).

3.1. Observation Definitions

We have a set of market measurements $\left\{{\omega }_{\beta }^{\alpha }\right\}$, where ${\omega }_{\beta }^{\alpha }$ is the quoted price of a derivative with maturity $T_{\alpha }$, $\alpha =1,\dots ,M_T$ and settings $S_{\beta }$, $\beta =1,\dots ,N$, assuming that $T_1\le \dots \le T_{M_T}$.

We minimize the following cost function:

${\displaystyle {\Gamma }_{\alpha }\left(\sigma \right)=\frac{1}{N}\sum _{\beta =1}^N{\left[v_{\beta }({\sigma }_{\alpha }\left({\tau }_{\alpha }\right);S_{\beta },T_{\alpha })-{\omega }_{\beta }^{\alpha }\right]}^2{\chi }_{\beta }^{\alpha },{\tau }_{\alpha }\in (0,T_{\alpha }],\alpha =1,\dots ,M_T,}$

where $v_{\beta }({\sigma }_{\alpha }\left({\tau }_{\alpha }\right);S_{\beta },T_{\alpha })$ is the numerical solution of (1) in setting $S_{\beta }$ and expiry time $T_{\alpha }$.

Now, the algorithms to solve inverse problems (1)–(5), (IP1), follow.

3.2. Algorithm

Step 1.1 We find ${\mu }_1$, which minimizes the cost function

${\displaystyle {\Gamma }_1\left(\sigma \right)=\frac{1}{N}\sum _{\beta =1}^N{\left[v_{\beta }({\sigma }_1\left({\tau }_1\right);S_{\beta },T_1)-{\omega }_{\beta }^1\right]}^2{\chi }_{\beta }^1,{\tau }_1\in (0,T_1].}$

Step 1.2 We assume that the volatility function on $(0,T_1]$ is constant, defined as $\left\{{\sigma }_1\left(\tau \right)\right\}={\mu }_1$. Then, we have

$\begin{array}{ccc}\hfill \sigma \left(\tau \right)& =& {\sigma }_1\left(\tau \right)\hfill \end{array}\mathrm{for}\tau \in [0,T_1].$

Step 2.1 We find ${\mu }_2$, which minimizes the cost function

${\displaystyle {\Gamma }_2\left(\sigma \right)=\frac{1}{N}\sum _{\beta =1}^N{\left[v_{\beta }({\sigma }_2\left({\tau }_2\right);S_{\beta },T_2)-{\omega }_{\beta }^2\right]}^2{\chi }_{\beta }^2,{\tau }_2\in (0,T_2].}$

Step 2.2 We assume that the volatility function for $(0,T_1]$ is linear, defined as ${\sigma }_2\left(\tau \right)=a\tau +b$ (Figure 1). Then, we have

$\begin{array}{ccc}\hfill \sigma \left(\tau \right)& =& {\sigma }_2\left(\tau \right)\hfill \end{array}\mathrm{for}\tau \in [0,T_2].$

Step 3 is repeated from $\alpha =3$ to $\alpha =M_T$.

Step 3.1 We find ${\mu }_{\alpha }:=\left\{{\sigma }_{\alpha }\right\}$, which minimizes the cost function

${\displaystyle {\Gamma }_{\alpha }\left(\sigma \right)=\frac{1}{N}\sum _{\beta =1}^N{\left[v_{\beta }({\sigma }_{\alpha }\left({\tau }_{\alpha }\right);S_{\beta },T_{\alpha })-{\omega }_{\beta }^{\alpha }\right]}^2{\chi }_{\beta }^{\alpha },{\tau }_{\alpha }\in (0,T_{\alpha }].}$

Step 3.2 We define the linear volatility function ${\sigma }_{\alpha }\left(\tau \right)$ on $[T_{\alpha -{\displaystyle \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}},T_{\alpha }]$ as (Figure 2)

${\sigma }_{\alpha }\left(\tau \right)={\displaystyle \frac{{\mu }_{\alpha }-{\mu }_{\alpha -1}}{T_{\alpha }-T_{\alpha -\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}}(\tau -T_{\alpha })+{\mu }_{\alpha }.$

If $\alpha =2$, then the volatility function is linear:

$\begin{array}{ccc}\hfill \sigma \left(\tau \right)& =& {\sigma }_2\left(\tau \right)\hfill \end{array}\mathrm{for}\tau \in [0,T_2].$

If $\alpha \ge 3$, the volatility function is piecewise linear:

$\sigma \left(\tau \right)=\left\{\begin{array}{cc}\hfill {\sigma }_2\left(\tau \right)& \mathrm{for}\tau \in [0,T_{{\displaystyle \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}],\hfill \\ \hfill {\sigma }_j\left(\tau \right)& \mathrm{for}\tau \in [T_{j-{\displaystyle \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}},T_{j-{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}]\mathrm{for}2<j<\alpha ,\hfill \\ \hfill {\sigma }_{\alpha }\left(\tau \right)& \mathrm{for}\tau \in [T_{\alpha -{\displaystyle \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}},T_{\alpha }],\hfill \end{array}\right.$

Finally, we arrive at the recovered volatility function $\sigma \left(\tau \right)$ for $\tau \in (0,T_{M_T}]$.

The algorithm could be slightly modified to solve the inverse problem (IP2).

3.3. Decomposition

The discrete inverse problem (12), given $\left\{{\omega }_{\beta }^{\alpha }\right\}$, which requires the determination of the pair $\left({\Sigma }^{j+1},{\mathrm{D}}^2{v}^{j+1}\right)$, where ${\Sigma }^{j+1}:={\left({\sigma }^{j+1}\right)}^2$, is (quadratically) nonlinear. We use the following approximation of the expression $a^{j+1}{b}^{j+1}$ at the point ${\tau }^{j+{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$:

$a^{j+1}{b}^{j+1}=a^{j+1}{b}^j+a^j{b}^{j+1}-a^j{b}^j+\mathcal{O}\left({(▵\tau )}^2\right).$

In case the nonlinearity is cubic, then the expression $a^{j+1}{b}^{j+1}{c}^{j+1}$ is linearized as follows:

$a^{j+1}{b}^{j+1}{c}^{j+1}=a^{j+1}{b}^j{c}^j+a^j{b}^{j+1}{c}^j+a^j{b}^j{c}^{j+1}-2a^j{b}^j{c}^j+\mathcal{O}\left({(▵\tau )}^2\right).$

Following [18], to find an approximate solution to (IP1) $\left\{{\Sigma }^{j+1},v^{j+1}\right\}$ at the new time layer, we place the decomposition

$v^{j+1}=z^{j+1}+{\Sigma }^{j+1}{w}^{j+1}$

in the linearized counterpart of (12). As a result, we obtain two linear systems for the auxiliary functions $z^{j+1}$ and $w^{j+1}$.

We integrate the decomposition at the last time layer with the cost function minimization approach as follows.

We apply the decomposition only on the last time step before each maturity $T_{\alpha }$. If the respective number of time layers is $j_{\alpha }$ (since we have reversed the time, it appears that ${\tau }^{j_{\alpha }}=t^1$ is the first time layer after the initial condition in the natural time t), then

(14)${\displaystyle {\Gamma }_{\alpha }(\Sigma )=\frac{1}{N}\sum _{\beta =1}^N{\left[\left(z_{\beta ;i_0}^{j_{\alpha }}+{\Sigma }^{j_{\alpha }}{w}_{\beta ;i_0}^{j_{\alpha }}\right)-{\omega }_{\beta }^{\alpha }\right]}^2{\chi }_{\beta }^{\alpha }.}$

Since we want the minimum of ${\Gamma }_{\alpha }$, the candidates are described using

${\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial \Sigma }}=0$

and after rearrangement of the terms in (14) we arrive at

${\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial \Sigma }}={\displaystyle \frac{2}{N}}\left[\underset{A}{\underbrace{\left({\displaystyle \sum _{\beta =1}^N(z_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }){\chi }_{\beta }^{\alpha }{w}_{\beta ,i_0}^{j_{\alpha }}}\right)}}+\underset{B}{\underbrace{\left({\displaystyle \sum _{\beta =1}^N{\left(w_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }}\right)}}{\Sigma }^{j_{\alpha }}\right].$

The root of the latter is ${\sigma }^{j_{\alpha }}=\pm \sqrt{-{\displaystyle \raisebox{1ex}{$A$}\!\left/ \!\raisebox{-1ex}{$B$}\right.}}$, but only the positive ones are financially meaningful, so eventually

${\displaystyle {\sigma }^{\mathrm{imp}}:={\sigma }^{j_{\alpha }}=\sqrt{-\frac{{\displaystyle \sum _{\beta =1}^N(z_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }){\chi }_{\beta }^{\alpha }{w}_{\beta ,i_0}^{j_{\alpha }}}}{{\displaystyle \sum _{\beta =1}^N{\left(w_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }}}}.}$

In the case of three parameters’ reconstruction, to decompose the solution to (IP2) $\left\{{\Sigma }^{j+1},{\eta }^{j+1},{\gamma }^{j+1},v^{j+1}\right\}$ at the new time layer, we use

(15)$v^{j+1}=z^{j+1}+{\Sigma }^{j+1}{w}^{j+1}+{\eta }^{j+1}{x}^{j+1}+{\gamma }^{j+1}{y}^{j+1}$

in the linearized counterpart of (12). Consequently, we obtain four linear systems for the auxiliary functions $z^{j+1}$, $w^{j+1}$, $x^{j+1}$ and $y^{j+1}$.

This decomposition (15) is again applied only on the last time step before each maturity $T_{\alpha }$, layer number $j_{\alpha }$. Then,

${\displaystyle {\Gamma }_{\alpha }(\Sigma ,\eta ,\gamma )=\frac{1}{N}\sum _{\beta =1}^N{\left[\left(z_{\beta ;i_0}^{j_{\alpha }}+{\Sigma }^{j_{\alpha }}{w}_{\beta ;i_0}^{j_{\alpha }}+{\eta }^{j_{\alpha }}{x}_{\beta ;i_0}^{j_{\alpha }}+{\gamma }^{j_{\alpha }}{y}_{\beta ;i_0}^{j_{\alpha }}\right)-{\omega }_{\beta }^{\alpha }\right]}^2{\chi }_{\beta }^{\alpha }.}$

Since we want the minimum of ${\Gamma }_{\alpha }$, the candidates are described using

$\left|\begin{array}{c}{\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial \Sigma }=0}\hfill \\ {\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial \eta }=0}\hfill \\ {\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial \gamma }=0}\hfill \end{array}\right..$

The system is written in the form

$\left|\begin{array}{c}{\displaystyle \underset{A}{\underbrace{\left(\sum _{\beta =1}^N{\left(w_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }^{j_{\alpha }}+\underset{B}{\underbrace{\left(\sum _{\beta =1}^N{w}_{\beta ,i_0}^{j_{\alpha }}{x}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\eta }^{j_{\alpha }}+\underset{C}{\underbrace{\left(\sum _{\beta =1}^N{w}_{\beta ,i_0}^{j_{\alpha }}{y}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\gamma }^{j_{\alpha }}=\underset{G}{\underbrace{-\sum _{\beta =1}^N\left(z_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }\right)w_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }}}}\hfill \\ {\displaystyle \underset{B}{\underbrace{\left(\sum _{\beta =1}^N{w}_{\beta ,i_0}^{j_{\alpha }}{x}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }^{j_{\alpha }}+\underset{D}{\underbrace{\left(\sum _{\beta =1}^N{\left(x_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }\right)}}{\eta }^{j_{\alpha }}+\underset{E}{\underbrace{\left(\sum _{\beta =1}^N{x}_{\beta ,i_0}^{j_{\alpha }}{y}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\gamma }^{j_{\alpha }}=\underset{H}{\underbrace{-\sum _{\beta =1}^N\left(z_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }\right)x_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }}}}\hfill \\ {\displaystyle \underset{C}{\underbrace{\left(\sum _{\beta =1}^N{w}_{\beta ,i_0}^{j_{\alpha }}{y}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }^{j_{\alpha }}+\underset{E}{\underbrace{\left(\sum _{\beta =1}^N{x}_{\beta ,i_0}^{j_{\alpha }}{y}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\eta }^{j_{\alpha }}+\underset{F}{\underbrace{\left(\sum _{\beta =1}^N{\left(y_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }\right)}}{\gamma }^{j_{\alpha }}=\underset{H}{\underbrace{-\sum _{\beta =1}^N\left(z_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }\right)x_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }}}}\hfill \end{array}\right.$

If $\mathrm{DENOM}:=F{B}^2-2BCE+D{C}^2+A{E}^2-ADF\ne 0$, then the solution to the system is given by

$\left|\begin{array}{c}{\displaystyle {\Sigma }^{j_{\alpha }}=\frac{E^{2}G+BFH-CEH-BEI+CDI-DFG}{\mathrm{DENOM}}}\hfill \\ {\displaystyle {\eta }^{j_{\alpha }}=\frac{C^{2}H-BCI-AFH+BFG-CEG+AEI}{\mathrm{DENOM}}}\hfill \\ {\displaystyle {\gamma }^{j_{\alpha }}=\frac{B^{2}I-BCH+AEH-BEG+CDG-ADI}{\mathrm{DENOM}}}\hfill \end{array}\right..$

When solving the 2D problem, we decompose the solution $\left\{{\Sigma }_1^{j+1},{\Sigma }_2^{j+1},{\mathbf{v}}^{j+1}\right\}$ at the new time layer and employ the decomposition

(16)${\mathbf{v}}^{j+1}={\mathbf{z}}^{j+1}+{\Sigma }_1^{j+1}{\mathbf{w}}^{j+1}+{\Sigma }_2^{j+1}{\mathbf{x}}^{j+1}.$

Obviously, we obtain three linear systems for the auxiliary functions ${\mathbf{z}}^{j+1}$, ${\mathbf{w}}^{j+1}$ and ${\mathbf{x}}^{j+1}$.

The decomposition (16) is again applied only on the last time step before each maturity $T_{\alpha }$, layer number $j_{\alpha }$. Then,

${\displaystyle {\Gamma }_{\alpha }({\Sigma }_1,{\Sigma }_2)=\frac{1}{N}\sum _{\beta =1}^N{\left[\left({\mathbf{z}}_{\beta ;i_0}^{j_{\alpha }}+{\Sigma }_1^{j_{\alpha }}{\mathbf{w}}_{\beta ;i_0}^{j_{\alpha }}+{\Sigma }_2^{j_{\alpha }}{\mathbf{x}}_{\beta ;i_0}^{j_{\alpha }}\right)-{\omega }_{\beta }^{\alpha }\right]}^2{\chi }_{\beta }^{\alpha }.}$

Since we want the minimum of ${\Gamma }_{\alpha }$, the candidates are described using

$\left|\begin{array}{c}{\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial {\Sigma }_1}=0}\hfill \\ {\displaystyle \frac{\partial {\Gamma }_{\alpha }}{\partial {\Sigma }_2}=0}\hfill \end{array}\right..$

The system is written in the form

$\left|\begin{array}{c}{\displaystyle \underset{A}{\underbrace{\left(\sum _{\beta =1}^N{\left({\mathbf{w}}_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }_1^{j_{\alpha }}+\underset{B}{\underbrace{\left(\sum _{\beta =1}^N{\mathbf{w}}_{\beta ,i_0}^{j_{\alpha }}{\mathbf{x}}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }_2^{j_{\alpha }}=\underset{D}{\underbrace{-\sum _{\beta =1}^N\left({\mathbf{z}}_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }\right){\mathbf{w}}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }}}}\hfill \\ {\displaystyle \underset{B}{\underbrace{\left(\sum _{\beta =1}^N{\mathbf{w}}_{\beta ,i_0}^{j_{\alpha }}{\mathbf{x}}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }_1^{j_{\alpha }}+\underset{C}{\underbrace{\left(\sum _{\beta =1}^N{\left({\mathbf{x}}_{\beta ,i_0}^{j_{\alpha }}\right)}^2{\chi }_{\beta }^{\alpha }\right)}}{\Sigma }_2^{j_{\alpha }}=\underset{E}{\underbrace{-\sum _{\beta =1}^N\left({\mathbf{z}}_{\beta ,i_0}^{j_{\alpha }}-{\omega }_{\beta }^{\alpha }\right){\mathbf{x}}_{\beta ,i_0}^{j_{\alpha }}{\chi }_{\beta }^{\alpha }}}}\hfill \end{array}\right.$

If $\mathrm{DENOM}:=AC-B^2\ne 0$, then the solution to the system is given by

$\left|\begin{array}{c}{\displaystyle {\Sigma }_1^{j_{\alpha }}=\frac{BE-CD}{\mathrm{DENOM}}}\hfill \\ {\displaystyle {\Sigma }_2^{j_{\alpha }}=\frac{BD-AE}{\mathrm{DENOM}}}\hfill \end{array}\right..$

4. Computational Experiments

In this section, we will provide ample numerical simulations in order to demonstrate the application of the algorithm. First, we will solve the direct problem, and we will use its solution to gather observations, which in turn will be used to solve the inverse problems.

4.1. Direct Problem

For our synthetic data test, we take

• $r_{min}=-0.13$;

• $r_{max}=0.13$;

• $T=\overline{1,5}$ years;

• $\eta \left(\tau \right)=0.005-0.001\tau$;

• ${\displaystyle \gamma \left(\tau \right)=\frac{1+e^{-\tau }}{20}}$;

• $\sigma \left(\tau \right)=0.005+0.004ln(\tau +{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.})$.

When solving the direct problem, we take $I=161$ with $i_0=87$, i.e., $r_0=0.00975$. The temporal step is $▵\tau ={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$365$}\right.}$.

Without loss of generality, we consider a bond with the following setting—puttable and coupon bond. The coupon rate is fixed on 2.5% and the put dates are annual, starting one year after issuing and ending one year before the maturation of the bond. The same informational model could be used if the bond is non-puttable or zero-coupon.

Figure 3 shows the solution to (1), (3) and (5) by means of (12) and (13).

4.2. Implied Volatility

In this subsection, we aim to solve (IP1), i.e., to recover the unknown volatility function $\sigma \left(\tau \right)$ given market measurements $\left\{{\omega }_{\beta }^{\alpha }\right\}$, following the algorithm in Section 3.2.

In Figure 4, we can observe the real and reconstructed volatility functions. The proposed volatility is very close to the true one, with a discrepancy for the small values of $\tau$ due to the curvature of the real function, on the one hand, and the linear part of the recovered volatility, on the other hand.

To evaluate how measurement errors influence the volatility estimates, we perform a test using noisy observations. Specifically, we revisit (IP1), introducing Gaussian noise into the option price data $\left\{{\omega }_{\beta }^{\alpha }\right\}$. The perturbed observations are defined as

${\omega }_{\beta }^{\alpha }={\omega }_{\beta }^{\alpha }+{\epsilon }_{\beta }^{\alpha },\mathrm{where}{\epsilon }_{\beta }^{\alpha }\sim \mathcal{N}\left(0,{\left({\sigma }_{\beta }^{\alpha }\right)}^2\right).$

This approach assumes that, with $95\%$ confidence, the deviation in a given measurement does not exceed $p\%$ of its true value. For our simulations, we set $p=1$. The outcomes of these tests are illustrated in Figure 5.

The recovery in the ‘perturbed’ case is not as good as that in the ‘exact’ observations case, but it is capable of catching and following the trend in the volatility as a function of time.

4.3. Implied Parameters

In this subsection we solve (IP2) with the same values of the parameters, but now our goal is to recover not only the volatility, but also the reversion level and the reversion speed, all as functions of time. The results follow in Figure 6, Figure 7 and Figure 8.

All the functions, $\sigma \left(\tau \right)$, $\eta \left(\tau \right)$ and $\gamma \left(\tau \right)$ are accurately recovered. There are little inaccuracies in the proposed functions for small values, $\tau$, but again this is due to the curved nature of the real functions, which cannot be closely followed by the linear parts in the beginning.

4.4. Two-Dimensional Problem

In this subsection, we will solve the two-dimensional problem in (6)–(8), aiming to identify the two volatility functions.

For this problem, we set the following:

• $r_{min}=-0.13$;

• $r_{max}=0.13$;

• $u_{min}=-0.1$;

• $u_{max}=0.1$;

• $T=\overline{1,5}$ years;

• $a=0.1$;

• $b=0.08$;

• $\theta =0.05$;

• ${\sigma }_1\left(\tau \right)=0.05+0.04ln(\tau +{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.})$;

• ${\displaystyle {\sigma }_2\left(\tau \right)=\frac{1+e^{-\tau }}{10}}$;

• $\rho =-0.3$.

We take $I=53$ with $i_0=28$, i.e., $r_0=0.005$, and $K=51$ with $k_0=26$, i.e., $u_0=0$. The temporal step is $▵\tau ={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$365$}\right.}$. Again, the coupon rate is fixed at 2.5% and the put dates are annual, starting one year after the issuing and ending one year before the maturation of the bond.

The solution to the direct problem is given in Figure 9.

The solution to the 2D inverse problem is shown in Figure 10.

Again, the same implications could be made—the recovered functions ${\sigma }_1\left(\tau \right)$ and ${\sigma }_2\left(\tau \right)$ follow the real ones closely, identifying the trend and the changes in time.

4.5. Discussion

The numerical experiments conducted demonstrate the robustness and efficiency of the proposed algorithm for solving the inverse problems. For the one-dimensional case (IP1), the recovered volatility function closely matches the true function, successfully capturing the temporal trend despite minor discrepancies for small values of $\tau$. These deviations arise from the curved nature of the real volatility function, which cannot be perfectly approximated by the linear segments used in the reconstruction process. Nonetheless, the algorithm reliably identifies the underlying structure. The observed small inaccuracies are entirely approach-induced, and cannot be dealt with in the current setting.

When measurement noise is introduced, the algorithm exhibits remarkable resilience. While the recovery in the noisy case is less precise compared to that in the exact data scenario, the reconstructed volatility function still follows the overall trend accurately. This indicates that the method is robust to realistic perturbations in observed option prices, which is an essential feature for practical applications in financial markets.

The inverse problem (IP2), which involves recovering the reversion level $\eta \left(\tau \right)$ and reversion speed $\gamma \left(\tau \right)$ alongside volatility $\sigma \left(\tau \right)$, demonstrates the algorithm’s capability to handle more complex parameter identification tasks. The reconstructed functions align well with the true parameters, even in cases where their behavior is nonlinear. Similar to the volatility recovery, minor inaccuracies are observed for small values of $\tau$, which are attributable to the initial linear assumptions. However, the overall accuracy underscores the algorithm’s adaptability and reliability in reconstructing multiple time-dependent parameters.

For the two-dimensional problem, where two volatility functions ${\sigma }_1\left(\tau \right)$ and ${\sigma }_2\left(\tau \right)$ are identified, the algorithm continues to perform effectively. The reconstructed functions track the real ones closely, capturing the temporal trends and changes with high fidelity. This highlights the method’s scalability and applicability to higher-dimensional settings, providing a tool for analyzing more complex financial instruments.

With such levels of reconstruction, the identified price functions (these are the solution to the direct problem using the recovered values of the parameters) practically coincide with the observed ones as the differences are indistinguishable from zero for every $\tau$.

5. Conclusions

In this study, we suggested a robust and efficient algorithm for recovering the volatility term structure in the Hull–White model, demonstrating its applicability to real market scenarios. The predictor–corrector mechanism employed in our approach is central to its performance, starting with an initial assumption of constant volatility and iteratively refining it through linear forward steps. By correcting the volatility at a half-backward time level for all but the final step, the method constructs a piecewise linear volatility function that aligns well with observed market data. The simplicity of the algorithm is another key advantage, as it avoids the need to invert complex formulas or equations, instead solving for scalar parameters at each step, ensuring robustness and computational efficiency. Overall, the proposed algorithm lays a foundation for addressing complex volatility reconstruction problems in finance. Its robustness, efficiency, and adaptability to real market data highlight its practical relevance, while paving the way for future advancements in the calibration and application of sophisticated financial models.

Future work can extend this methodology in several directions to further enhance its applicability and accuracy. One potential avenue is the integration of regularization techniques, such as Tikhonov regularization, to improve the stability of the inverse problem, especially in the presence of noisy data. What is more, greater noise levels and/or types of noise other than Gaussian are interesting to test the algorithm with. Another promising extension involves adapting the algorithm to other models, such as the Heston stochastic volatility model [21], to address broader market conditions and applications. Additionally, applying the method to the Dupire equation [5] could allow for a more comprehensive analysis of local volatility surfaces. Expanding the scope to include American and exotic options is also of great interest, as it would enable the valuation of a wider range of financial instruments with satisfying precision, as well as other approaches in economics, e.g., the Evans model with two-scale price dynamics [22,23].

Author Contributions

Conceptualization, S.G.G. and L.G.V.; methodology, S.G.G. and L.G.V.; software, S.G.G.; validation, S.G.G. and L.G.V.; formal analysis, L.G.V.; investigation, S.G.G. and L.G.V.; resources, S.G.G. and L.G.V.; data curation, S.G.G. and L.G.V.; writing—original draft preparation, S.G.G. and L.G.V.; writing—review and editing, S.G.G. and L.G.V.; visualization, S.G.G.; supervision, L.G.V.; project administration, S.G.G.; funding acquisition, S.G.G. and L.G.V. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The authors confirm that the data supporting the findings of this study are available within the article.

Acknowledgments

The authors are grateful to the anonymous referees for the useful suggestions and comments.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 2. [Forumla omitted. See PDF.].

Figure 3. Bond price with coupon [Forumla omitted. See PDF.] and annual put dates.

Figure 4. True and recovered volatility in case of reconstructing one parameter. The circles denote the half-time layers (same as the following figures).

Figure 5. True and recovered volatility in case of reconstructing one parameter with perturbed observations.

Figure 6. True and recovered volatility in case of reconstructing three parameters.

Figure 7. True and recovered reversion level in case of reconstructing three parameters.

Figure 8. True and recovered reversion speed in case of reconstructing three parameters.

Figure 9. Bond price with coupon [Forumla omitted. See PDF.] and annual put dates in 2D.

Figure 10. True and recovered volatilities in case of reconstructing two parameters in 2D.

References

1. Hull, J.; White, A. Pricing interest-rate-derivative securities. Rev. Fin. Stud.; 1990; 3, pp. 573-592. [DOI: https://dx.doi.org/10.1093/rfs/3.4.573]

2. Black, F.; Scholes, M. The pricing of options and corporate liabilities. J. Polit. Econ.; 1973; 81, pp. 637-654. [DOI: https://dx.doi.org/10.1086/260062]

3. Merton, R.C. Theory of Rational Option Pricing. Bell J. Econ. Manag. Sci.; 1973; 4, pp. 141-183. [DOI: https://dx.doi.org/10.2307/3003143]

4. Derman, E.; Kani, I. Riding on a smile. Risk; 1994; 7, pp. 32-39.

5. Dupire, B. Pricing with a smile. Risk; 1994; 7, pp. 18-20.

6. MathConsult GmbH. Technical Paper: Calibration of Interest Rate Models. 2009; Available online: https://www.unrisk.com/wp-content/uploads/Calibration-of-Interest-Rate-Models.pdf (accessed on 15 December 2024).

7. Georgiev, S.; Vulkov, L. Fast reconstruction of time-dependent market volatility for European options. Comput. Appl. Math.; 2021; 40, 30. [DOI: https://dx.doi.org/10.1007/s40314-021-01422-9]

8. Jin, Y.; Wang, J.; Kim, S.; Heo, Y.; Yoo, C.; Kim, Y.; Kim, J.; Jeong, D. Reconstruction of the time-dependent volatility function using the Black–Scholes model. Discr. Dyn. Nat. Soc.; 2018; 2018, 3093708. [DOI: https://dx.doi.org/10.1155/2018/3093708]

9. Achdou, Y.; Pironneau, O. Computational Methods for Option Pricing; SIAM Frontiers in Applied Mathematics: Philadelphia, PA, USA, 2005.

10. Aichinger, M.; Binder, A. A Workout in Computational Finance; John Wiley & Sons Ltd.: Hoboken, NJ, USA, 2013.

11. Halgašová, J.; Stehlíková, B.; Bučková, Z. Estimating the short rate from the term structures in the Vasicek model. Tatra Mount. Math. Publ.; 2014; 61, pp. 87-103. [DOI: https://dx.doi.org/10.2478/tmmp-2014-0029]

12. Košútová, L.; Stehlíková, B. Calibration of the Ueno’s shadow rate model of interest rates. Mathematics; 2024; 12, 3564. [DOI: https://dx.doi.org/10.3390/math12223564]

13. Guterding, D.; Boenkost, W. The Heston stochastic volatility model with piecewise constant parameters—Efficient calibration and pricing of window barrier options. J. Comp. Appl. Math.; 2018; 343, pp. 353-362. [DOI: https://dx.doi.org/10.1016/j.cam.2018.04.054]

14. van der Zwaard, T.; Grzelak, L.A.; Oosterlee, C.W. On the Hull-White model with volatility smile for Valuation Adjustments. arXiv; 2024; arXiv: 2403.14841

15. Binder, A.; Jadhav, O.; Mehrmann, V. Model order reduction for parametric high dimensional models in the analysis of financial risk. arXiv; 2020; arXiv: 2002.11976

16. Russo, V.; Torri, G. Calibration of one-factor and two-factor Hull-White models using swaptions. Comput. Manag. Sci.; 2019; 16, pp. 275-295. [DOI: https://dx.doi.org/10.1007/s10287-018-0323-z]

17. Lee, Y.; Yang, K. Finite difference method for the Hull-White partial differential equations. Mathematics; 2020; 8, 1719. [DOI: https://dx.doi.org/10.3390/math8101719]

18. Georgiev, S.; Vulkov, L. Computational recovery of time-dependent volatility from integral observations in option pricing. J. Comput. Sci.; 2020; 39, 101054. [DOI: https://dx.doi.org/10.1016/j.jocs.2019.101054]

19. Itkin, A.; Muravey, D. Semi-closed form prices of barrier options in the Hull-White model. arXiv; 2020; arXiv: 2004.09591

20. Wong, H.-Y.; Zhao, J. An artificial boundary condition method for the Hull–White model of Americal interest rate derivatives. Appl. Math. Comp.; 2011; 217, pp. 4627-4643. [DOI: https://dx.doi.org/10.1016/j.amc.2010.11.015]

21. Vasilev, G.S. Time-dependent Heston model. arXiv; 2014; arXiv: 1402.5679

22. He, J.-H.; He, C.-H.; Sedighi, H.-M. Evans model for dynamic economics revised. AIMS Math.; 2021; 6, pp. 9194-9206. [DOI: https://dx.doi.org/10.3934/math.2021534]

23. Zabolotnii, S.; Mogilei, S. Modifications of Evans price equilibrium model. Inform. Autom. Pomiary Gospod. Ochr. Środowiska; 2023; 13, pp. 58-63. [DOI: https://dx.doi.org/10.35784/iapgos.3507]