D.E.A de Probabilités Appliquées - Filiére Finances Université Pierre et Marie Curie - Paris VI

Finite dimensional approximation for the Musiela model in the Gaussian case

Tiziano Vargiolu

sous la direction de

Nicole El Karoui, Franco Flandoli et Jerzy Zabczyk

Année Académique 1996/’97

Contents

Introduction iv

1 The Musiela model 1 1.1 Presentation of the Musiela model . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Formulation of the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Preliminaries about Langevin equations in Hilbert spaces . . . . . . . . . . 4 1.4 The space H_{γ}^{1} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 The space H^{1} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 The interest rates 8 2.1 The Musiela model in the Gaussian case . . . . . . . . . . . . . . . . . . . . 8 2.2 From the Musiela model to the Vasicek model . . . . . . . . . . . . . . . . . 8 2.3 From the Musiela model to a state variables model . . . . . . . . . . . . . . 9 2.4 Considerations on the Hilbert space H . . . . . . . . . . . . . . . . . . . . . 11 2.5 Structuring the Musiela model with respect to state variables . . . . . . . . 11

2.5.1 The Vasicek model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5.2 The general state variables model . . . . . . . . . . . . . . . . . . . . 12

2.6 Approximating the Musiela model with state variables . . . . . . . . . . . . 12 2.7 Open problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Calibrating the model 16 3.1 An “implied” approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Two “historical” approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2.1 With less datas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.2 With more datas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3 Karhunen–Loeve expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.1 The representation problem . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.2 The “finite energy” case: Karhunen–Loeve expansion . . . . . . . . . 21 3.3.3 Calculation of the eigenvalues of R . . . . . . . . . . . . . . . . . . . 21 3.3.4 Rate of convergence of the series . . . . . . . . . . . . . . . . . . . . 21 3.3.5 Application: the Musiela model . . . . . . . . . . . . . . . . . . . . . 22 4 Robustness of the Musiela model 24

4.1 A functional differential equation for interest rates derivatives . . . . . . . . 24 4.2 Robustness of the Musiela model . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3 Sensibility of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.4 Open problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

A Positivity of rates 32 A.1 The Vasicek model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 A.2 The d-dimensional Gaussian model . . . . . . . . . . . . . . . . . . . . . . . 37 A.3 Open problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Introduction

The aim of this work is to analyze the Musiela model for the forward rates presented in [13], and in particular the Gaussian version of this model. The Musiela model is a reparametriziation of the HJM model [10], so it is essentially a set of relations which an interest rate model must satisfy in order not to have arbitrage opportunities in the market.

This work is intended to be an attempt to analyze this model, and in particular what happens if we try to approximate the forward rate dynamics with a finite or infinite dimensional Gaussian model. Various topics associated to this problem are touched, as the reduction of the model to a finite dimensional one, the calibration of the model to hystorical market datas in order to get the most possibly precise shape of the forward rate curve, and the so called “robustness” of the model, that is a kind of good dependence of the model from the extimated datas, in such a way that even if an agent extimates the coefficients of the model making errors, he (she) succeeds in obtaining a superhedging strategy, and that a little error on the extimation of the curve gives a little error on the price and the replication of the derivative assets. In this work we wanted to present ideas related to this topics, giving also open problems, without insisting too much on technicalities; this also because, while some parts (like the reduction of the model to a finite dimensional one) don’t present technical difficulties, other ones (in particular the part on the robustness of the model) would require heavy mathematical tools, whose complexity would be parhaps sufficient for a Ph. D. thesis! Finally there is an appendix devoted to the aspect that is considered on of the main drawbacks of the Gaussian models, that is the possibility of negative rates; also in this appendix we try to present the main ideas and don’t insist too much on the technical side.

In the work we proceed in this way. First we present the Musiela model: its core is the stochastic differential equation

dr(t, x) = (\frac{∂}{∂x}r(t, x) + ∑_{n=1}^{∞} τ_{n}(t, x) ∫_{0}^{x} τ_{n}(t, u) du) dt + ∑_{n=1}^{∞} τ_{n}(t, x) dw_{t}^{n}, t ≥ 0

r(0, ·) ∈ W_{loc}^{1,1}(ℝ^{+}) given

where r(t, ·) ∈ W_{loc}^{1,1}(ℝ^{+}) ℚ-a.s. ∀t, and (w^{n})_{n} are independent real Brownian motions; we often will write the equation in a shorter way:

