Next Article in Journal
Free Vibration of Graphene Nanoplatelet-Reinforced Porous Double-Curved Shells of Revolution with a General Radius of Curvature Based on a Semi-Analytical Method
Previous Article in Journal
The Mathematical Simulation of South Korea’s Financial and Economic Impacts from Real Estate Bubbles: Lessons from the China Evergrande Collapse
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Calibration of the Kennedy Model

by
Dalma Tóth-Lakits
*,† and
Miklós Arató
Department of Probability Theory and Statistics, Eötvös Loránd University, 1117 Budapest, Hungary
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2024, 12(19), 3059; https://doi.org/10.3390/math12193059
Submission received: 31 August 2024 / Revised: 26 September 2024 / Accepted: 27 September 2024 / Published: 29 September 2024

Abstract

:
The Kennedy model offers a robust framework for modeling forward rates, leveraging Gaussian random fields to accommodate emerging phenomena such as negative rates. In our study, we employ maximum likelihood estimations to determine the parameters of the Kennedy field, utilizing Radon–Nikodym derivatives for enhanced accuracy. We introduce an efficient simulation method for the Kennedy field and develop a Black–Scholes-like analytical pricing formula for diverse financial assets. Additionally, we present a novel parameter estimation algorithm grounded in numerical extreme value optimization, enabling the recalibration of parameters based on observed financial product prices. To validate the efficacy of our approach, we assess its performance using real-world par swap rates in the latter part of this article.

1. Introduction

In the 2010s, a new phenomenon, negative rates, appeared in the financial markets, which brought extreme uncertainty to the world, resulting in the mathematical models used to describe the dynamics of interest rates being reconsidered. The model, defined by Kennedy in the 1990s, describes the dynamics of the forward rates with Gaussian random fields [1,2]. This approach has several advantages; for example, it offers a solution to handle negative rates naturally and can be connected to the industry standard Heath–Jarrow–Morton (HJM) framework [3]. Additionally, maximum likelihood estimations of the parameters and analytical Black–Scholes-like pricing formulas for different financial assets can be derived due to the standard distribution properties of the Gaussian random fields.
This article summarizes the most critical issues related to using the Kennedy model in the financial world. In Section 2, we present the term structure model for describing forward interest rates based on the Gaussian random fields proposed by Kennedy. Among other things, we present the condition for the martingale property of the discounted bond price and show which cases coincide with the Gaussian Heath–Jarrow–Morton framework. Section 3 introduces the theoretical background of the parameter estimations. The results of the Radon–Nikodym derivative of Gaussian measures with different means are shown. This section derives the maximum likelihood and probability one estimation for the parameters in the Kennedy field. Section 4 shows a practical, simple, and fast way to simulate the Kennedy field with the help of the Brownian sheet. The following Section 5 contains the analytical, fair price of various financial products (caplet, floorlet, and swap). Section 6 summarizes the calibration methods for different financial products, in our case, the optimization algorithm, which is based on a numerical extreme value search to estimate the parameters of the field. Finally, use of the previously presented calibration algorithm on real swap par rate data can be found in Section 7.

2. Kennedy Model

The development of forward rates in the model proposed by Kennedy is described in the following equation.
F ( s , t ) = α ( s , t ) + X ( s , t )
where X ( s , t ) is a centered Gaussian random field with covariance structure specified by
c o v [ X ( s 1 , t 1 ) , X ( s 2 , t 2 ) ] = c ( s 1 s 2 , t 1 , t 2 ) , 0 s i t i , i = 1 , 2 .
The function c is given and satisfies c ( 0 , t 1 , t 2 ) = 0 . We assume that the drift function α ( s , t ) is deterministic and continuous for 0 s t , and the initial term structure of α ( 0 , t ) , ( where t 0 ) is specified. Additionally, we also have E F ( 0 , t ) = α ( 0 , t ) for t 0 . The covariance function c ( s 1 s 2 , t 1 , t 2 ) is symmetric in t 1 and t 2 , and it is non-negative definite in pairs ( s 1 , t 1 ) and ( s 2 , t 2 ) . The dependence on s 1 s 2 ensures that the Gaussian random field X ( s , t ) exhibits independent increments.
A sufficient condition for the drift surface is established to guarantee that the discounted zero-coupon bond prices are martingales. Therefore, the model can be used to price financial products in the future.
First, let us introduce the following notations, where 0 s t .
R ( t ) = F ( t , t )
F Δ ( s , t ) = 1 Δ t t + Δ F ( s , u ) d u
P ( s , t ) = e s t F ( s , u ) d u
Z ( s , t ) = e 0 s R ( u ) d u P ( s , t )
F ( s ) = σ { F ( u , v ) , 0 u s , u v }
where R ( t ) denotes the spot rate at time t, P ( s , t ) represents the price at time s of a bond paying one unit at time t s . Z ( s , t ) defines the discounted price of the previously defined bond at time 0, with the information available at time s captured in the F ( s ) σ -algebra, indicating that the entire yield curve is observable at each time point. We also introduce a new notation, F Δ ( s , t ) , for the continuously compounded forward rate over the interval [ t , t + Δ ] , ( where Δ > 0 ), which can be interpreted as an average of the forward rate for the current period, at time s.
An important theorem is emphasized in Kennedy’s article, which states the following [2].
Theorem 1
(Kennedy (1997) [2]). In the independent-increments case the following statements are equivalent:
(a) 
The discounted bond-price process { Z ( s , t ) , F ( s ) , ( 0 s t ) } is a martingale for each t 0 ;
(b) 
P ( s , t ) = E e s t R ( u ) d u | F ( s ) , for all ( s , t ) , ( 0 s t ) ; and
(c) 
α ( s , t ) = α ( 0 , t ) + 0 t [ c ( s v , v , t ) c ( 0 , v , t ) ] d v for all ( s , t ) , ( 0 s t ) .
The proof of the theorem is accessible in the original article written by Kennedy [1]. Furthermore, a different derivation of the theorem can be found in Appendix A.1. To complete the proof, it is necessary to include an additional statement, which formulates an equivalent form of defining the drift term with the covariance function.
Remark 1.
The two statements for the drift term in the Kennedy model are equivalent. For all 0 s t
α ( s , t ) = α ( t , t ) + s t [ c ( s , v , t ) c ( v , v , t ) ] d v if and only if
α ( s , t ) = α ( 0 , t ) + 0 t [ c ( s v , v , t ) c ( 0 , v , t ) ] d v .
Proof of Remark 1. 
The proof is given by showing that both directions are correct.
α ( s , t ) α ( 0 , t ) = α ( t , t ) + s t [ c ( s , v , t ) c ( v , v , t ) ] d v α ( t , t ) 0 t [ c ( 0 , v , t ) c ( v , v , t ) ] d v
= s t [ c ( s , v , t ) c ( 0 , v , t ) ] d v + 0 s [ c ( v , v , t ) c ( 0 , v , t ) ] d v = 0 t [ c ( s v , v , t ) c ( 0 , v , t ) ] d v
α ( s , t ) α ( t , t ) = α ( 0 , t ) + 0 t [ c ( s v , v , t ) c ( 0 , v , t ) ] d v α ( 0 , t ) 0 t [ c ( t v , v , t ) c ( 0 , v , t ) ] d v
= 0 t [ c ( s v , v , t ) c ( v , v , t ) ] d v = s t [ c ( s , v , t ) c ( v , v , t ) ] d v

Connection Between HJM and the Kennedy Model

