1]SAIL Department, University of Applied Sciences Upper Austria 2]Google LLC

# Good Vibes! Towards Phone-to-User Authentication Through Wristwatch
Vibrations

Jakob Dittrich<EMAIL_ADDRESS>Rainhard Dieter Findling
<EMAIL_ADDRESS>[  [

###### Abstract

While mobile devices frequently require users to authenticate to prevent
unauthorized access, mobile devices typically do not authenticate to their
users. This leaves room for users to unwittingly interact with different
mobile devices. We present GoodVibes authentication, a variant of mobile
device-to-user authentication, where the user’s phone authenticates to the
user through their wristwatch vibrating in their pre-selected authentication
vibration pattern. We implement GoodVibes authentication as an Android
prototype, evaluate different authentication scenarios with 30 participants,
and find users to be able to well recognize and distinguish their
authentication vibration pattern from different patters, from unrelated
vibrations, and from the pattern being absent.

###### keywords:

Mobile device-to-user authentication, vibration pattern, hardware phishing
attack, phone, wristwatch

##  1 Introduction

Authentication with mobile devices is usually done from users towards devices
(user-to-device authentication), using methods such as PINs, passwords, unlock
patterns, or biometrics. That devices typically do not authenticate towards
their users leaves room for accidentally interacting with other devices, as
well as for targeted hardware phishing attacks [1]. Device-to-user
authentication [2], through which users can verify the authenticity of their
device, is one way to reduce the risk of such errors and attacks.

In this paper we present GoodVibes authentication, which is a variant of
device-to-user authentication with multiple mobile devices. When users start
to interact with their phone, thereby wake the device screen, the phone
authenticates to them by making the paired wristwatch vibrate in their pre-
selected authentication vibration pattern. We implement GoodVibes
authentication as an Android prototype and evaluate it with 30 participants to
answer the following research questions: RQ1: How well can users recognize
when a phone authenticates to them through a wristwatch vibration pattern, and
when this vibration pattern is absent, incorrect, or occurs without them
interacting with the phone? RQ2: How do users perceive the usability of a
phone authenticating to them through a wristwatch vibration pattern?

##  2 Related Work

Much prior work focuses on mobile user-to-device authentication through
knowledge, tokens, or biometrics [3, 4]. Such work does not address devices
authenticating to users, hence leaves room for users to accidentally interact
with other devices, and for targeted hardware phishing attacks [1]. Similar to
website-based phishing, hardware phishing attacks try to trick users to
interact with phishing devices that appear identical to the user owned
devices. When users attempt to unlock such a device with their PIN, password,
unlock pattern, or biometrics, they unwittingly disclose this authentication
secret to the attackers. As attackers are already in possession of the user
device from previously replacing it with the phishing device, they can
immediately use that secret to access the user device and the private data on
it.

Device-to-user authentication [2] aims to make it easier for users to
recognize if they would interact with another device than their own. Prior
research investigates phones authenticating to their users through the phone
itself vibrating in an authentication pattern when the user wakes the device
screen [1]. Findings indicate that users can recognize and distinguish such
vibrations (median success rate of 97.5%). However, physical access to the
user phone also gives attackers access to the vibration pattern. This allows
them to configure the phishing device with the same vibration. GoodVibes
authentication addresses this shortcoming, as attackers in control of the
locked user phone have no access to the authentication vibration pattern.

Other vibration-based authentication includes Vibrate-to-Unlock [5], which
uses vibration patterns on mobile phones to authenticate users to RFID tags.
While this approach uses vibrations, it does so to perform user-to-device
authentication, hence does not address the need for device-to-user
authentication. Prior work on website-to-user authentication investigates
watermarking [6] and embedding secrets in website content [7], which aims to
make it challenging for attackers to correctly mimic a legitimate website.
This requires users to visually inspect the website for the watermark or
embedded secret in order to detect phishing attempts.

##  3 GoodVibes Authentication

To enroll, the user installs the GoodVibes application on their paired phone
and wristwatch, and in the application selects the authentication vibration
pattern in which the wristwatch should vibrate to indicate that their phone is
authenticating to them.

Figure 1: Overview of the GoodVibes authentication process.

The GoodVibes authentication process works as follows (fig. 1): when the user
starts to interact with their phone, thereby wakes the device screen (1), the
phone sends a ping to the paired wristwatch over the secure channel that
exists as a result of the devices being paired (2). The wristwatch receives
the ping and vibrates in the user-chosen authentication vibration pattern (3).
The user feels the wristwatch vibrating in their authentication vibration
pattern, which informs them that they are interacting with their phone.

If the wristwatch would not vibrate, or vibrate in a different pattern, then
this would indicate to the user that they might be interacting with a phone
that is not theirs. Also, if the user would feel their authentication
vibration pattern without interacting with their phone, then this would
indicate to them that someone else might be interacting with their phone.

##  4 Evaluation

We implement GoodVibes authentication as a two-part prototype app111The
GoodVibes authentication prototype app is publicly available at
https://github.com/Jakob-Dittrich/GoodVibesAuth. for a paired Android phone
and Android wristwatch, which can communicate securely with each other over
their paired Bluetooth channel. The prototype app can also be controlled
remotely by the study supervisor to inject or suppress vibrations, to simulate
and evaluate different situational authentication scenarios (sec. 4.2). As
authentication vibration patterns we use the ”2” and the ”1 3” patterns from
[1], with vibration bursts of 60 ms, vibration pauses between bursts of 200
ms, and a total duration of 180 ms and 560 ms, respectively.

###  4.1 Evaluation Procedure

A supervisor provided instructions and monitored the session. Participants
answered demographic and technical knowledge questions. They then took a seat
at a desk in an empty office room where only supervisor and participant were
present, with the phone on the desk in front of them, and attached the
wristwatch to their preferred hand. 50% of participants then chose their
preferred authentication vibration pattern, while the other 50% were randomly
assigned an authentication vibration pattern. They were instructed to indicate
any case of the wristwatch not vibrating in their personal authentication
vibration pattern when they started to interact with the phone, and of the
wristwatch vibrating in their personal authentication pattern without them
having started to interact with the phone. They were then engaged in various
tasks and objectives, such as browsing a news website, playing a game, reading
a page from a book, or writing a message. This simulated situations of
concentration and exposure to distraction, and included points of purposefully
starting and stopping to interact with the phone. Participants did those tasks
for roughly 35 minutes, during which they were exposed to 5 different
situational authentication scenarios (sec. 4.2), 24 times in total (9, 6, 3,
3, and 3 times each), in a fixed order. They then completed a questionnaire
with 5 Likert-scale questions to evaluate the usability of GoodVibes
authentication. Participants were recruited for our evaluation study through
word of mouth among students of the University of Applied Sciences Upper
Austria and their relatives and friends. Participants received no compensation
for their participation. 30 participants (mean age 34.1 years, age std 12.3;
17 female, 13 male, 0 other) completed the study, resulting in 720
authentication results and 60 questionnaire answers.

###  4.2 Evaluation Scenarios

We evaluate GoodVibes authentication in 5 scenarios that cover different cases
of users possibly incorrectly identifying vibration patterns or their absence.

In scenario 1, when the user starts to interact with their phone, the
wristwatch vibrates in their personal authentication vibration pattern. We
measure how often users incorrectly recognize this as the vibration being
absent or as not being theirs.

In scenario 2, without the user starting to interact with their phone, the
wristwatch vibrates in a pattern that is not their personal authentication
vibration pattern, simulating unrelated notifications, alerts, or alike. We
measure how often users incorrectly recognize this as their personal
authentication vibration pattern.

In scenario 3, without the user starting to interact with their phone, the
wristwatch vibrates in their personal authentication vibration pattern. We
measure how often users correctly identify this as their personal
authentication vibration pattern.

In scenario 4, when the user starts to interact with their phone, the
wristwatch does not vibrate. We measure how often users correctly identify
this lack of vibration.

In scenario 5, when the user starts to interact with their phone, the
wristwatch vibrates in a different pattern than their personal authentication
vibration pattern. We measure how often users correctly identify this
mismatch.

##  5 Results and Discussion

Overall, participants recognized vibrations and their absence well. In
scenario 1, participants incorrectly thought the vibration was absent or not
their authentication vibration pattern 1% of the time. This indicates that
users can reasonably well recognize their authentication vibration pattern in
normal authentication situations, and that GoodVibes authentication seems to
cause a reasonably low false alert rate.

In scenario 2, participants recognized 97% of the time that the vibration was
unrelated to their personal authentication vibration pattern. Similar to the
results of scenario 1, this too indicates that users can reasonably well
distinguish between GoodVibes authentication vibrations and unrelated
vibrations. In the remaining cases (3%) users could easily quickly check their
phone’s status to ensure that it is actually secure. This seems a sufficiently
low false alert rate to not overly annoy GoodVibes authentication users, which
questionnaire findings seem to confirm (see below).

In scenario 3, participants recognized 98% of the time that the vibration was
their personal authentication pattern without them having started to interact
with their phone. This indicates that GoodVibes authentication would likely
help users to recognize that someone else started to use their phone, e.g.
another person mistaking the phone for theirs, or someone having stolen the
phone – and in turn help the user to take timely countermeasures.

In scenario 4 and 5, participants recognized 91% of the time if their
wristwatch did not vibrate, and 94% of the time if their wristwatch vibrated
in a pattern that was not their personal authentication vibration pattern.
Though those error rates are slightly higher than in the scenarios above, this
indicates that GoodVibes authentication would likely also help users to
recognize cases of them mistaking another phone as theirs, cases of hardware
phishing attacks, and cases of their wristwatch coincidentally vibrating for a
different reason at the same time. This indicates that GoodVibes
authentication could help users to not continue their user-to-device
authentication in those moments, thereby protecting their authentication
secret towards their phone.

Participants with prior smartwatch experience more often correctly recognized
vibrations or their absence than participants without prior smartwatch
experience. Over all scenarios combined, 9 participants who use a smartwatch
daily achieved 97% correct vibration recognition, and 8 participants who had
sometimes used a smartwatch achieved 99% – while 13 participants without prior
smartwatch experience achieved 89%. This indicates that users seem able to
learn to recognize and distinguish wristwatch vibrations, which in turn
indicates that GoodVibes authentication users without prior smartwatch
experience should be able to over time learn to distinguish their
authentication vibration pattern from other vibrations. Also, over all
scenarios combined, participants who could choose their personal
authentication vibration pattern achieved better results over participants who
were assigned a pattern (2% and 5% overall error rate, respectively). This
indicates that it is beneficial to let users choose which vibration pattern
they prefer to use with their GoodVibes authentication.

Questionnaire results indicate that participants found GoodVibes
authentication easy to use (mean 4.9, std 0.3, from 1 ”very difficult” to 5
”very easy”), fast to use (mean 5, std 0, from 1 ”very slow” to 5 ”very
fast”), and easy to adapt to (mean 4.9, std 0.3, from 1 ”very hard” to 5 ”very
easy”). Slightly more users considered it likely than unlikely that they would
use GoodVibes authentication if it were available today (mean 3.4, std 0.9,
from 1 ”very unlikely” to 5 ”very likely”).

In summary we can answer the posed research questions as follows: RQ1: users
seem to well recognize wristwatch vibration patterns through which their phone
would authenticate to them (99%), and also to well distinguish them from
unrelated vibrations (97%). Users furthermore seem to recognize well if their
authentication vibration pattern occurs without them interacting with their
phone (98%), as well as when they expect the pattern to occur but it instead
is absent (91%) or different (94%). RQ2: users seem to find their phone
authenticating to them through wristwatch vibrations easy to use (4.9/5), fast
to use (5/5), easy to adapt to (4.9/5), and slightly more users considered it
likely than unlikely that they would use such an approach (3.4/5).

## Conclusion

We presented GoodVibes authentication, where phones authenticate to their
users through vibrations of their paired wristwatch. We evaluated how well
users can recognize and distinguish such authentication vibration patterns, or
their absence, in 5 different authentication scenarios, and how users perceive
the usability of such authentication vibrations. Findings include that users
seem to be able to well recognize and distinguish their authentication
vibration pattern from other vibrations, their pattern occurring without them
interacting with their phone, as well as the pattern being absent or
different. Users seem to perceive the usability of their phone authenticating
to them over vibrations of their wristwatch as good.

Our work is limited by accounting only for a limited threat model without
sophisticated attackers, and by assessing vibration pattern recognizably and
usability with users who only used GoodVibes authentication for a short time
in a lab study environment. Future research could expand our work by
considering a wider threat model with different attackers, and by assessing
outside-lab vibration pattern recognizably and usability in a long-term real-
world study where participants are exposed to diverse distractions and
environmental factors that occur in their everyday lives.

## References

  *   * Findling and Mayrhofer [2015] Findling, R.D., Mayrhofer, R.: Towards device-to-user authentication: Protecting against phishing hardware by ensuring mobile device authenticity using vibration patterns. In: MUM 2015, pp. 131–135. ACM, New York, NY, USA (2015). https://doi.org/10.1145/2836041.2836053
  * Mayrhofer and Sigg [2021] Mayrhofer, R., Sigg, S.: Adversary models for mobile device authentication. ACM Computing Surveys 54(9) (2021) https://doi.org/10.1145/3477601
  * Meng et al. [2015] Meng, W., Wong, D.S., Furnell, S., Zhou, J.: Surveying the development of biometric user authentication on mobile phones. IEEE Communications Surveys & Tutorials 17(3), 1268–1293 (2015) https://doi.org/10.1109/COMST.2014.2386915
  * Shah and Kanhere [2019] Shah, S.W., Kanhere, S.S.: Recent trends in user authentication – a survey. IEEE Access 7, 112505–112519 (2019) https://doi.org/10.1109/ACCESS.2019.2932400
  * Saxena et al. [2011] Saxena, N., Uddin, M.B., Voris, J., Asokan, N.: Vibrate-to-unlock: Mobile phone assisted user authentication to multiple personal RFID tags. In: PerCom 2011, pp. 181–188 (2011). https://doi.org/10.1109/PERCOM.2011.5767583
  * Singh et al. [2011] Singh, A.P., Kumar, V., Sengar, S.S., Wairiya, M.: Detection and prevention of phishing attack using dynamic watermarking. In: AIM 2011, pp. 132–137. Springer, Berlin, Heidelberg (2011). https://doi.org/10.1007/978-3-642-20573-6_21
  * Varshney et al. [2012] Varshney, G., Joshi, R.C., Sardana, A.: Personal secret information based authentication towards preventing phishing attacks. In: ACITY 2012, pp. 31–42. Springer, Berlin, Heidelberg (2012). https://doi.org/10.1007/978-3-642-31513-8_4


[C] left=1.3in, right=1.3in, top=1in,bottom=1in, includefoot,
headheight=13.6pt

# Randomization of Short-Rate Models, Analytic Pricing and Flexibility in
Controlling Implied Volatilities

LECH A. GRZELAK<EMAIL_ADDRESS>Mathematical Institute, Utrecht University,
Utrecht, the Netherlands  Financial Engineering, Rabobank, Utrecht, the
Netherlands

###### Abstract

We focus on extending existing short-rate models, enabling control of the
generated implied volatility while preserving analyticity. We achieve this
goal by applying the Randomized Affine Diffusion (RAnD) method [11] to the
class of short-rate processes under the Heath-Jarrow-Morton framework. Under
arbitrage-free conditions, the model parameters can be exogenously stochastic,
thus facilitating additional degrees of freedom that enhance the calibration
procedure. We show that with the randomized short-rate models, the shapes of
implied volatility can be controlled and significantly improve the quality of
the model calibration, even for standard 1D variants. In particular, we
illustrate that randomization applied to the Hull-White model leads to
dynamics of the local volatility type, with the prices for standard
volatility-sensitive derivatives explicitly available. The randomized Hull-
White (rHW) model offers an almost perfect calibration fit to the swaption
implied volatilities.

###### keywords:

Randomization, RAnD Method, Randomized Hull-White Model (rHW), Stochastic
Parameters, Stochastic Collocation, Quantization, Short-Rate Models

††footnotetext: The views expressed in this paper are the author’s personal
views and do not necessarily reflect the views or policies of his current or
past employers. The author has no competing interests.

##  1 Introduction

Choosing a pricing model, especially in risk management and for Valuation
Adjustments (xVA), is a balancing act between accuracy and pricing speed.
Often, the calculation speed is critical in determining the model of choice,
especially when handling large portfolios. The modern approach for modelling
interest rates is based on the Heath, Jarrow, and Morton (HJM) [13] framework
with an arbitrary term-structure of volatility and covariance of forward rates
across maturities. The framework is generic, and essentially any sensible term
structure dynamics can be modelled by the HJM technique. Unfortunately, the
number of models giving rise to analytical formulae for option pricing or even
zero-coupon bonds is minimal. Therefore, models that can generate realistic
implied volatilities while allowing for efficient derivative pricing are in
high demand, especially when the portfolios need to be evaluated multiple
times and for different market scenarios due to regulatory pressure.

A family of the HJM models that has gained particular acceptance from
practitioners and academia is the family of affine models [7]. Over many
years, various affine models have been studied, leading to significant
contributions regarding efficient pricing, simulation and calibration (for the
overview, see [4]). In particular, the classic affine short-rate models, like
the Hull-White [14] model, are popularised due to the closed formula for zero-
coupon bonds and semi-analytic swaption pricing. On the downside, they lack
sufficient flexibility to calibrate to market-implied volatilities. Additional
model flexibility without compromising numerical efficiency is highly desired.
A straightforward approach to incorporate and control the implied volatility
smile and skew into the short-rate models is to define a stochastic or a local
volatility process, for example, in [6, 8, 9]. Although such models give us a
closed-form expression for the zero-coupon bond, pricing more advanced
derivatives like options often involves advanced numerical techniques like
Fourier inversion or numerical integration, which is undesired for the pricing
of high-volume derivatives. A survey of the most popular stochastic volatility
models in the interest rate world can be found in [1].

This article focuses on extending existing interest rate models with model
parameters to be random via the so-called randomization method (RAnD) [11].
Randomization by making the model parameters stochastic is a technique for
more flexibility in the stochastic model-thus, improving the calibration
quality while preserving the analytic properties of the models considered. The
concept of randomization represents the uncertainty of potentially hidden
states that are not sufficiently captured by deterministic parameters.
Conventionally, under the affine framework, an extension by a stochastic
parameter would require the model to meet linearity conditions - this does not
have to be the case under RAnD. The method requires an affinity to hold, given
a specific realization of the stochastic parameter. In other words, it builds
an outer layer over the affine models and allows a stochastic parameter
setting in that layer. The quadrature rule handles the numerical complexity
associated with an infinite number of parameter values and reduces it to only
a few critical parameter realizations. The selection of these points is based
on the moments of the stochastic parameter.

The RAnD method is generic and can be applied to any pricing model. For the
sake of simplicity, however, we focus on one of the most popular interest rate
models, namely the Hull-White (HW) model [14]. In particular, we will show
that applying the RAnD method to the HW model results in a local-volatility
type of dynamics. The local volatility dynamics may seem like a “deal-
breaker”; typically, the local volatility models, especially in interest
rates, are notoriously tricky to operate, i.e., closed-form option pricing is
hardly possible. However, although the randomized Hull-White (rHW) model is of
the local volatility type, the option pricing is as complex as the standard,
unrandomized model. Furthermore, the presented model offers fast and accurate
calibration as the local volatility dynamics provide more flexibility in
generating volatility smiles and skews for various market conditions.

Our method also addresses a challenging problem stated by Brigo and Mercurio
in [3], where the mixture of lognormals was studied. In their paper, it was
shown that great flexibility to standard models like Black-Scholes was
introduced by considering a convex combination of the associated probability
density functions:
$\omega_{1}f_{1}(x;\sigma_{1})+\dots+\omega_{N}f_{N}(x;\sigma_{N})$ with
$\omega_{1}+\dots+\omega_{N}=1$ and some $\sigma_{i}>0$, $i=1,\dots,N.$
However, the problem encountered was the large number of the model parameters
that are difficult to interpret and relate to the corresponding implied
volatilities. As stated in their work: “the absence of bounds on the parameter
$N$ implies that a virtually unlimited number of parameters can be introduced
in the dynamics so as to be used for a better calibration to market data. (…),
one has to find the correct tradeoff between model flexibility and number of
parameters so as to avoid both poor calibration and over-parametrization
issues”. However, this problem is resolved with the RAnD method. The method
transitions from a continuous randomizer $\vartheta$ to a unique sequence of
weights $\omega_{i}$ and associated values $\sigma_{i}$. In this setting, the
model calibration is performed by varying the parameters of the randomizer,
$\vartheta$; therefore, the number of free parameters is manageable. Moreover,
as presented in this work, the varying parameters of $\vartheta$ enable a
transparent control of the implied volatilities.

This article comprises five sections. First, the details of the HJM framework
and pricing equations using the RAnD method are provided in Section 2, and a
discussion on the arbitrage-free conditions is covered in Section 2.1. In
Section 2.2, we consider the rHW model and derive the dynamics of the
randomized stochastic processes. Section 2.3 discusses different variants of
the rHW model parameter choices, namely the randomization using univariate or
bivariate distributions. Section 3 focuses on derivative pricing; in
particular, the pricing equations for options on zero-coupon bonds and
swaptions are provided. Numerical results are presented in Section 4, where in
Section 4.1, the swaption implied volatility surfaces for HW and rHW models
are compared. Section 4.2 analyzes the randomizers’ impact on controlling
implied volatilities for swaptions. Next, the calibration results are
presented in Section 4.3, and in Section 4.4, the extension to the bivariate
randomization case is established. Finally, the convergence results are
illustrated in Section 4.5 and Section 5 concludes.

##  2 The Heath-Jarrow-Morton models, affinity and randomization

Here, we will give an informal presentation of the Heath-Jarrow-Morton (HJM)
framework and its relation to arbitrage-free short-rate models. We define the
probability space $\left(\Omega,\mathcal{F}(t),\mathbb{Q}\right)$ on which the
HJM arbitrage-free condition requires that the instantaneous forward rates,
$f(t,T)$, are modelled by the following arbitrage-free SDE:

${\rm d}f(t,T)=\alpha(t,T){\rm d}t+\gamma(t,T){\rm
d}W(t),\;\;\;\alpha(t,T)=\gamma(t,T)\int_{t}^{T}\gamma(t,z){\rm d}z.$ (2.1)

The instantaneous forward rate, $f(t,T)$, is determined by the volatility
driver $\gamma(t,T)$, which can be defined as a deterministic function,
stochastic process or by a random variable. The case where $\gamma(t,T)$ is
defined as a random variable, thus time-invariant, will be associated with
parameter randomization.

We focus here on the short-rate dynamics under the HJM framework. Suppose that
$f(0,t)$, $\alpha(t,T)$ and $\gamma(t,T)$ are differentiable in their second
argument, with $\int_{0}^{T}\left|\frac{\partial}{\partial
t}f(0,t)\right|<\infty,$ then the short-rate process under the HJM framework
is given by,

$\displaystyle{\rm d}r(t)=\zeta(t){\rm d}t+\gamma(t,t){\rm d}W(t),$ (2.2)

with

$\zeta(t)=\alpha(t,t)+\frac{\partial}{\partial
t}f(0,t)+\int_{0}^{t}\frac{\partial}{\partial t}\alpha(z,t){\rm
d}z+\int_{0}^{t}\frac{\partial}{\partial t}\gamma(z,t){\rm d}W(z).$

A particular class of popular models are the models belonging to the Affine-
Diffusion interest rate models. It is shown [7] that in this class, for
process ${{r}}(t)$, the discounted characteristic function is of the following
form:

$\phi_{{r}}({{u}};t,T)=\mathbb{E}_{t}\left[\exp\left(-\int_{t}^{T}r(s){\rm
d}s+i{u}{{r}}(T)\right)\right]=\exp\left(A({u};\tau)+{{B}}({u};\tau){{r}}(t)\right),$
(2.3)

with the expectation under risk-neutral measure $\mathbb{Q}$ for $\tau=T-t$,
where $A({u};\tau)$ and ${{B}}({u};\tau)$ form the fundamental solution to the
corresponding pricing PDE and satisfying the complex-valued ODEs, see the work
by Duffie-Pan-Singleton [7]. The affinity conditions are rather strict as they
impose a linear structure on the model, i.e., the model dynamics and its
covariance structure needs to be linear in the state variables. The
representation in (2.3) is of particular importance as it gives an explicit
closed-form expression for the integrated short-rate process and therefore
determines the price of a zero-coupon bond (ZCB). Models not leading to an
analytic ZCB price are rarely used for pricing purposes. Unfortunately,
because of the affinity constraints, the model extensions for some model
parameters to be stochastic are minimal. In this article, we focus on relaxing
this constraint and letting the model parameters be random while benefiting
from the closed-form solutions for the ZCBs and pricing, as presented in
(2.3).

We consider a vector $\Theta=[{\vartheta}_{1},\dots,{\vartheta}_{d}]^{\rm T}$,
with $d\in\mathbb{N}$ representing the number of randomized parameters, where
each ${\vartheta}_{i}$ is a possibly correlated, time-invariant, random
variable. A realization of ${\vartheta_{i}}$ is indicated by $\theta_{i}$,
$\vartheta_{i}(\omega)=\theta_{i}$, and consequently the realization for
${\Theta}$ is indicated by $\theta=[\theta_{1},\dots,\theta_{d}]^{\rm T}$.
Under the RAnD framework, parameters are not driven by a stochastic process
but only by a random variable; however, an extension to piece-wise stochastic
parameters is also possible (see [11] for details).

In a nutshell, the RAnD method relies on quadrature integration and
conditional expectation. Therefore, the technique can be applied wherever we
deal with an expectation. In principle, randomization can be applied to
Equation (2.3), which would provide us with a randomized ZCB; on the other
hand, the randomization can be applied at the valuation level, where the
expectation of the expected payoff is computed. We prefer the latter approach
to construct the randomized variant of any interest rate model by simply
evaluating the model on a particular realization of the model parameters. This
will be particularly important when pricing non-standard derivatives.

In the randomized HJM framework, the volatility in (2.1) may be dependent on
random parameters $\vartheta_{1},\dots,\vartheta_{d}$, i.e.,
$\gamma(t,T;\\{\vartheta_{1},\dots,\vartheta_{d}\\}).$ Let us now consider a
stochastic model under the HJM framework for derivative pricing with the price
denoted as $V(t,r(t;\theta))$ with $\theta$ indicating a vector-valued
realization of the model parameter, $\Theta$. $V(t,r(t;\theta))$ may
correspond to a model of Vašiček, Hull-White, Black-Karasiński or any other
model that belongs to the class of HJM [4]. Then, the pricing under randomized
parameter space will be given by the following relation:

$\displaystyle
V(t,r(t;\Theta)):=\mathbb{E}_{t}[V(t,r(t;\Theta=\theta))]=\int_{\mathbb{R}^{d}}V(t,r(t;\Theta=\theta))f_{\Theta}(\theta){\rm
d}\theta,$ (2.4)

where the integration takes place over the parameter space
$\Theta\in\mathbb{R}^{d}$ and $f_{\Theta}(\theta)$ corresponds to the
associated $d-$dimensional probability density function. Equation (2.4) shows
that a random model parameter can be considered averaging over possible
parameter realizations. From a numerical perspective, the integration in (2.4)
is expensive.

The randomization enables for randomizing $d-$dimensional space of parameters;
however, for simplicity, we will mainly focus on single or bivariate parameter
randomization, $\vartheta_{1}(\omega)=\theta_{1}$ and
$\vartheta_{2}(\omega)=\theta_{2}$, thus also letting for the correlation
between stochastic parameters.

Our aim is twofold: we aim to provide numerically efficient methods for the
computation of the pricing Equation in (2.4) and to determine the stochastic
differential equation for the randomized short-rate processes,
$r(t;\vartheta)$ where $\vartheta$ is the chosen random parameter. The insight
into the dynamics of the randomized stochastic process will help us categorise
the randomized processes and the structure of their drift and volatility
coefficients.

######  Theorem 2.1 (Pricing formula for the randomized model).

Consider a random variable $\vartheta$, defined on some finite domain
$D_{\vartheta}:=[a,b]$, with its PDF, $f_{\vartheta}(x)$, CDF,
$F_{\vartheta}(x)$ and a realization $\theta$, $\vartheta(\omega)=\theta$ such
that for some $N\in\mathbb{N}$ the moments are finite,
$\mathbb{E}[\vartheta^{2N}]<\infty$. Let $V(t,r(t;\vartheta))$ be the price of
a financial derivative depending on the short-rate $r(t;\theta)$ with
parameter $\theta$, then the randomized price, $V(t,r(t;\vartheta))$ is given
by:

$\displaystyle V(t,r(t;\vartheta))=\int_{[a,b]}V(t,r(t;\theta)){\rm
d}F_{\vartheta}(\theta)=\sum_{n=1}^{N}\omega_{n}V(t,r(t;\theta_{n}))+\epsilon_{N},$
(2.5)

where the pairs $\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$ are the Gauss-
quadrature weights and the nodes based on the parameter cumulative
distribution function, $F_{\vartheta}(\cdot),$ determined by
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$
defined in A.2 and where the error $\epsilon_{N}$ is defined as:

$\epsilon_{N}=\frac{1}{(2N)!}\frac{\partial^{2N}}{\partial\xi^{2N}}V(t;r(t,\vartheta=\xi)),\;\;\;a<\xi<b.$
(2.6)

###### Proof.

Define a payoff depending on the short-rate $r(t)$, $0<t<T$, by $\chi(r,T)$,
then, for a stochastic parameter $\vartheta$ with the help of conditional
expectation, we find:

$\displaystyle
V(t,r(t;\vartheta))=\mathbb{E}_{t}\left[\mathbb{E}_{t}\left[\frac{M(t)}{M(T)}\chi(r,T)\Big{|}\vartheta=\theta\right]\right]=:\mathbb{E}_{t}\left[V(t,r(t;\theta))\right],$
(2.7)

with $M(t)$ being the money-savings account. By application of the quadrature
rule, Equation (2.7) becomes:

$\displaystyle V(t,r(t;\vartheta))=\int_{[a,b]}V(t,r(t;\theta)){\rm
d}F_{\vartheta}(\theta)=\sum_{n=1}^{N}\omega_{n}V(t,r(t;\theta_{n}))+\epsilon_{N}.$

Here, the error is given by $\epsilon_{N}$, for $a<\xi<b$, and the pair
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$ forms the Gauss-quadrature weights and
the nodes based on the parameter distribution, as proven in [5] (p.180,
Theorem 3.6.24). ∎

Theorem 2.1 gives us an explicit relation between a continuous randomizer
$\vartheta$ and the discretization in terms of pairs
$\\{w_{n},\theta_{n}\\}_{n=1}^{N}$, leading to a simplification of the
integral in (2.4). In this paper, these pairs are based on moments of the
randomizer, $\vartheta,$ and, for completeness, the detailed computation
procedure is given in A.2.

The associated error $\epsilon_{N}$ in (2.6) can be interpreted as a cost when
moving from the continuous to a discretized random parameter,
$\vartheta\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$. To some extent,
the pricing with a discrete set of parameter realizations, $\theta_{n}$, and
the associated probabilities, $w_{n}$, resembles the so-called regime-switch
method where a finite set of states is defined. However, from a practical
perspective, it is challenging to deal with $N$ pairs of the possible variable
states and the associated probabilities, especially in the calibration
procedure.

The procedure described in A.2, however, simplifies this problem as we can
control the quadrature pairs employing the parameters of the randomizer
$\vartheta$. To illustrate this process, let us consider a stochastic
parameter with a randomizer $\vartheta\sim\mathcal{N}(\mu,\sigma^{2})$; thus,
the stochastic parameter, $\vartheta$, is driven by two parameters $\mu$ and
$\sigma^{2}$. By application of A.2, for any set of parameters $\mu$ and
$\sigma^{2}$, we compute the corresponding pairs
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$; therefore, in the calibration
procedure we will only vary $\mu$, and $\sigma^{2}$. This procedure
drastically reduces the number of associated model parameters, thus
facilitating the calibration. We also stress that the set of optimal pairs is
solely obtained based on the moments of the randomizing random variable,
$\vartheta$, implying that every random variable with, preferably, closed-form
moments may be used for randomization. In Table 2, a few selected random
variables and their moments are tabulated. In Section 4.3, we will also
classify the randomizers, $\vartheta$, for which the computations can be
greatly reduced or even tabulated.

To emphasize the flexibility of the RAnD method, we focus on the application
in the Hull-White model, where we will apply randomization to each model
parameter and consider an extension to a bivariate case. In the following
sections, we consider a finite number of realizations of $\vartheta$. We will
denote them as $\theta_{1},\dots,\theta_{N}$, for some $N\in\mathbb{N}$. These
“specific” realizations we will interchangeably call either “collocation” [12]
or “quadrature” points. Finally, by $\vartheta(\hat{a},\hat{b})$, we denote
that $\vartheta$ is a random variable driven by parameters: $\hat{a}$,
$\hat{b}$.

###  2.1 Arbitrage-free conditions under the RAnD method

Before we analyze the specifics of the randomized Hull-White (rHW) model, let
us review the implications of randomization on the pricing, under the affine-
diffusion framework, of ZCB and discuss arbitrage-free aspects of the RAnD
method. By setting $u=0$ in (2.3), the risk-free pricing formula for a
randomized ZCB, $P(t,T;\vartheta)$ is given by:

$\displaystyle P(t,T;\vartheta)=\int_{[a,b]}P(t,T;\theta){\rm
d}F_{\vartheta}(\theta)=\int_{[a,b]}{\rm
e}^{A(\tau;\theta)+{{B}}(\tau;\theta){{r}}(t;\theta)}{\rm
d}F_{\vartheta}(\theta),\;\;a<b,$ (2.8)

with $A(\tau;\theta):=A(0,\tau;\theta)$ and $B(\tau;\theta):=A(0,\tau;\theta)$
in (2.3) for $\tau=T-t,$ and $r(t;\theta)$ indicates a short-rate model with
constant parameter realizations $\theta.$ By application of Theorem 2.1, the
randomized ZCB is known explicitly, and it is presented in Corollary 2.1.

######  Corollary 2.1 (ZCBs under Randomized Affine Jump Diffusion Processes).

Let ${r}(t;\theta)$ represent an affine short-rate process with some constant
parameter $\theta$. Assuming that the corresponding ChF,
$\phi_{r(T)|\vartheta=\theta}(\cdot)$, is well defined and $2N$ times
differentiable w.r.t. $\theta$, the randomized ZCB is given by,

$\displaystyle P(t,T;\vartheta)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\omega_{n}P(t,T;\theta_{n})+\epsilon_{N}=\sum_{n=1}^{N}\omega_{n}{\rm
e}^{A(\tau;\theta_{n})+{{B}}(\tau;\theta_{n})r(t;\theta_{n})}+\epsilon_{N}$
(2.9) $\displaystyle=:$ $\displaystyle
P(t,T;\\{\theta_{n}\\}_{n=1}^{N})+\epsilon_{N},$

with ${A}(\tau;\theta_{n})$, ${{B}}(\tau;\theta_{n})$ being the real-valued
functions obtained from Riccati-type of ODEs available for affine models [16].

###### Proof.

The proof is a consequence of Theorem 2.1. ∎

Representation (2.9) illustrates that the ZCB, $P(t,T;\vartheta)$, randomized
with stochastic parameter $\vartheta$ can be expressed as a weighted sum of
the ZCBs evaluated at a given realization of $\vartheta$ and where
$\omega_{1}+\dots+\omega_{N}=1$. Note that under the HJM framework, each ZCB,
$P(t,T;\theta_{n})$ in (2.9) is arbitrage-free; therefore, the convex
combination is too. This can be shown by checking whether
$P(S,T;\\{\theta_{n}\\}_{n=1}^{N})/M(S)$ is indeed a martingale for some
$t<S<T$:

$\displaystyle\mathbb{E}_{t}\left[\frac{P(S,T;\\{\theta_{n}\\}_{n=1}^{N})}{M(S)}\right]=\mathbb{E}_{t}\left[\frac{\sum_{n=1}^{N}\omega_{n}P(S,T;\theta_{n})}{M(S)}\right]=\sum_{n=1}^{N}\omega_{n}\mathbb{E}_{t}\left[\frac{P(S,T;\theta_{n})}{M(S)}\right],$

since every discounted ZCB, $P(S,T;\theta_{n})/M(S),$ is a martingale under
the $\mathbb{Q}$ measure we have,

$\displaystyle\mathbb{E}_{t}\left[\frac{P(S,T;\\{\theta_{n}\\}_{n=1}^{N})}{M(S)}\right]=\sum_{n=1}^{N}\omega_{n}\frac{P(t,T;\theta_{n})}{M(t)}=\frac{P(t,T;\\{\theta_{n}\\}_{n=1}^{N})}{M(t)}.$

This result shows that the application of the RAnD method enables the
transition from a continuous random variable to a discrete one; the cost of
this transition is expressed by $\epsilon_{N}$. Moreover, as shown above, the
error $\epsilon_{N}$ in (2.9) does not impact the arbitrage conditions-the
pricing under the RAnD method is simply a weighted average of arbitrage-free
prices. However, from a practical perspective, keeping this error as small as
possible is still essential. This is particularly important when calibrating
the randomized model to market data, i.e., the calibration procedure is
associated with varying parameters of the variable $\vartheta$. It is,
therefore, important that the discretized version resembles the continuous
version as closely as possible. In the numerical section, Section 4.5, we will
investigate the convergence aspects of different randomizers and an optimal
number of expansion terms, $N$.

###  2.2 The randomized Hull-White (rHW) short-rate model

This section provides a specification leading us to the randomized version of
the famous HW model [14]-a single-factor, no-arbitrage yield curve model, in
which an extended Ornstein-Uhlenbeck drives the short-term interest rate mean-
reverting process. Under the HJM framework and the arbitrage-free condition
for the drift in (2.1), the rHW model is specified by:

$\gamma(t,T)=\eta\cdot{\rm e}^{-\lambda(T-t)},\;\;\;t<T,$ (2.10)

where we consider three different randomization cases: the randomization of
the volatility parameter, $\eta$, the mean-reversion, $\lambda$, or the
randomization of both parameters using bivariate distribution:

$\eta\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta_{1},\;\;\text{or}\;\;\;\lambda\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta_{2},\;\;\;\text{or}\;\;\;\lambda|\eta\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta_{2}|\vartheta_{1}.$ (2.11)

To simplify the notation, we consider one parameter
$\theta\in\\{\eta,\lambda\\}$ with the corresponding random variable
$\vartheta$. Particular choices of the parameters will be given explicitly.

Given the HJM volatility in (2.10), under the risk-free measure, the dynamics
of the HW model read:

${\rm d}{r}(t)=\lambda(\psi(t)-{r}(t)){\rm d}t+\eta{\rm
d}W(t),\;\;\;r_{0}\equiv f(0,0),$ (2.12)

with

$\psi(t)=f(0,t)+\frac{1}{\lambda}f(0,t)+\frac{\eta^{2}}{2\lambda^{2}}\left(1-{\rm
e}^{-2\lambda t}\right),\;\;\;f(0,t)=-\frac{\partial\log P(0,t)}{\partial t},$
(2.13)

where $\psi(t)$ is a time-dependent drift term, which ensures the fit of the
calibration to a yield curve observed in the market, $W(t)$ is the Brownian
motion under measure $\mathbb{Q}$, and $f(0,t)$ indicates the instantaneous
forward rate computed from the yield curve, defined in terms of ZCBs,
$P(0,t)$. Parameter $\eta$ determines the overall level of the volatility, and
$\lambda$ is the reversion rate parameter. A large value of $\lambda$ causes
short-term rate movements to dampen out rapidly, reducing the long-term
volatility. Because of this interdependence between the model parameters,
often, in practical applications, $\lambda$ is fixed and set to a constant. In
contrast, $\eta$ is often set to be piece-wise constant and calibrated such
that the ATM volatilities are calibrated. We will show that under the RAnD
method, such a strategy is inadequate.

The short-rate $r(t)$ in (2.12) is thus normally distributed with
$r(t)\sim\mathcal{N}(\mu_{r(t)},\sigma^{2}_{r(t)}),$ with

$\displaystyle\mu_{r(t)}:=r_{0}{\rm e}^{-\lambda
t}+\lambda\int_{0}^{t}{\psi(z){\rm e}^{-\lambda(t-z)}{\rm
d}z},\;\;\;\sigma^{2}_{r(t)}:=\frac{\eta^{2}}{2\lambda}\left(1-{\rm
e}^{-2\lambda t}\right).$ (2.14)

Moreover, for $\psi(t)$ constant, i.e., $\psi(t)\equiv\psi$ (in this case we
deal with the Vašiček model [18]), we have
$\lim_{t\rightarrow\infty}\mathbb{E}_{t_{0}}\left[r(t)\right]=\psi.$ This
means that the first moment of the process converges to the mean-reverting
level $\psi$, for large values of $t$.

As a first step, we check the impact of randomization on the probability
density function, PDF, of the interest rate process $r(t)$ in (2.12). Since
the HW model belongs to the affine class of processes, we can benefit from the
available ChF. Given a stochastic parameter ${\bf\vartheta}$, the ChF is given
by:

$\displaystyle\phi_{{r}}({u};t,T):=\mathbb{E}_{t}\Big{[}{\rm
e}^{-\int_{t}^{T}r(s){\rm
d}s+i{u}{r}(T)}\Big{]}=\mathbb{E}_{t}\left[\mathbb{E}_{t}\Big{[}{\rm
e}^{-\int_{t}^{T}r(s){\rm
d}s+i{u}{r}(T)}\big{|}{\vartheta}={\theta}\Big{]}\right].$

By application of the quadrature rule in Theorem 2.1 and A.2, we find:

$\displaystyle\phi_{{r}}({u};t,T)=\int_{[a,b]}\phi_{{r}|{\vartheta}=\theta}({u};t,T){\rm
d}F_{\vartheta}(\theta)=\sum_{n=1}^{N}\omega_{n}\phi_{{{r}(T)}|\vartheta=\theta_{n}}({u};t,T)+\epsilon_{N}^{F},$
(2.15)

and by utilizing the Fourier transform, the PDF of the rHW model reads,

$\displaystyle f_{{{r}(T)}}(x)=\frac{1}{2\pi}\int_{\mathbb{R}}{\rm
e}^{-iux}\sum_{n=1}^{N}\omega_{n}\phi_{{{r}}|\vartheta=\theta_{n}}({u};t,T){\rm
d}u+\epsilon_{N}^{F}=\sum_{n=1}^{N}\omega_{n}f_{{r(T;\theta_{n})}}(x)+\epsilon_{N}^{F},$

where $r(t)|\vartheta=\theta_{n}$ indicates the HW model with a particular
realization, $\theta_{n}$, of the randomized variable. The representation
above presents the relation between the densities of the rHW process, $r(t)$,
as a convex combination of the HW processes, $r(t)|\theta_{n}$ with a
particular parameter $\theta_{n}$. We also highlight that the error
$\epsilon_{N}^{F}$ in (2.15) may differ from the quadrature error in (2.6).

Because of different randomization choices, we distinguish the following
randomization types and their associated error:

$f_{{{r}(T)}}(y)=\sum_{n=1}^{N}\omega_{n}f_{{r(T;\theta_{n})}}(y)+\epsilon_{N}^{F}=:\left\\{\begin{array}[]{ccc}f_{\overline{r}(T)}(y)+\overline{\epsilon}_{N},\;\;\;\text{for}\;\;\;\eta\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta,\\\
f_{\widetilde{r}(T)}(y)+\widetilde{\epsilon}_{N},\;\;\;\text{for}\;\;\;\lambda\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta,\\\
f_{\widehat{r}(T)}(y)+\widehat{\epsilon}_{N},\;\;\;\text{for}\;\;\;\eta\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta_{1}\;\;\&\;\;\lambda|\eta\stackrel{{\scriptstyle{\rm
d}}}{{=}}\vartheta_{2}.\end{array}\right.$ (2.16)

Under the RAnD method, the PDF, $f_{{{r}(T)}}(y)$ in (2.16) is given as a
linear combination of normal densities that depend on different parameter
realizations $\theta_{n}$, resembling a similar problem as presented in [3].

In the next section, we associate the PDF given in (2.16) and find the
corresponding SDE for a 1D stochastic representation.

###  2.3 Dynamics of the rHW model

The combination of Equations (2.11) and (2.12) shows us that the simulation of
the randomized short-rate models is explicit using the Monte Carlo technique.
However, by an explicit form of the corresponding SDE, we will get more
insight into the randomization and its impact on the dynamics of the driving
process.

As a start, under the rHW model, consider the randomization of the volatility
parameter, $\theta:=\eta$, with $\vartheta(\omega)=\theta$, and where the
corresponding quadrature pairs are given by
$\\{\omega_{n},\eta_{n}\\}_{n=1}^{N}$. For such a setting, under each
realization $\eta_{n}$, we have the associated HW model,
$\overline{r}_{n}(t):=\overline{r}_{n}(t;\eta_{n})$, with the dynamics given
by:

$\displaystyle{\rm
d}\overline{r}_{n}(t)=\lambda(\overline{\psi}_{n}(t)-\overline{r}_{n}(t)){\rm
d}t+\eta_{n}{\rm
d}W(t),\;\;\;\overline{r}_{n}(t_{0})=f(0,0),\;\;\;n=1,\dots,N,$ (2.17)

with common Brownian motion, $W(t)$, and initial value
$\overline{r}_{n}(t_{0})=f(0,0)$ in (2.13), where $\lambda$ is constant and
equal for all the $\overline{r}_{n}(t)$, and where $\overline{\psi}_{n}(t)$ is
defined in (2.13) and is given as a function of $\lambda$ and $\eta_{n}$.

Given the sequence of HW model processes,
$\overline{r}_{1}(t),\dots,\overline{r}_{N}(t)$, in (2.17) and the probability
density relation in (2.2), we consider the problem of finding the
corresponding SDE for the rHW process, $\overline{r}(t)$. Formally, we seek an
SDE, with the solution given in (2.16) and where each of the constituent
processes, $\overline{r}_{n}(t)$, is driven by (2.17). Thus, we consider the
following process,

$\displaystyle{\rm d}\overline{r}(t)=\overline{\lambda}(t,\overline{r}(t)){\rm
d}t+\overline{\eta}(t,\overline{r}(t)){\rm
d}W(t),\;\;\;\overline{r}(t_{0})=f(0,0),$ (2.18)

with some state-dependent drift, $\overline{\lambda}(t,\overline{r}(t))$, and
volatility, $\overline{\eta}(t,\overline{r}(t))$, and where Brownian motion
$W(t)$ is as in (2.17). As indicated in [2], the problem of finding an SDE
when marginal distributions and the corresponding weights are given is the
reverse of finding the marginal density function of the solution of an SDE
when the coefficients are known. Moreover, the problem of a mixture of normal
or lognormal processes is known [2]. Here, in Proposition 2.1, we adapt those
techniques and provide an explicit form for the SDE in (2.18).

######  Proposition 2.1 (Local volatility process for the HW model with
randomized volatility parameter, $\eta$).

Let us assume a sequence of positive constants $\eta_{n}$, $n=1,\dots,N.$
Then, the SDE

${\rm d}\overline{r}(t)=\overline{\lambda}(t,\overline{r}(t)){\rm
d}t+\overline{\eta}(t,\overline{r}(t)){\rm
d}W(t),\;\;\;\overline{r}(t_{0})=f(0,0),$ (2.19)

with

$\displaystyle\overline{\lambda}(t,y)=\sum_{n=1}^{N}\overline{\Lambda}_{n}(t,y)\lambda(\overline{\psi}_{n}(t)-y),\;\;\;\overline{\eta}^{2}(t,y)=\sum_{n=1}^{N}\eta_{n}^{2}\overline{\Lambda}_{n}(t,y),$
(2.20)

where:

$\overline{\Lambda}_{n}(t,y)=\frac{\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)}{\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)},$

