Calibration of the Gaussian Musiela model using the Karhunen-Loeve expansion

Tiziano Vargiolu Scuola Normale Superiore Piazza dei Cavalieri 7 - 56100 Pisa (PI)

                      email: vargiolu@cibs.sns.it

May 12, 2000

Abstract

In this paper we calibrate the stationary Gaussian Musiela model to time series of market data using the Karhunen-Loeve expansion in order to get an ortonormal basis (classically known as EOF, empirical orthonormal functions) in a separable Hilbert space. The basis found is optimal for representing the covariance of the invariant measure of the forward rates’ process.

1 Introduction

In this paper we apply the Karhunen-Loeve expansion in order to estimate from time series the functional coefficients of the stationary Gaussian Musiela model for the term structure of interest rates. The core of this model is the stochastic partial differential equation

\{ \begin{matrix}dr_{t}(x) = (\frac{∂}{∂x}r_{t}(x) + τ^{∗}(x) ∫_{0}^{x} τ(u) du) dt + τ^{∗}(x) dW_{t} \\ r_{0} ∈ W_{loc}^{1,1}(ℝ^{+}) \end{matrix}

where W is a k-dimensional standard Brownian motion, τ is a deterministic function with values in ℝ^{k}, and the symbol ∗ indicates the transposition in ℝ^{k}. In the equation, r_{t}(x) represents the spot forward rate, that is the instantaneous rate prevailing in t for a forward contract at time t + x. We study the case when the process takes values in the separable Hilbert space H_{γ}^{1} (that is the Sobolev space of functions in ℝ^{+} that, together with their first weak derivative, are square integrable with respect to the measure e^{−γx} dx) and we present some results already present in literature in [17] about existence and uniqueness of mild solutions of the equation and existence of Gaussian invariant measures N(b(·), Q_{∞}) in H_{γ}^{1}. Our aim is to find an orthonormal basis (e_{n})_{n} in H_{γ}^{1} and real numbers (μ_{n})_{n} such that if we substitute the (τ_{n})_{n} with (\sqrt{μ_{n}}e_{n})_{n} we make the least possible error. In order to do this, we present the Karhunen-Loeve expansion that is used to diagonalize the operator Q_{∞} in the form:

Q_{∞} = ∑_{n=1}^{∞} λ_{n}k_{n} ⊗ k_{n}

where λ_{n} ↘ 0 and (k_{n})_{n} is an orthonormal basis in H_{γ}^{1}. Following [12], we present an iterative method to implement the Karhunen-Loeve expansion and we apply it to the Musiela model, finding also a stopping criterion to find only the significant components of Q_{∞}. Unfortunately, the problem of passing from the (k_{n})_{n} to the (τ_{n})_{n} (that is to the (e_{n})_{n}) is rather difficult, so we present another approach: we prove that Q_{∞} satisfies the operator relation AQ_{∞} + Q_{∞}A^{∗} + ττ^{∗} = 0, where A = \frac{∂}{∂x} and ττ^{∗} = ∑_{n=1}^{∞} τ_{n} ⊗ τ_{n}. This implies that, under the assumption that Tr ττ^{∗} < +∞, we can apply Karhunen-Loeve expansion and the itherative algorithm to ττ^{∗} and get the (e_{n})_{n} we were looking for.

The Musiela model is a reparametrization of the better known model of Heath-Jarrow-Morton (HJM) (see [8]), in which f(t, T) represents the forward rate prevailing in t for the time T. In other words, we have that f(t, t + x) = r(t, x). The Musiela parametrization allows to consider the forward curve r(t, ·) as a Markov process in a suitable function space (while in general in the HJM model this is not possible) and is coherent with other forward rates models (see [2], [14] and [16]). There are two possible ways of solving the equation and of using the Musiela model. The first one is to solve it in its generality. This means solving a stochastic differential equation in a locally convex space. This is a rather difficult task for which there are very few existing references in literature. The second one is to choose a particular separable Hilbert space H ⊆ W_{loc}^{1,1}(ℝ^{+}) and to solve the equation in H. This is much easier and one can find many works in literature. In particular, throughout this work we refer to the book [5] which gives a nearly complete treatment of stochastic integration and stochastic differential equations in separable Hilbert spaces, even if there are other possible approaches (see for example [13] or [18]). This strategy leads us to two questions: can one find a separable Hilbert space that is the “right” one for the equation? Second, does there exists an orthonormal basis (i.e. a complete orthonormal system) in this space which minimizes the error when one approximates the dynamics with the first n components of this basis? The answer to the first question seems very arbitrary and depends on subjective tastes: it seems natural to consider some Sobolev space with at least the first derivative defined in some sense, but some choices (finite or infinite time horizon, weight to put in this horizon) remain. Once the space H is fixed, the second question is very interesting. From a theoretical point of view, it gives a good interpretation of the behaviour of r and from a pratical point of view one can choose to approximate the dynamics of r (which a priori is infinite dimensional) with the first n components of the “optimal” orthonormal basis, thus making the least possible error. In this way it would be possible to implement and calibrate the model efficiently. To this purpose, we notice that the usual procedure in market practise and in the literature is to fix a priori both the number of the sources of noise (i.e. the dimension of the Brownian motion) as the shape of the functional parameters τ_{i}, i = 1, …, n, regardless of the fact that these choices are efficient (see for example [1], [2], [3] for some ways to do this). The idea of finding an “efficient” basis appears for the first time (to our knowledge) in [6]: in that work the existence of optimal (τ_{n})_{n} to calibrate the model is proved in a generic separable Hilbert space H of continuous functions, but no operative way to obtain these (τ_{n})_{n} is presented. In this work, in order to calibrate the model, we chose to use time series because it gives much more information than an implied calibration does (used for example in [2] and [3]). In particular, we apply the method known as “Karhunen-Loeve expansion”. This method is classically used in fluid dynamics in order to find an orthonormal basis, whose vectors are known as “empirical orthonormal functions”, which diagonalize the covariance operator R of an infinite dimensional random variable, which is a compact selfadjoint operator. It turns out that this method gives good numerical results, in the sense that only a few components of the basis are sufficient to describe much of the evolution of the system, and there doesn’t exist another orthonormal basis that can improve the result. We present the method following [12] and we present an iterative algorithm to obtain the eigenvalues and the eigenvectors of R as a possible estimation of the number N of components necessary to describe sufficiently R. Then we present a way to apply the method to the Musiela model, considering a time series as a sample of the infinite dimensional random variable r_{t} having as law an invariant measure of the equation. Here there is a hidden detail: we always studied the Musiela model under the risk-neutral probability ℚ, while we can obtain from the time series an estimation of the historical probability ℙ. This doesn’t create problems: since ℚ and ℙ are equivalent, the Gaussian random variable we are studying has the same covariance under both ℙ and ℚ, so we can proceed safely. We show the steps to implement the model numerically, and we give an idea of how to pass from the covariance of the invariant measure

R = ∑_{n=1}^{∞} ∫_{0}^{∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} du = ∑_{n=1}^{∞} λ_{n}k_{n} ⊗ k_{n}

to the (τ_{n})_{n}.

The work is organized as follows. In section 2 the Gaussian Musiela model used in this work is presented. In section 3 some existing results on the Musiela equation are cited, namely about the existence and uniqueness of the solution, and the existence of Gaussian invariant measures for the forward rate curve. In section 4 we present Italian market of bonds, and we give a method to get samples of the forward curve starting from time series based on Italian bonds. In section 5 we present a theoretical framework for Karhunen-Loeve expansion. In section 6 we present an iterative method used in order to get the eigenvalues and the eigenvectors of the covariance operator in an efficient way. In section 7 we apply the results to the Musiela model.

The author wishes to thank Nicole El Karoui for having given him the idea of finding the “right” ortonormal basis for the Hilbert space used in the Musiela model that inspired his work, and Franco Flandoli for his invaluable advice and suggestion to use Karhunen-Loeve expansion.

2 Presentation of the Gaussian Musiela model

Now let’s introduce the Musiela model for the term structure of interest rates. We have a probability space (Ω, F, ℙ), endowed with the filtration (F_{t})_{t≥0}. We suppose that the price at time t of a bond expiring at time T is given by the process (B(t, T))_{t∈[0,T]}, and we suppose to have a random field (r(t, x))_{t,x≥0}, in which r(t, x) is called spot forward rate at time t for the maturity t + x and represents the rate at which at time t one can enter a forward contract at time t + x for a short (infinitesimal) period of time. Then the price at time t of a bond expiring at time T is given by:

B(t, T) = exp (− ∫_{0}^{T −t} r(t, u) du) .

We call spot rate the quantity r(t, 0); it represents the rate at which at time t one can enter a contract expiring immediately after. We also call price progress of the savings account the process (β(t))_{t}, given by:

β(t) = exp (∫_{0}^{t} r(u, 0) du) .

The actualized price at time t of a bond expiring at time T is given by:

B̃(t, T) = \frac{B(t, T)}{β(t)} = exp (∫_{0}^{t} r(u, 0) du − ∫_{0}^{T −t} r(t, u) du) .  (1)

Now we add the hypotheses that there exists a k-dimensional standard Brownian motion (W_{t})_{t≥0} adapted to (F_{t≥0})_{t}, and there exist two progressively measurable random fields (α(t, x))_{t,x≥0}, (τ(t, x))_{t,x≥0}, that can also depend on (r(t, x))_{t,x}, such that ∀x ≥ 0, (α(t, x))_{t} has trajectories in L^{1}(ℝ^{+}; ℝ), (τ^{∗}(t, x))_{t} has trajectories in L^{2}(ℝ^{+}; ℝ^{k}) and such that the forward rate satisfies the following stochastic differential equation:

\{ \begin{matrix}dr(t, x) = α(t, x) dt + τ^{∗}(t, x) dW_{t} \\ r(0, x) given. \end{matrix}  (2)

where τ^{∗}(t, x) dW_{t} stands for ∑_{i=1}^{k} τ_{i}(t, x) dW_{t}^{i}. A natural thing to ask is that the actualized bond price process (B̃(t, T))_{t∈[0,T]} is a martingale under a probability measure ℚ equivalent to ℙ, which is called risk-neutral probability or equivalent martingale measure. To this purpose, we cite the:

Theorem 1 If r(t, x) follows an equation of the form (2), then the following facts are equivalent:

1. ∀T > 0 the process (B̃(t, T))_{t∈[0,T]} is a local martingale with respect to ℚ.

2. (r(t, x))_{t,x} is such that ∀t ∈ ℝ^{+} the application x → r(t, x) ℚ-a.s. belongs to AC(ℝ^{+}) and Eq. (2) becomes

\{ \begin{matrix}dr_{t} = (Ar_{t} + c(t, ·)) dt + τ^{∗}(t, ·) dW_{t} \\ r_{0} ∈ AC(ℝ^{+}) \end{matrix}  (3)

where

A = \frac{∂}{∂x}, c(t, x) = ∑_{n=1}^{k} τ_{n}(t, x) ∫_{0}^{x} τ_{n}(t, u) du

If one of the conditions above is verified, then for all T > 0 the dynamics of B(t, T) is given by:

\{ \begin{matrix}dB(t, T) & = & B(t, T) (r(t, 0) dt + Γ(t, T) dW_{t}) ,  t ≤ T \\ B(T, T) & = & 1 \end{matrix}

where Γ(t, T) = − ∫_{0}^{T −t} τ(t, u) du.

Proof. See [15]. ⋄

If we make the further hypotheses that (F_{t})_{t} is the completion of the natural filtration of the Brownian motion (W_{t})_{t}, and that |τ_{t}(x)| ≤ M(T) ℙ-a.s. ∀t, x such that t + x ≤ T, then the process (B̃(t, T))_{t} is a martingale; then the hypotheses of the martingale representation theorem hold, and so every square integrable random variable (that is every contingent claim) can be represented as the sum of his expectation, which is the arbitrage free price of the claim, and of a stochastic integral with respect to (W_{t})_{t}. Besides, we can build a self-financing portfolio strategy based on B(·, T) and on β(·), which simulates the claim (see [15]).

Now we make some hypotheses under which this equation will have explicit solutions. We suppose that the process (τ_{t})_{t} is identically equal to a real k-valued deterministic function τ(x) belonging to W_{loc}^{1,1}(ℝ^{+}). Then the equation becomes:

\{ \begin{matrix}dr_{t} = (Ar_{t} + c) dt + τ^{∗} dW_{t} \\ r_{0} ∈ L^{2}(Ω, F_{0}, ℙ; H) \end{matrix}  (4)

This is a Langevin equation which, under right hypotheses, has an explicit solution.

We notice that the general model is in the space W_{loc}^{1,1}(ℝ^{+}), which is a locally convex space. In this kind of space, the theory of stochastic integration is rather difficult, while it becomes much simpler in Hilbert space. Though such a problem can be treated in its generality, the standard method is to choose a separable Hilbert space H ⊆ W_{loc}^{1,1}(ℝ^{+}) and to solve the equation in that space (see [3], [7], [17]). Natural choices for H can be H^{1}(ℝ^{+}) and H_{γ}^{1}(ℝ^{+}) with γ > 0, that is the Sobolev space of L^{2} functions with weak derivative also in L^{2}, and the same kind of space with respect to the exponential measure γe^{−γx} dx. The behaviour of the equation in the two kind of spaces is quite the same: in fact the formal expressions for the explicit solution and the law of the process are the same. However, the space H^{1} doesn’t contain constant functions and this, besides having a bad economic interpretation (a nonzero flat forward rate curve, which is the easiest interest rate model, is not allowed!), is the main reason why in H^{1} there is a unique invariant measure, while in H_{γ}^{1} there are infinitely many Gaussian ones (see [17]).

We now present some results of a study of the equation in H_{γ}^{1}, which is contained in [17], where we send the interested reader to. We use the theory of the stochastic integration in separable Hilbert spaces contained in [5] which is nearly complete.

3 The Musiela equation in the space H_{γ}^{1}

Theorem 2 . If τ_{i} ∈ H_{γ}^{1} ∩ H^{1} ∩ L_{γ}^{4} ∀n ∈ ℕ, ∑_{n=1}^{∞} ∥τ_{n}∥_{H_{γ}^{1}}^{2} < +∞, ∑_{n=1}^{∞} ∥τ_{n}∥_{H^{1}}^{4} < +∞, ∑_{n=1}^{∞} ∥τ_{n}∥_{L_{γ}^{4}}^{4} < +∞, and r_{0} ∈ L^{2}(Ω, F_{0}, ℚ; H_{γ}^{1}), then there exists a unique mild solution of Equation (4), given by

r_{t}(x) = r_{0}(x + t) + ∑_{i=1}^{k} ∫_{0}^{t} τ_{i}(x + t − u) (∫_{0}^{x+t−u} τ_{i}(v) dv) du +

+ ∫_{0}^{t} τ_{i}(x + t − u) dW_{u}^{i} =

= r_{0}(x + t) + \frac{1}{2} ∑_{i=1}^{k} ((∫_{0}^{x+t} τ_{i}(u) du)^{2} − (∫_{0}^{x} τ_{i}(u) du)^{2}) +

+ ∑_{i=1}^{k} ∫_{0}^{t} τ_{i}(x + t − u) dW_{u}^{i}  (5)

If r_{0} is a Gaussian random variable, then the solution is a Gaussian process with functional mean:

𝔼[r_{t}] = S_{t}𝔼[r_{0}] + ∫_{0}^{t} S_{t−u}c du  (6)

and functional covariance (if t ≤ v):

Cov (r_{t}, r_{v}) = S_{v}Cov (r_{0}, r_{0})S_{t}^{∗} + ∑_{i=1}^{k} ∫_{0}^{t} S_{t−u}τ_{i} ⊗ S_{v−u}^{∗}τ_{i} du  (7)

Since H_{γ}^{1} is a space of continuous functions, this means that ∀t, v, x, y ≥ 0:

𝔼[r_{t}(x)] = 𝔼[r_{0}(x + t)] + \frac{1}{2} ∑_{i=1}^{k} ((∫_{0}^{x+t} τ_{i}(u) du)^{2} − (∫_{0}^{x} τ_{i}(u) du)^{2})

Cov (r_{t}(x), r_{v}(y)) = Cov (r_{0}(t + x), r_{0}(v + y)) + ∑_{i=1}^{k} ∫_{0}^{t} τ_{i}(x + t − u)τ_{i}(y + v − u) du

Theorem 3 . Given Equation (4) in the Hilbert space H_{γ}^{1}, necessary and sufficient conditions to have an invariant measure are τ_{n} ∈ H_{γ}^{1}(ℝ^{+}) ∩ H^{1}(ℝ^{+}) ∩ L_{γ}^{4}(ℝ^{+}) ∀n ∈ ℕ, ∑_{n=1}^{∞} ∥τ_{n}∥_{H^{1}}^{2} < +∞, ∑_{n=1}^{∞} ∥τ_{n}∥_{L_{γ}^{4}}^{4} < +∞, and there exist infinitely many invariant measures. In particular, the measures of the kind:

δ_{b^{∗}(·)+b_{0}} ∗ N(0, Q_{∞}), b_{0} ∈ ℝ

are invariant measures, where:

Q_{∞} = ∑_{i=1}^{k} ∫_{0}^{+∞} τ_{i}(· + u) ⊗ τ_{i}(· + u) du  (8)

b^{∗}(x) = − ∑_{i=1}^{k} ∫_{0}^{x} τ_{i}(u) ∫_{0}^{u} τ_{i}(v) dv du = −\frac{1}{2} ∑_{i=1}^{k} (∫_{0}^{x} τ_{i}(u) du)^{2}  (9)

4 Forward rates in the Italian market

In this section we try to calibrate the model to the Italian market data. We chose to make the calibration using a historical approach for two reasons: the first is that in the Italian market there are very few derivative assets; the second is that a historical calibration can give much more information on the shape of the forward curve than an implicit calibration can give. In order to calibrate the model, we suppose that (r_{t})_{t} is a stationary solution of the Musiela equation (4) in the space H_{γ}^{1}. Then the marginal law of r_{t} under ℙ is

ℙ(r_{t}) = N(b, Q_{∞}) ∀t ≥ 0

with b ∈ H_{γ}^{1} and

Q_{∞} = ∫_{0}^{∞} ∑_{n=1}^{∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} du,

Since the process is stationary, in order to know its law it is sufficient to know Q_{∞}; in fact, once we know Q_{∞}, we are able to get the τ_{n} and to get c from them.

Here there is a hidden detail: the Musiela model is formulated under the risk-neutral probability ℚ, while from the time series we can obtain an estimation of the historical probability ℙ. This doesn’t create problems, because since ℚ and ℙ are equivalent, the Gaussian random variable we are studying has the same covariance both under ℙ as under ℚ, so we can proceed without difficulty.

Now we suppose to have N̄ observations r_{t_{i}}(x), i = 1, …, N̄ of the forward curve for x ∈ [0, T], where T is a fixed time horizon. Since the process r is stationary, we can consider N̄ observations r_{t_{i}} = r_{t_{i}}(ω), i = 1, …, N̄ as N̄ realizations of the random variable r_{t} for a given t ≥ 0, which has law N(b, Q_{∞}).

From now on, we indicate r_{i} = r_{t_{i}} for brevity. Then we can approximate:

𝔼[r_{t}] ≃ Ẽ[r_{t}] = \frac{1}{N̄} ∑_{i=1}^{N̄} r_{i}  (10)

Cov (r_{t}, r_{t}) = Q_{∞} ≃ \frac{1}{N̄} ∑_{i=1}^{N̄} r_{i} ⊗ r_{i}  (11)

Now let’s see how we can get the observations r_{i}. In the Italian market the bonds issued by the national bank (Banca d’Italia) can be assumed not to have risk of default. The Banca d’Italia issues BOT (Buoni Ordinari del Tesoro), which are zero-coupon bonds with maturity 3 months, 6 months or 1 year issued each month; CTZ (Certificati del Tesoro Zero-coupon), which are zero-coupon bonds with maturity 18 or 24 months; BTP (Buoni del Tesoro Pluriennali), which are fixed coupon bonds with maturity 3, 5, 7 or 10 years, issued every 3 months. We can get the price of zero-coupon bonds with typical maturities T_{j} = j · 3 months, j = 1, …, K (for example K = 40, that is 10 years). For maturities of up to 2 years we can look directly to the prices of the BOTs and CTZ, and for higher maturities we can decompose the prices of BTP in their different coupons to get the price of the single zero-coupons; this can be done for example using a “strip” technique (see [9]). We suppose we have data in the form B(t, T_{j}), T_{j} = j · 3 months, j = 1, …, K, t = i days, i = 1, …, N̄. We can get the forward rates in this way: since

B(t, T) = exp (− ∫_{0}^{T −t} r(t, u) du)

then

r(t, x) = −\frac{∂}{∂T} log B(t, T)_{|T =t+x}

Here there is a problem. Since T is a variable that can assume only the discrete values T_{j}, we cannot perform a real derivation. Instead, we can approximate the derivative with:

r(t, x) ≃ −\frac{log \frac{B(t,t+x+3 months)}{B(t,t+x)}}{3 months} = \frac{log \frac{B(t,t+x)}{B(t,t+x+3 months)}}{0.25} = 4 log \frac{B(t, t + x)}{B(t, t + x + 3 months)}

If we put x = T − t, then for each day t the values of (B(t, t + x))_{x} have a 1 day shift to the left, so if we have N̄ samples we have a discrete scheme with K · N̄ values. We can interpolate each curve and we can obtain N̄ samples of forward rates’ curves, and we can build the operator Q_{∞} by Eq. (11). Now let’s see how we can get the (τ_{n})_{n} in an efficient way using Karhunen-Loeve expansion.

5 Karhunen-Loeve expansion

In this part we follow the approach of [12]. The Karhunen-Loeve expansion is a tool that is classically used in fluid dynamics giving nice numerical results, and we believe that it also succeeds in treating our problem well.

We suppose we have a random variable u with values in a separable Hilbert space H and we indicate its scalar product with ⟨·, ·⟩ and the relative norm with ∥ · ∥. From now on, we suppose that 𝔼[u] = 0 (we can always reduce the problem to this case by defining u′ = u − 𝔼[u]).

We consider the bilinear functional defined by the variance of u in H:

⟨Rh, k⟩ = 𝔼[⟨u, h⟩⟨u, k⟩] ∀h, k ∈ H

Lemma 4 If u ∈ L^{2}(Ω, F, ℙ; H), then R is a compact selfadjoint operator.

Proof. We calculate the trace of R: if (e_{n})_{n} is an orthonormal basis of H, we have:

Tr R = ∑_{n=1}^{∞}⟨Re_{n}, e_{n}⟩ = ∑_{n=1}^{∞} 𝔼[⟨u, e_{n}⟩⟨u, e_{n}⟩] = 𝔼 [∑_{n=1}^{∞}⟨u, e_{n}⟩^{2}] = 𝔼[∥u∥^{2}] < +∞

this means that R is a trace class operator; this implies that R is compact. The fact that R is selfadjoint follows trivially from its definition. ⋄ Corollary 5 There exists a complete orthonormal set (k_{n})_{n} in H such that Rk_{n} = λ_{n}k_{n}, with λ_{n} ↘ 0, so we can represent R as

R = ∑_{m=1}^{∞} λ_{m}k_{m} ⊗ k_{m}  (12)

Proof. The result is a classical result of functional analysis, and it derives from the fact that R is compact and selfadjoint in H. For a proof, see for example [4] or [10]. ⋄

We notice that in order to know the operator R, it is sufficient to know the sets (λ_{n})_{n} and (k_{n})_{n}. The (k_{n})_{n} has also another interpretation. In fact the following holds:

Lemma 6 The (k_{n})_{n} are extremals of the functional

k → \frac{𝔼[⟨u, k⟩^{2}]}{∥k∥^{2}} = \frac{⟨Rk, k⟩}{∥k, k∥^{2}}

Proof. We calculate the first variation of the functional in the point k̄: ∀ε > 0, k ∈ H_{γ}^{1}, we calculate the functional in k̄ + εk and we obtain:

F (ε, k) = \frac{⟨Rk̄, k̄⟩ + 2ε⟨Rk̄, k⟩ + ε^{2}⟨Rk, k⟩}{∥k̄∥^{2} + 2ε⟨k̄, k⟩ + ε^{2}∥k∥^{2}}

Then k̄ is an extremal if and only if \frac{∂F}{∂ε} (ε, k)_{|ε=0} = 0 ∀k ∈ H_{γ}^{1}. This leads to

2⟨Rk̄, k⟩∥k̄∥^{2} = 2⟨k̄, k⟩⟨Rk̄, k̄⟩ ∀k ∈ H_{γ}^{1}

so all the eigenvectors of R are extremal points and the result follows. ⋄

From the lemma it follows that k_{1} maximizes the function above. In other words, the knowledge of the (k_{n})_{n} gives us the precise behaviour of R; this analysis is analogous to the principal components analysis used in statistics.

6 Calculation of the eigenvalues of R: an iterative method

Now we present an iterative method to find the first eigenvalue λ_{1} of R and the corresponding unitary eigenvector k_{1}.

We fix k^{0} ∈ H, and we define:

k^{n+1} = Rk^{n}  (13)

Lemma 7 We have that

\frac{∥k^{n+1}∥}{∥k^{n}∥} → λ_{1}, \frac{k^{n}}{λ_{1}^{n}} ⟶^{H_{γ}^{1}} ⟨k_{1}, k^{0}⟩k_{1}

Proof. From (12), it follows that k^{n} = ∑_{m=1}^{∞} λ_{m}^{n}k_{m}⟨k_{m}, k^{0}⟩. By induction:

k^{1} = Rk^{0} = ∑_{m=1}^{∞} λ_{m}k_{m}⟨k_{m}, k^{0}⟩

k^{n+1} = Rk^{n} = ∑_{m=1}^{∞} λ_{m}k_{m} ⟨k_{m}, ∑_{m′=1}^{∞} λ_{m′}^{n}k_{m′}⟨k_{m′}, k^{0}⟩⟩ = ∑_{m=1}^{∞} λ_{m}k_{m}λ_{m}^{n}⟨k_{m}, k^{0}⟩

so

\frac{k^{n}}{λ_{1}^{n}} − ⟨k_{1}, k^{0}⟩k_{1} = k_{1}⟨k_{1}, k^{0}⟩ + ∑_{m=1}^{∞} \frac{λ_{m}^{n}}{λ_{1}^{n}} k_{m}⟨k_{m}, k^{0}⟩ − ⟨k_{1}, k^{0}⟩k_{1} = ∑_{m=2}^{∞} (\frac{λ_{m}}{λ_{1}} )^{n} k_{m}⟨k_{m}, k^{0}⟩

and we have the second statement of the lemma. Furthermore,

\frac{∥k^{n+1}∥}{∥k^{n}∥} = \frac{∥k^{n+1}∥}{λ_{1}^{n+1}} \frac{λ_{1}^{n}}{∥k^{n}∥} · λ_{1} = λ_{1} \frac{∥k_{1}⟨k_{1}, k^{0}⟩ + ∑_{m=2}^{∞} (\frac{λ_{m}}{λ_{1}} )^{n+1} k_{m}⟨k_{m}, k^{0}⟩∥}{∥k_{1}⟨k_{1}, k^{0}⟩ + ∑_{m=2}^{∞} (\frac{λ_{m}}{λ_{1}} )^{n} k_{m}⟨k_{m}, k^{0}⟩∥} → λ_{1}

giving the first statement of the lemma. ⋄

In order to find the other eigenvalues and eigenvectors, we can apply the method to the operator R_{1} = R − λk_{1} ⊗ k_{1} and so on.

We search for a number N that would be in some way an index of how much of the terms of the series are “significant”; we define

N = \frac{∑_{n=1}^{∞} λ_{n}}{λ_{1}}

This estimation represents the fact that at least N terms of the series are needed for the sum to be significant. If we use this method, having obtained λ_{1} from the algorithm (13), we get N, since the sum of the eigenvalues of R is given by:

∑_{m=1}^{∞} λ_{m} = ∑_{m=1}^{∞}⟨Rk_{m}, k_{m}⟩ = 𝔼 [∑_{m=1}^{∞}⟨u, k_{m}⟩^{2}] =

= 𝔼 [⟨∑_{m=1}^{∞}⟨u, k_{m}⟩k_{m}, ∑_{m′=1}^{∞}⟨u, k_{m′}⟩k_{m′}⟩] = 𝔼[⟨u, u⟩]

7 Application to the Musiela model

Now we want to apply the Karhunen-Loeve expansion to the process (r_{t})_{t}. We suppose that (r_{t})_{t} is a stationary solution of the Musiela equation (4) in the space H_{γ}^{1} and that we have N̄ observations r_{i}, i = 1, …, N̄ of the forward curve for x ∈ [0, T]. As seen before, we can approximate mean and covariance of the Gaussian random variable r_{t} using Equations (10) and (11). The most intuitive approach would be to apply directly the Karhunen-Loeve expansion to the samples r_{i}, but this would only allow us to get a representation of Q_{∞}:

Q_{∞} = ∑_{n=1}^{∞} ∫_{0}^{∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} du = ∑_{n=1}^{∞} λ_{n}k_{n} ⊗ k_{n}

but, once we obtain the (k_{n})_{n}, we still have to obtain the (τ_{n})_{n}. This is a non trivial problem. One can be tempted to solve the problems

∫_{0}^{∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} du = λ_{n}k_{n} ⊗ k_{n} ∀n ∈ ℕ  (14)

but this is a problem that is solvable if and only if τ_{n} is an eigenvector of A, that is if and only if τ_{n}(x) = Ce^{−αx}. In fact, if τ_{n} is not an eigenvector of A, then the operator on the left hand side has rank greater than 1 and it can’t be represented by an operator of the form λ_{n}k_{n} ⊗ k_{n}. More generally, a reduction of the kind:

∑_{i=1}^{n} ∫_{0}^{∞} S_{u}τ_{i} ⊗ S_{u}τ_{i} du = ∑_{i=1}^{n} λ_{i}k_{i} ⊗ k_{i}

is possible if and only if the vector subspace [τ_{i}]_{i=1,...,n} in H_{γ}^{1} is stable under A. This leads us to take in consideration a family of curves of the kind [p(x)e^{−αx}] with α ∈ ℝ^{+} where p is a polynomial of degree ≤ n (see [1]), but this puts strong limitations on the form of the τ_{n}, thus making our work useless.

Example 8 The easiest case is to search for a τ_{n} which is an eigenvector of S_{u}, as we see in this case: if k_{1}(x) = Ce^{−αx}, then we search for a τ_{1} having the form τ_{1}(x) = τe^{−αx}. For a generic ϕ ∈ H_{γ}^{1}, we have

C^{2}λe^{−αx}⟨ϕ(x), e^{−αx}⟩ = ∫_{0}^{∞} τe^{−α(x+u)}⟨ϕ(x), τe^{−α(x+u)}⟩ du =
= ∫_{0}^{∞} τ^{2}e^{−α(x+u)}e^{−αu}⟨ϕ(x), e^{−αx}⟩ du

Then, simplifying on both sides the terms e^{−αx}⟨ϕ(x), e^{−αx}⟩ we obtain:

C^{2}λ = ∫_{0}^{∞} τ^{2}e^{−2αu} du = τ^{2} [−\frac{1}{2α}e^{−2αu}]_{0}^{+∞} = \frac{τ^{2}}{2αC^{2}}

and we find τ = \sqrt{2αλ}/C. Anyway the problem of inverting the application from the τ_{n} to the k_{n} in the general case does not seem to have an easy solution.

Let us try another approach.

Lemma 9 If S_{u}τ_{n} ⟶^{H_{γ}^{1}} 0 ∀n, then Q_{∞} satisfies the equation

AQ + QA^{∗} + ττ^{∗} = 0  (15)

where ττ^{∗} = ∑_{n=1}^{∞} τ_{n} ⊗ τ_{n}. Proof. We have:

AQ_{∞} + Q_{∞}A^{∗} = A ∑_{n=1}^{∞} ∫_{0}^{+∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} du + ∑_{n=1}^{∞} ∫_{0}^{+∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} duA^{∗} =
= ∑_{n=1}^{∞} ∫_{0}^{+∞} (A(S_{u}τ_{n} ⊗ S_{u}τ_{n}) + (S_{u}τ_{n} ⊗ S_{u}τ_{n})A^{∗}) du =
= ∑_{n=1}^{∞} ∫_{0}^{+∞} (S_{u}τ_{n} ⊗ (AS_{u}τ_{n}) + (AS_{u}τ_{n}) ⊗ S_{u}τ_{n}) du =
= ∑_{n=1}^{∞}[S_{u}τ_{n} ⊗ S_{u}τ_{n}]_{0}^{+∞} = − ∑_{n=1}^{∞} τ_{n} ⊗ τ_{n}

and the result follows. ⋄

Lemma 10 If τ ∈ H^{1}, then S_{u}τ ⟶^{H_{γ}^{1}} 0. Proof. We have:

∥S_{u}τ∥_{H_{γ}^{1}}^{2} ≤ ∥S_{u}τ∥_{H^{1}}^{2} = ∫_{0}^{+∞} (τ^{2}(x + u) + τ^{′2}(x + u)) dx =
= ∫_{u}^{+∞} (τ^{2}(x) + τ^{′2}(x)) dx → 0

since τ ∈ H^{1}, and we have the result. ⋄

Corollary 11 Q_{∞} satisfies Equation (15). Proof. Since Q_{∞} is an invariant measure then, from section 3, τ_{n} ∈ H^{1} ∀n, and the result follows from lemma 9 and 10. ⋄

Let’s see which consequences such a relation will have in our case. If Q = λk ⊗ k, then:

AQf = λ⟨k, f⟩k′ QA^{∗}f = λ⟨k′, f⟩k

so AQ and QA^{∗} are operators going from H_{γ}^{1} to two one dimensional subspaces generated respectively by k′ and k. In particular we have:

AQ = λk ⊗ k′, QA^{∗} = λk′ ⊗ k If Q is a (possibly infinite) sum, then A and A^{∗} apply to each component of Q in the way seen above. This means that we can implement the Karhunen-Loeve expansion in this way: instead of decomposing Q_{∞} and going back to ττ^{∗}, we can decompose ττ^{∗} itself! In order to do this, we only need this:

Assumption 12 Tr ττ^{∗} < +∞.

The assumption seems quite natural when dealing with an invariant measure. In fact requiring that Tr Q_{∞} < +∞ is equivalent to requiring that ∫_{0}^{+∞} ∑_{n=1}^{∞} ∥S_{u}τ_{n}∥^{2} du < +∞. This implies ∑_{n=1}^{∞} ∥S_{u}τ_{n}∥^{2} < +∞ du-a.s.; but (S_{u})_{u} is a C^{0} semigroup in H_{γ}^{1} such that ∥S_{u}∥ ≤ e^{\frac{γ}{2} u} (see [17]), then the function u → e^{− \frac{γ}{2} u}∥S_{u}τ∥ is decreasing ∀τ ∈ H_{γ}^{1}, so e^{− \frac{γ}{2} u} ∑_{n=1}^{∞} ∥S_{u}τ_{n}∥^{2} < +∞ ∀u > 0. Finally Tr Q_{∞} < +∞ implies ∑_{n=1}^{∞} ∥S_{u}τ_{n}∥^{2} < +∞ ∀u > 0. Besides, if Assumption 12 is not true, then it is not possible to fix a stopping criterion based on a number

N = \frac{Tr ττ^{∗}}{μ_{1}}

where μ_{1} is the first eigenvalue of ττ^{∗}, because N turns out to be infinity. Assumption 12 implies the following

Corollary 13 There exists a complete orthonormal set (e_{n})_{n} in H_{γ}^{1} such that ττ^{∗}e_{n} = μ_{n}e_{n} with μ_{n} ↘ 0, and we can represent ττ^{∗} as

ττ^{∗} = ∑_{m=1}^{∞} μ_{n}e_{n} ⊗ e_{n}

If Assumption 12 holds, then:

∑_{m=1}^{∞} μ_{m} = Tr ττ^{∗} = −Tr (AQ_{∞} + Q_{∞}A^{∗}) = − ∑_{m=1}^{∞}⟨(AQ_{∞} + Q_{∞}A^{∗})k_{m}, k_{m}⟩ =

= − ∑_{m=1}^{∞}⟨AQ_{∞}k_{m}, k_{m}⟩ − ∑_{m=1}^{∞}⟨Q_{∞}A^{∗}k_{m}, k_{m}⟩ =

= −𝔼 [∑_{m=1}^{∞}⟨r̄_{t}, k_{m}⟩⟨r̄_{t}, A^{∗}k_{m}⟩] − 𝔼 [∑_{m=1}^{∞}⟨r̄_{t}, A^{∗}k_{m}⟩⟨r̄_{t}, k_{m}⟩] =

= −2𝔼 [∑_{m=1}^{∞}⟨r̄_{t}, k_{m}⟩⟨Ar̄_{t}, k_{m}⟩] = −2𝔼[⟨r̄′_{t}, r̄_{t}⟩]

At last, we can write the Musiela model in this way:

\{ \begin{matrix}dr_{t}(x) = (\frac{∂}{∂x}r_{t}(x) + ∑_{n=1}^{∞} τ_{n}(x) ∫_{0}^{x} τ_{n}(u) du) dt + ∑_{n=1}^{∞} τ_{n}(x) dW_{t}^{n} \\ r_{0} ∈ H_{γ}^{1}(ℝ^{+}) \end{matrix}

with τ_{n} = \sqrt{μ_{n}}e_{n} ∀n ∈ ℕ^{∗}. The main advantage of this representation, as already outlined in [6], is the following: if we define r_{t}^{N} as the solution of the equation obtained by substituting in the previous one the infinite sum with a sum until N, then r^{N} has the diffusion term that best approximates the one of r in the class of H-valued processes that satisfy a stochastic differential equation driven by N Brownian motions.

Now let’s see how we can implement the Karhunen-Loeve expansion on ττ^{∗}. As before, we suppose we have N̄ observations r_{i}(x), i = 1, …, N̄ of the forward curve for x ∈ [0, T]. Since the process is stationary, we can consider a realization of r_{t} consisting of the N̄ observations r_{i}, i = 1, …, N̄.

1st step — we center r, passing to r̄_{t} = r_{t} − 𝔼[r_{t}]:

r̄_{j} = r_{j} − \frac{1}{N̄} ∑_{i=1}^{N̄} r_{i}

then r̄_{t} ∼ N(0, Q_{∞}), and Cov (r_{t}, r_{t}) = Q_{∞} is given by (11).

2nd step — we build ττ^{∗} = AQ_{∞} + Q_{∞}A^{∗}:

ττ^{∗} = −\frac{1}{N̄} ∑_{i=1}^{N̄}(r̄′_{i} ⊗ r̄_{i} + r̄_{i} ⊗ r̄′_{i})

3rd step — we apply the algorithm for ττ^{∗}: we fix an initial k^{0} and we define:

e^{n+1} = −\frac{1}{N̄} ∑_{i=1}^{N̄}(⟨r̄′_{i}, e^{n}⟩r̄_{i} + ⟨r̄_{i}, e^{n}⟩r̄′_{i})

and we find μ_{1} and e_{1}.

4th step — we define R_{1}:

R_{1} = −\frac{1}{N̄} ∑_{i=1}^{N̄}(r̄′_{i} ⊗ r̄_{i} + r̄_{i} ⊗ r̄′_{i}) − μ_{1}e_{1} ⊗ e_{1}

and we come back to the 3rd step. We can stop to the suitable N, obtained as:

N = \frac{∑_{m=1}^{∞} μ_{m}}{μ_{1}} = \frac{−2𝔼[⟨r̄′, r̄⟩]}{μ_{1}} ≃ − \frac{\frac{2}{N̄} ∑_{i=1}^{N̄}⟨r′_{i}, r_{i}⟩}{μ_{1}}

References

[1] T. Bjork, and B. J. Christensen, Forward rate models and invariant manifolds, preprint

[2] A. Brace, D. Gatarek, and M. Musiela, The market model of interest rate dynamics, Mathematical Finance, vol. 7, no. 2, 1997, 127-155.

[3] A. Brace, and M. Musiela, A multifactor Gauss-Markov implementation of Heath, Jarrow, and Morton, Mathematical Finance, vol. 4, no. 3, 1994, 259–283

[4] H. Brezis, Analyse fonctionnelle — Théorie et applications, Masson Editeur, 1983 [5] G. Da Prato and J. Zabczyk, Stochastic equations in infinite dimensions, Cambridge

University Press, 1992

[6] R. Douady, Yield curve smoothing and residual variance of fixed income positions, preprint

[7] B. Goldys, and M. Musiela, On partial differential equations related to term structure models, preprint

[8] D. Heath, R. Jarrow and A. Morton, Bond pricing and the term structure of interest rates: a new methodology, Econometrica 61, 1992, 77–105.

[9] J. Hull, Options, futures and other derivative securities, Prentice Hall, 1996

[10] T. Kato, Perturbation theory for linear operators, Springer Verlag, 1980

[11] D. Lamberton, and B. Lapeyre, Introduction au calcul stochastique appliquée à la Finance, Société de Mathématiques Appliquées et Industrielles, 1991

[12] J. L. Lumley, Stochastic tools in turbulence, New York Academic Press, 1970

[13] M. Metivier, Stochastic partial differential equations in infinite dimentional space, Scuola Normale Superiore di Pisa, 1988

[14] K. R. Miltersen, K. Sandmann, and D. Sondermann, Closed form solutions for term structure derivatives with log-normal interest rates, University of Bonn Discussion Paper B–308

[15] M. Musiela, Stochastic PDEs and term structure models, Journees Internationales de Finance, IGR-AFFI, La Baule, 1993

[16] K. Sandmann, and D. Sondermann, A term structure model and the pricing of interest rate derivatives, The Review of Futures Markets, 12 (2), 392–423

[17] T. Vargiolu, Invariant measures for the Musiela equation with deterministic diffusion term, to appear in Finance and Stochastics

[18] M. Yor, Sur les intégrales stochastiques à valeur dans un espace de Banach, Annales I. H. P., vol. X, n. 1, 1974, pp. 31–36