The Heath–Jarrow–Morton framework is a widely used model considered an industry standard [3]. This is also a term structure model, which creates a connection between bonds with different maturities. The HJM model is an infinite-dimensional framework. Therefore, the whole yield curve evolves in forward time instead of at a specific point.
This framework’s key point is to recognize an explicit relationship between drift and volatility parameters of the forward rate dynamics in a no-arbitrage world [4]. The critical assumption of the Heath–Jarrow–Morton model is that there is an elementary bond for each maturity. Overall, in an arbitrage-free term structure model, the forward rates must evolve like the following stochastic differential equation.
d F ( s , t ) = σ ( s , t ) s t σ ( s , u ) d u d s + σ ( s , t ) d W ( s )
Hereinafter, the main statement of the Heath–Jarrow–Morton model is that if there are no arbitrage opportunities, the forward rate drift is driven by volatility, known as the HJM no-arbitrage drift condition. This drift condition arises from the fact that the discounted process must be martingale [3].
Kennedy stated in his article that the Kennedy model encompasses the Heath–Jarrow–Morton framework in scenarios where the coefficients β ( s , t ) and σ i ( s , t ) in the underlying stochastic differential equations are deterministic, resulting in the rates F ( s , t ) being Gaussian [2]. Therefore, in this section, we will show precisely in which cases the two models can correspond to each other.
The notations of the HJM model are written consistently with those found in the book by Shreve ([4] on pages 423–435). We examine the case when a single Wiener process drives the forward interest rates, and β ( s , t ) and σ ( s , t ) are deterministic processes; then, the dynamics can be written as follows:
F ( s , t ) = F ( 0 , t ) + 0 s β ( u , t ) d u + 0 s σ ( u , t ) d W ( u )
where F ( 0 , t ) refers to the initial forward year curve known at time 0, W ( u ) is a Wiener process under the actual measure, and β ( s , t ) and σ ( s , t ) are deterministic processes in the variable s. Let us denote { ξ ( t ) } t 0 = { F ( 0 , t ) } t 0 , which is independent from process { W ( t ) } t 0 and is a Gaussian process.
The expected value and the covariance function of the Heath–Jarrow–Morton framework and the key Kennedy field conditions can be written in the following form.
  • The expected value function from the Heath–Jarrow–Morton model can be calculated in the following way.
    α ( s , t ) = E F ( s , t ) = E ξ ( t ) + 0 s β ( u , t ) d u = m ( t ) + 0 s β ( u , t ) d u
  • Similarly to the previously calculated expected value function, the covariance function is calculated as follows. Let us denote the covariance function between ξ ( t 1 ) , ξ ( t 2 ) with c o v ( ξ ( t 1 ) , ξ ( t 2 ) ) = r ( t 1 , t 2 )
    c ( s 1 , s 2 , t 1 , t 2 ) = c o v ( F ( s 1 , t 1 ) , F ( s 2 , t 2 ) ) = c o v ( ξ ( t 1 ) , ξ ( t 2 ) ) + 0 min ( s 1 , s 2 ) σ ( u , t 1 ) σ ( u , t 2 ) d u
    The covariance function is specified as a function of s 1 s 2 . This ensures that the Gaussian random field X ( s , t ) has independent increments in time s, which is also fulfilled due to point 2 in the HJM framework. This confirms that all Gaussian HJM models (where the drift and the volatility terms are deterministic) are the well-known Kennedy model.
  • By adding the martingale property into the Kennedy model (like in point (c) in Theorem 1), this guarantees that the conditional expected value of the discounted bond-price process is a martingale under the risk-neutral measure. As a result, the model is arbitrage-free. Then, by matching the equations of the expected values to each other, we obtain the famous condition of the HJM model, according to which the drift term can be obtained in the form below.
    α ( 0 , t ) + 0 s β ( u , t ) d u = α ( 0 , t ) + 0 t [ c ( s v , v , t ) c ( 0 , v , t ) ] d v
    α ( 0 , t ) + 0 s β ( u , t ) d u = α ( 0 , t ) + 0 t [ r ( v , t ) + 0 m i n ( s , v ) σ ( u , v ) σ ( u , t ) d u r ( v , t ) ] d v
    0 s β ( u , t ) d u = 0 t 0 m i n ( s , v ) σ ( u , v ) σ ( u , t ) d u d v
    β ( s , t ) = σ ( s , t ) s t σ ( s , v ) d v
    where the last equation equals the famous no-arbitrage HJM condition.
  • By adding the Markov property to the previous conditions, where the discounted bond price process is martingale, we obtain an even narrower class of models. We first define the following concepts based on Kennedy’s article for a random field { F ( s , t ) : 0 s t } [2].
    Definition 1
    (first Markov property). F satisfies the first Markov property if for all 0 s 1 s 2 < s 3 , s 1 t 1 and s 3 t 2 the following holds: F ( s 1 , t 1 ) F ( s 3 , t 2 ) | F ( s 2 , t 2 ) .
    Definition 2
    (second Markov property). F satisfies the second Markov property if for all 0 s 1 < s 2 and for any t 1 , t 2 with s 2 t 1 t 2 the following condition holds: F ( s 1 , t 1 ) F ( s 2 , t 2 ) | F ( s 2 , t 1 ) .
    Definition 3
    (Markov property). F is considered Markovian if it satisfies both the first and second Markov properties.
    Definition 4
    (Markov in t-direction). F is said to be Markovian in the t-direction, meaning in the maturity-time coordinate, if for all s t 1 t 2 t 3 the following condition holds
    F ( s , t 1 ) F ( s , t 3 ) | F ( s , t 2 )
    Definition 5
    (strict Markov property). F is considered strictly Markovian if it is both Markov and Markovian in the t-direction.
    Kennedy stated (in theorem 3.1 in [2]) that if a random field of forward rates is Markovian and satisfies the independent increments property, then the covariance function can be expressed in the following form.
    c ( s , t 1 , t 2 ) = f ( s ) g ( t 1 , t 2 ) ,
    where f is a monotone increasing, and g is a symmetric and positive, semidefinite function. This property can be written as follows for the HJM model.
    r ( t 1 , t 2 ) + 0 s σ ( u , t 1 ) σ ( u , t 2 ) d u = f ( s ) g ( t 1 , t 2 )
    Then, by deriving (23) according to the variable s we obtain
    σ ( s , t 1 ) σ ( s , t 2 ) = f ( s ) g ( t 1 , t 2 )
    By setting t 1 and t 2 equal to each other ( t 1 = t 2 = t ) , we obtain the following equality
    σ 2 ( s , t ) = f ( s ) g ( t , t )
    Consequently,
    σ ( s , t ) = b ( s ) g ( t ) ,
    where b ( s ) = f ( s ) and g ( t ) = g ( t , t ) . Therefore, it is shown that if the HJM model is Markovian, then the σ ( s , t ) function appears in the form of Equation (26). We thus obtained that in the Markovian case, the volatility function must be separable in the time parameters. Hence,
    σ ( s , t 1 ) σ ( s , t 2 ) = b 2 ( s ) g ( t 1 ) g ( t 2 )
    For s = 0 , Equation (23) can be written in the following form
    r ( t 1 , t 2 ) = f ( 0 ) · g ( t 1 , t 2 ) .
    From Equations (23), (28), and (39) it can be stated that
    r ( t 1 , t 2 ) = f ( s ) g ( t 1 , t 2 ) g ( t 1 ) g ( t 2 ) 0 s b 2 ( u ) d u
    f ( 0 ) g ( t 1 , t 2 ) = f ( s ) g ( t 1 , t 2 ) g ( t 1 ) g ( t 2 ) 0 s b 2 ( u ) d u
    g ( t 1 , t 2 ) = g ( t 1 ) g ( t 2 ) 0 s f ( u ) d u
    If the function f ( s ) is constant, then we obtain the trivial case when σ ( s , t ) = 0 for all ( s , t ) . In the non-trivial case, we obtain from Equation (31) that f ( s ) is not constant. Therefore, we obtain
    r ( t 1 , t 2 ) = c g ( t 1 ) g ( t 2 )
    Hence, we have shown that if the HJM model is Markovian, then functions σ and r occur in the previously derived form. Now, we show the opposite direction: if our covariance function has this shape, then the HJM model will be Markovian.
    c ( s , t 1 , t 2 ) = c g ( t 1 ) g ( t 2 ) + g ( t 1 ) g ( t 2 ) 0 s b 2 ( u ) d u
    = g ( t 1 ) g ( t 2 ) g ( t 1 , t 2 ) c + 0 s b 2 ( u ) d u f ( s )
    = f ( s ) g ( t 1 , t 2 )
    which is exactly the necessary condition (22).
    In 1992, Cheyette published an article in which a restriction was applied to the Heath–Jarrow–Morton model, which formed a subset of the original HJM models to make the model Markovian. This so-called Cheyette model is an arbitrage-free term structure model that is Markovian in a finite number of state variables and is consistent with any arbitrary initial term structure. Due to these favorable properties, the Cheyette model quickly spread throughout the industry and became widely used [5].
    In this case, the volatility function has to be separable into time- and maturity-dependent factors given by the following structure [6].
    σ ( s , t ) = α ( t ) β ( s ) α ( s )
    However, this condition is completely identical to the previously derived condition for the volatility term in the Markov case in the Kennedy model.
  • Kennedy further narrowed the model class by requiring stationarity in addition to the Markov property and the independent increments property (stated in Theorem 3.2 in [2]).
    Definition 6
    (stationary). F is stationary if, for each t > 0 , the joint distributions of { F ( s , t ) : 0 s t } are identical to those of { F ( s + u , t + u ) : 0 s t } for any fixed u > 0 .
    Therefore, the covariance function takes the form below:
    c ( s , t 1 , t 2 ) = e λ ( s ( t 1 t 2 ) ) · h ( | t 1 t 2 | )
    where λ 0 ,   | H ( x ) | h ( 0 ) e λ x 2 and x 0 .
    For the HJM framework, it was shown that r ( t 1 , t 2 ) = 0 , hence, according to point 5,
    c ( s , t 1 , t 2 ) = f ( s ) g ( t 1 , t 2 ) = g ( t 1 ) g ( t 2 ) c + 0 s b 2 ( u ) d u
    = e λ ( s ( t 1 t 2 ) ) · h ( | t 1 t 2 | )
    For s = 0 and t 1 = t 2 = t , it can be written
    c g 2 ( t ) = e λ t h ( 0 )
    g ( t ) = h ( 0 ) c · e λ t 2
    Returning to Equation (38)
    c + 0 s b 2 ( u ) d u = 1 g ( t 1 ) g ( t 2 ) e λ ( s ( t 1 t 2 ) ) · h ( | t 1 t 2 | )
    c + 0 s b 2 ( u ) d u = c h ( 0 ) e λ t 1 2 e λ t 2 2 e λ ( s ( t 1 t 2 ) ) · h ( | t 1 t 2 | )
    Now substituting s = 0
    c · h ( 0 ) = c exp λ 2 ( t 1 + t 2 ) λ ( t 1 t 2 ) h ( | t 1 t 2 | )
    h ( 0 ) = c exp λ 2 | t 1 t 2 | h ( | t 1 t 2 | )
    h ( u ) = h ( 0 ) exp λ 2 u
    Returning again to Equation (39), while substituting Equation (41)
    h ( 0 ) c exp λ 2 ( t 1 + t 2 ) c + 0 s b 2 ( u ) d u = e λ ( s ( t 1 t 2 ) ) · h ( | t 1 t 2 | )
    = e λ ( s ( t 1 t 2 ) ) · h ( 0 ) exp λ 2 | t 1 t 2 |
    1 c c + 0 s b 2 ( u ) d u = e λ s
    0 s b 2 ( u ) d u = c e λ s c
    By deriving the integral equation according to the variable s, we obtain the following solution
    b 2 ( u ) = c λ e λ u
    Therefore, the covariance function of the forward rates { F ( s , t ) : 0 s t } , when the rates are stationary, strictly Markov, and satisfy the independent-increments property, can be described with the following set of four parameters { σ , λ 0 , μ λ 2 , ν } and is of the form
    c o v [ F ( s 1 , t 1 ) , F ( s 2 , t 2 ) ] = σ 2 e λ min ( s 1 , s 2 ) + ( 2 μ λ ) min ( t 1 , t 2 ) μ ( t 1 + t 2 )
    The function of the expected value of the Gaussian random field can be easily derived from the covariance function.
    α ( s , t ) = ν σ 2 1 μ e μ ( t s ) 1 μ + 1 λ μ + e λ ( t s ) 1 λ μ

3. Parameter Estimation

In finance, where uncertainty reigns supreme and decisions are often made with incomplete information, accurate modeling of interest rate dynamics is paramount. This is where parameter estimation comes into play as a fundamental aspect of financial modeling, particularly in the context of Kennedy-type term structure models. While calibration is a widely adopted practice in finance, parameter estimation also holds significant importance.
The central assumption of these models is that interest rates follow stochastic processes, the parameters of which govern their behavior over time. These parameters determine the shape of the yield curve and influence the pricing of various financial instruments, such as bonds, options, and derivatives. Therefore, obtaining reliable estimates of these parameters is essential for making informed investment decisions, managing risk, and accurately pricing financial products.
Parameter estimation techniques enable practitioners to calibrate these models to observed market data, such as bond prices or interest rate derivatives. Among the most commonly used methods are maximum likelihood estimation (MLE), estimation with probability 1, and Radon–Nikodym derivatives, which allow determining parameter values that maximize the likelihood of observing the given market data under the model assumptions. Through rigorous statistical inference, these techniques provide a systematic framework for extracting information from observed market prices and estimating the underlying dynamics of interest rates.