has a strong solution whose marginal density is given by the mixture of normal
probability density functions:

$\displaystyle
f_{{\overline{r}(t)}}(y)=\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y),$
(2.21)

where $\sum_{n=1}^{N}\omega_{n}=1$ for $\omega_{n}\geq 0$, $n=1,\dots,N$ with
$f_{\overline{r}(t;\eta_{n})}(x)$ the PDF of the HW model with dynamics, given
by:

${\rm
d}\overline{r}_{n}(t)=\lambda(\overline{\psi}_{n}(t)-\overline{r}_{n}(t)){\rm
d}t+\eta_{n}{\rm d}W(t),\;\;\;\overline{r}_{n}(t_{0})=f(0,0),$ (2.22)

where $\overline{r}_{n}(t):=\overline{r}_{n}(t;\eta_{n})$ with
$\overline{\psi}_{n}(t)=f(0,t)+\frac{1}{\lambda}f(0,t)+\frac{\eta_{n}^{2}}{2\lambda^{2}}\left(1-{\rm
e}^{-2\lambda t}\right).$

###### Proof.

The problem we address is the derivation of the drift,
$\overline{\lambda}(t,y)$, and the volatility function,
$\overline{\eta}(t,\overline{r}(t))$, in the SDE:

${\rm d}\overline{r}(t)=\overline{\lambda}(t,\overline{r}(t)){\rm
d}t+\overline{\eta}(t,\overline{r}(t)){\rm d}W(t),$ (2.23)