\{ \begin{matrix}dr_{t} = (Ar_{t} + c_{t}) dt + \{τ_{t}, dW_{t}\}, t ≥ 0 \\ r_{0} ∈ W_{loc}^{1,1}(ℝ^{+}) given \end{matrix}

We notice that the choice of τ determines entirely the evolution of the forward rate; in particular we can choose to consider a finite or infinite number of driving Brownian motions, or we can choose to consider a deterministic τ, thus obtaining a Gaussian model. We also notice that the general model is in the space W_{loc}^{1,1}(ℝ^{+}), which is a locally convex space; in this kind of spaces the theory of stochastic integration is rather difficult, while it becomes much simpler in Hilbert spaces. Though such a problem can be treated in its generality, here we present a simpler problem, that is to analyze the model supposing that the forward rates take values in some Hilbert spaces contained in W_{loc}^{1,1}(ℝ^{+}). We present some results from an existing paper [15], regarding the existence and uniqueness of a so called mild solution of the Musiela equation in the simplest case, that is when τ is deterministic and there is a finite number of driving Brownian motions; the problems is solved in two different Hilbert spaces, namely H_{γ}^{1}(ℝ^{+}) and H_{γ}^{1}(ℝ^{+}), in which the equation has different behaviour. We present also some results linked to the existence and (possibly) the uniqueness of an invariant measure for the equation. This last problem is linked to the calibration of the model: we will find out that an implicit calibration gives very few informations on the model, while a particular hystorical calibration succeeds in giving much more information. But each time that a calibration based on time series is used (as in our case), there is always the implicit hypothesis that the process is stationary; this is equivalent to the existence of an invariant measure for the driving equation.

This study introduces us to the main topic of this work: even if we restrict the model to the Gaussian case, the natural framework for it is W_{loc}^{1,1}(ℝ^{+}), and the choices of the two Hilbert spaces appears arbitrary; moreover, there is no apparent reason to impose that there is only a finite number of driving Brownian motions. The main question is then: is it possible to find a “natural” Hilbert space in which the equation “behaves well”? More formally, does it exist a Hilbert space and an orthonormal basis (e_{n})_{n} such that, if we project the solution of the equation on the first N components, we make a small error. The simplest case in which this happens is when the model is itself finite dimensional, so we first analyze how the Musiela equation becomes when we have the Vasicek model (the simplest Gaussian model) or more generally a finite dimensional Gaussian model in which the state variables are particular reference rates. Starting from this, we make some qualitative considerations on which could be the “right” Hilbert space for the equation, but these considerations are neither rigorous nor fully satisfactory! We finish this session with a first attempt to calculate the error we make when we approximate the “real” dynamics of the model with a finite dimensional approximation. This leaves a first series of open problems, linked to the choices of the space, of the basis, and of the criterium to use to calculate the error.

Now we come to the part about the calibration of the model: we present an implicit approach, that is classically used in the market, but doesn’t give us much information; then we present two hystorical approaches, and we develop with more detail the first one; in particular, we apply the method known as “Karhunen–Loeve expansion”. This method is classically used in fluidodynamics in order to find an orthonormal basis, whose vectors are known as “empirical orthonormal functions”, which diagonalize the covariance operator of an infinite dimensional random variable; it turns out that this method gives nice numerical results, in the sense that 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 give a way to apply the method to our situation, considering a time series as a sample of an infinite dimensional random variable having as law an invariant measure of the equation. Here there is a hidden detail: in fact we always studied the Musiela model under the risk-neutral probability ℚ, while we can obtain from the time series an estimation of the hystorical 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 safely.

Finally we arrive to the part devoted to the robustness of the model. First of all we present a functional differential equation that gives, under no arbitrage assumptions, the price of an interest rate derivative under the hypothesis that the forward rates follow a (even non linear) Musiela model. We do this for two reason: the first is to make a link between this case and the classical finite dimensional cases (see [4], [16]) in which the price of interest rate derivatives can be derived via a PDE; besides, we can introduce the mathematical tools that we need later and an equation that is simpler, but similar, to the one known as Black–Scholes–Barenblatt equation, that we’ll meet later. Then we consider the real robustness problem: we suppose that the dynamics of the forward rates is given by the Musiela equation, but that we aren’t able to determine precisely the process τ, but we only know that τ takes values in some compact set Σ ⊆ H^{d}. This model is a quite general one, and it is an analogous of the stochastic volatility models for stock prices analyzed in literature in the last years (see for example [8] and [1]). Now we consider an operator in the market that wants to create and sell an European interest rate derivative asset with payoff F (rT̄) depending on the value in the maturity T̄ of all the forward rates curve rT̄. In many relevant cases, the payoff F will depend not on all the curve rT̄, but only on the values of N particular bonds with maturity T_{i} > T̄, so it will be of the form

F (B(T̄, T_{1}), …, B(T̄, T_{N})), T̄ < T_{1} < … < T_{N}

From now on, we suppose that F is of the form above. In order to hedge this asset, the operator chooses to use the homogeneous Gaussian Musiela model described in Chapter 1. Then he (she) builds a self-financing strategy in this way: he (she) fixes the price C_{t} of the option at time t as

C_{t} = u(t, r_{t})

where u is the solution of a nonlinear functional differential equation, namely a Hamilton-Jacobi-Bellmann type equation;

\{ \begin{matrix}\frac{∂u}{∂t}(t, r) + sup_{τ∈Σ} Lu(t, r) − r(0)u(t, r) = 0, (t, r) ∈ [0, T̄] × H \\ u(T̄, r) = F (r), r ∈ H \end{matrix}

where L is a particular second order operator acting on functionals on the Hilbert space H. We remark that the finite-dimensional versions of this equation appears for the first time in one dimension in [1], and it is generalized to the finite-dimensional case in [14]; as in the previous works the Equation above will be called from now on the Black-Scholes-Barenblatt equation associated to the contingent claim F. To our knowledge, this is the first time that such an approach is used in hedging interest rates derivatives, so it is also the first time that such an equation is analyzed in this context. In order to hedge the claim F, the agent builds the hedging portfolio holding the quantity Δ_{t}^{i} of the zero-coupon bond B(t, T_{i}) in the hedging portfolio Π defined in the classical delta-hedging way as the derivative of the solution of the equation with respect to the corresponding asset. Finally we define the tracking error e_{t} = Π_{t} − C_{t}, with e_{0} = 0. Intuitively, the tracking error gives the error made by the operator in his (her) estimation, or better the difference between the hedging portfolio held by the agent and the option he (she) sold. We find out that if C is a C^{1,2}([0, T) × R^{n}) solution of a suitable finite dimensional version of the BSB equation, then (e_{t})_{t} is an increasing process, that is the delta-hedging gives a superstrategy for our operator. This is the part where the open problems are more deep, not only because of the intrinsic difficulty of the study of the BSB equation in a Hilbert space, but also because all the cases of interest in Finance are such that the final payoff is not even C^{1}; it seems so that it would be necessary to develop some tools analogous to what viscosity solution are in finite dimensions.

The appendix is devoted to one of the main drawbacks of the Gaussian models, that is the possibility of negative rates with positive probability. First of all we concentrate on the one-dimensional Vasicek model, and calculate the probability that at some time there will be a negative forward rate; this probability being 1, we calculate the probability that there will be a negative forward rate before a fixed instant T. While we succeed to calculate exactly the first probability in terms of the solution of a ODE, we can’t calculate explicitly the second probability, because it is given in terms of the solution to a PDE, that doesn’t appear to have an explicit solution. Then we give two extimations of this probability, one using Kolmogorov’s maximal inequality and the fact that the process (e^{at}r_{t})_{t} is a martingale, and the other one making a time change, seeing the stochastic integral appearing in the solution r_{t} as a Brownian motion with respect to the new time scale, and using a sort of Cameron-Martin theorem. Then, we try to do a similar work for the d-dimensional Gaussian model, even if the work in this direction is far to be complete. Anyway, the final aim of this appendix would be to extend all or some of the results found in the case of Vasicek’s model to the general Gaussian Musiela model in a Hilbert space: from a mathematical point of view, we feel that this kind of problem is interesting and far from being trivial.

The work is organized as follows. In Chapter 1 there is an introduction to the Musiela model, with the main technical results concerning the equation in the spaces H_{γ}^{1}(ℝ^{+}) and H^{1}(ℝ^{+}). In Chapter 2 we want to evaluate when the homogeneous Gaussian Musiela model can be reduced to a finite factor model; moreover, we want to find a Hilbert space H which behaves well with respect to this interpretation. When this is not possible, we want to evaluate the error we do when we approximate a term structure having a Musiela dynamics with a finite factor model. In Chapter 3 we concentrate on the derivation of the forward curve from economical datas; in particular, we present an “implicit” approach, that derives bonds’ volatility from the prices of derivatives in the market, and two “historical” ones, that derive bonds’ volatility from the prices of bonds in the past. In particular, we develop one hystorical approach using the Karhunen–Loeve expansion. In Chapter 4 we present the functional differential equation satisfied by the bonds’ prices under no arbitrage assumptions. Then we present the Black–Scholes–Barenblatt equation, and a finite dimensional version of this equation, and prove that in certain cases the strategy built starting from the BSB equation is a superhedging one. Finally we present the concept of extended duration and show the possible links of this concept with our model. At last in the Appendix 1 we concentrate on the positivity of interest rates; first we calculate the probability to have a negative rate for the Vasicek model, then we make the same calculations in the case of a Gaussian d-dimensional model.

The author wishes to thank Nicole El Karoui for having given the idea of finding a “right” Hilbert space for the Musiela model, which is the idea that inspired this work; Franco Flandoli for his precious advices and his suggestion to use Karhunen–Loeve expansion; the Laboratory of Probability of the University of Paris VI, which has hosted him during the revision of this work; Vincent Lacoste for interesting talkings about robustness; Marek Musiela, which in a talk suggested the idea of the robustness of his model; root@linuz.sns.it (that is, Massimiliano Ghilardi and Michele Mader), for having allowed him to finish this work in spite of various accidents happened to the Calculus Center of the Scuola Normale Superiore in the late June 1997; Jerzy Zabczyk for having assisted him in the part on the positivity of rates.

Chapter 1

The Musiela model

1.1 Presentation of the Musiela model

Let’s introduce now the Musiela interest rate term structure model. We have a probability space (Ω, F, ℙ), with the filtration (F_{t})_{t≥0}. We suppose to have only a single bond in the market, and 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}. We suppose to have a two parameters process (r(t, x))_{t,x≥0}, called forward rate process, in which r(t, x) represents the rate at which at time t one can enter a contract (to borrow or lend money) at time T = t + x, for a short period of time. Then the price at time t of a bond expiring at time T = t + x is given by:

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

We call spot rate process the process (r(t, 0))_{t}; 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 + x is given by:

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

Now we add the hypotheses that there exists a k-dimensional standard Brownian motion (W_{t})_{t} adapted to (F_{t})_{t}, and there exist two two-parameter progressively measurable processes (α(t, x))_{t,x≥0}, (τ(t, x))_{t,x≥0} such that ∀x ≥ 0, (α(t, x))_{t} has trajectories in L^{1}([0, +∞); ℝ), (τ^{∗}(t, x))_{t} has trajectories in L^{2}([0, +∞); ℝ^{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}  (1.1)

A natural thing to ask is that the actualized bond price process (P̃(t, T))_{t} is a martingale under a probability measure ℚ equivalent to ℙ, which is called risk-neutral probability. To this purpose, we cite the:

Theorem 1.1 : if we model the forward rate with an equation of the form (1.1), then the following facts are equivalent:

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

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

α(t, x) = \frac{∂}{∂x} (r(t, x) + \frac{1}{2} |∫_{0}^{x} τ(t, u) du|^{2})

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 [13]. ⋄

The forward rate model obtained in this way is then the following:

\{ \begin{matrix}dr(t, x) = (\frac{∂}{∂x}r(t, x) + τ^{∗}(t, x) ∫_{0}^{x} τ(t, u) du) dt + τ^{∗}(t, x) dW_{t} \\ r(0, x) given \end{matrix}  (1.2)

We notice that this model is very general, and for particular choices of the process τ we can obtain models already studied in literature: for example if τ(x) = σe^{−ax}, we obtain the Vasicek model, if τ(x) = σ, we obtain the Ho-Lee model, and if τ(x) = σc(x)\sqrt{r(t, 0)}, for a particular deterministic function c, we obtain the CIR model; the first two models belong to the class of Gaussian models, which we will treat with more detail; we notice that it is possible to obtain also non-Gaussian models, as the CIR or others.

Now we will study this model using more powerful instruments. We now view the process (r(t, x))_{t,x} no more as a real two parameter process, but as a one parameter process taking values in the absolutely continuous functions space on ℝ^{+}, which will be identified with the Sobolev space W_{loc}^{1,1}(ℝ^{+}), and we write the equation (1.2) as a stochastic differential equation taking values in W_{loc}^{1,1}(ℝ^{+}).

Taking:

A = \frac{∂}{∂x}, τ_{t}(x) = τ(t, x), c_{t}(x) = τ_{t}^{∗}(x) ∫_{0}^{x} τ_{t}(u) du

then equation (1.2) becomes the following equation taking values in W_{loc}^{1,1}(ℝ^{+}):

\{ \begin{matrix}dr_{t} = (Ar_{t} + c_{t}) dt + τ_{t}^{∗} dW_{t} \\ r_{0} ∈ W_{loc}^{1,1}(ℝ^{+}) \end{matrix}  (1.3)

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) ∀t, x such that t + x ≤ T ℚ-a.s., then the process (B̃(t, T))_{t} is a martingale; then the hypotheses of the martingale representation theorem hold, and so every 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 construct a self-financing portfolio strategy based on B(·, T) and on β(·), which simulates the claim (see [13]).

1.2 Formulation of the problem

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}  (1.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 spaces the theory of stochastic integration is rather difficult, while it becomes much simpler in Hilbert spaces. Though such a problem can be treated in its generality, here we present a simpler problem, that is to analyze the model supposing that the forward rates take values in some Hilbert spaces contained in W_{loc}^{1,1}(ℝ^{+}).

We now present a study of the equation in two particular separable Hilbert spaces contained in W_{loc}^{1,1}(ℝ^{+}); this study is contained in ([15]), and here we present only the main results. We use the theory of the stochastic integration in separable Hilbert spaces contained in [5] which is somewhat complete. In particular, we will study the equation in the spaces:

H_{γ}^{1}(ℝ^{+}) = W_{γ}^{1,2}(ℝ^{+}) = {u : ℝ^{+} → ℝ | u ∈ L_{γ}^{2}(ℝ^{+}), u′ ∈ L_{γ}^{2}(ℝ^{+})}
H^{1}(ℝ^{+}) = W^{1,2}(ℝ^{+}) = {u : ℝ^{+} → ℝ | u ∈ L^{2}(ℝ^{+}), u′ ∈ L^{2}(ℝ^{+})}

where u′ is the weak derivative of u with respect to x, and:

L_{γ}^{2}(ℝ^{+}) = {u : ℝ^{+} → ℝ | u meas. and s.t. ∫_{0}^{+∞} u^{2}(x)e^{−γx} dx < +∞}

L^{2}(ℝ^{+}) = {u : ℝ^{+} → ℝ | u meas. and s.t. ∫_{0}^{+∞} u^{2}(x) dx < +∞}

H_{γ}^{1} is a Hilbert space with the scalar product:

⟨f, g⟩_{H_{γ}^{1}} = γ ∫_{0}^{+∞} f(x)g(x)e^{−γx} dx + γ ∫_{0}^{+∞} f′(x)g′(x)e^{−γx} dx  (1.5)

and H^{1} is a Hilbert space with the scalar product:

⟨f, g⟩_{H^{1}} = ∫_{0}^{+∞} f(x)g(x) dx + ∫_{0}^{+∞} f′(x)g′(x) dx  (1.6)

In all the spaces cited above we will examine the following problems:

• to find an explicit mild solution

• to find the marginal law of the solutions r_{t} ∀t > 0 and the joint law of the process (r_{t})_{t>0}

• existence, uniqueness (when it holds) and characterization of the invariant measures for the solution process

1.3 Preliminaries about Langevin equations in Hilbert spaces

Now we give a bit of general theory about the Langevin equation on separable Hilbert spaces.

We suppose we have a separable Hilbert space H, such that r_{0} ∈ L^{2}(Ω, F_{0}, ℙ; H), c ∈ H, that A : D(A) → H is the infinitesimal generator of a C^{0} semigroup of linear continuous operators (see appendix), that (W_{t})_{t} is a Q-Wiener process in H, that is a (F_{t})_{t}-stochastic process such that:

1. W_{0} = 0

2. (W_{t})_{t} has continuous trajectories

3. ∀s < t, W_{t} − W_{s} is independent of F_{s}

4. ∀s < t, W_{t} − W_{s} has gaussian law equal to N(0, (t − s)Q), where Q : H → H is linear, bounded, symmetric, positive and such that Tr Q < +∞

and that, if we call Q_{t} = S_{t}QS_{t}^{∗}, then ∫_{0}^{t} Tr Q_{u} du < +∞. Then the equation

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

has exactly one mild solution, which is the Ornstein-Uhlenbeck process in the Hilbert space H:

r_{t} = S_{t}r_{0} + ∫_{0}^{t} S_{t−u}c du + ∫_{0}^{t} S_{t−u} dW_{u}

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  (1.8)

and functional covariance (if t ≤ v):

Cov (r_{t}, r_{v}) = S_{v}Cov (r_{0}, r_{0})S_{t}^{∗} + ∫_{0}^{t} S_{v−u}QS_{t−u}^{∗} du  (1.9)

This means that ∀f, g ∈ H:

𝔼[⟨f, r_{t}⟩_{H}] = ⟨f, S_{t}𝔼[r_{0}]⟩_{H} + ⟨f, ∫_{0}^{t} S_{t−u}c du⟩_{H}

Cov (⟨f, r_{t}⟩_{H}, ⟨g, r_{v}⟩_{H}) = ⟨(S_{v}Cov (r_{0}, r_{0})S_{t}^{∗})f, g⟩_{H} + ∫_{0}^{t}⟨S_{v−u}QS_{t−u}^{∗}f, g⟩_{H} =
= ⟨Cov (r_{0}, r_{0})S_{t}^{∗}f, S_{v}^{∗}g⟩_{H} + ∫_{0}^{t}⟨QS_{t−u}^{∗}f, S_{v−u}^{∗}g⟩_{H} du

For a proof of these facts, see [5].

Applying these results to our equation, since A = \frac{∂}{∂x}, the natural semigroup is the translation semigroup (S_{t})_{t} such that:

S_{t}f(x) = f(t + x), t, x ≥ 0

Moreover, the Wiener process is of the kind (∑_{i=1}^{k} τ_{i}W_{t}^{i})_{t}, where τ ∈ H^{k} and (W_{t})_{t} is a real k-dimenional Wiener process; then Q = ∑_{i=1}^{k} τ_{i} ⊗ τ_{i}, in the sense that:

⟨Qf, g⟩ = ∑_{i=1}^{k}⟨τ_{i}, f⟩_{H}⟨τ_{i}, g⟩_{H} ∀f ∈ H

Then Tr Q = ∑_{i=1}^{k} ∥τ_{i}∥_{H}^{2} < +∞, and Tr Q_{t} = ∑_{i=1}^{k} ∥S_{t}τ_{i}∥^{2}. Since, if (S_{t})_{t} is a C^{0}- semigroup in H, then it exist M, α ∈ ℝ^{+} such that ∥S_{t}∥_{L(H,H)} ≤ M𝔼^{αt}, then the condition

∫_{0}^{t} Tr Q_{u} du = ∑_{i=1}^{k} ∫_{0}^{t} ∥S_{u}τ_{i}∥^{2} du ≤ ∑_{i=1}^{k} ∥τ_{i}∥^{2} ∫_{0}^{t} M𝔼^{αu} du < +∞

is always verified. The solution becomes:

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

= 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}  (1.10)

The solution is a Gaussian process with functional mean (1.8) and functional covariance (1.9); this means that ∀f, g ∈ H:

𝔼[⟨r_{t}, f⟩_{H}] = ⟨r_{0}(· + t), f⟩_{H} + \frac{1}{2} ∑_{i=1}^{k} ⟨((∫_{0}^{·+t} τ_{i}(u) du)^{2} − (∫_{0}^{·} τ_{i}(u) du)^{2}) , f⟩_{H}

Cov (⟨r_{t}, f⟩_{H}, ⟨r_{v}, g⟩_{H}) =

= ⟨S_{v}Cov (r_{0}, r_{0})S_{t}^{∗}f, g⟩_{H} + ∑_{i=1}^{k} ⟨(∫_{0}^{s} τ_{i}(· + t − u) ⊗ τ_{i}(· + v − u) du) f, g⟩_{H} =

= ⟨Cov (r_{0}, r_{0})S_{t}^{∗}f, S_{v}^{∗}g⟩_{H} + ∑_{i=1}^{k} ∫_{0}^{t} ⟨τ_{i}(· + t − u), f⟩_{H}⟨τ_{i}(· + v − u), g⟩_{H} du

In some applications it can be useful to know only mean and covariance of particular rates and not necessarily all the forward rates curve (for example, to know the spot rate at time t, it is sufficient to know only r_{t}(0)). This can be done with the theory given before: if the Hilbert space H is a subset of the space of the continuous functions C^{0}([0, +∞)), then ∀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

1.4 The space H_{γ}^{1}

Theorem 1.2 : if τ_{i} ∈ H_{γ}^{1} ∩ H^{1} ∩ L_{γ}^{4} ∀i = 1, …, k and r_{0} ∈ L_{γ}^{2}, then it exists a unique mild solution of the equation (1.4), given by Eq. (1.10). The solution is a Gaussian process with functional mean (1.8) and functional covariance (1.9).

Theorem 1.3 : given the equation (1.4) in the Hilbert space H_{γ}^{1}, a necessary and sufficient condition to have an invariant measure is τ_{i} ∈ H_{γ}^{1}(ℝ^{+}) ∩ H^{1}(ℝ^{+}) ∀i = 1, …, k, and there exist infinitely many invariant measures. In particular, the measures of the kind:

δ_{b^{∗}(·)+b_{0}} ∗ N(0, Q_{∞}), b_{0} ∈ ℝ
μ_{f} ∗ N(0, Q_{∞}), f ∈ H_{γ}^{1} periodical

are invariant measures, where:

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

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}  (1.12)

f ∈ H_{γ}^{1} is periodical with period T and μ_{f} is the measure induced by the uniform density on [0, T] via the application t → f(t + ·).

Lemma 1.4 If τ ∈ H_{γ}^{1} ∩ H^{1}, and the initial data is a deterministic r_{0} ∈ H_{γ}^{1}, such that lim_{x→+∞} r_{0}(x) = r_{0}(∞) ∈ ℝ, lim_{x→+∞} r′_{0}(x) = 0, then the marginal law of r_{t} converges for t → +∞ to the invariant measure N (r_{0}(∞) + b_{0} + b(x), Q_{∞}), where b and Q_{∞} are given respectively by (1.12) and (1.11), and b_{0} is given by

b_{0} = \frac{1}{2} ∑_{i=1}^{k} (∫_{0}^{+∞} τ_{i}(u) du)^{2} .  (1.13)

1.5 The space H^{1}

Theorem 1.5 : if ∀i = 1, …, k, τ_{i}(x) ∈ H^{1}(ℝ^{+}) ∩ L^{4}(ℝ^{+}), \sqrt{x}τ_{i}(x) ∈ L^{2}(ℝ^{+}), \sqrt{x}τ′_{i}(x) ∈ L^{2}(ℝ^{+}) and r_{0} ∈ H^{1}(ℝ^{+}), then the equation (1.4) has an unique solution in H^{1}(ℝ^{+}), given by (1.10), which is a gaussian process with functional mean (1.8) and functional covariance (1.9).

Theorem 1.6 : given the equation (1.4) in the Hilbert space H^{1}, necessary and sufficient conditions to have an invariant measure are:

∀i = 1, …, k : \{ \begin{matrix}τ_{i}(x) ∈ H^{1}(ℝ^{+}) ∩ L^{1}(ℝ^{+}) \\ \sqrt{x}τ_{i}(x) ∈ L^{2}(ℝ^{+}), \sqrt{x}τ′_{i}(x) ∈ L^{2}(ℝ^{+}) \\ ∫_{x}^{+∞} τ_{i}(u) du ∈ L^{2}(ℝ^{+}) \end{matrix}

in these hypotheses only the following invariant measure exists:

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

where b^{∗} is given by (1.12), b_{0} is given by (1.13), and Q_{∞} is given by (1.11).

The interest rates

2.1 The Musiela model in the Gaussian case

We suppose that the term structure of interest rates has the dynamics described from the following SPDE:

\{ \begin{matrix}dr_{t}(x) = (Ar_{t}(x) + c_{t}(x)) dt + τ_{t}(x) dW_{t} \\ r_{0}(·) ∈ H \end{matrix}  (2.1)

where A = \frac{∂}{∂x}, (W_{t})_{t} is a d-dimensional Brownian motion, H is a separable Hilbert space of functions from ℝ^{+} to ℝ (for example, L^{2}(ℝ^{+}) or H^{1}(ℝ^{+})), τ_{t}(x) = (τ_{t}^{1}(x), …, τ_{t}^{d}(x)) is a H^{d}-valued adapted process, and c(x) = ∑_{i=1}^{d} τ_{t}^{i}(x) ∫_{0}^{x} τ_{t}^{i}(u) du is a H-valued process. A solution for Eq. (2.1) has to be intended in the mild sense, as presented in [5].

If we make the hypothesis that τ is a deterministic function of x, and it belongs to H, then Eq. (2.1) has the unique mild solution:

r_{t}(x) = r_{0}(t + x) + ∫_{0}^{t} c(x + t − s) ds + ∫_{0}^{t} τ(x + t − s) dW_{s}  (2.2)

2.2 From the Musiela model to the Vasicek model

We suppose that the interest rates follow a Vasicek dynamics: this means that the spot rate r_{t}(0) is Markov and follows the dynamics:

dr_{t}(0) = a(b − r_{t}(0)) dt + σ dW_{t}  (2.3)

where (W_{t})_{t} is a real Brownian motion, and a, b, σ ∈ ℝ^{+}.

Theorem 2.1 In order for the spot rate to follow a Vasicek dynamic, it is necessary and sufficient that τ(·) and r_{0}(·) are of the form:

τ(x) = e^{−ax}σ, r_{0}(t) = e^{−at}r_{0}(0) + (1 − e^{−at})b  (2.4)

Proof. If the spot rate r_{t}(0) follows a Vasicek dynamics, then we have two possible expressions for it:

r_{t}(0) = e^{−at}r_{0}(0) + ∫_{0}^{t} e^{−a(t−s)}ab ds + ∫_{0}^{t} e^{−a(t−s)}σ dW_{s} =  (2.5)
= r_{0}(t) + ∫_{0}^{t} τ(t − s) dW_{s}

In order for the two expressions to give the same quantity, we must impose that τ(x) = e^{−ax}σ; moreover, r_{0}(·) can assume only the particular form given in the thesis. In fact:

r_{0}(t) = e^{−at}r_{0}(0) + ∫_{0}^{t} e^{−a(t−s)}ab ds = e^{−at}r_{0}(0) + [−be^{−as}]_{0}^{t} =
= e^{−at}r_{0}(0) + b(1 − e^{−at})

Conversely, if Eqs. (2.4) hold, then r_{t}(0) is given by (2.5), which is the solution of (2.3). Since the knowledge of the spot rate is sufficient to derive the entire term structure, whe have the thesis. ⋄

2.3 From the Musiela model to a state variables model

Now we suppose that the interest rates follow a particular state variables model: we suppose that all the instantaneous forward rates can be written in the form:

r_{t}(θ) = ∑_{i=1}^{n} α_{i}(θ) · r_{t}(θ_{i−1}) + β(θ) = α(θ) · R_{t} + β(θ)  (2.6)

where R_{t} = (r_{t}(0), r_{t}(θ_{1}), …, r_{t}(θ_{n−1})) are particular reference rates, with θ_{i} = iθ_{1}, and where R_{t} has the dynamics:

dR_{t} = (C − AR_{t}) dt + σ dW_{t}  (2.7)

where C ∈ ℝ^{n}, and A, σ ∈ ℝ^{n×n}. This corresponds to say that the only relevant quantities to know are the rates in R_{t}, and all the other rates can be derived from these through the weight coefficients α_{i}(θ). It can be proved (see [6]) that the matrix A satisfies:

e^{−Aθ_{1}} = ( \begin{matrix}0 & 1 & 0 & … & 0 \\ ⋮ & ... & ... & ... & ⋮ \\ ⋮ & & ... & ... & 0 \\ 0 & … & … & 0 & 1 \\ μ_{1} & … & … & … & μ_{n} \end{matrix} )

Theorem 2.2 In order for the rates (R_{t})_{t} to be state variables as in the model (2.6), it is necessary and sufficient that τ(·) and r_{0}(·) are of the form:

τ_{j}(x) = (e^{−Ax}σ)_{1j}  (2.8)
r_{0}(t) = e^{−At}R_{0}|_{1} + ∫_{0}^{t} e^{−A(t−s)}C|_{1} ds − ∫_{0}^{t} c(t − s) ds  (2.9)

Besides, the coefficients α(·) must be of the form:

α(θ) = (1, 0, …, 0)e^{−Aθ}

Proof. If the rates in R are state variables which satisfy (2.6) and (2.7), then we have two expressions for the solution:

r_{t}(θ_{i}) = e^{−At}R_{0}|_{i+1} + ∫_{0}^{t}(e^{−A(t−s)}C)|_{i+1} ds + (∫_{0}^{t} e^{−A(t−s)}σ dW_{s}) |_{i−1} = (2.10)
= r_{0}(θ_{i} + t) + ∫_{0}^{t} c(θ_{i} + t − s) ds + ∫_{0}^{t} τ(θ_{i} + t − s) dW_{s}

where the symbol x|_{i} indicates the i-th component of the vector x ∈ ℝ^{n}. In order for the two expressions to give the same quantity, we must impose that τ(θ_{i} + x) = (e^{−ax}σ)|_{i}, that is:

τ_{j}(x) = (e^{−Ax}σ)_{1j}

Now, let’s consider the rates with a maturity different from θ_{i}: we have

r_{t}(θ) = r_{0}(θ + t) + ∫_{0}^{t} c(θ + t − s) ds + ∫_{0}^{t} τ(θ + t − s) dW_{s} =
= α(θ)e^{−At}R_{0} + ∫_{0}^{t} α(θ)e^{−A(t−s)}C ds + ∫_{0}^{t} α(θ)e^{−A(t−s)}σ dW_{s} + β(θ)

From this expression, we can derive the right weight coefficients:

τ(θ + t − s) = α(θ)e^{−A(t−s)}σ = e^{−A(θ+t−s)}σ|_{1}

so:

α(θ) = (1, 0, …, 0)e^{−Aθ}

Conversely, if τ and r_{0} satisfy (2.8) and (2.9), then R is given by (2.10), which is the solution of (2.7). Once again, since the knowledge of the spot rate is sufficient to derive the entire term structure, we have the thesis. ⋄

We notice that in this model, as in the Vasicek model, r_{0}(·) can assume only particular forms, though if it can have a more complex structure than in the Vasicek model.

2.4 Considerations on the Hilbert space H

If we want to use the Musiela model in a state variable framework, then from the comparison of the models above we can derive qualitative considerations on the Hilbert space H and on its norm:

• the best kind of norm is an “homogeneous” one, in which no particular maturity has a particular weight (for example, the norms defined by the measures dx or γe^{−γx} dx, with γ > 0). In fact, while for the state variables model the relevant maturities are θ_{i}, for the Musiela model all the maturities are relevant, because of the particular semigroup involved in the equation (the translation semigroup S_{t}f = f(t + ·)): if the norm were not homogeneous, the semigroup would move this weight in the Musiela model in a way that would be incoherent with the way the state variables model behave.

• if we have a infinite horizon (that is, if we search a solution for all t ∈ ℝ^{+}), then we need all the forward curve r_{0}(t), t ∈ ℝ^{+}; we notice that a flat structure of interest rates r_{0}(t) = r ∀t doesn’t belong to L^{2} (or H^{1}), while it belongs to L_{γ}^{2} (or H_{γ}^{1}). Moreover, we have to choose well γ: in fact, for the equation to have sense, we must assume this condition on γ:

γ > ∥e^{−Aθ_{1}}∥_{∞} = max{ℜλ | λ eigenvalue of e^{−Aθ_{1}}}

• if we have a finite horizon (t ∈ [0, T]), then we only need the forward curve r_{0}(t), t ∈ [0, 2T], and the space L^{2} (as H^{1}) allows a flat structure of interest rates (such as the spaces L_{γ}^{2}, H_{γ}^{1}), so we don’t have particular limitations.

2.5 Structuring the Musiela model with respect to state variables

Now we want to impose on the Hilbert space a structure that gives a good interpretation in terms of the state variables models seen above.

2.5.1 The Vasicek model

First of all we notice that the particular τ found for the Musiela model is stable with respect to a translation; this tells us that the noise in the Hilbert space will be always in the component [e^{−au}], where the symbol [·] indicates the linear subspace generated by the argument in H.

We can reduce the Musiela model to the Vasicek model in this way: we can find a complete orthonormal basis (e_{i})_{i} in H, such that [e_{1}(u)] = [e^{−au}], and [e_{1}(u), e_{2}(u)] = [e^{−au}, r_{0}]: so with the first two components of our Hilbert space we can represent all the structure given by the Vasicek model. Then we can define a Hilbert seminorm on H that gives weight 1 to e_{1} and e_{2}, and weight 0 to all the other components:

||| · ||| = ⟨e_{1}, ·⟩ + ⟨e_{2}, ·⟩

If we quotientate H with respect to the seminorm ||| · |||, we obtain a new Hilbert space H̄, endowed with the norm ||| · |||, which is isomorphic to ℝ^{2}; then the initial data r_{0} can be decomposed in its component in H̄.

2.5.2 The general state variables model

We notice that, if A is diagonalizable, then the particular τ_{j} found for the Musiela model in this case are linear combinations of exponential functions:

τ_{j}(u) = ∑_{i=1}^{N} a_{ij}e^{−λ_{i}u}, where (λ_{i})_{i} are the eigenvalues of A

so the τ_{j} are stable with respect to a translation; this tells us that the noise in the Hilbert space will be always in the subspace [e^{−λ_{i}u}]_{i=1,...,N}.

As before, we can reduce the Musiela model to the state variables model in this way: we can find a complete orthonormal basis (e_{i})_{i} in H, such that

[e_{i}(u)]_{i=1,...,n} = [e^{−λ_{i}u}]_{i=1,...,n} ∀n ≤ N
[e_{i}(u)]_{i=1,...,N+k} = [e^{−λ_{1}u}, …, e^{−λ_{N} u}, e^{−2λ_{1}u}, …, e^{−2λ_{k}u}] ∀k = 1, …, N
[e_{i}(u)]_{i=1,...,2N+1} = [e^{−λ_{1}u}, …, e^{−λ_{N} u}, e^{−2λ_{1}u}, …, e^{−2λ_{N} u}, r]

We add the functions e^{−2λ_{i}u} in order to be able to represent the vector c(u).

We define a Hilbert seminorm on H that gives weight 1 to e_{i}, i = 1, …, 2N + 1, and weight 0 to all the other components:

||| · ||| = ∑_{i=1}^{2N+1}⟨e_{i}, ·⟩

If we quotientate H with respect to the seminorm ||| · |||, we obtain a new Hilbert space H̄, endowed with the norm ||| · |||, which is isomorphic to ℝ^{2N+1}; then the initial data r_{0} can be decomposed in its components in H̄. Moreover, the covariance σ of the rates R_{t} influences τ and the a_{i}.

2.6 Approximating the Musiela model with state variables

Now we suppose that the term structure of interest rates has the dynamics described from the following SPDE:

\{ \begin{matrix}dr_{t}(x) = (Ar_{t}(x) + c(x)) dt + τ(x) dW_{t} \\ r_{0}(·) ∈ H (= L_{(γ)}^{2}, H_{(γ)}^{1}) \end{matrix}  (2.11)

where:

• A = \frac{∂}{∂x},

• the object τ(x) dW_{t} represents the integration with respect to

∑_{n=1}^{∞} τ_{n}e_{n}(·)w_{t}^{n}

where ∑_{n=1}^{∞} τ_{n}^{2} < +∞, {(w_{t}^{n})_{t}}_{n} are real independent Brownian motions, (e_{n})_{n} is a complete orthonormal basis of H such that:

[e_{i}(u)]_{i=1,...,n} = [e^{−λ_{i}u}]_{i=1,...,n} ∀n ∈ ℕ

where the λ_{i} are particular numbers (see section 7),

• c is given by:

c(x) = ∑_{n=1}^{∞} τ_{n}e_{n}(x) ∫_{0}^{x} τ_{n}e_{n}(u) du

then the initial data r_{0} can be decomposed on the ONB (e_{n})_{n}:

r_{0}(x) = ∑_{n=1}^{∞} r_{n}e_{n}(x)

If r_{k} = τ_{k} = 0 ∀k ≥ n, then we obtain the state variables model seen before.

Now we suppose to modelize (making a mistake) the forward rates with the model given by Eq. (2.7). Then we try to approximate the term structure given by the Musiela model by a term structure (r̃_{t})_{t} in this way:

r̃_{t}(θ_{i}) = R_{t}^{i+1} ∀i = 0, …, n − 1
r̃_{t}(θ) = ∑_{i=1}^{n} α(θ)r̃_{t}(θ_{i}) + β(θ) = α(θ) · R_{t} + β(θ) ∀θ ≠ θ_{i}

Now we want to compare the “real” rates (r_{t})_{t} with the approximating ones (r̃_{t})_{t}.

Theorem 2.3 If we build σ, A, R_{0} and C such that

σ_{ij} = τ_{j}(θ_{i})  (2.12)
−AR_{0} + C = (r′_{0}(θ_{i}) + c(θ_{i}))_{i}  (2.13)

then the error r_{t}(θ) − r̃_{t}(θ) goes to 0 as N → ∞ ℚ-a.s. ∀t, θ.

Proof. The relations above are obtained “inverting” the relations of Theorem 2.2; in fact, relation (2.8) is equivalent to τ_{j}(θ_{i} + x) = (e^{−Ax}σ)_{ij}, and letting x = 0, we obtain relation (2.12). Relation (2.9) is equivalent to

e^{−At}R_{0}|_{i+1} + ∫_{0}^{t} e^{−As}C|_{i+1} ds = r_{0}(θ_{i} + t) + ∫_{0}^{t} c(θ_{i} + s) ds

Deriving these relations with respect to t we obtain:

−Ae^{−At}R_{0} + e^{−At}C = (r′_{0}(θ_{i} + t) + c(θ_{i} + t))_{i}

which, letting t = 0, gives us relation (2.13). Under these conditions, the error we make on the rates with maturity θ_{i} is:

r_{t}(θ_{i}) − r̃_{t}(θ_{i}) = r_{0}(θ_{i} + t) + ∫_{0}^{t} c(θ_{i} + t − s) ds + ∫_{0}^{t} τ(θ_{i} + t − s) dW_{s} +

−e^{−At}R_{0}|_{i−1} − ∫_{0}^{t}(e^{−A(t−s)}C)|_{i−1} ds − (∫_{0}^{t} e^{−A(t−s)}σ dW_{s}) |_{i−1} =

= ∑_{n=N+1}^{∞} r_{n}e_{n}(θ_{i} + t) + ∫_{0}^{t} ∑_{n=N+1}^{∞} τ_{n}^{2}e_{n}(θ_{i} + t − s) ∫_{0}^{θ_{i}+t−s} e_{n}(u) du +

+ ∫_{0}^{t} ∑_{n=N+1}^{∞} τ_{n}e_{n}(θ_{i} + t − s) dw_{s}^{n}

and for the others is:

r_{t}(θ) − r̃_{t}(θ) = r_{0}(θ + t) + ∫_{0}^{t} c(θ + t − s) ds + ∫_{0}^{t} τ(θ + t − s) dW_{s} − β(θ)+

− ∑_{i=1}^{N} α_{i}(θ) (e^{−At}R_{0}|_{i−1} − ∫_{0}^{t}(e^{−A(t−s)}C)|_{i−1} ds − (∫_{0}^{t} e^{−A(t−s)}σ dW_{s}) |_{i−1}) =

= ∑_{n=N+1}^{∞} r_{n}e_{n}(θ_{i} + t) + ∫_{0}^{t} ∑_{n=N+1}^{∞} τ_{n}^{2}e_{n}(θ_{i} + t − s) ∫_{0}^{θ_{i}+t−s} e_{n}(u) du +

+ ∫_{0}^{t} ∑_{n=N+1}^{∞} τ_{n}e_{n}(θ_{i} + t − s) dw_{s}^{n}

this doesn’t surprise us, because we can see the state variables model as a particular Musiela model driven by a finite dimensional Brownian motion. The consequence is that the error |r_{t}(θ) − r̃_{t}(θ)| goes to 0 as N goes to infinity, and the speed of convergence depends from the coefficients (r_{n})_{n} and (τ_{n})_{n}. ⋄

2.7 Open problems

• how to choose the “right” numbers λ_{i}? In the state variables model, they were the eigenvalues of A, so a natural choice would be that the λ_{i} are the eigenvalues of the matrix A in the “wrong” model.

• for a given orthonormal basis e_{1}, e_{2}, … in a given Hilbert space H (for example L^{2}(ℝ^{+}), L_{γ}^{2}(ℝ^{+}), H^{1}(ℝ^{+}), H_{γ}^{1}(ℝ^{+}), H^{1}([0, T]), etc.), analyse the “distance” between r_{t}, t ≤ T and r_{t}^{N}, t ≤ T, where

r_{t}^{N} = ∑_{n=1}^{N}⟨r_{t}, e_{n}⟩e_{n}, t ≤ T

The distance can be for example

1. sup_{t≤T} 𝔼 [∥r_{t} − r_{t}^{N}∥_{H}^{2}]
2. 𝔼 [sup_{t≤T} ∥r_{t} − r_{t}^{N}∥_{H}^{2}]
3. sup_{t≤T} ℙ {∥r_{t} − r_{t}^{N}∥_{H}^{2} ≥ δ}

4. the prices of some interest rate derivatives (such as call options on zero coupon

bonds, caps, swaptions or futures contracts) The distance will depend on the choice of the basis e_{1}, e_{2}, … and on the coefficients of the equation c, τ, that is, due to the particular structure of the model, only on τ.

• analyse the distance with respect to the choice of the basis; consider specific basis, such as trigonometric basis, Haar basis, Karhunen-Loeve basis or other ones.

• choose basis which is “good” for the volatility function τ with respect to observed

datas. Still are to be chosen:

• T > 0 (even if we could fix T = 10, 20 or 30 years)

• the basis e_{1}, e_{2}, …

• the number N of finite components (a “good” number would be N = 3 or 4) Instead of a separable Hilbert space H, one could consider Schauder basis in the space C([0, T]) and do a similar analysis.

Calibrating the model

In this chapter we try to calibrate the model to the market datas, that is to the assets traded in the market. Tipically market datas are the prices of some coupon bonds, from which one can derive prices of zero-coupon bonds with typical maturities T = j · 3 months, j = 1, …, N (for example N = 40, that is 10 years), and prices of some interest rate derivatives, such as futures or options on some bonds, caps, floors and swaptions.

3.1 An “implied” approach

Now we try to find the diffusion function τ implicitly from the prices of the caps in the market. We will find some general finite dimensional constraints for τ, which it must satisfy in order to price the caps.

A cap is the sum of several caplets, each one having a different maturity T_{i} = i·θ, i = 1, …, N, where θ usually is equal to 0.25 or 0.5 (three or six months), and the same exercise rate.

A caplet is a claim giving the right to exchange the linear rate R_{L}(T_{i−1}, θ), defined by

1 + θR_{L}(T_{i−1}, θ) = \frac{1}{B(T_{i−1}, T_{i})}

with the fixed rate R^{f}, multiplied by θ, in the maturity T_{i}; hence the caplet gives the payoff

(θR_{L}(T_{i−1}, θ) − θR^{f})^{+}

at time T_{i}, so it has the price:

Caplet_{i}(t) = 𝔼_{ℚ}^{t} [e^{∫_{t}^{T_{i}} r_{u} du} (θR_{L}(T_{i−1}, θ) − θR^{f})^{+}]

where we indicate with 𝔼_{ℚ}^{t} the conditional expectation with respect to the probability ℚ and the σ-algebra F_{t}. We notice that the caplet will be exercised if and only if it gives a positive payoff, that is if and only if we are in the so called exercise set

E_{i} = \{\frac{B(T_{i−1}, T_{i})}{B(T_{i}, T_{i})} ≤ \frac{1}{1 + θR^{f}} \}

We also notice that the set E_{i} is F_{T_{i−1}}-measurable, so one can know in T_{i−1} whether to exercise the cap in T_{i} or not. Now, introducing the forward neutral probabilities ℚ_{T_{i−1}} and ℚ_{T_{i}} (see [7]), we calculate explicitly the price of a caplet:

Caplet_{i}(t) = 𝔼_{ℚ}^{t} [e^{∫_{t}^{T_{i}} r_{u} du} (\frac{1}{B(T_{i−1}, T_{i})} − (1 + θR^{f}))^{+}] =
= 𝔼_{ℚ}^{t} [e^{∫_{t}^{T_{i}} r_{u} du} \frac{1}{B(T_{i−1}, T_{i})}1_{E_{i}} − (1 + θR^{f})1_{E_{i}}] =
= 𝔼_{ℚ}^{t} [e^{∫_{t}^{T_{i−1}} r_{u} du}𝔼_{ℚ}^{T_{i−1}} [e^{∫_{T_{i−1}}^{T_{i}} r_{u} du}] \frac{1}{B(T_{i−1}, T_{i})}1_{E_{i}}] −
−(1 + θR^{f})𝔼_{ℚ}^{t} [e^{∫_{t}^{T_{i}} r_{u} du}1_{E_{i}}] =
= 𝔼_{ℚ}^{t} [e^{∫_{t}^{T_{i−1}} r_{u} du}1_{E_{i}}] − (1 + θR^{f})B(t, T_{i})ℚ_{T_{i}}^{t}(E_{i}) =
= B(t, T_{i−1})ℚ_{T_{i−1}}^{t}(E_{i}) − (1 + θR^{f})B(t, T_{i})ℚ_{T_{i}}^{t}(E_{i})

Now we suppose that

\frac{dB(t, T)}{B(t, T)} = r_{t} dt + Γ(t, T) dW_{t}

for all T = T_{i}, where W is a d-dimensional (infinite dimensional?) Brownian motion under ℚ, and Γ is deterministic; this is equivalent to suppose that

dr_{t}(x) = (Ar_{t}(x) + c(t, x)) dt + τ(t, x) dW_{t}

with τ deterministic function of t and x such that

τ(t, x) = −\frac{∂}{∂T} Γ(t, t + x), that is Γ(t, T) = − ∫_{0}^{T −t} τ(t, u) du  (3.1)

Then the price of a caplet is given by:

Caplet_{i}(t) = B(t, T_{i−1})N(d_{0,i}) − B(t, T_{i})(1 + θR^{f})N(d_{1,i})

where

d_{0,i} = \frac{1}{σ_{t,T_{i}}} log (\frac{B(t, T_{i−1})}{(1 + θR^{f})B(t, T_{i})}) + \frac{1}{2}σ_{t,T_{i}}, d_{1,i} = d_{0,i} − σ_{t,T_{i}}

and

σ_{t,T_{i}}^{2} = ∫_{t}^{T_{i−1}} ∥Γ(s, T_{i}) − Γ(s, T_{i−1})∥^{2} ds

where ∥ · ∥ is the norm on ℝ^{d}.

We notice that this formula is very similar to the Black-Scholes formula; in particular, Caplet_{i}(t) it is an increasing function of σ_{t,T_{i}}. This means that if one knows the price of a caplet and all the variables in the price, except σ_{t,T_{i}}, one can derive it inverting the function.

The problem is that in the market caplets are not priced, but one can read only the price of caps, which are given by

Cap(t) = ∑_{i=1}^{N} Caplet_{i}(t)

so it is not well clear how to derive the σ_{t,T_{i}} for the single caplet. One method could be to ricave a volatility for each cap priced in the market, and then to make a linear interpolation for the volatility of each caplet; this method is presented in [3] as being the market practise.

Once that we have the quantities σ_{t,T_{i}}, we obtain a constraint for the function τ:

\frac{∂}{∂t}σ_{t,T_{i}}^{2} = −∥Γ(t, T_{i}) − Γ(t, T_{i−1})∥^{2} = − ∥∫_{T_{i−1}−t}^{T_{i}−t} τ(t, u) du∥^{2}

however, these are only finite dimensional constraints, and it is still to be specified:

• the “dimension” of the model, that is the number d of Brownian motions to use

• the shape of the functions τ: for example, in [3] the authors use cubic splines, but say that the market practise is linear interpolation.

3.2 Two “historical” approaches

Now we present two historical approaches, each one based on the observation of the prices of zero-coupon bonds, which we suppose to have dynamics

\frac{dB(t, T)}{B(t, T)} = r_{t}(0) dt + Γ(t, T) dW_{t}

To this purpose we define R(t, T) and R̃(t, T) in this way:

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

then, since R̃(t, T) = log B(t, T), we have:

dR̃(t, T) = (r_{t}(0) + 1/2Γ^{2}(t, T)) dt + Γ(t, T) dW_{t}

and since R(t, T) = \frac{R̃(t,T)}{T −t}, we have:

dR(t, T) = −\frac{R̃(t, T)}{(T − t)^{2}} dt + \frac{dR̃(t, T)}{T − t} \frac{d⟨R̃(·, T)⟩_{t}}{2(T − t)} =
= −\frac{R(t, T)}{T − t} dt + \frac{r_{t}(0) + 1/2Γ^{2}(t, T)}{T − t} dt + \frac{Γ(t, T)}{T − t} dW_{t}

3.2.1 With less datas

In this approach we suppose to have datas in the form B(t, T), T = j·3 months, j = 1, …, N (for example N = 40, that is 10 years), t = i days, i = 1, …, 90. We can get the rates in this way:

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

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

(we can get only approximate values)

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 90 samples we have a discrete scheme with 90 · N values. If we interpolate each curve (which way? linear? spline? exponential?), then we could have 90 samples of forward rates’ curves. Then we can use Karhunen-Loewe method to find a “good” basis, and try to remake the interpolation using Karhunen-Loewe basis...

3.2.2 With more datas

In this approach we suppose to have datas in the form B(t, T_{i}), T_{i} = i·3 months, i = 1, …, N (for example N = 40, that is 10 years), and datas in t with a difference in t very near (for example, if we have tick by tick datas: this datas are taken in each transaction made in the stock market, and tipically the difference between two t’s near is 5-10 seconds). From these datas, we can get the volatility Γ of the bonds

dB(t, T) = B(t, T)(r_{t} dt + Γ(t, T) dW_{t})

where Γ is a d-dimensional vector.

A possible way is the following: we fix the j-th day, and let t vary only into the j-th day. We suppose that Γ(t, T_{i}) ≃ Γ(T_{i}) for t ∈ [i day, i+1 day], and we pose X_{t}^{i} = log B(t, T_{i}), so we have

X_{t}^{i} = X_{0}^{i} + (r_{t} − 1/2Γ^{2}(T_{i}))t + Γ(T_{i})W_{t}

Now we discretize [i day, i + 1 day] with a step h = 1 day/n, and we pose:

σ̂_{n}^{lm} = ∑_{i=0}^{n−1}(X_{(i+1)h}^{l} − X_{ih}^{l})(X_{(i+1)h}^{m} − X_{ih}^{m})

Then we have:

σ̂_{n}^{lm} = ∑_{i=0}^{n−1} ((r_{t} − 1/2Γ^{2}(T_{l}))h + Γ(T_{l})(W_{(i+1)h} − W_{ih})) ·
· ((r_{t} − 1/2Γ^{2}(T_{m}))h + Γ(T_{m})(W_{(i+1)h} − W_{ih})) =
= (r_{t} − 1/2Γ^{2}(T_{l}))(r_{t} − 1/2Γ^{2}(T_{m})) ∑_{i=0}^{n−1} h^{2} + Γ(T_{l})Γ(T_{m}) ∑_{i=0}^{n−1}(W_{(i+1)h} − W_{ih})^{2} +

+ ((r_{t} − 1/2Γ^{2}(T_{l}))Γ(T_{m}) + (r_{t} − 1/2Γ^{2}(T_{m}))Γ(T_{l})) h ∑_{i=0}^{n−1}(W_{(i+1)h} − W_{ih}) =

= (r_{t} − 1/2Γ^{2}(T_{l}))(r_{t} − 1/2Γ^{2}(T_{m}))\frac{t^{2}}{n} + Γ(T_{l})Γ(T_{m}) ∑_{i=0}^{n−1}(W_{(i+1)h} − W_{ih})^{2} +

+(r_{t} − 1/2Γ(T_{l})Γ(T_{m}))(Γ(T_{l}) + Γ(T_{m}))\frac{t}{n}W_{t}

As n tends to infinity, the first and the third addend tend to zero, while the second tends to Γ(T_{l})Γ(T_{m})t; so the estimator σ̂_{n} = (σ̂_{n}^{lm})_{lm} tends to Γ · Γ^{t}, where Γ is the matrix having as rows Γ(T_{l}), l = 1, …, d. We can then get Γ making the square root of the resulting matrix.

Once we have Γ(t, T), we get:

τ(t, x) = −\frac{∂}{∂T} Γ(t, t + x)

(as before, we can get only approximate values)

3.3 Karhunen–Loeve expansion

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

We suppose we have a random variable u with values in a function space H, dual of another space K.

Example 3.1 In our case, we could take H = K = L^{2}(ℝ^{+}), H^{1}(ℝ^{+}), L^{2}([0, T]) or H^{1}([0, T]).

From now on, we suppose that 𝔼[u] = 0 (we can always reduce the problem to this case by defining u′ = u − 𝔼[u]). We suppose then that the bilinear function defined by the variance R of u from H^{2} into ℝ^{+}

R(k, k) = 𝔼[⟨u, k⟩^{2}] ∀k ∈ K

is continuous; then, for the kernel representation theorem, we can find a function R(t, t′) such that

R(k, k) = ∫ R(t, t′)k(t)k(t′) dt dt′

then we can interpret R as an operator going from K to H.

3.3.1 The representation problem

Now we want to maximize the expression

λ = \frac{𝔼[⟨u, k⟩^{2}]}{⟨k, k⟩} = \frac{R(k, k)}{⟨k, k⟩}

in k ∈ K ∩ H. From a standard procedure by calculus of variations, the extremal points are the k such that

R(k, k′) = λ⟨k, k′⟩ ∀k′ ∈ K ∩ H

that is, the extremal k are the eigenvectors of R:

Rk = λk

3.3.2 The “finite energy” case: Karhunen–Loeve expansion

Now we suppose that (f_{n})_{n} is a complete (non necessary orthonormal) system in K (in L^{2}), decreasing exponentially to infinity, and that

∑_{n=1}^{∞} 𝔼[⟨u, f_{n}⟩^{2}] < +∞

this is equivalent to ∑_{n=1}^{∞}⟨u, f_{n}⟩^{2} < +∞ ℙ-a.s.; then u ∈ L^{2} ℙ-a.s.; besides, R is symmetric and nonnegative-definite, so it exists a complete orthonormal set (k_{n})_{n} such that Rk_{n} = λ_{n}k_{n}, with λ_{n} ≥ 0.

3.3.3 Calculation of the eigenvalues of R

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} ∈ K, and we define:

k^{n+1} = Rk^{n}

then we have that:

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

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

3.3.4 Rate of convergence of the series

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

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

this means that at least N terms of the series are needed in order that the sum is significative. If we suppose that L is large enough to include in [0, L]^{2} “most” of the support of the function R(·, ·), then we have:

λ_{1} ≥ \frac{1}{L} ∫_{0}^{L} ∫_{0}^{L} R(t, t′) dt dt′

and we have the estimation:

N ≤ L\frac{∫_{0}^{∞} R(t, t) dt}{∫_{0}^{L} ∫_{0}^{L} R(t, t′) dt dt′}

3.3.5 Application: the Musiela model

Now we want to apply this expansion to the process (r_{t})_{t} in this way: we suppose that (r_{t})_{t} is a stationary solution of the Musiela equation:

\{ \begin{matrix}dr_{t}(x) = (Ar_{t}(x) + c_{t}(x)) dt + τ_{t}(x) dW_{t} \\ r_{0}(·) ∈ L^{2}(Ω, F, ℙ; H) \end{matrix}

where τ is a priori an element of H^{∞} such that ∑_{n=1}^{∞} ∥τ_{n}∥^{2} < +∞. Then the marginal law of r_{t} is

ℙ(r_{t}) = N(b^{∗}(·) + b_{0}, Q_{∞}) ∀t ≥ 0

with b_{0} ∈ ℝ and

Q_{∞} = ∫_{0}^{∞} ∑_{n=1}^{∞} S_{u}τ_{n} ⊗ S_{u}τ_{n} du, b^{∗}(·) = − ∫_{0}^{x} c(u) 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: in fact we always studied the Musiela model under the risk-neutral probability ℚ, while we can obtain from law of the time series an estimation of the hystorical 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 troubles.

Now we suppose to have N̄ observations r_{t_{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 N̄ observations r_{t_{i}}, i = 1, …, N̄ as N̄ realizations of the random variable r_{t} for a given t ≥ 0.

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

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

In order to implement the algorithm, we follow these procedure:

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_{∞}), that is R = Q_{∞}.

2nd step we apply the algorithm: we fix an initial k^{0} and we define:

k^{n+1} = \frac{1}{N̄} ∑_{i=1}^{N̄}⟨r_{i}, k^{n}⟩r_{i}

and we find λ_{1} and k_{1}.

3rd step we define R_{1}:

R_{1} = \frac{1}{N̄} ∑_{i=1}^{N̄} r_{i} ⊗ r_{i} − λ_{1}k_{1} ⊗ k_{1}

and we come back to the 2nd step (we can stop to the suitable N).

Once we obtain the (k_{n})_{n}, we can obtain the (τ_{n})_{n} in this way: we have two expressions of the covariance Q_{∞}:

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

each term in the two sums corresponding to a different component of the Brownian motion W. Since the components of W are independent, then we have:

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

so we have:

∫_{0}^{∞} τ_{n}(x + u) du = \sqrt{λ_{n}}k_{n}(x) ∀x ∈ ℝ^{+}, n ∈ ℕ

Deriving under the integral sign, we obtain

\sqrt{λ_{n}}k′_{n}(x) = ∫_{0}^{∞} τ′_{n}(x + u) du = [τ_{n}(x + u)]_{0}^{+∞} = −τ_{n}(x) ∀x ∈ ℝ^{+}, n ∈ ℕ

because lim_{x→+∞} τ_{n}(x) = 0; so we have that

τ_{n} = −\sqrt{λ_{n}}k′_{n}

Robustness of the Musiela model

4.1 A functional differential equation for interest rates derivatives

We start from the Musiela equation

\{ \begin{matrix}dr_{t}(x) = (Ar_{t}(x) + c(t, r_{t}(·))(x) dt + τ(t, r_{t}(·))(x) dW_{t}, t ∈ [0, T] \\ r_{0}(·) ∈ H \end{matrix}  (4.1)

where A = \frac{∂}{∂x}, (W_{t})_{t} is a d-dimensional Brownian motion, H is a separable Hilbert space of functions from ℝ^{+} to ℝ (for example, L^{2}(ℝ^{+}) or H^{1}(ℝ^{+})), τ : ℝ^{+} × H → H^{d} and c : ℝ^{+} × H → H, defined by c(t, r(·))(x) = ∑_{i=1}^{d} τ_{i}(t, r(·))(x) ∫_{0}^{x} τ_{i}(t, r(·))(u) du, satisfy the following conditions:

(A1) the mapping c : R_{+} × H → H is measurable and

∥c(t, r) − c(t, r′)∥ ≤ K∥r − r′∥ ∀r, r′ ∈ H, ∀t ∈ [0, T]

(A2) the mapping τ : R_{+} × H → H^{d} is measurable and

∥τ(t, r) − τ(t, r′)∥ ≤ K∥r − r′∥ ∀r, r′ ∈ H, ∀t ∈ [0, T]

These conditions comprehend the stationary Gaussian case we analyzed in Chapter 1 as the particular case in which τ is a constant function.

Proposition 4.1 If (A1) and (A2) hold then for every H-valued random variable r_{0} independent of the Wiener process W there exists a unique continuous solution (r_{t})_{t} to (1.1) which is a strong Markov process on H. Moreover, if 𝔼[∥r_{0}∥^{p}] < ∞ for a certain p ≥ 2 then for any T > 0:

sup_{t≤T} 𝔼[∥r_{t}∥^{p}] ≤ C_{T,p} (1 + 𝔼[∥r_{0}∥^{p}]) .

Moreover, if τ and c are Lipschitz transformations on H^{1}(0, ∞) and r_{0} ∈ H^{1}(0, ∞) then r_{t} ∈ H^{1}(0, ∞) ℚ-almost surely for all t ≥ 0.

The proof is an immediate consequence of the general theory developed in [5].

From now on we assume that our system is homogeneous in time, that is c and τ don’t depend explicitly on time, but only on the rate r ∈ H. We know that the bonds’ prices are deterministic functions of the spot forward rates:

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

More in general, we can consider a zero coupon bond as the simplest case of an European contingent claim. We say that F : H → ℝ defines an European interest rate derivative or European interest rate contingent claim with maturity T if 𝔼[F (r_{T})] < +∞; the claim gives to his owner the final payoff F (r_{T}). In particular the case F ≡ 1 corresponds to the zero-coupon bond B(t, T) with maturity T. Due to the Markov property of the process r the price C_{T} of the claim F at time t ≤ T is

C_{T} (t) = 𝔼 [exp (− ∫_{t}^{T} r_{u}(0) du) F (r(T))| F_{t}]
= 𝔼 [exp (− ∫_{t}^{T} r_{u}(0) du) F (r(T))| r(t)] .  (4.2)

As in the finite dimensional case, we have the following result, that links the European contingent claim F with the equation

\{ \begin{matrix}\frac{∂u}{∂t}(t, r) + Lu(t, r) − r(0)u(t, r) = 0 \\ u(T, r) = F (r). \end{matrix}  (4.3)

where L is defined on D(A) by

LF (r) = \frac{1}{2}⟨D^{2}F (r)τ(r), τ(r)⟩ + ⟨Ar + c(r), DF (r)⟩  (4.4)

Equation (4.3) is a Feynman-Kac equation in the Hilbert space H. To make our statements more precise, for m ≥ 0 we define the space

C_{m}(H) = \{F : H → R | \frac{F (·)}{1 + ∥ · ∥^{m}} bounded and uniformly continuous on H\}

The space C_{m}(H) is endowed with the norm

∥F ∥_{m} = sup_{r∈H} (\frac{|F (r)|}{1 + ∥r∥^{m}} ) .

Let C_{m}^{n}(H) denote the subspace of C_{m}(H) containing all functions F which are n times Fréchet differentiable on H and

sup_{r∈H} (\frac{∥D^{k}F (r)∥}{1 + ∥r∥^{m}}) < ∞, k = 1, …, n,

where D^{k}F (ϕ) denotes the k-th Fréchet derivative of F. If F ∈ C_{m}^{1}(H) then the derivative DF (ϕ) of F at ϕ ∈ H can be identified with an element of H and therefore DF : H → H is a continuous mapping. If F ∈ C_{m}^{2}(H) then the second derivative D^{2}F (ϕ) is a symmetric linear operator in H and we can define a continuous mapping D^{2}F : H → L(H) from H into the algebra of bounded operators on H. Finally by C^{n}(H) we denote the vector space of all real-valued continuous functions on H with n continuous Fréchet derivatives. Let ψ_{1}, …, ψ_{n} ∈ dom (A^{∗}) and let P_{n} denote the orthogonal projection on the linear span H_{n} of the vectors ψ_{1}, …, ψ_{n}. Moreover, let

D_{0} = {F ∈ C_{α}(H) : F = f ∘ P_{n}, f ∈ C_{m}^{2}(R^{n}), n = 1, …}.

Note that every F ∈ C_{m}(H) is a pointwise limit of elements of D_{0}.

Theorem 4.2 If τ and c are twice differentiable on H and such that

∥τ(ϕ)∥ + ∥c(ϕ)∥ ≤ C(1 + ∥ϕ∥).

then for every F ∈ C_{m}^{2}(H) the function C_{T} is a unique solution of the Feynman-Kac equation (4.3) in the following sense.

i) The function C_{T} : R_{+} × H → R is bounded and continuous with respect to each variable.

ii) For every t ≥ 0 we have C_{T} (t, ·) ∈ C^{2}(H).

iii) C_{T} ∈ C^{1}([0, T], H^{1}).

iv) Equation (4.3) holds for every ϕ ∈ dom (A) and t ≥ 0. Moreover, C_{T} is given by (4.2).

For a proof, see [9].

This theorem is a first result, but it is not enough: in fact the principal examples of European interest rates derivatives (calls, puts, swaptions) are of class C^{0} but not of class C^{1}.

4.2 Robustness of the Musiela model

We suppose that the dynamics of the forward rates is given by Equation

\{ \begin{matrix}dr_{t}(x) = (Ar_{t}(x) + c̄_{t}(x)) dt + τ̄_{t}(x) dW_{t} \\ r_{0}(·) ∈ H \end{matrix}  (4.5)

where A = \frac{∂}{∂x}, (W_{t})_{t} is a d-dimensional Brownian motion, H is a separable Hilbert space of functions from ℝ^{+} to ℝ (for example, L^{2}(ℝ^{+}) or H^{1}(ℝ^{+})), τ̄_{t}(x) is a H^{d}-valued adapted process, and c̄(x) = ∑_{i=1}^{d} τ̄_{t}^{i}(x) ∫_{0}^{x} τ̄_{t}^{i}(u) du is a H-valued process; moreover, the process τ is supposed to take values in some compact set Σ ∈ H^{d}. This model is a quite general one, and it is an analogous of stochastic volatility models for stock prices analyzed in literature in the last years (see for example [8] and [1]).

Now we consider an operator in the market that wants to create and sell an European interest rate derivative asset with payoff F (rT̄) depending on the value in the maturity T̄ of all the forward rates curve rT̄. In many relevant cases, the payoff F will depend not on all the curve rT̄, but only on the values of N particular bonds with maturity T_{i} > T̄, so it will be of the form

F (B(T̄, T_{1}), …, B(T̄, T_{N})), T̄ < T_{1} < … < T_{N}  (4.6)

From now on, we suppose that F is of the form (4.6). In order to hedge this asset, the operator chooses to use the homogeneous Gaussian Musiela model described in Chapter 1. Then he (she) builds a self-financing strategy in this way: he (she) fixes the price C_{t} of the option at time t as

C_{t} = u(t, r_{t})

where u is the solution of the following nonlinear functional differential equation:

\{ \begin{matrix}\frac{∂u}{∂t}(t, r) + sup_{τ∈Σ} Lu(t, r) − r(0)u(t, r) = 0, (t, r) ∈ [0, T̄] × H \\ u(T̄, r) = F (r), r ∈ H \end{matrix}  (4.7)

where L is the operator (4.4). Equation (4.7) is a Hamilton-Jacobi-Bellmann type equation; we notice that its finite-dimensional versions appear for the first time in one dimension in [1], and it is generalized to the finite-dimensional case in [14]; as in the previous works Equation (4.7) will be called from now on the Black-Scholes-Barenblatt equation associated to the contingent claim F.

Lemma 4.3 If F is of the form (4.6), and the function u(t, r) is a solution of Equation (4.7) and is of the form

u(t, r) = C(t, r_{t}(0), B(t, T_{1}), …, B(t, T_{N}))

where C : [0, T̄] × ℝ × ℝ^{N}, r_{t}(0) is the spot rate at time t and B(t, T_{i}) are the values in t of the zero-coupon bonds maturing in T_{i}, i = 1, …, N, then C is solution of the following equation

\frac{∂C}{∂t} (t, r_{0}, B) + r_{0}⟨DC(t, r_{0}, B), B⟩ + \frac{1}{2} sup_{τ∈Σ} tr (D^{2}C(t, r_{0}, B)(B̄Γ)(B̄Γ)^{∗}) = r_{0}C(t, B),

C(T̄, r_{0}, B) = F (B), r_{0} ∈ ℝ, B ∈ ℝ^{N}

(4.8)

where DC(t, r_{0}, B) = (\frac{∂C}{∂B_{1}} (t, r_{0}, B), …, \frac{∂C}{∂B_{n}} (t, r_{0}, B)), and D^{2}C(t, r_{0}, B) = (\frac{∂^{2}C}{∂B_{i}B_{j}} (t, r_{0}, B))_{ij}, Γ(t, T) is a d × N matrix which the i-th row is the volatility Γ(t, T_{i}) of the bond B(t, T_{i}), given by Equation (3.1), and

B̄ = diag (B) = ( \begin{matrix}B_{1} & 0 & … & 0 \\ 0 & ... & ... & ⋮ \\ ⋮ & ... & ... & 0 \\ 0 & … & 0 & B_{n} \end{matrix} )

Proof. We notice that

\frac{∂u}{∂r} = \frac{∂C}{∂B} \frac{∂B}{∂r} , \frac{∂^{2}u}{∂r^{2}} = (\frac{∂B}{∂r} )^{∗} \frac{∂^{2}C}{∂B^{2}} (\frac{∂B}{∂r} ) + \frac{∂C}{∂B} \frac{∂^{2}B}{∂r^{2}}

and that for all maturities T_{i} we have:

\frac{∂B(t, T_{i})}{∂r} = −1_{[0,T_{i}−t]}B, \frac{∂^{2}B(t, T_{i})}{∂r^{2}} = 1_{[0,T_{i}−t]}^{∗}B1_{[0,T_{i}−t]}

then from Equation (4.7) we get:

r(0)u(t, r) = \frac{∂u}{∂t}(t, r) + max_{τ∈Σ} Lu(t, r) =

= \frac{∂u}{∂t}(t, r) + max_{τ∈Σ} (⟨\frac{∂u}{∂r}(t, r), Ar + c(r)⟩ + \frac{1}{2} ⟨\frac{∂^{2}u}{∂r^{2}}(t, r)τ(r), τ(r)⟩) =

= \frac{∂C}{∂t} (t, B) + max_{τ∈Σ} (− ∑_{i=1}^{N} \frac{∂C}{∂B_{i}} (t, B)B(t, T_{i})⟨1_{[0,T_{i}−t]}, Ar + c(r)⟩+

+\frac{1}{2} ∑_{i=1}^{N} τ(r) (1_{[0,T_{i}−t]}^{∗}B(t, T_{i})\frac{∂^{2}C}{∂B^{2}} (t, B)B(t, T_{i})1_{[0,T_{i}−t]}+

+\frac{∂C}{∂B} 1_{[0,T_{i}−t]}^{∗}B_{i}1_{[0,T_{i}−t]}) τ(r)) =

= \frac{∂C}{∂t} (t, B) + max_{τ∈Σ} (− ∑_{i=1}^{N} \frac{∂C}{∂B_{i}} (t, B)B(t, T_{i})(r(T_{i} − t) − r(0) + ⟨1_{[0,T_{i}−t]}, c(r)⟩)+

+\frac{1}{2}Tr (∑_{i=1}^{N} Γ(t, T_{i})^{∗}B(t, T_{i})\frac{∂^{2}C}{∂B^{2}} (t, B)B(t, T_{i})Γ(t, T_{i})+

+ Γ(t, T_{i})^{∗} \frac{∂C}{∂B_{i}} B(t, T_{i})Γ(t, T_{i}))) =

= \frac{∂C}{∂t} (t, B) + max_{τ∈Σ} (∑_{i=1}^{N} \frac{∂C}{∂B_{i}} (t, B)B(t, T_{i})r(0)+

+ \frac{1}{2}Tr (∑ i = 1^{N}Γ(t, T_{i})^{∗}B(t, T_{i})\frac{∂^{2}C}{∂B^{2}} (t, B)B(t, T_{i})Γ(t, T_{i})))

(− ∑_{i=1}^{N} \frac{∂C}{∂B_{i}} (t, B)B(t, T_{i})(r(T_{i} − t)))

and Equation (4.8) follows. ⋄

Equation (4.8) will be called from now on the reduced Black-Scholes-Barenblatt equation for F.

In order to hedge the claim F, the agent builds the hedging portfolio holding the quantity Δ_{t}^{i} of the zero-coupon bond B(t, T_{i}) in the hedging portfolio Π as

Δ_{t}^{i} = \frac{∂C}{∂B_{i}} (t, B(t, T))

where we indicated as B(t, T) the vector (B(t, T_{1}), …, B(t, T_{N})); Π follows the dynamics:

\{ \begin{matrix}dΠ_{t} = r_{t}(0)(Π_{t} − ⟨Δ_{t}, B(t, T)⟩) dt + ⟨Δ_{t}, dB(t, T)⟩ \\ Π_{0} = P_{0} \end{matrix}

(for a proof, see [11]).

We define the tracking error:

e_{t} = Π_{t} − C_{t}

with e_{0} = 0. Intuitively, the tracking error gives the error made by the operator in his (her) estimation, or better the difference between the hedging portfolio held by the agent and the option he (she) sold.

We want to find a superstrategy where e_{t} is increasing ∀t Q-a.s.

Theorem 4.4 If C is a C^{1,2}([0, T) × R^{n}) solution of Equation (4.8), then (e_{t})_{t} is an increasing process.

Proof. The evolution of (e_{t})_{t} is the following:

de_{t} = r_{t}(0)(Π_{t} − ⟨Δ_{t}, B(t, T)⟩) dt + ⟨Δ_{t}, dB(t, T)⟩ − \frac{∂C}{∂t} (t, B(t, T)) dt −

−⟨DC(t, B(t, T)), dB(t, T)⟩ −

−\frac{1}{2}tr (D^{2}C(t, B(t, T))(B̄(t, T)Γ(t, T))(B̄(t, T)Γ(t, T))^{∗}) dt =

= r_{t}(0)(Π_{t} − ⟨Δ_{t}, B(t, T)⟩) dt + r_{t}(0)⟨B(t, T), DC(t, B(t, T))⟩ dt − r_{t}(0)C_{t} dt +

+\frac{1}{2} max_{τ∈Σ} tr (D^{2}C(t, B(t, T))(B̄(t, T)Γ(t, T))(B̄(t, T)Γ(t, T))^{∗}) dt −

−\frac{1}{2}tr (D^{2}C(t, B(t, T))(B̄(t, T)Γ(t, T))(B̄(t, T)Γ(t, T))^{∗}) dt =

= r_{t}(0)e_{t} dt + \frac{1}{2} max_{τ∈Σ} tr (D^{2}C(t, B(t, T))(B̄(t, T)Γ(t, T))(B̄(t, T)Γ(t, T))^{∗}) −

−\frac{1}{2}tr (D^{2}C(t, B(t, T))(B̄(t, T)σ(t, T))(B̄(t, T)σ(t, T))^{∗}) =

= r_{t}(0)e_{t} dt + R_{t} dt

where

R_{t} = \frac{1}{2} max_{τ∈Σ} tr (D^{2}C(t, B(t, T))(B̄(t, T)Γ(t, T))(B̄(t, T)Γ(t, T))^{∗}) −

−\frac{1}{2}tr (D^{2}C(t, B(t, T))(B̄(t, T)σ(t, T))(B̄(t, T)σ(t, T))^{∗}) ≥ 0

We call the quantity R_{t} consumption rate. We have:

e_{t} = M_{t} ∫_{0}^{t} \frac{1}{M_{u}} R_{u} du

where M_{t} = exp(∫_{0}^{t} r_{u}(0) du). Then, to have an increasing e_{t} is equivalent to have R_{t} ≥ 0. Since this is verified, we have the thesis. ⋄

4.3 Sensibility of the model

First of all we suppose that the rates evolve following the Vasicek model:

dr_{t} = a(b − r_{t}) dt + σ dW_{t}

The zero coupon bond is

B(t, T) = e^{−(T −t)Y (t,T)} = 𝔼 [e^{− ∫_{t}^{T} r_{u} du}|F_{t}]

The explicit solution for the spot rate is:

r_{u} = e^{a(u−t)} + ∫_{0}^{u−t} e^{a(u−t−v)}ab dv + ∫_{0}^{u−t} e^{a(u−t−v)}σ dW_{v}

that is, r_{u} depends linearly on r_{t}, on a deterministic term and on a (Gaussian) random variable independent on F_{t}, so we can represent the zero coupon bond in this way:

B(t, T) = exp (−r_{t}(T − t)β(T − t) − (T − t)d(T − t))

where

β(T − t) = \frac{1}{T − t} ∫_{t}^{T} e^{a(u−t)} du = \frac{e^{a(T −t)} − 1}{a(T − t)}

Then we have:

Y (t, T) = r_{t}β(T − t) + d(T − t)
dY (t, T) = β(T − t) dr_{t} − β′(T − t)r_{t} dt − d′(T − t) dt

By applying Itô formula, we obtain that the volatility of B(t, T) is

Γ(t, T) = (T − t)β(T − t)σ

Now we take a portfolio having the deterministic payoffs F_{i} in the dates T_{i}, i = 1, …, n; then its evolution is given by:

P_{t} = ∑_{i=1}^{n} F_{i}B(t, T_{i})

dP_{t} = ∑_{i=1}^{n} F_{i}dB(t, T_{i}) = r_{t}P_{t} dt + ∑_{i=1}^{n} F_{i}B(t, T_{i})Γ(t, T_{i}) dW_{t}

Then the volatility of P_{t} is

\frac{∑_{i=1}^{n} F_{i}B(t, T_{i})Γ(t, T_{i})}{∑_{i=1}^{n} F_{i}B(t, T_{i})} = \frac{∑_{i=1}^{n} F_{i}B(t, T_{i})(T_{i} − t)β(T_{i} − t)σ}{P_{t}}

Now we define the duration of the portfolio P with respect to the direction f in this way:

D_{f}P = \frac{∑_{i=1}^{n} F_{i}B(t, T_{i})(T_{i} − t)f(T_{i} − t)}{P_{t}}

then the volatility of P_{t} can be interpreted as the duration of the portfolio with respect to the direction f(x) = σβ(x).

There is a link with the classical definition of duration: in fact we have that

\frac{1}{P_{t}} \frac{∂P_{t}}{∂r_{t}} = \frac{1}{P_{t}} ∑_{i=1}^{n} F_{i} \frac{∂B(t, T_{i})}{∂r_{t}} = \frac{∑_{i=1}^{n} F_{i}B(t, T_{i})(T_{i} − t)β(T_{i} − t)σ}{P_{t}}

If more in general we have that Y (t, T) is a deterministic function of k reference quantities (which can be rates or other variables)

Y (t, T) = β_{1}(t, T)R_{t}^{1} + … + β_{k}(t, T)R_{t}^{k} + α(t, T)

then:

\frac{1}{P_{t}} \frac{∂P_{t}}{∂R_{t}^{j}} = \frac{1}{P_{t}} ∑_{i=1}^{n} F_{i} \frac{∂B(t, T_{i})}{∂R_{t}^{j}} = \frac{∑_{i=1}^{n} F_{i}B(t, T_{i})(T_{i} − t)β_{j}(t, T_{i})}{P_{t}}

that can be interpreted as the sensibility of P_{t} in the j-th direction.

4.4 Open problems

• existence and uniqueness of the solution of the BSB equation (4.7); while in the case in which the final condition F ∈ C^{2}(H) the proof would probably be quite standard, in the case F ∈ C^{0} is not clear; probably one has to introduce a notion like viscosity solutions in infinite dimensions.

• is it possible to reduce the sup operator in both the equations (4.7) and (4.8) to a max? Moreover, when is it possible to project the Σ ∈ H^{∞} in a compact set in ℝ^{N}? Finally, is it possible to reduce the reducted BSB equation (4.8) to a partial differential equation in ℝ^{N+2} (the N + 2 variables are N bonds, the short rate and the time)?

• if Y (t, T) = G(r_{t}(·)) is a linear functional of all the forward rate curve r_{t}(·), is it possible to write something like

\frac{1}{P_{t}} \frac{∂P_{t}}{∂f} = \frac{∑_{i=1}^{n} F_{i}B(t, T_{i})(T_{i} − t)\frac{∂G}{∂f}}{P_{t}}

where \frac{∂G}{∂f} is the Frechet derivative of G with respect to F?

Appendix A

Positivity of rates

One of the main weaknesses of Gaussian interest rates models is the possibility of negative rates with positive probability. This appendix partially investigates this point. First of all we concentrate on the one-dimensional Vasicek model, and calculate the probability that there will be a negative forward rate at a certain time; this probability being 1, we calculate the probability that there will be a negative forward rate before a fixed instant T. Later, we do all this in the more general finite dimensional Gaussian model.

A.1 The Vasicek model

We take R ≥ 0, and we define

τ_{R} = inf{t | r_{t} < 0 or r_{t} > R}, τ = inf{t | r_{t} < 0}

We notice that τ_{R} → τ a.s. as R tends to +∞.

Lemma A.1 τ < +∞ a.s.

Proof. First of all we search for a function h solution of Ah = 0, where A is the second order operator associated to (r_{t})_{t}, that is:

Ah(x) = a(b − x)h′(x) + \frac{1}{2}σ^{2}h′′(x)  (A.1)

Lemma A.2 A solution of Ah = 0 is of the form

h(r) = D + C ∫_{0}^{r} e^{\frac{a}{σ^{2}} (y−b)^{2}} dy

where C, D are arbitrary real numbers.

Proof. We define g = h′; then the equation Ah = 0 becomes:

g′(x) = \frac{2a(x − b)}{σ^{2}}g(x)

then

g(x) = Ce^{∫ \frac{2a(y−b)}{σ^{2}} dy} = Ce^{\frac{a(x−b)^{2}}{σ^{2}}}

since h = ∫ v(x) dx, we have the thesis. ⋄

If h is of the above form, then M_{t} = h(r_{t∧τ_{R}}) is a bounded martingale (because 0 ≤ r_{t∧τ_{R}} ≤ R and h is continuous), so it converges a.s.; but, if τ_{R} = +∞, then h(r_{t∧τ_{R}}) = h(r_{t}) cannot converge, so τ_{r} < +∞ a.s.; besides, since (M_{t})_{t} is a bounded martingale, we have 𝔼[M_{t}] = 𝔼[M_{0}] ∀t ≥ 0, so:

h(x) = 𝔼[h(r_{0})] = lim_{t→+∞} 𝔼[h(r_{t∧τ_{R}})] = 𝔼 [ lim_{t→+∞} h(r_{t∧τ_{R}})] = 𝔼[h(r_{τ_{R}})]

where the third equality follows from Lebesgue theorem, so:

h(x) = h(0)ℙ{r_{τ_{R}} = 0} + h(R)ℙ{r_{τ_{R}} = R} = h(0)ℙ{r_{τ_{R}} = 0} + h(R)(1 − ℙ{r_{τ_{R}} = 0})

and we obtain:

ℙ\{r_{τ_{R}} = 0\} = \frac{h(x) − h(R)}{h(0) − h(R)} = \frac{∫_{x}^{R} e^{\frac{a}{σ^{2}} (y−b)^{2}} dy}{∫_{0}^{R} e^{\frac{a}{σ^{2}} (y−b)^{2}} dy}

We call v(x, R) = ℙ{r_{τ_{R}} = 0} under the condition r_{0} = x. Then we have the intuitive conditions v(0, R) = 1 and v(R, R) = 0. Moreover we have that

{τ < +∞} = {∃R such that r_{τ_{R}} = 0} = ⋃_{R>0}{r_{τ_{R}} = 0}

so

ℙ{τ < +∞} = lim_{R→+∞}{r_{τ_{R}} = 0} = lim_{R→+∞} v(x, R) = 1

⋄

We have seen that ℙ{τ < +∞} = 1, that is there exists almost surely a T such that r_{T} < 0. Now, let’s fix T > 0, and calculate ℙ{τ ≤ T }, that is the probability that the spot rate is positive up to the time T. First of all we give the result in terms of the solution of a parabolic PDE with boundary condition.

Lemma A.3 If there exist a solution v(t, x) ∈ C_{b}^{1} of the parabolic partial differential equation

\{ \begin{matrix}\frac{∂v}{∂t} + Av = 0 & (t, x) ∈ [0, T) × (0, +∞) \\ v(T, x) = 0 & x ∈ (0, +∞) \\ v(t, 0) = 1 & t ∈ [0, T) \end{matrix}

where A is given by Eq. (A.1), then ℙ{τ ≤ T } = v(0, r_{0}).

Proof. We apply Itô formula to the function v:

dv(t, r_{t}) = (\frac{∂v}{∂t} + Av) (t, r_{t}) dt + \frac{∂v}{∂x}(t, r_{t})σ dW_{t}

then, because of the PDE, we have:

v(T ∧ τ, r_{T ∧τ}) = v(0, r_{0}) + ∫_{0}^{T ∧τ} \frac{∂v}{∂x}(t, r_{t})σ dW_{t}

making the expectation, since v ∈ C_{b}^{1}, then the expectation of the stochastic integral is 0, and we get:

v(0, r_{0}) = 𝔼[v(T ∧ τ, r_{T ∧τ})] = 𝔼[v(T, r_{T})1_{τ>T} + v(τ, r_{τ})1_{τ≤T}] =
= 𝔼[0 · 1_{τ>T} + 1 · 1_{τ≤T}] = ℙ{τ ≤ T }
⋄

The problem with this approach is that there is not an explicit solution to this PDE, so we have no idea of the behaviour of the solution. The following proposition partially answers this question, giving an estimation from above to the searched probability.

Proposition A.4

ℙ\{τ ≤ T \} ≤ \frac{σ\sqrt{e^{2aT} − 1}}{2x\sqrt{πa}}

Proof. We notice that

ℙ{τ ≤ T } = ℙ { min_{t∈[0,T]} r_{t} ≤ 0} = ℙ { max_{t∈[0,T]}(−r_{t}) ≥ 0} = ℙ { max_{t∈[0,T]}(−e^{at}r_{t}) ≥ 0} =
= ℙ { max_{t∈[0,T]} (−b(e^{at} − 1) − ∫_{0}^{t} e^{as}σ dW_{s}) ≥ r_{0}} ≤
≤ ℙ { max_{t∈[0,T]} (− ∫_{0}^{t} e^{as}σ dW_{s}) ≥ r_{0}}

the process in the last expression is a martingale, so we can apply Kolmogorov’s maximal inequality, and we get:

ℙ\{τ ≤ T \} ≤ \frac{1}{r_{0}} 𝔼 [(∫_{0}^{T} e^{as} dW_{s})^{+}]

Now we need the

Lemma A.5 If X is a random variable with law N(0, 1), then 𝔼[X^{+}] = (2π)^{−1/2}.

Proof.

E[X^{+}] = \frac{1}{\sqrt{2π}} ∫_{0}^{+∞} xe^{− \frac{1}{2} x^{2}} dx = \frac{1}{\sqrt{2π}} [−e^{− \frac{1}{2} x^{2}}]_{0}^{+∞} = \frac{1}{\sqrt{2π}}

⋄

Since ∫_{0}^{T} e^{as} dW_{s} has law N(0, ∫_{0}^{T} e^{2as} ds), that is N(0, \frac{σ^{2}}{2a} (e^{2aT} − 1)), then we can calculate the last expectation, and finally we get:

ℙ\{τ ≤ T \} ≤ \frac{1}{r_{0}} \sqrt{\frac{σ^{2}}{2a}(e^{2aT} − 1)} = \frac{σ\sqrt{e^{2aT} − 1}}{2r_{0}\sqrt{πa}}

⋄

Now we try a different approach. We see the stochastic integral in the solution as a Brownian motion with respect to another time scale, that is the one of the local times.

We define:

c(t) = ∫_{0}^{t}(e^{as}σ)^{2} ds = \frac{σ^{2}}{2a}(e^{2at} − 1)

and the local time:

τ(u) = c^{−1}(u) = \frac{1}{2a} log (1 + \frac{2a}{σ^{2}}u)

then τ(u) are (deterministic) stopping times such that

B(u) = ∫_{0}^{τ(u)} e^{as}σ dW_{s}

is a Brownian motion with respect to the filtration (F_{τ(u)})_{t}; then we can write:

b(e^{aτ(u)} − 1) = b (exp (a\frac{1}{2a} log (1 + \frac{2a}{σ^{2}}u)) − 1) = b (\sqrt{1 + \frac{2a}{σ^{2}}u} − 1)

and we finally obtain:

e^{aτ(u)}r_{τ(u)} = r_{0} + b (\sqrt{1 + \frac{2a}{σ^{2}}u} − 1) + B(u), 0 ≤ u ≤ c(T) = \frac{σ^{2}}{2a}(e^{2at} − 1)

Now we define:

X_{u} = r_{0} + b (\sqrt{1 + \frac{2a}{σ^{2}}u} − 1) + B(u) = r_{0} + ϕ(u) + B(u),
Y_{u} = r_{0} + B(u)

then the processes X and Y induce the laws ℙ_{X} and ℙ_{Y} on C_{0}([0, c(T)]), which are equivalent and which mutual density is

\frac{dℙ_{X}}{dℙ_{Y}} = exp (∫_{0}^{c(T)} ϕ′(t) dY_{t} − \frac{1}{2} ∫_{0}^{c(T)} ϕ′(t)^{2} dt) =

= exp (∫_{0}^{c(T)} ϕ′(t) dB(t) − \frac{1}{2} ∫_{0}^{c(T)} ϕ′(t)^{2} dt)

Applying Hölder’s inequality, we have the following extimation:

ℙ{τ ≤ T } = ℙ { min_{u∈[0,c(T)]} X_{u} ≤ 0} = ∫_{Ω} 1_{min_{u∈[0,c(T)]} X_{u}≤0} dℙ_{X} =

= ∫_{Ω} 1_{min_{u∈[0,c(T)]} Y_{u}≤0} \frac{dℙ_{X}}{dℙ_{Y}} dℙ_{Y} ≤ 𝔼[1_{min_{u∈[0,c(T)]} Y_{u}≤0}^{p}]^{1/p}𝔼 [(\frac{dℙ_{X}}{dℙ_{Y}} )^{q}]^{1/q} =
= (ℙ \{ min_{u∈[0,c(T)]} Y_{u} ≤ 0\})^{1/p} 𝔼 [(\frac{dℙ_{X}}{dℙ_{Y}} )^{q}]^{1/q}

for all p, q > 1 such that 1/p + 1/q = 1. We have that:

𝔼 [(\frac{dℙ_{X}}{dℙ_{Y}} )^{q}] = 𝔼 [exp (q ∫_{0}^{c(T)} ϕ′(t) dB(t) − \frac{q}{2} ∫_{0}^{c(T)} ϕ′(t)^{2} dt)] =

= exp (−\frac{q}{2} ∫_{0}^{c(T)} ϕ′(t)^{2} dt) 𝔼 [exp (q ∫_{0}^{c(T)} ϕ′(t) dB(t))] =

= exp ((\frac{q^{2}}{2} − \frac{q}{2}) ∫_{0}^{c(T)} ϕ′(t)^{2} dt)

and:

∫_{0}^{c(T)} ϕ′(t)^{2} dt = ∫_{0}^{c(T)} \frac{b^{2}a^{2} du}{σ^{4}(1 + \frac{2a}{σ^{2}} u)} = ∫_{0}^{c(T)} \frac{b^{2} du}{\frac{2σ^{2}}{a} (\frac{σ^{2}}{2a} + u)} =

= [\frac{b^{2}a}{2σ^{2}} log (u + \frac{σ^{2}}{2a})]_{0}^{c(T)} = \frac{b^{2}a}{2σ^{2}} log ( \frac{\frac{σ^{2}}{2a} e^{2aT}}{\frac{σ^{2}}{2a}} ) = \frac{b^{2}a^{2}T}{σ^{2}}

Moreover, if we call τ_{a} = inf{t | B(t) = a} the stopping time of the first passage of the Brownian motion from a, then

ℙ { min_{u∈[0,c(T)]} Y_{u} ≤ 0} = ℙ { min_{u∈[0,c(T)]} B(u) ≤ −r_{0}} = ℙ{τ_{−r_{0}} ≤ c(T)} =
= \frac{2}{\sqrt{2πc(T)}} ∫_{−r_{0}}^{+∞} e^{− \frac{x^{2}}{2c(T)}} dx = \frac{2}{\sqrt{2π}} ∫−∞^{\frac{r_{0}}{c(T)}} e^{− \frac{y^{2}}{2}} dy

then we arrive to:

ℙ\{τ ≤ T \} = (2N (\frac{r_{0}}{c(T)}))^{1/p} exp ((q − 1)\frac{b^{2}a^{2}T}{2σ^{2}} ) =
= (2N (\frac{r_{0}}{c(T)}))^{1/p} exp (\frac{1}{p − 1} \frac{b^{2}a^{2}T}{2σ^{2}} )

for all p > 1, where N is the distribution function of the normal law N(0, 1):

N(d) = \frac{1}{\sqrt{2π}} ∫_{−∞}^{d} e^{− \frac{x^{2}}{2}} dx

We can try to find the optimal exponent p_{1}, that is the exponent p_{1} > 1 such that the quantity on the right hand side is minimal.

Lemma A.6 The optimal exponent is

p_{1} = \frac{log A − \sqrt{−B log A}}{B + log A}

where

A = 2N (\frac{r_{0}}{c(T)}) , B = \frac{b^{2}a^{2}T}{2σ^{2}}

Proof. We have to minimize the function

f(p) = A^{1/p} exp (\frac{B}{p − 1}) = exp (\frac{log A}{p} + \frac{B}{p − 1})

We notice that minimizing f is equivalent to minimize its logarithm, then we calculate (log f)′:

\frac{d}{dp} (\frac{log A}{p} + \frac{B}{p − 1}) = −\frac{log A}{p^{2}} − \frac{B}{(p − 1)^{2}}

Imposing (log f)′ ≥ 0, we find:

p^{2}(log A + B) − 2p log A + log A ≤ 0

(log f)′ is decreasing from 1 to p_{1}, then it is increasing until a point p_{2}, then it is decreasing, so p_{1} is the unique local minimum. ⋄

A.2 The d-dimensional Gaussian model

As before, we take L ≥ 0, and we define

τ_{L} = inf{t | R_{t} ∉ [0, L]^{d}}, τ = inf{t | R_{t} ∉ (ℝ^{+})^{d}}

We notice that τ_{L} → τ a.s. as L tends to +∞. First of all we search for a function h solution of Ah = 0, where A is the second order operator associated to (R_{t})_{t}, that is:

Ah(x) = ⟨A(c − x), Dh(x)⟩ + \frac{1}{2}Tr (σσ^{∗}D^{2}h(x))  (A.2)

where Dh is the gradient of h and D^{2}h indicates the Hessian matrix of h. If h is such a function, then M_{t} = h(R_{t∧τ_{L}}) is a bounded martingale (because R_{t∧τ_{L}} ∈ [0, L]^{d} and h is continuous), so it converges a.s.; but, if τ_{L} = +∞, then h(R_{t∧τ_{L}}) = h(R_{t}) cannot converge, so τ_{L} < +∞ a.s.; besides, since (M_{t})_{t} is a bounded martingale, we have 𝔼[M_{t}] = 𝔼[M_{0}] ∀t ≥ 0, so:

h(x) = 𝔼[h(R_{0})] = lim_{t→+∞} 𝔼[h(R_{t∧τ_{L}})] = 𝔼 [ lim_{t→+∞} h(R_{t∧τ_{L}})] = 𝔼[h(R_{τ_{L}})]

where the third equality follows from Lebesgue theorem. We call v(x, L) = ℙ{R_{τ_{L}}|_{i} = 0 for some i} under the condition R_{0} = x. Then we have the intuitive conditions v(0, L) = 1 and v(L, L) = 0. Moreover we have that

{τ < +∞} = {∃L such that R_{τ_{L}}|_{i} = 0 for some i} = ⋃_{L>0}{R_{τ_{L}}|_{i} = 0 for some i}

so

ℙ{τ < +∞} = lim_{L→+∞}{R_{τ_{L}}|_{i} = 0 for some i} = lim_{L→+∞} v(x, L) = 1

We have seen that ℙ{τ < +∞} = 1, that is there exists almost surely a T such that r_{T} < 0. Now, let’s fix T > 0, and calculate ℙ{τ ≤ T }, that is the probability that the rates are positive up to the time T. First of all we give the result in terms of the solution of a parabolic PDE with boundary condition.

Lemma A.7 If there exist a solution v(t, x) ∈ C_{b}^{1} of the parabolic partial differential equation

\{ \begin{matrix}\frac{∂v}{∂t} + Av = 0 & (t, x) ∈ [0, T) × (0, +∞)^{d} \\ v(T, x) = 0 & x ∈ (0, +∞)^{d} \\ v(t, x) = 1 & t ∈ [0, T), x ∈ ∂(0, +∞)^{d} \end{matrix}

where A is given by Eq. (A.2), then ℙ{τ ≤ T } = v(0, R_{0}).

Proof. We apply Itô formula to the function v:

dv(t, R_{t}) = (\frac{∂v}{∂t} + Av) (t, R_{t}) dt + \frac{∂v}{∂x}(t, R_{t})σ dW_{t}

then, because of the PDE, we have:

v(T ∧ τ, R_{T ∧τ}) = v(0, r_{0}) + ∫_{0}^{T ∧τ} \frac{∂v}{∂x}(t, R_{t})σ dW_{t}

making the expectation, since v ∈ C_{b}^{1}, then the expectation of the stochastic integral is 0, and we get:

v(0, R_{0}) = 𝔼[v(T ∧ τ, R_{T ∧τ})] = 𝔼[v(T, R_{T})1_{τ>T} + v(τ, R_{τ})1_{τ≤T}] =
= 𝔼[0 · 1_{τ>T} + 1 · 1_{τ≤T}] = ℙ{τ ≤ T }
⋄

As before, the problem with this approach is that there is not an explicit solution to this PDE, so we have no idea of the behaviour of the solution.

A.3 Open problems

• generalize the two kinds of estimation of ℙ{τ ≤ T } made for the Vasicek model to the d-dimensional Gaussian model

• make a similar work for the general Gaussian Musiela model, which is infinite dimensional

Bibliography

[1] M. Avellaneda, A. Levy, and A. Parás, Pricing and hedging derivative securities in markets with uncertain volatilities, preprint

[2] A. Brace, D. Gatarek, and M. Musiela, The market model of interest rate dynamics, preprint

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

[4] J. Cox, J. Ingersoll, S. Ross, A theory of the term structure of the interest rates, Econometrica, 53 (1985), 385–407

[5] G. Da Prato and J. Zabczyk, Stochastic equations in infinite dimensions, Cambridge University Press, 1992

[6] N. El Karoui, H. Geman, and V. Lacoste, On the role of state variables in interest rates models, preprint

[7] N. El Karoui, H. Geman, and J. C. Rochet, Changes of numeraire, changes of probability measures and option pricing, Journal of Applied Probability 32, 1995, 443–458.

[8] N. El Karoui, M. Jeanblanc-Pique, and S. E. Shreve, Robustness of the Black and Scholes formula, submitted to Math. Finance, 1995

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

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

[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. Musiela, Stochastic PDEs and term structure models, Journees Internationales de Finance, IGR-AFFI, La Baule, 1993 [14] S. Romagnoli, and T. Vargiolu, Robustness of the Black–Scholes approach in the case

of options on several assets, preprint

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

[16] O. Vasicek, An equilibrium characterisation of the term structure, Journal of Financial Economics, 5 (1977), 177–188