3.1. Maximum Likelihood Estimations

This section presents the theoretical background of maximum likelihood estimations for Gaussian functionals, drawing on the work of Rozanov and Arató [7,8]. The following definitions, theorems, and theoretical results were published in a conference proceedings in Barcelona, where we gave a presentation [9].
Definition 7
(Gaussian functional). Let ( Ω , A , P ) be a probability space, and let T be a parameter set. In this context, ξ : Ω × T R is considered a Gaussian functional if, for any n N and c 1 , , c n R , t 1 , , t n T , the following expression is normally distributed:
i = 1 n c i ξ t i
In this case, P is referred to as a Gaussian measure on ( Ω , F ξ ) . For simplicity, we can assume that A = F ξ . The expected value and the covariance of the Gaussian functional are denoted as follows:
m ( t ) = E ξ ( t ) , B ( s , t ) = c o v [ ξ ( s ) , ξ ( t ) ]
It is well-established that two Gaussian measures are either equivalent or orthogonal.

3.1.1. The Case of Different Expected Values

Let ξ : Ω × T R be a Gaussian functional. Assume that the expected value of the Gaussian functional under the measure P is 0, while the expected value under the measure P 1 is m.
E P ξ ( t ) = 0 , E P 1 ξ ( t ) = m ( t ) , t T
Let U represent the linear space of variables structured in the following manner, while U ¯ is the Hilbert space obtained by closing U:
i = 1 n c i ξ t i , n N , c 1 , , c n R , t 1 , , t n T .
In addition, take the following scalar product:
< u , v > = Ω u v d P .
The upcoming theorem can be found in [8].
Theorem 2
(Rozanov). The measures P and P 1 are equivalent if and only if there exists an η U ¯ such that
m ( t ) = ω ξ ( t ) η ( t ) d P , t T .
The Radon–Nikodym derivative of the two measures in the case of equivalence is given by
d P 1 d P = e η < η , η > 2
A straightforward consequence of this theorem is the following statement made by Arató [7].
Theorem 3
(Arató). Let ξ : Ω × T R be a Gaussian functional. Assume that the expected value of the Gaussian functional under the measure P is 0 and under the measure P 1 is m · a ( t ) . The measures P and P 1 are equivalent if and only if there exists an η U ¯ such that
a ( t ) = ω ξ ( t ) η ( t ) d P , t T .
The Radon–Nikodym derivative of the two measures in the case of equivalence is given by
d P 1 d P = e m η m 2 < η , η > 2
Theorem 4
(Maximum likelihood estimation). Using the notations from the previous theorem, the maximum likelihood estimate of m is given by
m ^ = η < η , η > .
The estimation is normally distributed and unbiased, and the standard deviation is
D P 1 2 m ^ = 1 < η , η >
Proof of Theorem 4. 
The form of the estimation is directly obtained from the Radon–Nikodym derivative. To ascertain the expected value and the variance, we must calculate the following expected value if X N ( 0 , σ 2 ) .
E ( X k e m X ) = x k e m x 1 2 π σ e x 2 2 σ 2 d x = e m 2 σ 2 2 x k 1 2 π σ e ( x m σ 2 ) 2 2 σ 2 d x
we obtain the following values by calculating the first two moments:
  • if k = 1 , then → E ( X e m X ) = m σ 2 e m 2 σ 2 2
  • if k = 2 , then → E ( X 2 e m X ) = ( σ 2 + m 2 σ 4 ) e m 2 σ 2 2
For the expected value, we obtain the following:
E P m m ^ = 1 < η , η > Ω η d P m = 1 < η , η > Ω η e m η m 2 < η , η > 2 d P
= 1 < η , η > e m 2 < η , η > 2 m < η , η > e m 2 < η , η > 2 = m
Similarly to the first moment, we can derive the second moment.
E P m m ^ 2 = 1 < η , η > 2 Ω η 2 d P m = 1 < η , η > 2 Ω η 2 e m η m 2 < η , η > 2 d P
= 1 < η , η > 2 e m 2 < η , η > 2 ( < η , η > + m 2 < η , η > 2 ) e m 2 < η , η > 2 = 1 < η , η > + m 2 .
The standard deviation can be derived directly from these previous calculations. Given that these are Gaussian functionals, the normality is evident. □

3.1.2. The Case of Constant Expected Value