such that

$f_{{\overline{r}(t)}}(y)=\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y),$
(2.24)

where $f_{\overline{r}(t;\eta)}(y)$ is the PDF of the Hull-White process in
(2.17), with parameter $\eta$. Under the HW model dynamics, the so-called
linear-growth condition holds, i.e., $\eta_{n}^{2}\leq C_{n}(1+y^{2})$ for
$i=1,\dots,N.$ Assuming that the same non-explosion condition holds for
(2.23), i.e., $\overline{\eta}^{2}(t,y)\leq C(1+y^{2}),$ uniformly in $t,$ we
start by deriving the Fokker-Planck equation for $\overline{r}(t)$:

$\displaystyle\frac{\partial}{\partial
t}f_{\overline{r}(t)}(y)=\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}\left(\overline{\eta}^{2}(t,y))f_{\overline{r}(t)}(y)\right)-\frac{\partial}{\partial
y}\overline{\lambda}(t,y)f_{\overline{r}(t)}(y),$ (2.25)

while for each individual interest rate process, $r_{n}(t)$, in (2.22) we
have:

$\displaystyle\frac{\partial}{\partial
t}f_{\overline{r}(t;\eta_{n})}(y)=\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}\left(\eta_{n}^{2}f_{\overline{r}(t;\eta_{n})}(y)\right)-\frac{\partial}{\partial
y}\left(\lambda(\overline{\psi}_{n}(t)-y)\right)f_{\overline{r}(t;\eta_{n})}(y).$
(2.26)