In cases where the expected value of our Gaussian process is constant, a ( t ) = 1 for every t T . Let F = σ { ( ξ ( t ) ξ ( s ) ), s , t T } . Fix a point t 0 T and define h ( ξ ) = E P [ ξ ( t 0 ) F ) ] . Assuming that D 2 [ ξ ( t 0 ) h ( ξ ) ] > 0 , the maximum likelihood estimate of m is given by
m ˜ = ξ ( t 0 ) h ( ξ ) .
Proof. 
The proof uses the law of total expectation and the E p m ˜ ( ξ ( t 0 ) h ( ξ ) ) = D P 2 ( m ˜ ) statement.
E P ( m ˜ ξ ( t ) ) = E P ( m ˜ ( ξ ( t ) ξ ( t 0 ) + ξ ( t 0 ) h ( ξ ) + h ( ξ ) ) )
= E P ( m ˜ ( ξ ( t ) ξ ( t 0 ) + h ( ξ ) ) ) + D P 2 ( m ˜ )
= E P ( E P ( ( m ˜ ( ξ ( t ) ξ ( t 0 ) + h ( ξ ) ) F ) ) + D P 2 ( m ˜ ) = D P 2 ( m ˜ )
Based on the previous derivations, the maximum likelihood estimation is the following:
m ^ = m ˜ / D P 2 ( m ˜ ) D P 2 ( m ˜ / D P 2 ( m ˜ ) )

3.1.3. Some Simple Examples

The established results of various stochastic processes commonly used to model financial processes directly stem from the previous theorems.
For instance, consider a Gaussian process with an expected value of m and the same covariance as the Wiener process over the interval [ a , b ] , where a > 0 and a < b . In this case, the maximum likelihood estimation of the Wiener process is the value of the process at the start, m ˜ = ξ ( a ) .
Similarly, a stationary Ornstein–Uhlenbeck process can be observed over the interval [ 0 , T ] . In this case, the value of λ > 0 is known in advance, as well as the expected value and the covariance matrix of the process.
E P m [ ξ ( t ) ] = m
c o v P m [ ( ξ ( s ) , ξ ( t ) ) ] = σ 2 e λ | t s | , s , t [ 0 , T ]
Therefore, the following covariances can be easily determined:
E P [ ξ ( 0 ) ξ ( t ) ] = σ 2 e λ t , E P [ ξ ( t ) ξ ( T ) ] = σ 2 e λ ( T t )
E P 0 T ξ ( s ) d s · ξ ( t ) = σ 2 2 e λ t e λ ( T t ) λ .
By leveraging the fact that, in this case, the maximum likelihood estimate is unbiased, we obtain the well-known result from Grenander [10]:
m ^ = ξ ( 0 ) + ξ ( T ) + λ 0 T ξ ( s ) d s 2 + λ T .

3.2. Parameter Estimations of the Kennedy Field

From now on, we investigate the case when the random field of forward rates is strictly Markov, stationary, and satisfies the property of independent increments. Then, as we have seen previously, the covariance and expected value functions have the form (52) and (53), and these functions are defined by four parameters ( ν , μ , α and σ ). Therefore, the expected initial forward curve is easily obtained.
α ( 0 , t ) = ν σ 2 μ + σ 2 μ e μ t + σ 2 λ μ e μ t σ 2 λ μ e λ t
= ν + σ 2 μ e μ t 1 + σ 2 λ μ e μ t e λ t
In addition, from the equation above, it can be straightforwardly seen that the parameter ν refers to the expected value of the spot curve.
E F ( s , s ) = E R ( s ) = α ( s , s ) = ν σ 2 μ + σ 2 μ + σ 2 λ μ σ 2 λ μ = ν
It is evident that the F ( s , s + t ) field constitutes an Ornstein–Uhlenbeck process with respect to the variable s, therefore
c o v [ F ( s 1 , s 1 + t ) , F ( s 2 , s 2 + t ) ] = σ 2 e λ t e μ | s 1 s 2 |
This implies that if we can observe the F ( s , s + t ) process over an interval defined by s for some specific value of t, then σ 2 e λ t μ is determined with probability 1. If we can achieve this for two distinct values of t, then both σ 2 μ and λ are defined with probability 1.
Examining another covariance from the field,
c o v F log s 1 λ , t , F log s 2 λ , t = σ 2 e λ t min ( s 1 , s 2 )
This means that σ 2 e λ t is defined with probability 1, which in turn implies that both σ 2 and μ are also defined with probability 1.
In the following, we observe the field on a region marked with T. The following ξ ( s , t ) auxiliary random field is introduced, where the expected value under the measure P ν is ν .
ξ ( s , t ) = F ( s , t ) + σ 2 1 μ e μ ( t s ) 1 μ 1 λ μ + e λ ( t s ) 1 λ μ .
where
W ( x i , y j ) = k = 1 i l = 1 j ξ ( k , l )
We demonstrate that the following estimate is the maximum likelihood estimate of this parameter.
ν ^ = e λ b 1 μ ξ ( a , b 1 ) + e λ b 2 μ λ ξ ( a , b 2 ) + b 1 b 2 e λ ν ξ ( a , v ) d v e λ b 2 1 λ + 1 μ λ + e λ b 1 1 μ 1 λ .
First, we obtain that E P 0 ( ξ ( s , t ) ν ^ ) gives the same value for every ( s , t ) T . On the other hand, E P ν ( ν ^ ) = ν . Thus, based on Theorems 2 and 3, ν ^ is the maximum likelihood estimate.

4. Simulation of the Kennedy Field

This section aims to simulate the Kennedy field in n × m points. We can consider this as an n × m normally distributed vector whose expected value and covariance matrix are known. However, for sufficiently large n and m, simulating a multidimensional, normally distributed vector becomes exceedingly slow, due to the size of the covariance matrix. A much more effective, faster, and simpler way is to observe that if W ( x , y ) is a Brownian sheet, then α ( s , t ) + σ e μ t W ( e λ s , e ( 2 μ λ ) t ) forms a Kennedy field with the appropriate covariance structure.
The question is how can we most efficiently generate a Brownian sheet at the ( x i , y j ) points ( x 1 < < x n , y 1 < < y m ) , where the division is not necessarily equidistant. Let us take independent random variables with distributions N ( 0 , ( x i x i 1 ) ( y j y j 1 ) ) , ( where x 0 = y 0 = 0 ) and denote them as η ( i , j ) . Consequently, the Brownian sheet can be expressed in the following form:
W ( x i , y j ) = k = 1 i l = 1 j η ( k , l )
Therefore, the subsequent matrix operation should be implemented efficiently to obtain the desired results.
A B : B ( i , j ) = k = 1 i l = 1 j A ( k , l )
Fortunately, ready-made, fast algorithms exist for this double summation.
In Figure 1, two different realizations of the forward rate field (denoted as F ( s , t ) , 0 s t ) are depicted, generated using the simulation algorithm described earlier. The figures contain 10,000 simulated points ( 100 × 100 = 10,000 ). The higher rates are represented in light yellow, while as the values decrease, they transition into dark blue.
Since the rates have not yet been calibrated to market data, the field exhibits considerable volatility. However, it is noteworthy that the Kennedy model is capable of producing negative forward rate values. It is also important to highlight that the figures clearly show a general increase in forward rates for longer maturities, which is consistent with the empirical observations typically seen in the market.

5. Option Pricing

This section aims to show that fair prices for various financial assets can be derived analytically if we assume that the forward rates evolve according to a Gaussian random field.

5.1. European Caplet

In the case of options, compounded forward rates are used instead of instantaneous forward rates. Consequently, it is necessary to transition from the instantaneous forward rate described earlier by the Kennedy field to a discrete forward rate for a given time period, often denoted as L ( t , T i ) , concerning the LIBOR rate. Consistently with the following discretization scheme, the discretized version of the HJM framework, which is considered the industry standard, is the LIBOR market model (LMM).
1 + L ( s , t ) Δ = e t t + Δ F ( s , u ) d u = e Δ F Δ ( s , t )
This derivation is equivalent to the one derived by Kennedy but uses a different approach. We aim to calculate the price of an interest rate caplet with a strike K for the period from t to t + Δ . This can be interpreted as a European option on the forward rate given by F Δ ( s , t ) = 1 Δ t t + Δ F ( s , u ) d u which is exercised at time t if f Δ ( t , t ) > K , resulting in a payoff at time t + Δ . The payoff function of this transaction is shown below.
V ( t , K ) = [ ( e Δ F Δ ( t , t ) 1 ) ( e Δ K 1 ) ] + = [ e Δ F Δ ( t , t ) e Δ K ] +
The discount factor from time s to time t is defined as follows.
D ( s , t ) = e s t r ( u ) d u
A cap normally consists of a series of such options for successive time periods; however, it is sufficient in this context to consider only a single time period. The discounted payoff of the option at time s is given by the following expression:
D ( s , t + Δ ) V ( t , K ) = e s t + Δ r ( u ) d u ( e Δ F Δ ( t , t ) e Δ K ) +
The price of a financial asset is obtained by taking the expected value of the discounted payoff function. The definition of the drift term guarantees that the model is under a risk-neutral measure, just like in the Heath–Jarrow–Morton framework.
P c a p l e t ( s ) = E [ e s t + Δ r ( u ) d u ( e Δ F Δ ( t , t ) e Δ K ) + ]
For the sake of simplicity, two additional variables are introduced ( ξ ( s , t ) and η ( s , t ) ) to denote the time range over which the forward rate is integrated.
ξ ( s , t ) = s t r ( u ) d u = s t F ( u , u ) d u
η ( s , t ) = s t F ( t , u ) d u
Hence, in the case of a caplet, we deal with the following special case of ξ and η .
ξ ( s , t + Δ ) = s t + Δ r ( u ) d u = s t + Δ F ( u , u ) d u
η ( t , t + Δ ) = Δ F Δ ( t , t ) = t t + Δ F ( t , u ) d u
Due to the properties of the Gaussian random field, ( ξ ( s 1 , t 1 ) , η ( s 2 , t 2 ) ) follows a multivariate normal distribution. Henceforth, except for necessary cases, we omit the corresponding time indices to indicate the expected value, standard deviation, and correlation between ξ and η . Consequently, let us denote them with the following notations E ξ = μ 1 , D 2 ( ξ ) = σ 1 2 , E η = μ 2 , and D 2 ( η ) = σ 2 2 . From now on, the conditional normal distribution theorem can be used. As a result, the conditional distribution of ξ given η is the following:
ξ | η N μ 1 + ρ σ 1 η μ 2 σ 2 , σ 1 2 ( 1 ρ 2 )
where c o r r ( ξ ( s 1 , t 1 ) , η ( s 2 , t 2 ) ) = ρ ( s 1 , t 1 , s 2 , t 2 ) . Therefore, the fair price of the European option can be calculated as follows.
E [ e s t + Δ r ( u ) d u ( e Δ F Δ ( t , t ) e Δ K ) + ] = E [ e ξ ( e η e Δ K ) + ]
= E [ E ( e ξ ( e η e Δ K ) + | η ] = E [ ( e η e Δ K ) + · E ( ( e ξ ) | η ) ]
During the derivations, the law of total expectation and the fact that ( e η e Δ K ) + is measurable for η is used. As we can see, ξ N ( μ 1 , σ 1 ) is normally distributed; therefore, ξ N ( μ 1 , σ 1 ) , where c o r r ( ξ , η ) = ρ . Therefore, ξ | η N ( μ 1 ρ σ 1 η μ 2 σ 2 , σ 1 2 ( 1 ρ 2 ) ) . Since the conditional distribution of ξ given η is known, E [ e ξ | η ] can be calculated as the expectation of a lognormally distributed random variable.
E [ e ξ | η ] = e μ 1 ρ σ 1 η μ 2 σ 2 + 1 2 σ 1 2 ( 1 ρ 2 )
Returning to the pricing formula,
E [ ( e η e Δ K ) + · E [ e ξ | η ] ] = E ( e η e Δ K ) + · e μ 1 ρ σ 1 η μ 2 σ 2 + 1 2 σ 1 2 ( 1 ρ 2 )
= e μ 1 + 1 2 σ 1 2 ( 1 ρ 2 ) + ρ μ 2 σ 1 σ 2 E ( e η e Δ K ) + · e ρ η σ 1 σ 2
= e μ 1 + 1 2 σ 1 2 ( 1 ρ 2 ) + ρ μ 2 σ 1 σ 2 Δ K e x ( 1 ρ μ 2 σ 1 σ 2 ) e Δ K x ρ σ 1 σ 2 1 2 π σ 2 e ( x μ 2 ) 2 2 σ 2 2 d x
= e μ 2 μ 1 + σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 Δ K 1 2 π σ 2 e ( x μ 2 σ 2 2 + ρ σ 1 σ 2 ) 2 2 σ 2 2 d x
e Δ K μ 1 + 1 2 σ 1 2 Δ K 1 2 π σ 2 e ( x μ 2 + ρ σ 1 σ 2 ) 2 2 σ 2 2 d x
Finally, by subtracting the values of the two integrals from each other, we obtain the analytical pricing formula for the European call option in the case of the Kennedy fields.
P c a p l e t ( s ) = e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) Φ μ 2 + σ 2 2 ρ σ 1 σ 2 Δ K σ 2 e Δ K μ 1 + 1 2 σ 1 2 Φ μ 2 ρ σ 1 σ 2 Δ K σ 2

Expected Values and Variances

Based on the calculations in Appendix A.2 for the caplet pricing, the expected value of ξ and η , their standard deviation, and the correlation between them are as follows.
μ 1 = E ξ ( s , t + Δ ) = ν σ 2 μ ( t + Δ s )
μ 2 = E η ( t , t + Δ ) = ν σ 2 μ Δ σ 2 μ 2 e μ Δ 1 σ 2 μ ( λ μ ) e μ Δ 1 + σ 2 λ ( λ μ ) e λ Δ 1
σ 1 2 = D 2 ξ ( s , t + Δ ) = 2 σ 2 μ 2 ( t + Δ s ) μ + e μ ( t + Δ s ) 1
σ 2 2 = D 2 η ( t , t + Δ ) =
= σ 2 ( λ μ ) λ e λ Δ 1 + σ 2 μ ( μ λ ) e μ Δ 1 + σ 2 μ ( μ λ ) e μ Δ e λ Δ + σ 2 λ μ 1 e λ Δ
c o v = c o v ξ ( s , t + Δ ) , η ( t , t + Δ )
= σ 2 μ 2 1 e μ Δ e μ ( t s ) + e μ ( t + Δ s ) +
+ σ 2 λ ( μ λ ) + σ 2 λ μ e λ Δ 1 + 2 σ 2 μ ( μ λ ) e λ Δ μ Δ 1
ρ = c o r r ξ ( s , t + Δ ) , η ( t , t + Δ ) = c o v ξ ( s , t + Δ ) , η ( t , t + Δ ) D ξ ( s , t + Δ ) D η ( t , t + Δ ) = c o v σ 1 σ 2
c o v ( ξ , η ) = σ 2 μ 2 1 e μ Δ e μ ( t s ) + e μ ( t + Δ s ) +
+ σ 2 λ ( μ λ ) + σ 2 λ μ e λ Δ 1 + 2 σ 2 μ ( μ λ ) e λ Δ μ Δ 1

5.2. European Floorlet

Similarly to the previously derived European caplet, the price of an interest rate floorlet is now derived. Thereby, using the put-call parity, the pricing formula of the swap can be easily calculated. The payoff function of this transaction is shown below.
V ( t , K ) = [ ( e Δ K 1 ) ( e Δ F Δ ( t , t ) 1 ) ] + = [ e Δ K e Δ F Δ ( t , t ) ] +
The discount factor from time s to time t is defined as follows: D ( s , t ) = e s t r ( u ) d u . A floorlet typically consists of a series of such options for successive time periods; however, it is sufficient in this context to consider only a single time period. The following expression gives the discounted payoff of the option at time s:
D ( s , t + Δ ) V ( t , K ) = e s t + Δ r ( u ) d u ( e Δ K e Δ F Δ ( t , t ) ) +
The price of a financial asset is obtained by taking the expected value of the discounted payoff function. The definition of the drift term guarantees that the model is under a risk-neutral measure, just like in the Heath–Jarrow–Morton framework.
P f l o o r l e t ( s ) = E [ e s t + Δ r ( u ) d u ( e Δ K e Δ F Δ ( t , t ) ) + ]
The derivation is exactly the same as the price of the caplet product presented earlier and can be found in Appendix A.3. Therefore, the analytical pricing formula for the European floorlet option in the case of Kennedy fields is as follows.
P f l o o r l e t ( s ) = e μ 1 + 1 2 σ 1 2 ( 1 ρ 2 ) e Δ K + 1 2 ρ 2 σ 1 2 Φ Δ K μ 2 + ρ σ 1 σ 2 σ 2 e μ 2 + 1 2 ( σ 2 ρ σ 1 ) 2 Φ Δ K μ 2 σ 2 2 + ρ σ 1 σ 2 σ 2
= e Δ K μ 1 + 1 2 σ 1 2 Φ Δ K μ 2 + ρ σ 1 σ 2 σ 2 e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) Φ Δ K μ 2 σ 2 2 + ρ σ 1 σ 2 σ 2

5.3. Swap

An interest rate swap is a forward contract exchanging a floating and fixed rate for a predetermined period. The special financial asset in which the floating versus fixed rate exchange only applies to one period is called a swaplet. In this section, the fair price, and hence the conditional expected value of the discounted payoff function under the risk-neutral measure for one period, is derived.
We first examine the simplest case, a one-period swap. In this case, the interest rate exchange takes place at only one point in time, T in the swap product. In this case, the swap is similar to a caplet with an extreme cap value, where the product is definitely worth calling. The price of the swaplet product at time s is as follows.
P s w a p l e t ( s ) = E [ e s s + Δ F ( u , u ) d u ( e s s + Δ F ( s , u ) d u e Δ K ) ]
In this case, the previously introduced ξ and η are interpreted in the following time period.
ξ ( s , s + Δ ) = s s + Δ F ( u , u ) d u
η ( s , s + Δ ) = Δ F Δ ( s , s ) = s s + Δ F ( s , u ) d u
As we can see, the definition of η is unchanged; therefore, only the value of μ 1 , σ 1 , and the covariance change.
μ 1 = ν σ 2 μ Δ
σ 1 2 = D 2 ξ = 2 σ 2 μ 2 Δ μ + e μ Δ 1
c o v ( ξ , η ) = σ 2 ( λ μ ) λ e λ Δ 1 + σ 2 μ ( μ λ ) e μ Δ 1
+ σ 2 μ ( μ λ ) e μ Δ e λ Δ + σ 2 λ μ 1 e λ Δ
σ 2 2 = D 2 η = c o v ( ξ , η ) = ρ σ 1 σ 2
As we can see, in that case, the covariance of ξ and η equals the variance of η .
This calculation is easy to understand if we use the previously introduced multidimensional normally distributed random variables ( ξ , and η ) because, in this case, e ξ and e η random variables are lognormally distributed. It is well-known that the quotient of two lognormally distributed random variables with correlation ρ is also lognormally distributed with the following expected value and standard deviation.
ξ N ( μ 1 , σ 1 ) , η N ( μ 2 , σ 2 ) ( η ξ ) N ( μ 2 μ 1 , σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 )
Hence,
E ( e ξ ) = e μ 1 + 1 2 σ 1 2
E ( e ξ e η ) = E e η e ξ = E ( e η ξ ) = e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 )
Therefore, we can easily obtain the previously calculated result.
P s w a p l e t ( s ) = e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) e Δ K μ 1 + 1 2 σ 1 2
Furthermore, the price of a one-period long swap, the so-called swaplet at time s, can be easily obtained using the previously derived caplet and floorlet pricing formulas and the put-call parity. Therefore, this is he difference between the calculated fair price of the caplet and the floorlet option.
The fair price of a fixed vs. floating swap for more time periods at time 0 can be found in Appendix A.2.

5.4. Par Swap Rate

In the previous section, we derived the fair price of a one-period swap, the so-called swaplet.
P s w a p l e t ( 0 ) = e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) e Δ K μ 1 + 1 2 σ 1 2
Therefore, let us first adjust the previously defined ξ and η variables to the following time periods.
ξ ( s , s + Δ ) = s s + Δ F ( u , u ) d u
η ( s , s + Δ ) = Δ F Δ ( s , s ) = s s + Δ F ( s , u ) d u
The so-called swap quote can be easily expressed from that equality, which equals the par swap rate. The par rate is the fixed rate that results in the swap having a zero present value, meaning it is the rate that equates the value of both legs of the swap. The derivation of this rate is important because, in many cases, the financial data contain par swap rates instead of swap prices; in other words, this financial product is quoted using the par swap rate.
P s w a p l e t ( 0 ) = 0 = e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) e Δ K μ 1 + 1 2 σ 1 2
e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) = e Δ K μ 1 + 1 2 σ 1 2
μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) = Δ K μ 1 + 1 2 σ 1 2
μ 2 + 1 2 σ 2 2 ρ σ 1 σ 2 = Δ K
K = 1 Δ ( μ 2 + 1 2 σ 2 2 ρ σ 1 σ 2 )
After that, we want to express the par swap rate with the original parameters of the Kennedy field.
Δ K = μ 2 + 1 2 σ 2 2 ρ σ 1 σ 2 = μ 2 + 1 2 σ 2 2 σ 2 2 = μ 2 1 2 σ 2 2
= ν σ 2 μ Δ σ 2 μ 2 e μ Δ 1 σ 2 μ ( λ μ ) e μ Δ 1 + σ 2 λ ( λ μ ) e λ Δ 1
σ 2 2 ( λ μ ) λ e λ Δ 1 + σ 2 2 μ ( λ μ ) e μ Δ 1 + σ 2 2 μ ( λ μ ) e μ Δ e λ Δ + σ 2 2 λ μ e λ Δ 1
= ν σ 2 μ Δ + e μ Δ 1 σ 2 μ 2 σ 2 μ ( λ μ ) + σ 2 2 μ ( λ μ ) + σ 2 2 μ ( λ μ )
+ e λ Δ 1 σ 2 λ ( λ μ ) + σ 2 2 λ ( λ μ ) + σ 2 2 λ μ σ 2 2 μ ( λ μ )
= ν σ 2 μ Δ + σ 2 μ 2 1 e μ Δ σ 2 2 λ ( λ μ ) + σ 2 2 μ λ σ 2 2 μ ( λ μ ) 1 e λ Δ
= ν σ 2 μ Δ + σ 2 μ 2 1 e μ Δ
From here, the par swap rate can be easily written with the original parameters of the Kennedy field.
K = ν σ 2 μ + σ 2 μ 2 Δ 1 e μ Δ
Therefore, if the swap par rate can be observed for at least four different tenors, then three parameters of the original four ( ν , μ , and σ ) can be determined with a probability of 1. However, it is worth mentioning that the parameter λ is omitted from the description of the par swap rate.

6. Calibration on Simulated Data

Simulated financial caplet, floorlet, and swaplet prices were generated using the previously derived analytical pricing formulas with different maturities and strikes. Therefore, first of all, we just wanted to test the punctuality of the calibration engine.
Numerical calibration is an extreme value optimization problem. The method aimed to find the parameter set that minimized the squared deviation error between the previously generated financial caplet, floorlet, and swaplet data and the analytically calculated prices with the calibrated parameters. The calibration engine was based on the stochastic gradient descent method. The extreme value optimization was based on the article of Mikhaliov and Nögel, and the implementation in Python was based on the work of Emerick and Tatsat [11,12].
The Figure 2 shows that the prices calculated with the back-estimated parameters fit our synthetic dataset almost perfectly. The calibration returned the used parameters; the difference was negligible and could be considered a numerical error.

7. Calibration on Real Data

A time series of swap par rates was obtained with the help of the Bloomberg terminal. The financial dataset contains the par swap rates of the USD SOFR fixed versus floating interest rate swaplet from July 2018 to April 2023. The historical dataset includes par swap rates for 28 different maturities daily.
We calibrated the model daily for different maturities for par swap rates with the calibration algorithm going 100 days back. As we can see in Figure 3, the Kennedy field fit the dataset nicely; however, it slightly overestimated the values for shorter maturities, while slightly underestimating the par swap rates at long maturities.
In addition to the analytical results, it can also be seen that λ did not play a role in the numerical implementation, since this back-estimated parameter value is highly volatile. Meanwhile, the value of the other three parameters varied on a much smaller scale, and similar trends could be observed in Figure 4.
As a result, our guess is that the λ parameter, which is not included in the par swap rate, describes a temporal relationship; in other words, the term structure of the model; σ , greatly influences the standard deviation of the field; while ν is used to describe the level of the yield curves, since ν is the parameter that describes the expected value of the spot rate ( ν = E F ( s , s ) ) .
In the following, we plotted Kennedy fields (shown in Figure 5) for describing forward interest rates with the three parameters back-estimated from the par swap rate dataset ( ν = 0.05171817 , μ = 0.56028928 and σ = 0.11315586 ) and three different λ values, to see what rates would be generated in a realistic case.

8. Discussion

Kennedy introduced a model based on Gaussian random fields for modeling forward interest rates in the 1990s, which, due to its normal distribution, can generate negative interest rates. At that time, this scenario was not considered feasible; however, in the interest rate environment of the 2010s, negative interest rates emerged. Since this model naturally handles them, this underscored the relevance of the model. Calibration on actual par swap rates demonstrated that our model fits well with the current interest rate environment and effectively describes the market.
Moving forward, our primary objective is to go beyond analytical pricing formulas and utilize models based on artificial intelligence, including LSTM and neural networks, for parameter estimation. We aim to compare the accuracy of parameters recovered through calibration with these AI-based models. Additionally, we plan to investigate the temporal stability of parameters and compare them with industry-standard models such as SABR.

9. Conclusions

Our article focused on the mathematical model based on Gaussian random fields introduced by Kennedy to describe forward interest rates. Among other things, we provided a novel proof for the equivalence of conditions regarding the martingale property of discounted bond prices. We demonstrated the relationship between the Kennedy model and the HJM framework in special cases (Markov property, stationarity). Additionally, utilizing Radon–Nikodym derivatives, we derived maximum likelihood estimates and estimates with a probability of one for the original parameters of the field. We presented a new, efficient method, based on Brownian sheets, to simulate the Kennedy field.
Subsequently, we derived analytical pricing formulas resembling Black–Scholes for various financial products, including caplets, floorlets, swaplets, and swaps. Finally, we calibrated the field using a numerical extreme value search algorithm based on stochastic gradient descent on a simulated synthetic dataset to recover the original parameters. We then calibrated it on actual par swap rates to examine how our model performed in a market environment.
In summary, we can conclude that a new result has been derived in estimating the parameters of the Kennedy field, leading to the development of an effective calibrator. This serves as a strong foundation for further investigation and comparison with other, more complex models.
Additionally, our research aims to calibrate the parameters of negative interest rate models using machine learning algorithms, enabling us to compare these results with previously derived analytical estimates.