After substituting (2.24) into (2.25):

$\displaystyle\frac{\partial}{\partial
t}\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)=\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}\left[\overline{\eta}^{2}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)\right]-\frac{\partial}{\partial
y}\overline{\lambda}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y).$

Due to the linearity of the derivative operator, we find,

$\displaystyle\sum_{n=1}^{N}\omega_{n}\frac{\partial}{\partial
t}f_{\overline{r}(t;\eta_{n})}(y)=\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}\left[\overline{\eta}^{2}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)\right]-\frac{\partial}{\partial
y}\overline{\lambda}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y).$

By substitution of (2.26) and operator linearity,

$\displaystyle\small\frac{\partial^{2}}{\partial
y^{2}}\sum_{n=1}^{N}\omega_{n}\frac{1}{2}\eta_{n}^{2}f_{\overline{r}(t;\eta_{n})}(y)-\sum_{n=1}^{N}\omega_{n}\frac{\partial}{\partial
y}\lambda(\overline{\psi}_{n}(t)-y)f_{\overline{r}(t;\eta_{n})}(y)=$
$\displaystyle\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}\left[\overline{\eta}^{2}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)\right]-\frac{\partial}{\partial
y}\overline{\lambda}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y).$

Finally, by matching the appropriate terms and integration, we find:

$\displaystyle\sum_{n=1}^{N}\omega_{n}\eta_{n}^{2}f_{\overline{r}(t;\eta_{n})}(y)$
$\displaystyle=$
$\displaystyle\overline{\eta}^{2}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)+C_{1}(t)y+C_{2}(t),$
$\displaystyle\sum_{n=1}^{N}\omega_{n}\lambda(\overline{\psi}_{n}(t)-y)f_{\overline{r}(t;\eta_{n})}(y)$
$\displaystyle=$
$\displaystyle\overline{\lambda}(t,y)\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)+C_{3}(t),$