Author Contributions

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

Funding

This work was supported by the KDP-2021 program and the ELTE TKP 2021-NKTA-62 funding scheme of the Ministry of Innovation and Technology from the source of the National Research, Development and Innovation Fund.

Data Availability Statement

Data used in this study include historical par swap rates of the USD SOFR fixed vs. floating interest rate swaps (Bloomberg ticker: USOSFR1Z BGN Curncy) obtained from Bloomberg. The data spans from 20 April 2018 to 20 April 2023, with daily frequency in USD, sourced from BGN. Due to proprietary restrictions, these data are not publicly available. Access to the data is subject to Bloomberg’s terms and conditions and is not shared openly. For more information, please contact the authors.

Acknowledgments

We would like to take this opportunity to thank Csaba Kőrössy, for his help and dedication to the research. Even though he joined the research late on, he spent much time understanding it, and his comments and suggestions immensely helped this article. We want to thank Fáth Gábor, who involved us in Risklab and since then has regularly consulted, given ideas, and helped on this topic, in which he also has great expertise. Finally, we would like to thank András Ványolos for his help, interest, and the enthusiasm with which he joined our project; who is motivated simply by his love for mathematical research. We enjoy doing math with you.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

Appendix A.1

The proof of Theorem 1 can be found in the original article by Kennedy [1]. However, we present an alternative derivation to prove the statement. To facilitate this, we first introduce an auxiliary lemma, which will be utilized in the proof of the theorem.
Lemma A1.
{ X , ξ ( u ) , u U } is a Gaussian system, where D 2 ξ ( u ) > 0 for every u U . Then, the following statements are equivalent:
(a) 
E ( e X | ξ ( u ) , u U ) = 1 and
(b) 
X and ξ ( u ) are independent for all u U and E X + 1 2 D 2 X = 0 .
Proof of Lemma A1. 
The equivalence is proven by deriving both directions from each other.
  • ( a ) ( b )
    First, we prove the direction of the (a) to (b). If E ( e X | ξ ( u ) , u U ) = e 0 = 1 is true, then it can be stated that
    e E ( X | ξ ( v ) ) + 1 2 D 2 ( X | ξ ( v ) ) = E ( e X | ξ ( v ) ) = E ( E ( e X | ξ ( u ) , u U ) | ξ ( v ) ) = 1
    Thus,
    0 = E ( X | ξ ( v ) ) + 1 2 D 2 ( X | ξ ( v ) )
    = E X + corr ( X , ξ ( v ) ) D X D ξ ( v ) ( ξ ( v ) E ξ ( v ) ) + 1 2 D 2 X ( 1 corr ( X , ξ ( v ) ) 2 )
    Therefore, we can conclude that corr ( X , ξ ( v ) ) = 0 and E X + 1 2 D 2 X = 0 .
  • ( b ) ( a )
    Then, we deduce that if (b) is fulfilled, then statement (a) is also true. If X and ξ ( u ) are independent for all u U and E X + 1 2 D 2 X = 0 , then
    E ( e X | ξ ( u ) , u U ) = E ( e X ) = e E X + 1 2 D 2 X = 1
Proof of Theorem 1. 
It is obvious that R ( u ) , where 0 u s , and P ( s , t ) are F ( s ) measurable random variables, just as Z ( s , t ) . Therefore, it can be stated that there is no problem with the existence of expected values. The equivalence of the statements is proved circularly.
  • ( a ) ( b )
    Let us start with the statement that the discounted bond price is a martingale. Hence,
    Z ( s , t ) = E [ Z ( t , t ) | F ( s ) ] = E [ e 0 t R ( u ) d u | F ( s ) ]
    P ( s , t ) = e 0 s R ( u ) d u Z ( s , t ) = E [ e s t R ( u ) d u | F ( s ) ]
    From statement (a), we quickly deduced that the discount factor occurs in the given form.
  • ( b ) ( c )
    Henceforth, we derive the drift term from the discount factor
    E e s t ( F ( s , u ) R ( u ) ) d u | F ( s ) = 1
    According to Lemma A1, this is equivalent to the fact that ξ ( s , t ) and F ( v 1 , v 2 ) are independent and E ξ ( s , t ) + 1 2 D 2 ξ ( s , t ) = 0 , where v 1 s , v 1 v 2 , ξ ( s , t ) = s t ( F ( s , u ) R ( u ) ) d u . Since we are dealing with Gaussian variables, it is sufficient to examine the covariance.
    cov ( F ( s , u ) R ( u ) , F ( v 1 , v 2 ) ) = c ( s v 1 , u , v 2 ) c ( u v 1 , u , v 2 ) = c ( v 1 , u , v 2 ) c ( v 1 , u , v 2 ) = 0
    Since ξ ( s , t ) is equal to s t ( F ( s , u ) R ( u ) ) d u , therefore ξ ( s , t ) and F ( v 1 , v 2 ) are independent.
    E ξ ( s , t ) = s t ( α ( s , u ) α ( u , u ) ) d u
    The variance is a bit more complicated to calculate.
    D 2 ξ ( s , t ) = s t s t c ( u v , u , v ) d v d u + s t s t c ( s , u , v ) d v d u 2 s t s t c ( s , u , v ) d v d u
    = s t s t ( c ( u v , u , v ) c ( s , u , v ) ) d v d u
    Let us apply the Leibniz integral rule to the following function.
    f ( t ) = E ξ ( s , t ) + 1 2 D 2 ξ ( s , t ) = s t ( α ( s , u ) α ( u , u ) ) d u + 1 2 s t s t ( c ( u v , u , v ) c ( s , u , v ) ) d v d u
    d f d t = α ( s , t ) α ( t , t ) + 1 2 s t ( c ( t v , t , v ) c ( s , t , v ) ) d v + s t ( c ( u t , u , t ) c ( s , u , t ) ) d u ,
    from which
    d f d t = α ( s , t ) α ( t , t ) + s t ( c ( v , v , t ) c ( s , v , t ) ) d v .
    Since f ( s ) = 0 , thus
    E ξ ( s , t ) + 1 2 D 2 ξ ( s , t ) = 0 , 0 s t α ( s , t ) = α ( t , t ) + s t [ c ( s , v , t ) c ( v , v , t ) ] d v ,
    for all 0 s t . Using Remark 1, we have established the implication from (b) to (c).
  • ( c ) ( a )
    First, we show that part (b) of the theorem is satisfied by showing that the drift term has the form (c), and this is sufficient because part (b) immediately demonstrates that Z ( s , t ) is a regular martingale. It can be easily seen that in Lemma A1, using the previous notations, ξ ( s , t ) = s t ( F ( s , u ) R ( u ) ) d u and F ( v 1 , v 2 ) are independent. During the derivations, Remark 1 is also used.
    E ξ ( s , t ) = s t α ( s , u ) α ( u , u ) d u = s t α ( u , u ) s u c ( s , v , u ) c ( v , v , u ) d v α ( u , u ) d u
    = s t s u c ( s , v , u ) c ( v , v , u ) d v d u
    D 2 ξ ( s , t ) = s t s t ( c ( u v , u , v ) c ( s , u , v ) ) d v d u
    = s t s u ( c ( v , u , v ) c ( s , u , v ) ) d v d u + s t u t ( c ( u , u , v ) c ( s , u , v ) ) d v d u
Let us apply the Leibniz rule again for the following f ( t ) function using the fact that the covariance function c ( s 1 s 2 , t 1 , t 2 ) is symmetric in t 1 and t 2 .
f ( t ) = E ξ ( s , t ) + 1 2 D 2 ξ ( s , t ) = s t s u c ( s , v , u ) c ( v , v , u ) d v d u
+ 1 2 s t s u ( c ( v , u , v ) c ( s , u , v ) ) d v d u + 1 2 s t u t ( c ( u , u , v ) c ( s , u , v ) ) d v d u
= 1 2 s t s u c ( s , u , v ) c ( v , u , v ) d v d u + 1 2 s t u t ( c ( u , u , v ) c ( s , u , v ) ) d v d u
= 1 2 s t s t c ( s , u , v ) c ( v , u , v ) d v d u 1 2 s t s t ( c ( u , u , v ) c ( s , u , v ) ) d v d u
= 1 2 s t s t c ( s , u , v ) c ( v , u , v ) c ( u , u , v ) + c ( s , u , v ) d v d u = 0
Finally, the theorem is proved. □

Appendix A.2

This subsection calculates the expected value and standard deviation of the expressions previously marked with ξ ( s , t ) and η ( s , t ) and their correlation.
ξ ( s , t ) = s t F ( u , u ) d u = s t r ( u ) d u
μ 1 ( s , t ) = E ξ ( s , t ) = E s t F ( u , u ) d u = s t E F ( u , u ) d u
= s t ν σ 2 μ + σ 2 μ e μ ( u u ) + σ 2 λ μ e μ ( u u ) σ 2 λ μ e λ ( u u ) d u
= s t ν σ 2 μ d u = ν σ 2 μ ( t s )
Now, let us move on to the expected value of η ( s , t ) . The steps of the derivation are similar to what we have seen above.
η ( s , t ) = s t F ( s , u ) d u
μ 2 ( s , t ) = E η ( s , t ) = E s t F ( s , u ) d u = s t E F ( s , u ) d u
= s t ν σ 2 μ + σ 2 μ e μ ( u s ) + σ 2 λ μ e μ ( u s ) σ 2 λ μ e λ ( u s ) d u
= s t ν σ 2 μ + σ 2 μ e μ u + μ s + σ 2 λ μ e μ u + μ s σ 2 λ μ e λ u + λ s d u
= ν u σ 2 μ u σ 2 μ 2 e μ u + μ s σ 2 μ ( λ μ ) e μ u + μ s + σ 2 λ ( λ μ ) e λ u + λ s u = s u = t
= ν σ 2 μ ( t s ) σ 2 μ 2 e μ ( t s ) 1 σ 2 μ ( λ μ ) e μ ( t s ) 1 + σ 2 λ ( λ μ ) e λ ( t s ) 1
After the derivation of the expected values, the standard deviation of the expressions marked with ξ ( s , t ) and η ( s , t ) are calculated.
c o v [ F ( u , u ) , F ( v , v ) ] = σ 2 exp { 2 μ min ( u , v ) μ ( u + v ) } = σ 2 e μ | v u |
D 2 ξ ( s , t ) = s t s t σ 2 e μ | v u | d u d v
At first, let us deal with the inner integral, then move on to the outer integral.
s t e μ | v u | d u = s v e μ u e μ v d u + v t e μ v e μ u d u = e μ v e μ u μ u = s u = v + e μ v e μ u μ u = v u = t
= 1 μ 1 e μ ( v s ) e μ ( t v ) + 1 = 1 μ 2 e μ ( v s ) e μ ( t v )
σ 2 μ s t 2 e μ ( v s ) e μ ( t v ) d v = σ 2 μ 2 ( t s ) e μ s e μ v μ v = s v = t e μ t e μ v μ v = s v = t
= σ 2 μ 2 ( t s ) + 1 μ ( e μ ( t s ) 1 ) 1 μ ( 1 e μ ( t s ) ) = σ 2 μ 2 2 μ ( t s ) + 2 e μ ( t s ) 2
D 2 ξ ( s , t ) = 2 σ 2 μ 2 μ ( t s ) + e μ ( t s ) 1
The variance of η ( s , t ) can be similarly derived.
c o v [ F ( s , u ) , F ( s , v ) ] = σ 2 exp { λ s + ( 2 μ λ ) min ( u , v ) μ ( u + v ) }
D 2 η ( s , t ) = s t s t σ 2 e λ s + ( 2 μ λ ) min ( u , v ) μ ( u + v ) d u d v
Similarly to the previous calculation, the derivation starts by calculating the inner integral without the σ 2 multiplier.
s t e λ s + ( 2 μ λ ) min ( u , v ) μ ( u + v ) d u = s v e λ s + ( 2 μ λ ) u μ u μ v d u + v t e λ s + ( 2 μ λ ) v μ u μ v d u
= s v e λ s + ( μ λ ) u μ v d u + v t e λ s + ( μ λ ) v μ u d u = e λ s μ v s v e ( μ λ ) u d u + e λ s v t e μ v λ v μ u d u
= e λ s μ v e ( μ λ ) u μ λ u = s u = v + e λ s + μ v λ v e μ u μ u = v u = t
= e λ s μ v 1 μ λ e μ v λ v e μ s λ s + e λ s + μ v λ v 1 μ e μ v e μ t
= 1 μ λ e λ s λ v 1 μ λ e μ s μ v 1 μ e λ s + ( μ λ ) v μ t + 1 μ e λ s λ v
Now, let us move to the outer integral per term.
= 1 μ λ e λ s s t e λ v d v = e λ s μ λ e λ v λ v = s t = e λ s ( λ μ ) λ e λ t e λ s
= 1 ( λ μ ) λ e λ ( t s ) 1
= 1 μ λ e μ s s t e μ v d v = 1 μ λ e μ s e μ v μ u = s t = 1 μ ( μ λ ) e μ s e μ t e μ s
= 1 μ ( μ λ ) e μ ( t s ) 1
= s t 1 μ e λ s μ t e ( μ λ ) v d v = 1 μ ( μ λ ) e λ s μ t e ( μ λ ) s e ( μ λ ) t
= 1 μ ( μ λ ) e μ ( t s ) e λ ( t s )
= s t 1 μ e λ s λ v d v = 1 μ e λ s s t e λ v d v = 1 λ μ e λ s e λ t e λ s = 1 λ μ 1 e λ ( t s )
Therefore, by adding the σ 2 multiplier, we obtain the variance of η ( s , t ) .
D 2 η ( s , t ) = σ 2 ( λ μ ) λ e λ ( t s ) 1 + σ 2 μ ( μ λ ) 2 e μ ( t s ) e λ ( t s ) 1 + σ 2 λ μ 1 e λ ( t s )
The last variable to be calculated is ρ ( s 1 , s 2 , t 1 , t 2 ) , indicating the correlation between ξ ( s 1 , t 1 ) = s 1 t 1 F ( u , u ) d u and η ( s 2 , t 2 ) = s 2 t 2 F ( s 2 , v ) d v . However, for all financial products used in the article, the values of t 1 and t 2 were equal, denoted by t. Furthermore, we can assume that s 1 s 2 , since s 1 represents the time at which we discount, the time at which we want the value of the financial product, while s 2 is the starting time of the transaction, which can start now or even later. Thus, let us suppose that s 1 < s 2 < t .
c o v [ F ( u , u ) , F ( s 2 , v ) ] = σ 2 exp { λ min ( u , s 2 ) + ( 2 μ λ ) min ( u , v ) μ ( u + v ) }
c o v ( ξ ( s 1 , t ) , η ( s 2 , t ) ) = s 1 t s 2 t σ 2 exp { λ min ( u , s 2 ) + ( 2 μ λ ) min ( u , v ) μ ( u + v ) } d v d u
= s 1 s 2 s 2 t σ 2 exp { λ u + ( 2 μ λ ) u μ ( u + v ) } d v d u
+ s 2 t s 2 t σ 2 exp { λ t + ( 2 μ λ ) min ( u , v ) μ ( u + v ) } d v d u
The two terms of the summation are calculated separately.
= σ 2 s 1 s 2 s 2 t exp { μ u μ v } d v d u = σ 2 s 1 s 2 e μ u d u s 2 t e μ v d v
= σ 2 e μ u μ u = s 1 s 2 e μ v μ v = s 2 t = σ 2 μ 2 e μ s 2 e μ s 1 e μ s 2 e μ t
= σ 2 μ 2 1 e μ ( t s 2 ) e μ ( s 2 s 1 ) + e μ ( t s 1 )
= σ 2 s 2 t s 2 t exp { λ t + ( 2 μ λ ) min ( u , v ) μ ( u + v ) } d v d u
= σ 2 s 2 t s 2 u exp { λ t + ( 2 μ λ ) min ( u , v ) μ ( u + v ) } d v d u
+ σ 2 s 2 t u t exp { λ t + ( 2 μ λ ) min ( u , v ) μ ( u + v ) } d v d u
= σ 2 e λ t s 2 t e μ u s 2 u e μ v λ v d v d u + σ 2 e λ t s 2 t e μ u λ u u t e μ v d v d u
= σ 2 e λ t s 2 t e μ u e ( μ λ ) v μ λ v = s 2 u d u + σ 2 e λ t s 2 t e ( μ λ ) u e μ v μ v = u t d u
= σ 2 μ λ e λ t s 2 t e μ u e ( μ λ ) u e ( μ λ ) s 2 d u + σ 2 μ e λ t s 2 t e ( μ λ ) u e μ u e μ t d u
= σ 2 μ λ s 2 t e λ t λ u e λ ( t s 2 ) μ ( u s 2 ) d u + σ 2 μ s 2 t e λ ( t u ) e μ ( t u ) λ ( t u ) d u
= σ 2 μ λ e λ t e λ u λ u = s 2 t e μ u + μ s 2 λ s 2 μ u = s 2 t + σ 2 μ e λ t e λ u λ u = s 2 t e μ u λ u μ t μ λ u = s 2 t
= σ 2 λ ( μ λ ) e λ t e λ s 2 e λ t + σ 2 μ ( μ λ ) e λ t e μ t + μ s 2 λ s 2 e μ s 2 + μ s 2 λ s 2
+ σ 2 λ μ e λ t e λ s 2 e λ t + σ 2 μ ( μ λ ) e λ t e ( μ λ ) s 2 μ t e ( μ λ ) t μ t
= σ 2 λ ( μ λ ) e λ ( t s 2 ) 1 + σ 2 μ ( μ λ ) e λ ( t s 2 ) μ ( t s 2 ) 1 + σ 2 λ μ e λ ( t s 2 ) 1 + σ 2 μ ( μ λ ) e λ ( t s 2 ) μ ( t s 2 ) 1
= σ 2 λ ( μ λ ) + σ 2 λ μ e λ ( t s 2 ) 1 + 2 σ 2 μ ( μ λ ) e λ ( t s 2 ) μ ( t s 2 ) 1
By adding the calculated two terms back together, we obtain the covariance between ξ ( s 1 , t ) and η ( s 2 , t ) .
c o v ( ξ ( s 1 , t ) , η ( s 2 , t ) ) = σ 2 μ 2 1 e μ ( t s 2 ) e μ ( s 2 s 1 ) + e μ ( t s 1 )
+ σ 2 λ ( μ λ ) + σ 2 λ μ e λ ( t s 2 ) 1 + 2 σ 2 μ ( μ λ ) e λ ( t s 2 ) μ ( t s 2 ) 1
Therefore, the correlation between ξ and η is calculated as follows:
c o r r ( ξ , η ) = c o v ( ξ , η ) D ξ D η

Appendix A.3