for some time-dependent functions $C_{1}(t)$, $C_{2}(t)$ and $C_{3}(t).$ Since
the LHS and RHS of the equations above need to satisfy uniform convergence
requirements implying that for $y\rightarrow\infty$, they converge to $0$, the
following needs to hold: $C_{1}(t)=C_{2}(t)=C_{3}(t)=0$, $\forall t.$
Therefore, the expressions for $\overline{\lambda}(t,y)$ and
$\overline{\eta}(t,y)$ read:

$\displaystyle\overline{\eta}^{2}(t,y)=\frac{\sum_{n=1}^{N}\omega_{n}\eta_{n}^{2}f_{\overline{r}(t;\eta_{n})}(y)}{\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)},\;\;\;\;\;\overline{\lambda}(t,y)=\frac{\sum_{n=1}^{N}\omega_{n}\overline{\lambda}(\overline{\psi}_{n}(t)-y)f_{\overline{r}(t;\eta_{n})}(y)}{\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)}.$
(2.27)

By setting

$\overline{\Lambda}_{n}(t,y)=\frac{\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)}{\sum_{n=1}^{N}\omega_{n}f_{\overline{r}(t;\eta_{n})}(y)}\;\;\;\text{for}\;\;\;n=1,\dots,N,$
(2.28)

we can write

$\displaystyle\overline{\lambda}(t,y)=\sum_{n=1}^{N}\overline{\Lambda}_{n}(t,y)\lambda(\overline{\psi}_{n}(t)-y),\;\;\;\text{and}\;\;\;\overline{\eta}^{2}(t,y)=\sum_{n=1}^{N}\eta_{n}^{2}\overline{\Lambda}_{n}(t,y).$
(2.29)

Finally, by taking $\eta_{*}:=\max_{i=1,\dots,N}\eta_{n}$ and since
$\sum_{n=1}^{N}\overline{\Lambda}_{n}(t,y)=1,$ $\forall y,$ we have:

$\displaystyle\overline{\eta}^{2}(t,y)=\sum_{n=1}^{N}\eta_{n}^{2}\overline{\Lambda}_{n}(t,y)\leq\sum_{n=1}^{N}\eta_{*}^{2}\overline{\Lambda}_{n}(t,y)=\eta_{*}^{2}=C.$

Since the volatility parameter $\overline{\eta}^{2}(t,y)$ is bounded by a
constant, the uniform convergence criterion is satisfied. The uniqueness of
the strong solution follows from Theorem 12.1 in [17] while in [2] (Theorem
2.1), a proof for a generic case for a normal mixture is provided. ∎

Proposition 2.1 illustrates that under the randomized volatility parameter for
the HW model, the normal mixture dynamics resemble a one-dimensional local-
volatility-type diffusion process. The local volatility function is expressed
as a weighted volatility squared, $\eta_{n}^{2}$, of the constituent
processes, $\overline{r}_{n}(t)$, and where the weights are functions of the
quadrature coefficients, $\omega_{n}$, and the corresponding PDFs,
$f_{\overline{r}(t;\eta_{n})}(y)$.

Following the same strategy, we derive the dynamics of a process of the rHW
model with the randomized mean-reversion parameter, $\lambda$. Proposition 2.2
provides the details.

######  Proposition 2.2 (Dynamics of the HW model with randomized mean-
reversion parameter, $\lambda$).

Let us assume a sequence of positive constants $\lambda_{m}$, $m=1,\dots,M$,
then the SDE

${\rm d}\widetilde{r}(t)=\widetilde{\lambda}(t,\widetilde{r}(t)){\rm
d}t+\widetilde{\eta}(t,\widetilde{r}(t)){\rm
d}W(t),\;\;\;\widetilde{r}(t_{0})=f(0,0),$ (2.30)

with $\widetilde{\eta}(t,y)=\eta$, and

$\displaystyle\widetilde{\lambda}(t,y)=\sum_{m=1}^{M}\widetilde{\Lambda}_{m}(t,y)\lambda_{m}(\psi_{m}(t)-y),\;\;\;\text{where}\;\;\;\widetilde{\Lambda}_{m}(t,y)=\frac{\varpi_{m}f_{\widetilde{r}(t;\lambda_{m})}(y)}{\sum_{m=1}^{M}\varpi_{m}f_{\widetilde{r}(t;\lambda_{m})}(y)},$
(2.31)

has a strong solution whose marginal density is given by the mixture of normal
probability density functions:

$\displaystyle
f_{{\widetilde{r}(t)}}(x)=\sum_{m=1}^{M}\varpi_{m}f_{\widetilde{r}(t;\lambda_{m})}(x),$
(2.32)

where $\sum_{m=1}^{M}\varpi_{m}=1$ for $\varpi_{m}\geq 0$, $m=1,\dots,M$, with
$f_{\widetilde{r}(t;\lambda_{m})}(x)$ is the PDF the HW model whose dynamics
is given by:

${\rm d}\widetilde{r}_{m}(t)=\lambda_{m}(\psi_{m}(t)-\widetilde{r}_{m}(t)){\rm
d}t+\eta{\rm d}W(t),\;\;\;\widetilde{r}_{m}(t_{0})=f(0,0),$

where $\widetilde{r}_{n}(t):=\widetilde{r}_{n}(t;\lambda_{n})$ with
$\psi_{m}(t)=f(0,t)+\frac{1}{\lambda_{m}}f(0,t)+\frac{\eta^{2}}{2\lambda_{m}^{2}}\left(1-{\rm
e}^{-2\lambda_{m}t}\right).$

###### Proof.

The proof is analogous to Proposition 2.1. ∎

With the rHW model, the randomization of either model parameters enables us to
determine the corresponding 1D SDE. For the random volatility parameter,
$\eta$, we derived a local-volatility type short-rate process; however, in the
case of random $\lambda$, the corresponding 1D process has a different
structure, i.e., the volatility coefficient, $\widetilde{\eta}(t,y)=\eta$, in
(2.30), stays constant. When analyzing the dynamics of the corresponding ZCB
(see Corollary 2.1), which is a tradable asset, the dynamics will resemble a
local-volatility type process. To derive the corresponding SDE, one needs to
follow the same strategy as in Propositions 2.1 and 2.2 with the process for
the ZCB, instead of the short-rate process.

In the final part of this section, the case of joint randomization of both HW
model parameters will be discussed.

####  2.3.1 Dynamics of the rHW model with bivariate distribution for
$\lambda$ and $\eta$.

An extension of the RAnD method is to consider both HW model parameters random
and follow a bivariate distribution. Such an extension may benefit from the
interconnection between model parameters and possibly the correlation
coefficient; however, it would require that the conditional moments are known
explicitly. This is troublesome because only for a few random distributions
the moment functions are known in the closed form. However, if we stay, for
example, within the Gaussian world, such an extension to the 2D case is
possible.

######  Corollary 2.2 (Random parameters with bivariate distribution).

Under a bivariate distribution $\Theta=[\vartheta_{1},\vartheta_{2}]$ with
$\zeta(\vartheta_{1})=\\{\omega_{1,n},\theta_{1,n}\\}_{n=1}^{N}$ and
conditioned on
$\zeta(\vartheta_{2}|\vartheta_{1})=\\{\omega_{2,m},\theta_{2,m}\\}_{m=1}^{M}$,
the randomized ChF is given by:

$\displaystyle\phi_{{X}}({u};t,T)=\sum_{n=1}^{N}\omega_{n}\sum_{m=1}^{M}\omega_{m}\phi_{{X}|\vartheta_{1}=\theta_{n},\vartheta_{2}=\theta_{m}}({u};t,T)+\epsilon_{N,M}^{b},$
(2.33)

where $N$ and $M$ indicate the number of expansion terms for $\vartheta_{1}$
and $\vartheta_{2}|\vartheta_{1}$ respectively, $\vartheta_{2}|\vartheta_{1}$
indicates a conditional random variable, $\epsilon_{N,M}^{b}$ is the
corresponding aggregated error, and the remaining specification follows
Theorem 2.1.

The PDF computation for the bivariate case requires the sequential computation
of the associated weights and the corresponding points, which can be
established utilizing Propositions 2.1 and 2.2. In the first iteration step,
we compute the grid points associated with the volatility parameter, $\eta$,
and for each realization $\eta_{n}$, we establish the corresponding
conditional PDF. By summing over all possible pairs, the unconditional PDF for
the rHW model can be computed,

$\displaystyle
f_{r(T;\eta_{n})}(x)=\sum_{m=1}^{M}\varpi_{m,n}f_{{r(T;\eta_{n},\lambda_{m})}}(x),\;\;\;f_{{r(T)}}(x)=\sum_{n=1}^{N}\sum_{m=1}^{M}\omega_{n}\varpi_{m,n}f_{{r(T;\eta_{n},\lambda_{m})}}(x),$
(2.34)

where $\varpi_{m,n}$ is defined in (2.32), $\omega_{n}$ is given in (2.21) and
$f_{{r(T;\eta_{n},\lambda_{m})}}(x)$ indicates the PDF of the HW model with
parameters $\eta_{n}$ and $\lambda_{m}$, and is defined in (2.14). Proposition
2.3 provides the dynamics of the associated rHW model.

######  Proposition 2.3 (Local volatility process for the HW model with
randomized parameters).

Let us assume a sequence of positive constants $\eta_{n}$, $n=1,\dots,N$, and
$\lambda_{m,n}$, $m=1,\dots,M$, then the SDE:

${\rm d}\widehat{r}(t)=\widehat{\lambda}(t,\widehat{r}(t)){\rm
d}t+\widehat{\eta}(t,\widehat{r}(t)){\rm
d}W(t),\;\;\;\widehat{r}(t_{0})=f(0,0),$ (2.35)

with

$\displaystyle\widehat{\eta}^{2}(t,\widehat{r}(t))$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}\widehat{\Lambda}_{n,m}(t,y)\eta_{n}^{2},\;\;\;\widehat{\lambda}(t,\widehat{r}(t))=\sum_{n=1}^{N}\sum_{m=1}^{M}\widehat{\Lambda}_{n,m}(t,y)\lambda_{m,n}(\psi_{m,n}(t)-y),$

where

$\widehat{\Lambda}_{n,m}(t,y)=\frac{\omega_{n}\varpi_{m,n}f_{{r(t;\eta_{n},\lambda_{m})}}(y)}{\sum_{n=1}^{N}\sum_{m=1}^{M}\omega_{n}\varpi_{m,n}f_{r(t;\eta_{n},\lambda_{m,n})}(y)},$

has a strong solution whose marginal density is given by the following mixture
of normal probability density functions:

$\displaystyle
f_{{\widehat{r}(t)}}(y)=\sum_{n=1}^{N}\sum_{m=1}^{M}\omega_{n}\varpi_{m,n}f_{{r(t;\eta_{n},\lambda_{m,n})}}(y),$

where $\sum_{n=1}^{N}\omega_{n}=1$, $\sum_{m=1}^{M}\varpi_{m,n}=1$ for
$\omega_{n},\varpi_{m,n}\geq 0$, $n=1,\dots,N$, $m=1,\dots,M$ with
$f_{r(t;\eta_{n},\lambda_{m,n})}(y)$ being the PDF of the HW model whose
dynamics are given by:

${\rm d}r_{m,n}(t)=\lambda_{m,n}(\psi_{m,n}(t)-r_{m,n}(t)){\rm
d}t+\eta_{n}{\rm d}W(t),\;\;\;r_{m,n}(t_{0})=f(0,0),$

and $r_{n,m}(t):=r(t;\eta_{n},\lambda_{m})$ with
$\psi_{m,n}(t)=f(0,t)+\frac{1}{\lambda_{m,n}}f(0,t)+\frac{\eta_{n}^{2}}{2\lambda_{m,n}^{2}}\left(1-{\rm
e}^{-2\lambda_{m,n}t}\right).$

###### Proof.

The proof is analogous to Proposition 2.1. ∎

The concept of multivariate model parameters can be extended further. In the
case of the HW model with piece-wise constant parameters, it is possible to
randomize each piece-wise element, to mimic the parameter controlled by the
stochastic process. Unfortunately, as the number of ranges increases, the
number of terms in the summation in (2.34) will grow exponentially.

The impact of bivariate model parameters on implied volatilities will be
analyzed further in Section 4.

##  3 Pricing under the randomized Hull-White (rHW) model

This section focuses on the pricing under the rHW model. The presented pricing
equations utilize Theorem 2.1. Throughout the section, we denote by $V_{{\rm
HW}}(t,r(t;\theta))$ with $\theta\in\\{\eta,\lambda\\}$, the Hull-White price
of a derivative depending on the short-rate $r(t;\theta)$ with parameters
$\lambda$ and $\eta$, by $V_{{\rm rHW}}(t,r(t;\vartheta))$ we denote the
pricing under the rHW model.

###  3.1 Swaptions under rHW model

While in Equation (2.5), there is a relationship between HW and rHW prices,
here explicit pricing formulas will be given. The valuation of swaption
contracts in the rHW model can be carried out at the same complexity level as
in the standard HW model.

The critical ingredient for pricing swaptions under the rHW model is an
analytical expression for the options of the pure ZCB. Lemma 3.4 provides the
pricing equations under the rHW model.

######  Lemma 3.4 (Pricing of options on a ZCB under the rHW model).

Under the randomized HW model parameters, $\theta:=\eta$ or $\theta:=\lambda$,
with $\vartheta(\omega)=\theta$, the price of a European-style option on a
ZCB, $P(T,S)$, with option’s expiry, $T$, strike, $K$, and maturity of the
underlying bond $S$, is given by:

$\displaystyle
V^{\text{Z}}_{\chi}(t,T,S,K;\vartheta)=\sum_{n=1}^{N}\omega_{n}V^{\text{Z}}_{\chi}(t,T,S,K;\theta_{n}),$
(3.1)

where,

$\displaystyle V^{\text{Z}}_{\chi}(t,T,S,K;\theta_{n})$ $\displaystyle=$
$\displaystyle\mathbb{E}^{\mathbb{Q}}_{t}\left[\frac{M(t)}{M(T)}\max(\chi(P(T,S;\theta_{n})-K),0)\right]$
(3.2) $\displaystyle=$ $\displaystyle\chi P(t,S)F_{\mathcal{N}(0,1)}(\chi
d_{n})-\chi KP(t,T)F_{\mathcal{N}(0,1)}(\chi(d_{n}-\bar{\sigma}_{n})),$

with $P(t,S)$ and $P(t,T)$ the ZCB computed from the associated yield curve;
$P(T,S;\theta_{n})$ is the ZCB obtained from the model with the randomized
parameter, $\vartheta$; $\chi=1$ and $\chi=-1$ corresponds to call and put
options, respectively. Depending on the randomized parameter, we have:

$\displaystyle\text{for}\;\;\eta:\;\bar{\sigma}^{2}_{n}=\frac{\eta_{n}^{2}}{2\lambda}\left(1-{\rm
e}^{-2\lambda(T-t)}\right)B^{2}(T,S;\lambda),$ (3.3)
$\displaystyle\text{for}\;\;\lambda:\;\bar{\sigma}^{2}_{n}=\frac{\eta^{2}}{2\lambda_{n}}\left(1-{\rm
e}^{-2\lambda_{n}(T-t)}\right)B^{2}(T,S;\lambda_{n}),$ (3.4)

with

$d_{n}=\frac{1}{\bar{\sigma}_{n}}\log\frac{P(t,S)}{P(t,T)K}+\frac{\bar{\sigma}_{n}}{2},\;\;\;B(T,S;\lambda)=\frac{1}{\lambda}\left(1-{\rm
e}^{-\lambda(S-T)}\right),$

where the HW model is defined in (2.12) and $F_{\mathcal{N}(0,1)}(\cdot)$
corresponds the standard normal CDF.

###### Proof.

The proof is a direct consequence of combining Theorem 2.1 with the pricing of
options on a ZCB under the HW model (as given in [4, 16]). ∎

Now, we consider an interest rate swap with a fixed rate, $K$, and payment
times $\mathcal{T}=\\{T_{i},T_{i+1},\dots,T_{m}\\}$ and the corresponding
reset rates $\\{T_{i-1},T_{i+1},\dots,T_{m-1}\\}$ with the payoff given by:

$H^{Swap}_{\text{P/R}}(\mathcal{T},K)=\bar{\alpha}\sum_{k=i}^{m}\tau_{k}\left(\ell(T_{k-1};T_{k-1},T_{k})-K\right),$

for $\tau_{k}=T_{k}-T_{k-1},$ with $P$ indicating a swap payer for
$\bar{\alpha}=1$ and a swap receiver for $\bar{\alpha}=-1$ and where
$\ell(t;T_{k-1},T_{k})$ stands for the libor rate over the period
$[T_{k-1},T_{k}]$ observed at time $t.$

To determine today’s value of the swap, we evaluate the corresponding
expectation of the discounted future cash flows, i.e., each payment which
takes place at the time points, $T_{i},\dots,T_{m}$, needs to be discounted to
today,

$\displaystyle
V_{\text{P/R}}^{\text{Swap}}(t,K,\mathcal{T})=\bar{\alpha}\sum_{k=i}^{m}\tau_{k}\mathbb{E}^{\mathbb{Q}}_{t}\left[\frac{M(t)}{M(T_{k})}\big{(}\ell_{k}(T_{k-1})-K\big{)}\right]=\bar{\alpha}\sum_{k=i}^{m}\tau_{k}P(t,T_{k})\big{(}\ell_{k}(t)-K\big{)},$
(3.5)

with $\ell_{k}(t):=\ell(t;T_{k-1},T_{k}).$

We derive the valuation formula for a swaption contract with strike $K$ and
option expiry $T=T_{i-1}$. Using (3.5), the pricing equation reads:

$\displaystyle
V^{\text{Swpt}}_{\text{P/R}}(t,T,K,\mathcal{T})=\mathbb{E}^{\mathbb{Q}}_{t}\left[\frac{M(t)}{M(T)}\max\left(V_{\text{P/R}}^{\text{Swap}}(T,K,\mathcal{T}),0\right)\right]=P(t,T)\mathbb{E}^{T}_{t}\left[\max\left(V_{\text{P/R}}^{\text{Swap}}(T,K,\mathcal{T}),0\right)\right],$

where the expectation, $\mathbb{E}^{T}_{t}[\cdot]$, is taken under the
$T-$forward measure and where the price of a swap price at $T=T_{i-1}$ is
given by:

$\displaystyle
V_{\text{P/R}}^{\text{Swap}}(T,K,\mathcal{T})=\bar{\alpha}\left[1-P(T,T_{m})-K\sum_{k=i}^{m}\tau_{k}P(T_{i},T_{k})\right]=\bar{\alpha}-\bar{\alpha}\sum_{k=i}^{m}c_{k}P(T_{i},T_{k}),$

with $c_{k}=K\tau_{k}$ for $k=i,\dots,m-1$, $c_{m}=1+K\tau_{m}$, and
$\tau_{k}=T_{k}-T_{k-1}$.

So, the pricing equation becomes:

$\displaystyle
V^{\text{Swpt}}_{\text{P/R}}(t,T,K,\mathcal{T})=\bar{\alpha}P(t,T)\mathbb{E}^{T}_{t}\left[\max\left(1-\sum_{k=i}^{m}c_{k}P(T_{i},T_{k}),0\right)\right].$

Up to this point, the pricing equations do not depend on specific model
choices. Now, however, we consider pricing under the rHW model, which, via the
conditional expectation approach, as presented in Lemma 2.1, yields:

$\displaystyle V^{\text{Swpt}}_{\text{P/R}}(t,T,K,\mathcal{T};\vartheta)$
$\displaystyle=$
$\displaystyle\bar{\alpha}P(t,T)\mathbb{E}_{t}^{T}\left[\mathbb{E}^{T}_{t}\left[\max\Big{(}1-\sum_{k=i}^{m}c_{k}P(T_{i},T_{k}),0\Big{)}\Big{|}\vartheta=\theta\right]\right]$
$\displaystyle=$
$\displaystyle\bar{\alpha}P(t,T)\sum_{n=1}^{N}\omega_{n}\mathbb{E}^{T}_{t}\left[\max\Big{(}1-\sum_{k=i}^{m}c_{k}P(T_{i},T_{k};\theta_{n}),0\Big{)}\right]+\epsilon_{N}.$

Here, $\epsilon_{N}$ is the associated quadrature error defined in (2.6).
Since the inner expectation in (3.1) resembles the expression for a swaption
under the classical HW, we follow the standard procedure and apply the so-
called Jamshidian trick [15] to exchange the maximum operator and the
expectation. The resulting pricing equations for swaptions under the rHW model
are given in Lemma 3.5.

######  Lemma 3.5 (Pricing of Swaptions under randomized Hull-White model).

Consider the rHW model, with parameters $\\{\lambda,\eta\\}$ and the
randomizing random variable $\vartheta$, which randomizes either of the model
parameters. For a unit notional, a constant strike, $K$, option expiry
$T=T_{i-1}$ and a strip of swap payments
$\mathcal{T}=\\{T_{i},\dots,T_{m}\\}$, with $T_{i}>T_{i-1}$ and accruals
$\tau_{i}=T_{i}-T_{i-1}$, the prices of swaption payer and receiver,
$\text{P/R}:=\text{Payer/Receiver}$, are given by:

$\displaystyle
V_{\text{P/R}}^{\text{Swpt}}(t_{0},T,\mathcal{T},K;\vartheta)=\sum_{n=1}^{N}\omega_{n}\sum_{k=i}^{m}c_{k}V^{\text{Z}}_{\chi}(t_{0},T,T_{k},\hat{K}_{k}(\theta_{n});\theta_{n}),$
(3.7)

with a swaption payer, P, for $\chi=-1$, swaption receiver, R, with $\chi=1$,
where $V^{\text{Z}}_{\chi}(\cdot)$ is defined in (3.2) and where the strike
price
$\hat{K}_{k}(\theta_{n})=\exp\left(A(T,T_{k};\theta_{n})+B(T,T_{k};\theta_{n})r^{*}_{n}\right)$.
Here, $r^{*}_{n}$ is determined by solving, for each parameter realization
$\theta_{n}$, the following equation:

$1-\sum_{k=i}^{m}c_{k}\exp\Big{(}A(T,T_{k};\theta_{n})-B(T,T_{k};\theta_{n})r^{*}_{n}\Big{)}=0,\;\;\;n=1,\dots,N,$
(3.8)

where

$\displaystyle A(T,T_{k};\\{\lambda,\eta\\})$ $\displaystyle=$
$\displaystyle\log\frac{P(0,T_{k})}{P(0,T)}+B(T,T_{k};\\{\lambda,\eta\\})f(0,T)-\frac{\eta^{2}}{4\lambda}\left(1-{\rm
e}^{-2\lambda T}B^{2}(T,T_{k};\\{\lambda,\eta\\})\right),$ $\displaystyle
B(T,T_{k};\\{\lambda,\eta\\})$ $\displaystyle=$
$\displaystyle\frac{1}{\lambda}\left(1-{\rm e}^{-\lambda(T_{k}-T)}\right),$

with $c_{k}=K\tau_{k}$ for $k=i,\dots,m-1$, $c_{m}=1+K\tau_{m}$ and where the
pairs $\\{\omega_{n},\theta_{n}\\}$, $n=1,\dots,N$, are based on the
randomizer $\vartheta$ and computed using A.2.

Under the rHW model, the pricing Equation (3.7) for swaptions is a direct
application of Theorem 2.1. The pricing under randomized parameters is simply
an average of non-randomized prices, with varying model parameters,
accompanied by weights based on the randomizer.

From the computational perspective, the swaption pricing in (3.7) requires $N$
swaption prices for different realizations $\theta_{n}$, $n=1,\dots,N$;
therefore, $N$ optimization problems in (3.8) need to be solved. This
computational complexity can be greatly reduced by employing the multi-d
Newton-Raphson algorithm to determine optimal $r^{*}_{n}$ for $n=1,\dots,N.$

######  Remark 3.6 (Computation of sensitivities).

The computation of sensitivities under the RAnD method is straightforward,
i.e., because the pricing is expressed as a convex combination, the
sensitivities are expressed as a weighted sum of individual derivatives. For
example, in the case of the swaption pricing, the sensitivity to a particular
market quote $q$ is expressed as a sum of sensitivities of individual options
on ZCBs:

$\displaystyle\frac{\partial}{\partial
q}V_{\text{P/R}}^{\text{Swpt}}(t_{0},T,\mathcal{T},K;\vartheta)=\sum_{n=1}^{N}\omega_{n}\sum_{k=i}^{M}c_{k}\frac{\partial}{\partial
q}V^{\text{Z}}_{\chi}(t_{0},T,T_{k},\hat{K}_{k}(\theta_{n});\theta_{n}),$
(3.9)

with the specification as given in Lemma 3.5.

On the other hand, the sensitivity to the quadrature pairs,
$\\{\omega_{n},\theta_{n}\\}$, $\partial\theta_{n}/\partial\hat{a}$ and
$\partial\omega_{n}/\partial\hat{a}$ may be, for some specific cases, computed
analytically (see Section 4.3). Still, in a generic setting, it is recommended
to compute these derivatives numerically, with, for example, finite
differences:

$\frac{\zeta(\vartheta(\hat{a}+\delta_{\hat{a}}))-\zeta(\vartheta(\hat{a}-\delta_{\hat{a}}))}{2\delta_{\hat{a}}}\approx\left\\{\frac{\partial\omega_{n}}{\partial\hat{a}},\frac{\partial\theta_{n}}{\partial\hat{a}}\right\\},$

where $\vartheta(\hat{a})$ indicates the dependence of the random variable and
parameter $\hat{a}$, $\delta_{\hat{a}}$ is the shock size and
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$
is defined in A.2. Due to the applied finite difference shocks to $\hat{a}$,
an additional bias will be introduced. We expect, however, this error to be of
acceptable magnitude, as is commonly observed in current financial practice.

##  4 Numerical experiments

In this section, several pricing experiments will be performed. First, we
present a detailed analysis of the rHW model in realistic pricing scenarios.
In the first experiment, we explore the implied volatility smile evolution in
time, comparing the implied volatility surface from the HW and rHW models. As
the next step, the study of the parameter randomization on shapes of implied
volatilities will be illustrated. Finally, the calibration results with market
data will be presented in the conclusive experiment. This section will end
with numerical experiments involving bivariate distribution for the model
parameters.

###  4.1 Evolution of implied volatilities

We analyze the swaption implied volatilities observed in the market and
compare them to implied volatilities generated using the HW and rHW models. We
focus here on the flexibility in generating realistic implied volatility
shapes. Throughout the section, we will consider the shifted implied
volatilities computed by inverting Black’s formula, which, for
$\bar{\alpha}\in\\{1,-1\\}$ for call and put options, respectively, reads:

$\displaystyle V_{\text{B}}(T,K,F_{0},\sigma_{s},\bar{\alpha},s)$
$\displaystyle=$ $\displaystyle\bar{\alpha}\cdot(F_{0}+s)\cdot
F_{\mathcal{N}(0,1)}(\bar{\alpha}d_{1})-\bar{\alpha}(K+s)F_{\mathcal{N}(0,1)}(\bar{\alpha}d_{2}),$
(4.1) $\displaystyle d_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma\sqrt{T}}\left[\log(F_{0}+s)/(K+s)+1/2\sigma_{s}^{2}T\right],$
$\displaystyle d_{2}$ $\displaystyle=$ $\displaystyle
d_{1}-\sigma_{s}\sqrt{T},$

with shift parameter $s$, $F_{0}$ being the forward rate, $\sigma_{s}$ is the
corresponding volatility coefficient, $K$ being the strike, and $T$
corresponds to the time to option expiry.

Then, the shifted Black’s formula for swaptions is given by:

$\displaystyle
V^{\text{B,Swpt}}_{P/R}=V_{\text{B}}(T,K,S(t_{0}),\sigma_{s},\bar{\alpha},s)\sum_{i=1}^{m}\tau_{i}P(t_{0},T_{i}),\;\;\;S(t_{0})=\frac{P(t_{0},T_{i-1})-P(t_{0},T_{m})}{\sum_{i=1}^{m}\tau_{i}P(t_{0},T_{i})},$
(4.2)

with swaption payer, $P$, for $\bar{\alpha}=1$ and swaption receiver, $R$, for
$\bar{\alpha}=-1$ and where $S(t_{0})$ is the corresponding swap rate. The
shift parameter, $s$, is typically perceived as the lower bound for the
interest rates by the market participants and varies depending on the
currency. It is important to note that when inverting the Black’s formula in
(4.2), the corresponding implied volatility $\sigma_{s}$ is a function of the
shift parameter, i.e., for different choices of $s$, the implied volatilities
are different. In principle, this is not a problem as long as the implied
volatilities from the market and model volatilities are computed with the same
shift coefficient.

In the first experiment, we compare the implied volatility surfaces of the HW
model and the rHW model. In the rHW model, the speed of mean-reversion is
randomized by a uniform distribution,
$\lambda\sim\mathcal{U}([\hat{a},\hat{b}])$, on an interval
$[\hat{a},\hat{b}].$ The numerical results are illustrated in Figure 1. The
results demonstrate that the HW model can only generate implied volatility
skew. At the same time, the randomization of the mean-reversion parameter,
$\lambda$, shows implied volatility skew and smile. Although the randomization
is not time-dependent, i.e., the parameters are stochastic but stationary, we
observe a time evolution of the implied volatilities. The same phenomenon has
been observed in [3, 11], where the randomized Black-Scholes model was
considered.

Figure 1: Swaption volatility evolution for the HW and rHW models implied by
the shifted Black’s model. The simulation was performed for varying swaption
option expiry, $T$, and a fixed tenor of $1y$. The parameters specified in the
experiment are: for the HW model: $\eta=0.005$, $\lambda=0.001$ and for the
rHW model: $\eta=0.005$ and $\lambda~{}\sim\mathcal{U}([-0.15,0.6]).$ In the
experiment, the implied volatilities are computed with zero shift parameter,
$s=0.$

The reported numerical results are promising. By a stochastic mean-reversion
parameter, the HW model has an additional degree of freedom that improves the
model’s flexibility. In the following sections, we will analyze the
randomization effect of all model parameters and the different choices of
randomizing distributions.

###  4.2 Randomization and impact on implied volatilities

Here, we study the impact of different randomization choices on the implied
volatility shape. Understanding how parameters affect the generated
volatilities is crucial in model calibration. A clear relationship between the
model parameters and the associated skew, curvature and volatility level is
fundamental, enabling traders to react swiftly to market movements. It is also
the basis for setting up a calibration routine when a particular market
movement can be encapsulated in a parameter change.

In the first experiment, we analyze the randomization of volatility parameter
$\eta$. We consider two randomization cases where $\eta$ either follows a
uniform or a normal distribution. The details regarding the computation of the
associated moments and the corresponding quadrature points can be found in A.
The experiment is set up so that both approaches’ mean values,
$\mathbb{E}[\eta]$, are equal.

Figure 2: Impact of randomized volatility parameter $\eta$ on implied
volatilities with a fixed speed of mean-reversion parameter $\lambda=0.009.$
Left: $\eta\sim\mathcal{U}([\hat{a},\hat{b}])$. Right:
$\eta\sim\mathcal{N}(\hat{a},\hat{b}).$

The results are presented in Figure 2, where it is shown that the
randomization of the volatility parameter $\eta$ has a pronounced effect on
the level of implied volatility with minimal effect on the smile. The
curvature is visible for uniform randomization. It is essential to note that a
higher variance for the normal distribution may give rise to a higher
curvature but may also cause issues related to negative volatilities. Our
experiments have shown that even for distributions defined in the positive
domain, the impact on the smile is limited, even for fat-tailed random
variables.

A much richer spectrum of implied volatility shapes is obtained when the
randomization technique is applied to the mean-reversion parameter $\lambda$.
Figure 3 presents the randomization with either uniform or normal random
variables. A substantial amount of curvature can be generated by taking the
mean-reversion random. We also report that the curvature change affects the
overall volatility level, i.e., it is impossible to keep the level fixed and
only adjust the smile. However, the implied volatility level can be fixed by
adjusting the volatility parameter $\eta$. This strategy will be discussed
further in the context of model calibration.

Figure 3: Impact of randomized volatility parameter $\lambda$ on implied
volatilities with a fixed volatility parameter $\eta=0.0050.$ Left:
$\lambda\sim\mathcal{U}([\hat{a},\hat{b}])$. Right:
$\lambda\sim\mathcal{N}(\hat{a},\hat{b}).$

In the final experiment of this section, we consider the randomization of
$\lambda\sim\mathcal{U}([\hat{a},\hat{b}])$ using uniform distribution and
check how the parameters $\hat{a}$ and $\hat{b}$ affect the implied
volatilities. In Figure 4, the results show an interesting pattern: the
curvature level is mainly driven by the distance $|\hat{b}-\hat{a}|$, i.e.,
the larger the distance, the more implied volatility smile is generated.
Changes of either of the parameters affects the volatility level; therefore,
some of the volatility effect, $\eta$, can be offset by the interval
$[\hat{a},\hat{b}].$

Figure 4: Impact of randomized volatility parameter
$\lambda\sim\mathcal{U}([\hat{a},\hat{b}])$ on implied volatilities with a
fixed volatility parameter $\eta=0.0050$. Left: varying parameter $\hat{a}$.
Right: varying parameter $\hat{b}.$

Given the numerical results presented above, we conclude that the most
considerable impact on implied volatilities comes from the mean-reversion,
$\lambda$, and not from the volatility parameter, $\eta.$ This is particularly
interesting when we confront these results with Proposition 2.1 and
Proposition 2.2, where the dynamics of the corresponding short-rate processes
were derived and where it was shown that randomization of the mean-reversion
did not lead to the local volatility type of dynamics. In contrast, the
randomization of $\eta$ does, which can be explained by the nature of interest
rate derivatives, where derivative prices are driven by the dynamics of the
ZCBs, but not directly by the short-rate process. However, since under the HW
model, the volatility of the ZCBs is given by both parameters (see [16]), the
randomization of either of them will imply a local volatility type of dynamics
for the ZCBs.

###  4.3 Calibration of the randomized Hull-White model

Every model that may be considered successful must also show its ability to
calibrate to the market quotes. Moreover, the calibration and pricing process
needs to be computationally efficient, especially when considering the pricing
of large portfolios, for example, in the context of xVA. Therefore, this
section focuses on calibrating the rHW model to swaption implied volatilities.

Calibration of model parameters requires multiple iterations over the
parameter space until specific optimization criteria are met. This implies
that the pricing needs to be repeated at every choice of parameter candidate.
By randomization of the model parameters, the number of free parameters will
increase. The additional degrees of freedom will depend on the type of
randomizer chosen. Let us consider, for example, parameter randomization with
$\vartheta(\hat{a},\hat{b})$, which is driven by two model parameters
$\hat{a}$ and $\hat{b}$. Two parameters of the randomizer $\vartheta$
correspond to one additional degree of freedom that can be used for model
calibration.

From the computational perspective, using the RAnD method requires that for
every iteration step in the calibration procedure, the mapping between the
randomizing variable $\vartheta$ and the corresponding pairs
$\\{\omega_{n},\theta_{n}\\}$, $n=1,\dots,N$, needs to be employed. Although
the algorithm given in A.2 is straightforward (it mainly depends on the
computation of eigenvalues), it can be further simplified. When randomizing
variables can be expressed as a linear combination of some base random
variable, for example, normal or uniform, the calculations of the weights and
the corresponding nodes can be significantly simplified. In the case of a
normal randomizer, $\vartheta\sim\mathcal{N}(\hat{a},\hat{b}^{2})$, we can
benefit from the linearity of the normal distribution. For the standard normal
and its associated points, we have
$\theta_{n}=\hat{a}+\hat{b}\cdot\theta_{\mathcal{N},n}$, where
$\theta_{\mathcal{N},n}$ are the nodes corresponding to standard normal,
$\mathcal{N}(0,1).$ This implies that we can simply tabulate the results for
the standard normal and scale the points accordingly. In the case of the
weights $\omega_{n}$, they stay invariant to a linear transformation [12]. A
similar property holds for $\vartheta\sim\mathcal{U}([\hat{a},\hat{b}])$,
where the nodes can be computed for $\vartheta\sim\mathcal{U}([0,1])$ and
scaled appropriately.

In the calibration experiment, we consider the market data for the USD market
as of the 18th of August, 2022. In all calibration exercises, we consider a
fixed tenor of $1y$ and analyze the accuracy for varying expires and strikes.
The calibrated model parameters for both HW and rHW models are tabulated in
Table 1. The results are intriguing, i.e., in the randomization for $\lambda$,
we were able to calibrate all swaptions while keeping the mean of the
randomizer fixed at $0.1$. This shows that having more degrees of freedom is
not necessary to improve the calibration results. As in the standard HW model,
we have only used two parameters.

Table 1: Calibration of the HW and rHW model: parameters determined in swaption calibration.  | Hull-White | RAnD Hull-White
---|---|---
$T$, expiry | $\eta$ | $\lambda$ | $\eta$ | $\lambda$
$1y$ | 0.0094 | 0.0090 | 0.0091 | $\lambda\sim\mathcal{N}(0.1,0.45^{2})$
$2y$ | 0.0082 | 0.0035 | 0.0080 | $\lambda\sim\mathcal{N}(0.1,0.33^{2})$
$5y$ | 0.0069 | 0.0020 | 0.0079 | $\lambda\sim\mathcal{N}(0.1,0.16^{2})$
$8y$ | 0.0067 | 0.0095 | 0.0080 | $\lambda\sim\mathcal{N}(0.1,0.12^{2})$
$10y$ | 0.0067 | 0.0090 | 0.0082 | $\lambda\sim\mathcal{N}(0.1,0.11^{2})$
$15y$ | 0.0064 | 0.0080 | 0.0085 | $\lambda\sim\mathcal{N}(0.1,0.09^{2})$
$20y$ | 0.0060 | 0.0080 | 0.0086 | $\lambda\sim\mathcal{N}(0.1,0.08^{2})$

The calibration fit is presented in Figures 5, 6 and 7. We report an excellent
calibration fit for all considered option expiries, varying from 1y to 20y. We
have used two parameters in all the calibration cases, just as in the HW
model. The results confirm that the RAnD method has great potential for
improving existing pricing methods, even with the same number of degrees of
freedom.

Figure 5: Calibration results of the HW and the rHW models. The market implied
volatilities for swaptions were obtained on 18/08/2022 for the USD market.
Option expiry: $T=1y$ and $T=2y$ and the implied volatility shift: $s=1\%$.
Calibrated parameters are presented in Table 1.

Figure 6: Calibration results of the HW and the rHW models. The market implied
volatilities for swaptions were obtained on 18/08/2022 for the USD market.
Option expiry: $T=5y$ and $T=10y$ and the implied volatility shift: $s=1\%$.
Calibrated parameters are presented in Table 1.

Figure 7: Calibration results of the HW and the rHW models. The market implied
volatilities for swaptions were obtained on 18/08/2022 for the USD market.
Option expiry: $T=15y$ and $T=20y$ and the implied volatility shift: $s=1\%$.
Calibrated parameters are presented in Table 1.

###  4.4 Pricing under bivariate distributions for the model parameters

This section extends the rHW model and considers a bivariate distribution for
both model parameters. Under a bivariate distribution
$\Theta=[\vartheta_{1},\vartheta_{2}]$ with
$\zeta(\vartheta_{1})=\\{\omega_{1,n},\theta_{1,n}\\}_{n=1}^{N}$ and
conditioned on
$\zeta(\vartheta_{2}|\vartheta_{1})=\\{\omega_{2,m},\theta_{2,m}\\}_{m=1}^{M}$,
the randomized prices are:

$\displaystyle V_{{\rm
rHW}}(t,r(t);\vartheta_{1},\vartheta_{2})=\sum_{n=1}^{N}\omega_{n}\sum_{m=1}^{M}\omega_{m}V_{{\rm
HW}}(t,r(t);\eta_{n},\lambda_{m})+\epsilon_{N,M}^{b},$ (4.3)

where $N$ and $M$ indicate the number of expansion terms for $\vartheta_{1}$
and $\vartheta_{2}|\vartheta_{1}$, respectively, $\vartheta_{2}|\vartheta_{1}$
indicates a conditional random variable, $\epsilon_{N,M}^{b}$ is the
corresponding aggregated error. The remaining specification follows Theorem
2.1.

As an example, let us consider a bivariate normal distribution for the pair
$(\vartheta_{1},\vartheta_{2})$ with the corresponding realizations
$(\eta,\lambda)$ for which we have:

$\displaystyle\vartheta_{2}|\vartheta_{1}=\eta_{n}\sim\mathcal{N}\left(\mu_{\lambda}+\frac{\sigma_{\lambda}}{\sigma_{\eta}}\rho(\eta_{n}-\mu_{\eta}),(1-\rho^{2})\sigma_{\lambda}^{2}\right),$
(4.4)

where $\rho$ is the correlation coefficient between $\vartheta_{1}$ and
$\vartheta_{2}$. In Figure 8, we illustrate the impact of the correlation
coefficient $\rho$ on the swaption implied volatilities. We report that
similar to, e.g., the Heston model, the correlation controls the implied
volatility skew. Furthermore, we observe that a higher positive correlation
generates more skew, while a negative correlation generates more volatility
curvature (smile).

Figure 8: Impact of randomized volatility parameters $\eta$ and $\lambda$ on
implied volatilities driven by bivariate normal distribution. Results are
presented for varying correlation, $\rho$. Left:
$\eta\sim\mathcal{N}(0.008,0.002^{2})$ and
$\lambda\sim\mathcal{N}(0.5,0.05^{2})$. Right:
$\eta\sim\mathcal{N}(0.01,0.002^{2})$ and
$\lambda\sim\mathcal{N}(0.5,0.2^{2})$.

######  Remark 4.7 (Volatility term structure and feasible strategy for model
calibration).

Commonly, the volatility-parameter $\eta$ is piece-wise constant, so the ATM
volatilities are well calibrated. This strategy will also work with the
randomized mean-reversion parameter, i.e., the mean-reversion parameter can be
used for smile/skew calibration, while the volatility parameter, $\eta$, will
ensure a proper fit of the ATM level.

###  4.5 Convergence results

As presented in Theorem 2.1, applying the RAnD method produces a quadrature
error. In this section, we analyze the convergence of the error depending on
the number of expansion terms, $N$.

The convergence speed depends on the distribution and its parameters, as
presented in Figure 9. We report an excellent convergence factor; however, as
expected, it depends on the variance of the randomizing random variable. Both
cases will achieve satisfactory results for already $N=5$. Randomization using
different random variables showed equivalent patterns of convergence.

Figure 9: Convergence results for randomization of volatility parameter, and
mean reversion parameter $\eta$ and $\lambda$ under the HW model respectively.
The base parameters in the experiment were $\eta=0.00625$ and $\lambda=0.002$.

##  5 Conclusion

In this paper, we have applied the randomization technique to enhance the
flexibility of short-rate models for interest rates. We have shown that the
model parameters driven by a random variable, instead of being deterministic,
facilitate a practical extension of standard, well-popularized models. In
addition, this article points out how the normal mixture (a sum of the Hull-
White PDFs) can be expressed as a one-dimensional diffusion process with a
local volatility function.

We have illustrated that, for the randomized Hull-White model (rHW), one can
utilize the available closed-form pricing equations and benefit from
flexibility in controlling implied volatilities. In particular, we have shown
that the rHW model results in almost perfect swaption calibration.

Finally, the RAnD method is generic and is not limited to any particular
modelling choice; therefore, it opens many possibilities for improving
existing pricing frameworks.

## References

  * [1] J. Andreasen.  Stochastic volatility interest rate models.  Encyclopedia of Quantitative Finance, 2010.
  * [2] D. Brigo.  The general mixture-diffusion SDE and its relationship with an uncertain-volatility option model with volatility-asset decorrelation.  arXiv preprint arXiv:0812.4052, 2008.
  * [3] D. Brigo and F. Mercurio.  Lognormal-mixture dynamics and calibration to market volatility smiles.  International Journal of Theoretical and Applied Finance, 5(04):427–446, 2002.
  * [4] D. Brigo and F. Mercurio.  Interest rate models- theory and practice: with smile, inflation and credit.  Springer Finance, 2nd edition, 2007.
  * [5] R. Bulirsch and J. Stoer.  Introduction to numerical analysis, volume 3.  Springer, 2002.
  * [6] J. Casassus, P. Collin-Dufresne, and B. Goldstein.  Unspanned stochastic volatility and fixed income derivatives pricing.  Journal of Banking & Finance, 29(11):2723–2749, 2005.
  * [7] D. Duffie, J. Pan, and K. Singleton.  Transform analysis and asset pricing for affine jump-diffusions.  Econometrica, 68:1343–1376, 2000.
  * [8] H. G. Fong and O. A. Vašiček.  Fixed-income volatility management.  Journal of Portfolio Management, 17(4):41, 1991.
  * [9] D. Ga̧tarek and J. Jabłecki.  A nonparametric local volatility model for swaptions smile.  Journal of Computational Finance, 21(5):35–62, 2018.
  * [10] G. H. Golub and J. H. Welsch.  Calculation of Gauss quadrature rules.  Mathematics of Computation, 23(106):221–230, 1969.
  * [11] L. A. Grzelak.  On randomization of affine diffusion processes with application to pricing of options on VIX and S&P 500.  arXiv:2208.12518, 2022.
  * [12] L. A. Grzelak, J. Witteveen, M. Suarez-Taboada, and C. W. Oosterlee.  The stochastic collocation Monte Carlo sampler: highly efficient sampling from expensive distributions.  Quantitative Finance, 19(2):339–356, 2019.
  * [13] D. Heath, R. Jarrow, and A. Morton.  Bond pricing and the term structure of interest rates: A new methodology for contingent claims valuation.  Econometrica, 1(60):77–105, 1992.
  * [14] J. Hull and A. White.  Pricing interest-rate derivative securities.  The Review of Financial Studies, 3:573–592, 1990.
  * [15] F. Jamshidian.  An exact bond option formula.  The Journal of Finance, 44:205–209, 1989.
  * [16] C. W. Oosterlee and L. A. Grzelak.  Mathematical modeling and computation in finance.  World Scientific, 2019.
  * [17] L. C. G. Rogers and D. Williams.  Diffusions, Markov processes and martingales, volume 2 of Cambridge Mathematical Library.  Cambridge University Press, 2 edition, 2000.
  * [18] O. Vašiček.  An equilibrium characterization of the term structure.  Journal of Financial Economics, 5:177–188, 1977.

##  Appendix A Moments and optimal quadrature pairs,
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$

###  A.1 Moments and optimal quadrature pairs,
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$

Table 2: Selected distributions for the stochastic parameters. For the normal random variable with some mean, $\mu$, and variance, $\sigma^{2}$, it is sufficient to consider standard normal distribution, $\mathcal{N}(0,1)$, and properly scale the $\theta_{n}$ points, obtained from A.2. name | raw moment | domain
---|---|---
$\vartheta\sim\mathcal{U}([\hat{a},\hat{b}])$ | $\mathbb{E}[\vartheta^{n}]=\frac{\hat{b}^{n+1}-\hat{a}^{n+1}}{(n+1)(\hat{b}-\hat{a})}$ | $[\hat{a},\hat{b}]$
$\vartheta\sim\exp(\hat{a})$ | $\mathbb{E}[\vartheta^{n}]=\frac{n!}{\hat{a}^{n}}$ | $\mathbb{R}^{+}$
$\vartheta\sim\mathcal{N}(0,1)$ |  $\mathbb{E}[\vartheta^{n}]=(n-1)!!$ if $n$ even; $0$ otherwise | $\mathbb{R}$
$\vartheta\sim\Gamma(\hat{a},\hat{b})$ | $\mathbb{E}[\vartheta^{n}]=\hat{b}^{n}\Gamma(n+\hat{a})/\Gamma(\hat{b})$ | $\mathbb{R}^{+}$
$\vartheta\sim\chi^{2}(\hat{a},\hat{b})$ | $\mathbb{E}[\vartheta^{n}]=2^{n-1}(n-1)!(\hat{a}+n\hat{b})+\sum_{j=1}^{n-1}\frac{(n-1)!2^{j-1}}{(n-j)!}(\hat{a}+j\hat{b})\mathbb{E}[\vartheta^{n-j}]$ | $\mathbb{R}^{+}\cup\\{0\\}$

###  A.2 Optimal quadrature pairs, $\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$,
based on moments

Let $\\{p_{i}\\}_{i=0}^{N}$ with $\deg(p_{i})=i$ be a sequence of orthogonal
polynomials in $L^{2}$, with respect to PDF, $f_{\vartheta}(\vartheta)$, of
$\vartheta$, then the following holds,

$\mathbb{E}\left[p_{i}(\vartheta)p_{j}(\vartheta)\right]=\int_{\mathbb{R}}p_{i}(x)p_{j}(x)f_{\vartheta}(x){\rm
d}x=\delta_{i,j}\mathbb{E}\left[p_{i}^{2}(\vartheta)\right],\;\;\;i,j=0,\dots,N,$
(A.1)

with $\mathbb{R}$ the support of $\vartheta$, $\delta_{i,j}$ the Kronecker
delta. Orthogonal polynomials $\\{p_{i}\\}_{i=0}^{N}$ satisfy the following
recursion relation:

$\displaystyle{p}_{i+1}(x)$ $\displaystyle=$
$\displaystyle(x-\alpha_{i}){p}_{i}(x)-\beta_{i}{p}_{i-1},\;\;\;i=0,\dots,N-1,$
$\displaystyle{p}_{-1}(x)$ $\displaystyle\equiv$ $\displaystyle
0,\;\;{p}_{0}(x)\equiv 1,$

where $\alpha_{i}$ and $\beta_{i}$ are determined in terms of the moments of a
random variable $\vartheta$. Consider the monomials
$m_{i}(\vartheta)=\vartheta^{i}$, and define $\mu_{i,j}$ as
$\mu_{i,j}=\mathbb{E}\left[{m}_{i}(\vartheta){m}_{j}(\vartheta)\right]$, From
all moments $\mu_{i,j}$, we construct the Gram matrix
$M=\\{\mu_{i,j}\\}_{i,j=0}^{N}$, which is symmetric and contains all moments
$\\{1,\mathbb{E}[\vartheta^{1}],\dots,\mathbb{E}[\vartheta^{2N}]\\}$. Since
matrix $M$ is, by definition, positive definite [10], we decompose $M=R^{\rm
T}R$, by the Cholesky decomposition of $M.$

The following step relates the Cholesky upper-triangular matrix $R$ to the
orthogonal polynomials. This relationship has been established in [10] and is
given by,

$\displaystyle\alpha_{j}=\frac{r_{j,j+1}}{r_{j,j}}-\frac{r_{j-1,j}}{r_{j-1,j-1}},\;\;\;j=1,\dots,N,\;\;\;\text{and}\;\;\;\beta_{j}=\left(\frac{r_{j+1,j+1}}{r_{j,j}}\right)^{2},\;\;\;j=1,\dots,N-1,$
(A.2)

with $r_{0,0}=1$ and $r_{0,1}=0$ and where $r_{i,j}$ is the $(i,j)$-th element
of matrix $R$. Now, we determine the symmetric tridiagonal matrix, $J$,

$\displaystyle{J}:=\left[\begin{array}[]{ccccc}\alpha_{1}&\sqrt{\beta_{1}}&0&0&0\\\
\sqrt{\beta_{1}}&\alpha_{2}&\sqrt{\beta_{2}}&0&0\\\ &\ddots&\ddots&\ddots&\\\
0&0&\sqrt{\beta_{N-2}}&\alpha_{N-1}&\sqrt{\beta_{N-1}}\\\
0&0&0&\sqrt{\beta_{N-1}}&\alpha_{N}\\\
\end{array}\right]\in\mathbb{R}^{N\times N},$ (A.8)

with spectral factorization:

$J=W\Lambda W^{T},\;\;\;\Lambda=\text{\rm
diag}[\lambda_{1}^{*},\lambda_{2}^{*},\dots\lambda_{N}^{*}],\;\;\;WW^{\rm
T}=I,$

Due to the positivity of the off-diagonal entries, the eigenvalues
$\lambda_{j}^{*}$ are distinct, and all $W$’s first row entries are non-
vanishing. Moreover, it is well known [10] that the nodes and weights,
$\\{x_{i},w_{i}\\}_{i=1}^{N}$, of the Gauss rule are, for $i=1,\dots,N$, given
by $x_{i}=\lambda_{i}^{*}$ and $w_{i}=({\bf\epsilon}_{1}^{\rm
T}W{\bf\epsilon}_{i})^{2}$, where ${\bf\epsilon}_{j}$ is the $i$th axis vector
111The accompanying Python and MATLAB codes can be found at
https://github.com/LechGrzelak/Randomization..