In this subsection of the appendix, the analytical fair price of the European floorlet option is also derived, similarly to the previously derived European caplet option. As we have seen before, the fair price of the floorlet is the expected value of the payoff function under the risk-neutral measure.
P f l o o r l e t ( s ) = E [ e s t + Δ r ( u ) d u ( e Δ K e Δ F Δ ( t , t ) ) + ]
We use the previously introduced variables: ξ and η .
ξ ( s , t + Δ ) = s t + Δ r ( u ) d u = s t + Δ F ( u , u ) d u
η ( t , t + Δ ) = Δ F Δ ( t , t ) = t t + Δ F ( t , u ) d u
Similarly to the previously derived pricing formula for the caplet, the conditional expected value of ξ to η follows a normal distribution. Hence, the conditional standard distribution theorem can also be used with the previously defined parameters in this case. Therefore, the fair price of the European floorlet can be calculated as follows:
E [ e s t + Δ r ( u ) d u ( e Δ K e Δ F Δ ( t , t ) ) + ] = E [ e ξ ( e Δ K e η ) + ]
= E [ E ( e ξ ( e Δ K e η ) + | η ] = E [ ( e Δ K e η ) + · E ( e ξ | η ) ]
During the derivations, the law of total expectation and the fact that ( e Δ K e η ) + is measurable for η is used.
As we can see, ξ N ( μ 1 , σ 1 ) is normally distributed, therefore ξ N ( μ 1 , σ 1 ) , where c o r r ( ξ , η ) = ρ . Therefore, ξ | η N ( μ 1 ρ σ 1 η μ 2 σ 2 , σ 1 2 ( 1 ρ 2 ) ) . Since the conditional distribution of ξ given η is known, E [ e ξ | η ] can be calculated as the expected value of a lognormal distribution.
E [ e ξ | η ] = e μ 1 ρ σ 1 η μ 2 σ 2 + 1 2 σ 1 2 ( 1 ρ 2 )
The integral returns the expected value of a random variable that is lognormally distributed. Returning to the pricing formula
E [ ( e Δ K e η ) + · E [ e ξ | η ] ] = E ( e Δ K e η ) + · e μ 1 ρ σ 1 η μ 2 σ 2 + 1 2 σ 1 2 ( 1 ρ 2 )
= e μ 1 + 1 2 σ 1 2 ( 1 ρ 2 ) E ( e Δ K e η ) + · e ρ σ 1 η μ 2 σ 2
= e Δ K + 1 2 ρ 2 σ 1 2 Δ K 1 σ 2 2 π e ( x ( μ 2 ρ σ 1 σ 2 ) ) 2 2 σ 2 2 d x
e μ 2 + ( σ 2 2 ρ σ 1 σ 2 ) 2 2 σ 2 2 Δ K 1 σ 2 2 π e ( x ( μ 2 + σ 2 2 ρ σ 1 σ 2 ) ) 2 2 σ 2 2 d x
Therefore, the analytical pricing formula for the European floorlet option in Kennedy fields is as follows:
P f l o o r l e t ( s ) = e Δ K μ 1 + 1 2 σ 1 2 Φ Δ K μ 2 + ρ σ 1 σ 2 σ 2 e μ 2 μ 1 + 1 2 ( σ 1 2 + σ 2 2 2 ρ σ 1 σ 2 ) Φ Δ K μ 2 σ 2 2 + ρ σ 1 σ 2 σ 2

Appendix A.4

The fair price of a fixed vs. floating swap for several periods (k) at time s can be written using the formula below. Let us denote the time periods by s < T 0 < T 1 < < T k and τ j = T j T j 1 .
P s w a p ( s ) = E j = 1 k e s T j r ( u ) d u ( e Δ F Δ ( T j 1 , T j 1 ) e τ j K )
= E j = 1 k e s T j r ( u ) d u + Δ F Δ ( T j 1 , T j 1 ) e s T j r ( u ) d u + τ j K
= j = 1 k E [ e s T j r ( u ) d u + Δ F Δ ( T j 1 , T j 1 ) ] j = 1 k E [ e s T j r ( u ) d u + τ j K ]
= j = 1 k E [ e s T j r ( u ) d u + Δ F Δ ( T j 1 , T j 1 ) ] j = 1 k e τ j K E [ e s T j r ( u ) d u ]
As we have done previously, more additional variables are introduced, ξ j and η j , in the following way.
ξ j ( s , T j ) = s T j r ( u ) d u = s T j F ( u , u ) d u
η j ( T j 1 , T j ) = Δ F Δ ( T j 1 , T j 1 ) = T j 1 T j F ( T j 1 , u ) d u
Referring to previous calculations, we know the expected value and standard deviation of ξ and η . Therefore, the expected values of the variables are the following:
μ ξ j = E ξ j ( s , T j ) = E s T j F ( u , u ) d u = ν σ 2 μ ( T j s )
μ η j = E η j ( T j 1 , T j ) = E Δ F Δ ( T j 1 , T j 1 ) = T j 1 T j E F ( T j 1 , u ) d u
= ν σ 2 μ τ j σ 2 μ 2 e μ τ j 1 σ 2 μ ( λ μ ) e μ τ j 1 + σ 2 λ ( λ μ ) e λ τ j 1
and the covariance is
σ ξ j 2 = D 2 ξ j ( s , T j ) = 2 σ 2 μ 2 ( T j s ) μ + e μ ( T j s ) 1
σ η j 2 = D 2 η j ( T j 1 , T j ) = σ 2 ( λ μ ) λ e λ τ j 1 + σ 2 μ ( μ λ ) e μ τ j 1
+ σ 2 μ ( μ λ ) e μ τ j e λ τ j + σ 2 λ μ 1 e λ τ j
Due to the properties of the Gaussian random field, ( ξ j , η j ) always follows a multivariate normal distribution, with the covariance matrix shown before.
c o v ( ξ j , η j ) = σ 2 μ 2 1 e μ ( T j T j 1 ) e μ ( T j 1 s ) + e μ ( T j s )
+ σ 2 λ ( μ λ ) + σ 2 λ μ e λ ( T j T j 1 ) 1 + 2 σ 2 μ ( μ λ ) e λ ( T j T j 1 ) μ ( T j T j 1 ) 1
The correlation between ξ j and η j is the value of the covariance normalized with the standard deviation of ξ j and η j .
c o r r ( ξ j , η j ) = c o v ( ξ j , η j ) D ξ j D η j
Because of the properties of the normal distribution, the distribution of ξ j + η j is also normally distributed where the mean of the convolution is the sum of the means, and the variance is the following.
ξ j + η j N ( μ ξ j + μ η j , σ ξ j 2 + σ η j 2 2 ρ σ ξ j σ η j )
Therefore, the price of the interest rate swap can be easily calculated:
P s w a p ( s ) = j = 1 k E ( e ξ j + η j ) j = 1 k e τ j K E ( e ξ j )
= j = 1 k e μ ξ j + μ η j + 1 2 ( σ ξ j 2 + σ η j 2 2 ρ σ ξ j σ η j ) j = 1 k e τ j K e μ ξ j + 1 2 σ ξ j 2
= j = 1 k e μ ξ j + μ η j + 1 2 ( σ ξ j 2 + σ η j 2 2 ρ σ ξ j σ η j ) j = 1 k e τ j K μ ξ j + 1 2 σ ξ j 2

References

  1. Kennedy, D.P. The term structure of interest rates as a Gaussian random field. Math. Financ. 1994, 4, 247–258. Available online: https://ideas.repec.org/a/bla/mathfi/v4y1994i3p247-258.html (accessed on 15 September 2021). [CrossRef]
  2. Kennedy, D.P. Characterizing Gaussian models of the term structure of interest rates. Math. Financ. 1997, 7, 107–116. Available online: https://ideas.repec.org/a/bla/mathfi/v7y1997i2p107-118.html (accessed on 15 September 2021). [CrossRef]
  3. Heath, D.C.; Jarrow, R.A.; Morton, A. Bond pricing and term structure of interest rates: A new methodology for contigent claims valuation. Econometrica 1992, 60, 77–105. [Google Scholar] [CrossRef]
  4. Shreve, S.E. Stochastic Calculus for Finance I-II, 1st ed.; Springer Finance: Pittsburg, PA, USA, 2004. [Google Scholar]
  5. Cheyette, O. Markov representation of the Heath-Jarrow-Morton model. SSRN Electron. J. 2001. Available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=6073 (accessed on 12 February 2024). [CrossRef]
  6. Beyna, I.; Wystup, U. On the Calibration of the Cheyette Interest Rate Model; CPQF Working Paper Series No. 25; Frankfurt School of Finance and Management: Frankfurt am Main, Germany, 2010. [Google Scholar]
  7. Arató, N.M. Mean estimation of Brownian sheet. Comput. Math. Appl. 1997, 33, 12–25. Available online: https://core.ac.uk/download/pdf/82766418.pdf (accessed on 3 November 2021). [CrossRef]
  8. Rozanov, J. Infinite-dimensional Gaussian distributions: Proceedings of the Steklov Institute of Mathematics, No. 108 (1968), 3rd ed.; American Mathematical Society: Providence, Rhode Island, 1971. [Google Scholar]
  9. Arató, M.; Tóth-Lakits, D. Modeling Negative Rates. In Contributions to Risk Analysis: RISK 2022; Fundacion MAPFRE: Madrid, Spain, 2022; pp. 251–260. [Google Scholar]
  10. Grenander, U. Stochastic Processes and Statistical Inference. Ark. Mat. 1950, 1, 195–277. Available online: https://archive.ymsc.tsinghua.edu.cn/pacm_download/116/6848-11512_2007_Article_BF02590638.pdf (accessed on 2 February 2022). [CrossRef]
  11. Emerick, J.; Tatsat, H. Stochastic Volatility Models—Heston Model Calibration to option prices. QuantPy 2022. Available online: https://quantpy.com.au/stochastic-volatility-models/heston-model-calibration-to-option-prices/ (accessed on 10 January 2023).
  12. Mikhailov, S.; Nögel, U. Heston’s Stochastic Volatility Model Implementation, Calibration and Some Extensions. Fraunhofer Inst. Ind. Math. 2003, 74–79. Available online: https://www.maths.univ-evry.fr/pages_perso/crepey/Equities/051111_mikh%20heston.pdf (accessed on 15 January 2023).
Figure 1. Simulated Kennedy fields
Figure 1. Simulated Kennedy fields
Mathematics 12 03059 g001
Figure 2. Simulated Monte Carlo market prices (mesh) vs. calibrated Kennedy model prices (markers) for a caplet.
Figure 2. Simulated Monte Carlo market prices (mesh) vs. calibrated Kennedy model prices (markers) for a caplet.
Mathematics 12 03059 g002
Figure 3. Par swap rate market prices (mesh) vs. calibrated Kennedy model prices (markers).
Figure 3. Par swap rate market prices (mesh) vs. calibrated Kennedy model prices (markers).
Mathematics 12 03059 g003
Figure 4. Historical parameter estimations in time.
Figure 4. Historical parameter estimations in time.
Mathematics 12 03059 g004
Figure 5. Differently parameterized Kennedy fields for describing forward rates.
Figure 5. Differently parameterized Kennedy fields for describing forward rates.
Mathematics 12 03059 g005
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Tóth-Lakits, D.; Arató, M. On the Calibration of the Kennedy Model. Mathematics 2024, 12, 3059. https://doi.org/10.3390/math12193059

AMA Style

Tóth-Lakits D, Arató M. On the Calibration of the Kennedy Model. Mathematics. 2024; 12(19):3059. https://doi.org/10.3390/math12193059

Chicago/Turabian Style

Tóth-Lakits, Dalma, and Miklós Arató. 2024. "On the Calibration of the Kennedy Model" Mathematics 12, no. 19: 3059. https://doi.org/10.3390/math12193059

APA Style

Tóth-Lakits, D., & Arató, M. (2024). On the Calibration of the Kennedy Model. Mathematics, 12(19), 3059. https://doi.org/10.3390/math12193059

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop