format: text

text:
# Pseudo-diffusions and Quadratic term structure models

## 1. Introduction

To account for fat tails, skewness and excessive kurtosis of empirical probability distributions of returns in real Financial Markets, it has become increasingly popular to model the dynamics of market factors as a Lévy process. Lévy models are more realistic than Gaussian ones but the latter are much more tractable. Indeed, in the Gaussian framework, explicit pricing formulas are known for a wide range of options and other contingent claims both without and with early exercise features, whereas in the Lévy models, most of the pricing formulas have been obtained for contingent claims of the European type, with the deterministic life-span. There are some explicit analytic results for options with early exercise features: see Boyarchenko and Levendorskiǐ (2000, 2001, 2002a, b), Mordecki (2002) and the bibliography therein for pricing of perpetual American options, and Boyarchenko and Levendorskiǐ (2002b, c) for pricing of barrier options and first touch digitals. However, the pricing formulas are complicated and difficult for numerical implementation except for a rather special case of pricing of perpetual American options under exponential jump-diffusions or spectrally one-sided processes.

Another obstacle for non-Gaussian modelling arises when one considers more general Markov processes. The explicit pricing formulas in affine term structure models and certain Lévy-driven Ornstein-Uhlenbeck models are known in the case of contingent claims with the deterministic life span only – see Duffie et al. (2000, 2002), Chacko and Das (2002), and Barndorff-Nielsen and Shephard (2001b), Barndorff-Nielsen et al (2002), respectively; for non-Gaussian variants of the HJM-model, see Eberlein and Raible (1999). In the general case, the dependance on the state variable does not allow one to obtain explicit analytical answers.

The following observation helps to obtain efficient approximate solutions. As Barndorff-Nielsen and Levendorskiǐ (2001) notice, typically, a good fit to the data can be achieved with Lévy processes whose Lévy densities decay as $`\mathrm{exp}(\lambda |x|)`$ or faster, where $`\lambda >0`$ (the steepness parameter of the exponential Lévy process) is large. They used this property to derive an asymptotic pricing formula for European options under certain class of Feller processes. The same observation was used in Boyarchenko and Levendorskiǐ (2002a,b,d) and Kudryavzev and Levendorskiǐ (2002) to derive efficient approximate formulas for perpetual American and Bermudan options, and first-touch-digitals, respectively.

It was shown in Boyarchenko and Levendorskiǐ (2002a, b) that the simple approximate formula is of the same form as the corresponding formula in a Gaussian model even when the underlying Lévy process has no Gaussian component. It can be shown that the leading term of the approximate pricing formula in Barndorff-Nielsen and Levendorskiǐ (2001) can also be written as the pricing formula in a Gaussian model. These observations can serve as an analytical explanation of relatively good performance of Gaussian models in apparently non-Gaussian situations. Thus, as far as pricing formulas are concerned, Lévy processes with large steepness parameters behave almost as the Brownian motion, and Feller processes with large steepness parameters considered in Barndorff-Nielsen and Levendorskiǐ (2001) behave almost as Gaussian diffusions. It seems reasonable to use the nomer pseudo-diffusions for Lévy processes and more general Lévy-like Feller processes with large steepness parameters .

The modelling with pseudo-diffusions allows one to obtain an efficient approximation to the price; in some situations, the asymptotic expansion of the price can be obtained, of the form

(1.1)
$$P(x,t)P_0(x,t)(1+\lambda ^1P_1(x,t)+\lambda ^2P_2(x,t)+\mathrm{}),$$

where the leading term, $`P_0`$, is the price in the Gaussian model with the same instantaneous drift and variance. The first correction term takes into account the moments of order three as well (skewness), the second correction term accounts for moments of order four, etc. Notice that though the leading term looks as the pricing formula in the Gaussian model, the “drift” and “variance-covariance matrix” used in the formula for the leading term are not the same as the ones of the Gaussian component of the process unless it is purely Gaussian. Indeed, a Lévy process may have no diffusion component at all.

The aim of the paper is to apply the approximate pricing approach to quadratic term structure models (QTSM) when the stochastic factor follows a mean-reverting pseudo-diffusion process of the simplest form (it is unlikely that in the QTSM model, an explicit pricing formula can be obtained unless the process process is Gaussian), and derive a pricing formula of the form (1.1). For the discussion about advantages of the Gaussian QTSM model, see Ahn et al (2002, 2003) and Chen and Poor (2002). Cheng and Scaillet (2002) consider an affine-quadratic model, and allow for jumps but only in the dynamics of affine variables of the model. Notice that the use of jumps in QTSM models adds additional flexibility in joint modelling under the historic and a risk-neutral measures, and one may hope that the performance of QTSM models can be improved by introducing jumps.111The author is grateful to an anonymous referee for this suggestion Another improvement (and quite sizable one) is expected in pricing of out-of-the-money options near expiry, where the main contribution to the price comes from the jump part of the process. Near expiry, however, a different approximate formulas are needed, which use more detailed information about the jump part of the processes than the skewness and kurtosis. These formulas are similar to approximate formulas for out-of-the-money options on stocks developed in Levendorskiǐ (2003), and can be derived by the same reasoning.

### 1.1. Plan of the paper

In Section 2, we list families of exponential Lévy processes used in empirical and theoretical studies of financial markets. In Section 3, we formulate the pricing problem for an interest rate derivative of the European type, and by using the Feynman-Kac theorem, reduce the pricing problem to the boundary problem for an integro-differential equation. We also explain the scheme of the asymptotic pricing. In Section 4, we recall the solution of the bond pricing problem in the one-factor Gaussian case, and indicate the properties of the solution which are crucial for our asymptotic method. In Section 5, we demonstrate our method in the simplest case of the one-factor Lévy model for the bond price, and present numerical examples. In Section 6, we consider possible specifications of the market price of risk, the generalization for the multi-factor case, derive approximate formulas for interest rate derivatives near expiry, and suggest a procedure of parameter fitting based on the asymptotic expansions. In Section 7, we summarize our results, and compare the Lévy QTSM with multi-factor Gaussian QTSM. In the appendix, technical results are proven.

## 2. Lévy processes in financial modelling

As early as in 1963, Mandelbrot suggested to use stable Lévy processes. The modelling with stable Lévy processes is not quite realistic since the tails of Lévy stable distributions are too fat (polynomially decaying), whereas the tails of distributions of returns observed in real financial markets exhibit exponential decay. Moreover, the second moment of a Lévy stable distribution is infinite (unless it is a Gaussian one). This contradicts the observed convergence to the Gaussian distribution over a longer time scale, and even worse, the underlying stock itself should have the infinite price under the stable Lévy process, which makes the model inconsistent for pricing purposes. Starting with the beginning of the 90th, several families of Lévy processes with probability distributions having exponentially decaying tails have been used to describe the behavior of stock prices in real financial markets:

* Variance Gamma Processes (VGP) constructed and used by Madan and co-authors in a series of papers during 90th (see Madan et al. (1998) and the bibliography therein);
* Hyperbolic Processes (HP) were constructed and used by Eberlein and co-authors (see Eberlein et al. (1998), Eberlein and Prause (1999)); hyperbolic distributions were constructed by Barndorff-Nielsen (1977));
* Normal Inverse Gaussian Processes (NIG) were introduced by Barndorff-Nielsen (1998) and used to model German stocks by Barndorff-Nielsen and Jiang (1998);
* Truncated Lévy Processes (TLP) constructed by Koponen (1995) were used for modeling in real financial markets by Bouchaud and Potters (1997), Cont et al (1997) and Matacz (2001); the extended Koponen family was constructed in Boyarchenko and Levendorskiǐ (2000) (the generalization was needed since probability distribution of Koponen’s family have tails of the same rate of exponential decay whereas in real financial markets, the left tail is usually much fatter; in Carr et al (2002) and Boyarchenko and Levendorskiǐ (2002a,b), the extended Koponen family is called CGMY-model and KoBoL family, respectively).
* Normal Tempered Stable Lévy processes were constructed in Barndorff-Nielsen and Levendorskiǐ (2001) and Barndorff-Nielsen and Shephard (2001a); they contain NIG as a subclass.

In Boyarchenko and Levendorskiǐ (2000), a general class of Lévy processes, which contained all the classes listed above modulo certain reservation about VGP was introduced, under the name Generalized Truncated Lévy Processes. Later, in Barndorff-Nielsen and Levendorskiǐ (2001), the name: “Regular Lévy processes of exponential type” (RLPE) was suggested. For a more detailed exposition, see Boyarchenko and Levendorskiǐ (2002a, 2002b). In order to present examples, recall that a Lévy process can be completely specified by its characteristic exponent, $`\psi `$, definable from the equality $`E[e^{i\xi ,X(t)}]=e^{t\psi (\xi )}`$. The characteristic exponent is given by the Lévy-Khintchine formula

(2.1)
$$\psi (\xi )=ib,\xi +\frac{1}{2}A\xi ,\xi +_{𝐑^n}(1+i\xi ,y\mathrm{𝟏}_{||1}(y)e^{i\xi ,y})F(dy),$$

where $`A:=\mathrm{\Sigma }\mathrm{\Sigma }^T`$ is the variance-covariance matrix of the Gaussian component, $`b𝐑^n`$, and $`F(dx)`$ is the Lévy density (density of jumps), which satisfies

$$_{𝐑^n}\mathrm{min}\{|x|^2,1\}F(dx)<\mathrm{}.$$

Any generating triplet $`A,b,F(dx)`$ with these properties defines a Lévy process (see e.g. Sato (1999)). If $`\mathrm{\Sigma }=0`$, then we have a pure jump process.

Wide families of jump-diffusion processes are subclasses of the class of RLPE. In the first example, we introduce the family which is widely used in affine term structure models (see Duffie et al (2000) and Chacko and Das (2002)).

Example 2.1. Let $`X`$ be a Lévy process with the Lévy density

$$F(dx)=c_+\lambda _+e^{\lambda _+x}\mathrm{𝟏}_{(\mathrm{},0)}(x)dx+c_{}(\lambda _{})e^{\lambda _{}x}\mathrm{𝟏}_{(0,+\mathrm{})}(x)dx,$$

where $`\lambda _+>0,\lambda _{}<1`$ and $`c_\pm >0`$. Then

$$\psi (\xi )=\frac{\sigma ^2}{2}\xi ^2ib\xi +\frac{ic_+\xi }{\lambda _++i\xi }+\frac{ic_{}\xi }{\lambda _{}+i\xi },$$

where $`\sigma ^20`$ and $`b𝐑`$ are the variance and drift of the Gaussian component. The $`\psi (\xi )`$ is analytic in the strip $`\mathrm{}\xi (\lambda _{},\lambda _+)`$.

Example 2.2. The characteristic exponent of a process of KoBoL family in 1D is of the form

(2.2)
$$\psi (\xi )=i\mu \xi +c\mathrm{\Gamma }(\nu )[\lambda _+^\nu (\lambda _++i\xi )^\nu +(\lambda _{})^\nu (\lambda _{}i\xi )^\nu ],$$

where $`\nu (0,2),\nu 1,c>0,\lambda _{}<0<\lambda _+`$, and $`\mu 𝐑`$; it is analytic in a strip $`\mathrm{}\xi (\lambda _{},\lambda _+)`$, and (3.8)-(3.9) are satisfied in this strip.

Example 2.3. The characteristic exponent of a Normal Inverse Gaussian process in 1D is of the form

(2.3)
$$\psi (\xi )=i\mu \xi +\delta [(\alpha ^2(\beta +i\xi )^2)^{1/2}(\alpha ^2\beta ^2)^{1/2}],$$

where $`\nu (0,2),\delta >0,`$ and $`\alpha >|\beta |`$; it is analytic in the strip $`\mathrm{}\xi (\alpha +\beta ,\alpha +\beta )`$, and (3.8)-(3.9) are satisfied in this strip, with $`\nu =1`$.

Since the sum of the characteristic exponents of two RLPE’s is the characteristic exponent of an RLPE, the list of model examples can easily be expanded. For multi-dimensional examples, see Boyarchenko and Levendorskiǐ (2002b).

Examples 2.1–2.3 are examples of pseudo-diffusions if $`\lambda _+,|\lambda _{}|`$, and $`\alpha \pm \beta `$ are large. Typically, processes observed in empirical studies of financial markets (hyperbolic processes and variance gamma processes including) enjoy this property.

The majority of papers on Lévy models deal with asset pricing. Eberlein and Raible (1999) consider the HJM-model driven by a Lévy process (see also Eberlein and Özkan (2001)). For the usage of jump-diffusion processes and more general Lévy processes in affine term structure models of interest rates, see Duffie et al. (2000, 2002), Chacko and Das (2002) and the bibliography therein. Barndorff-Nielsen and Shephard (2001b) suggested to use Lévy-driven Ornstein-Uhlenbeck processes for interest rate modelling purposes. For the subsequent developments, see Barndorff-Nielsen et al (2002).

## 3. The model

### 3.1. Lévy-driven QTSM

In the Gaussian QTSM, the instantaneous interest rate is represented as a quadratic function of the state variables, and the latter are specified as diffusions. We assume that under an EMM chosen by the market, the SDE of the state variables can be written as

(3.1)
$$dX(t)=(\stackrel{~}{\theta }(t)\kappa X(t))dt+dZ(t),$$

where $`\{Z(t)\}`$ is an $`n`$-dimensional Lévy process, $`\stackrel{~}{\theta }:𝐑^n𝐑`$ is a continuous vector-function, and $`\kappa `$ is a constant $`n\times n`$ matrix, whose eigenvalues $`\lambda _j`$ satisfy the condition

(3.2)
$$\mathrm{}\lambda _j>0.$$

The interest rate is modelled as

(3.3)
$$r(X(t))=R_0+2R_1,X(t)+\mathrm{\Gamma }X(t),X(t),$$

where $`R_0𝐑`$, $`R_1𝐑^n`$ are constant scalar and vector, $`,`$ is the standard inner product in $`𝐑^n`$, and $`\mathrm{\Gamma }`$ is a positively definite symmetric matrix. The last condition ensures that

$$r(X(t))=\mathrm{\Gamma }(X(t)+\mathrm{\Gamma }^1R_1),X(t)+\mathrm{\Gamma }^1R_1+R_0\mathrm{\Gamma }^1R_1^2$$

is semi-bounded from below. By choosing $`R_0,R_1`$ and $`\mathrm{\Gamma }`$ appropriately, one can ensure any lower bound on $`r(X(t))`$. Notice that if one wishes to price a derivative of a stock whose dynamics is characterized by $`X`$, then one may allow $`r`$ to depend only on some of the factors $`X_j(t)`$, say, $`r=r(X_1(t),\mathrm{},X_m(t))`$, where $`m<n`$; in this case, in (3.3), $`R_1𝐑^m`$, and $`\mathrm{\Gamma }`$ is an $`m\times m`$ matrix.

If $`Z`$ has no jump component then the bond pricing problem reduces to a system of ODE (Riccati equations), which can easily be solved numerically, and in the one-factor case, even analytically. (In the multi-factor case, a system of Riccati equations can be reduced to a linear system; for the explicit realization in the framework of the Gaussian QTSM, see Kim (2003)). It seems unlikely that a reasonably simple exact solution exists for a general Lévy process but we manage to obtain an asymptotic solution if $`X`$ is a pseudo-diffusion, that is, the Lévy density of $`Z`$ decays exponentially, and the rate of decay is large. The leading term of the asymptotics is the price in the Gaussian model with the same instantaneous moments of order one and two, and the correction terms are polynomials in the factors with coefficients depending on the time to expiry. After the leading term is found, they can be calculated recursively, by using only integration procedures in 1D. Thus, the suggested method is relatively simple (though in multi-factor models, the number of additional integration procedures may be rather large; it is important that all the integrations remain one-dimensional, even in a multi-factor model). In the one-factor case, the first correction term is proportional to skewness, and the second one depends on the skewness and kurtosis; to be more precise, the first correction is proportional to skewness, and the second one is the sum of two terms, one of which is proportional to the square of the skewness, and the other to the kurtosis. In many cases, the contribution of the kurtosis is small relative to the other terms; if we omit the last term, then the pricing formula becomes a sum of the leading term which looks as the price in the Gaussian model, and the correction term, which is a quadratic polynomial w.r.t. to skewness.

Similar formula for the forward rate and numerical examples show that the first correction term has a pronounced upward hump, if the skewness is negative; in the result, the corrected forward rate curve can be hump shaped even when the Gaussian forward rate curve is not, and all parameters of the model are time-independent. By changing the parameters, various shapes of the forward rate curve can be obtained.

Empirical studies show that both skewness and kurtosis can be fairly large, and hence, the corrections to the Gaussian price quite sizable. Consider, for instance, the statistics for the daily change interest rates ($`dr`$) from Table 1 in Das (2002). (The table presents descriptive statistics for the Fed Funds rate over the period January 1988 to December 1997, and the unit is 1 percent). Mean: $`m=0.0005`$; standard deviation: $`\sigma =0.2899`$; skewness: $`\lambda _3=0.3950`$; excess kurtosis: $`k_4=19.8667`$. Recall that for probability distribution $`P(dx)`$,

$$m:=x:=_{\mathrm{}}^+\mathrm{}xP(dx),\sigma ^2:=(xm)^2,$$

$$\lambda :=(xm)^3/\sigma ^3,k_4:=(xm)^4/\sigma ^43,$$

and that if $`P(dx)=P_{\mathrm{\Delta }t}(dx)`$ is the probability distribution of a Lévy process with the characteristic exponent $`\psi `$, then

$$m(\mathrm{\Delta }t)/\mathrm{\Delta }t=i\psi ^{}(0);\sigma ^2(\mathrm{\Delta }t)/\mathrm{\Delta }t=\psi ^{\prime \prime }(0);$$

$$(xm)^3(\mathrm{\Delta }t)/\mathrm{\Delta }t=i\psi ^{(3)}(0);[(xm)^4(\mathrm{\Delta }t)3\sigma ^4(\mathrm{\Delta }t)]/\mathrm{\Delta }t=\psi ^{(4)}(0).$$

We see that the coefficients in the third and fourth terms in the Taylor series for $`\psi `$ around zero are smaller than the second one but non-negligible whereas in the Gaussian case all coefficients starting from the third one are zero.

The skewness and kurtosis of the process under an EMM can assume essentially arbitrary values provided they are small w.r.t. variance; in particular, one should expect that the skewness of the process under EMM is negative even when the one under historic measure is positive as in the empirical example above. This means that even the one-factor approximate non-Gaussian model has two free additional parameters (albeit small) which can be used to get a better fit to the data than in the Gaussian model. In multi-factor models, the number of additional free parameters is larger still.

### 3.2. Reduction of a pricing problem to a boundary problem

Consider a contingent claim with the maturity date $`T`$ and payoff $`g(X(T))`$. Its price at time $`t<T`$ is given by

(3.4)
$$f(X(t),t)=E_t\left[\mathrm{exp}\left(_t^Tr(X(s))𝑑s\right)g(X(T))\right].$$

(We consider the pricing under a risk-neutral measure chosen by the market). In applications, the payoff $`g`$ is measurable (usually, continuous), and it may grow at infinity. In the latter case, additional conditions on $`Z`$ may be needed. For instance, if $`g`$ grows not faster than an exponential:

(3.5)
$$|g(x)|Ce^{\omega |x|},$$

where $`C`$ and $`\omega >0`$ are independent of $`x`$, then it suffices to assume that there exists $`\lambda >\omega `$ such that for all $`\mu `$ in the ball $`|\mu |<\lambda `$, and some $`t>0`$,

(3.6)
$$E[e^{\mu ,Z(t)}]<\mathrm{},$$

which implies that the tails of probability densities of the process $`Z`$ decay exponentially: faster than $`\mathrm{exp}(\rho |x|)`$, for any $`\rho <\lambda `$.

It follows from (3.6), that for any $`\xi =\eta +i\tau 𝐂^n`$ from the tube domain $`𝐑^n+iU_\lambda :=\{\xi ||\mathrm{}\xi |=|\tau |<\lambda \}`$ (in the one-factor case, a tube domain is a strip), and any $`t>0`$,

(3.7)
$$E[e^{i\xi ,Z(t)}]<\mathrm{}.$$

(Instead of balls $`U_\lambda `$, one can use more general open sets containing the origin.) It is immediate from (3.7), that $`\psi (\xi )`$ and its derivatives w.r.t. the complex argument $`\xi `$ are well-defined in the same tube domain $`𝐑^n+iU_\lambda `$ (one says that $`\psi (\xi )`$ is analytic in $`𝐑^n+iU_\lambda `$), and we may use the latter condition on $`\psi `$ instead of the former condition (3.7). To justify the use of the Feynman-Kac formula, we assume that $`Z`$ is a regular Lévy process of exponential type (RLPE). This means that $`\psi `$ admits a representation

(3.8)
$$\psi (\xi )=i\mu ,\xi +\varphi (\xi ),$$

where $`\mu 𝐑^n`$, and $`\varphi `$ satisfies the following condition: there exist $`c>0,\nu (0,2]`$ and $`\nu _1<\nu `$ such that as $`\xi \mathrm{}`$ in the tube domain $`𝐑^n+iU_\lambda `$,

(3.9)
$$\varphi (\xi )=c|\xi |^\nu +O(|\xi |^{\nu _1})$$

(see Boyarchenko and Levendorskiǐ (2002b)). The $`\nu `$ and $`U_\lambda `$ are called the order and type of the process.

To simplify the justification of the use of the Feynman-Kac formula, we add unnecessary condition: for any multi-index $`\alpha =(\alpha _1,\mathrm{},\alpha _n)`$, there exists a constant $`C_\alpha `$ such that for all $`\xi `$ in the tube domain $`𝐑^n+iU_\lambda `$,

(3.10)
$$|^\alpha \varphi (\xi )|C_\alpha (1+|\xi |)^{\nu |\alpha |},$$

where $`|\alpha |=\alpha _1+\mathrm{}+\alpha _n`$. Notice that this condition holds for all model classes of RPPE’s.

In the appendix, by making use of the Feynman-Kac formula, we will prove the following theorem.

###### Theorem 3.1.

Let the stochastic factor satisfy (3.1), (3.2), (3.3), (3.6), (3.8), (3.9), and (3.10), let $`r`$ be given by (3.3), and let $`g`$ be a continuous function, which admits a bound (3.5).

Then a) the stochastic expression (3.4) defines a continuous function $`f`$, which admits an estimate

(3.11)
$$|f(x,t)|C_1e^{\omega |x|},$$

where $`C_1`$ is independent of $`x`$ and $`tT`$;

b) $`f`$ is a unique solution to the following problem:

(3.12) $`(_t+\stackrel{~}{\theta }(t)\kappa x,_x+Lr(x))f(x,t)`$ $`=`$ $`0,t<T,`$
(3.13) $`f(x,T)`$ $`=`$ $`g(x),`$
where $`L`$ is the infinitesimal generator of $`Z`$.

Recall that the infinitesimal generator of the Lévy process $`Z`$, $`L`$, can be represented in the form of a pseudo-differential operator (PDO) with the symbol $`\psi `$: $`L=\psi (D_x)`$. A PDO $`A=a(D)`$ with the symbol $`a`$ acts on sufficiently regular functions as follows:

$$(Au)(x)=(2\pi )^n_{𝐑^n}e^{ix,\xi }a(\xi )\widehat{u}(\xi )𝑑\xi ,$$

where $`\widehat{u}`$ is the Fourier transform of $`u`$:

$$\widehat{u}(\xi )=_{𝐑^n}e^{ix,\xi }u(x)𝑑x.$$

In particular, the partial derivative $`_x`$ is the PDO with the symbol $`i\xi `$.

### 3.3. Asymptotic pricing

The asymptotic pricing formulas will be derived under the following conditions. Assume that the characteristic exponent of the driving Lévy process depends on a small parameter $`ϵ>0`$: $`\psi (\xi )=\psi (ϵ,\xi )`$ and satisfies the following three conditions. First, we require that the $`\lambda `$ in the definition of the tube domain $`𝐑^n+iU_\lambda `$ satisfies $`\lambda >>ϵ^{1/2}`$. The next two conditions are formulated for $`\xi `$ in the tube domain $`𝐑^n+iU_\lambda `$:

1) in the region $`|\xi |>ϵ^{1/2}`$, $`\psi (ϵ,\xi )`$ admits an estimate

(3.14)
$$\mathrm{}\psi (ϵ,\xi )c|\xi |^\nu ,$$

where $`\nu (0,2]`$ and $`c>0`$ are independent of $`(ϵ,\xi )`$ in the region;

2) in the region $`|\xi |ϵ^{1/2}`$, $`\psi (ϵ,\xi )`$ admits an asymptotic expansion: in the one-factor case,

(3.15)
$$\psi (ϵ,\xi )=i\mu \xi +\frac{\sigma ^2}{2}\xi ^2\underset{j=3}{\overset{\mathrm{}}{}}ϵ^{j2}k_j(i\xi )^j,$$

where the coefficients $`k_j`$ are uniformly bounded:

(3.16)
$$|k_j|C\sigma ^2/2,$$

where $`C`$ is independent of $`j`$; in the multi-variate case, (3.15) is replaced with

(3.17)
$$\psi (ϵ,\xi )=i\mu ,\xi +\frac{1}{2}\mathrm{\Sigma }^T\xi ^2\underset{j=3}{\overset{\mathrm{}}{}}ϵ^{j2}k_j(i\xi ),$$

where $`k_j(\xi )`$ is a homogeneous polynomial of order $`j`$, which admits a bound

(3.18)
$$|k_j((\mathrm{\Sigma }^T)^1i\xi )|C|\xi |^j,$$

where $`C`$ is independent of $`j`$.

The asymptotic solution will be found in the following sections. Here we explain the main idea in the one-factor case. We look for the solution in the form

(3.19)
$$f=f_0+ϵf_1+ϵ^2f_2+\mathrm{}.$$

From (3.15), we can formally write

(3.20)
$$L=\mu _x+\frac{\sigma ^2}{2}_x^2+\underset{j=3}{\overset{\mathrm{}}{}}ϵ^{j2}k_j_x^j,$$

and by substituting (3.19) and (3.20) into (3.12), we obtain a formal equality

(3.21)
$$\left(L_0+\underset{j=1}{\overset{\mathrm{}}{}}ϵ^jL_j\right)\left(f_0+\underset{l=1}{\overset{\mathrm{}}{}}ϵ^lf_l\right)=0,t<T,$$

where

$$L_0=_t+(\stackrel{~}{\theta }(t)+\mu \kappa x)_x+\frac{\sigma ^2}{2}_x^2r(x)$$

is of the same form as the operator in the Gaussian model, and

$$L_l=k_{l+2}_x^{l+2},l=1,2,\mathrm{}.$$

By multiplying out in (3.21) and gathering terms of the same order in $`ϵ`$, we obtain the following series of problems. The leading term of the asymptotics is found from

$`L_0f_0(x,t)`$ $`=`$ $`0,t<T;`$
$`f_0(x,T)`$ $`=`$ $`g(x),`$
which is the pricing problem in the Gaussian model; and the following terms are found step by step, by solving problems

$`L_0f_l(x,t)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{l}{}}}k_{j+2}_x^{j+2}f_{lj}(x,t),t<T,`$
$`f_l(x,T)`$ $`=`$ $`0,`$
for $`l=1,2,\mathrm{}.`$ We believe that for practical purposes, it suffices to use an approximate formula (3.19) with terms up to order 2; this allows one to take into account the skewness and kurtosis. This approximate solution can be written as

(3.22)
$$ff_0+ϵk_3f_1+(ϵk_3)^2f_{21}+ϵ^2k_4f_{22},$$

where $`f_1`$, $`f_{21}`$ and $`f_{22}`$ solve equations

$`L_0f_1(x,t)`$ $`=`$ $`_x^3f_0(x,t),`$
$`L_0f_{21}(x,t)`$ $`=`$ $`_x^3f_1(x,t),`$
$`L_0f_{22}(x,t)`$ $`=`$ $`_x^4f_0(x,t)`$
in the half-space $`t<T`$, subject to zero boundary condition. The explicit formulas for the bond price can be found in Section 4 and Section 5. Formula (3.22) may seem somewhat inconvenient for practical applications since it depends on the small parameter $`ϵ`$, which is not explicitly specified. Notice, however, that

$$ϵk_3=i\psi ^{(3)}(0)/3!,ϵ^2k_4=\psi ^{(4)}(0)/4!,$$

and the derivatives of the characteristic exponent at 0 can be inferred from empirical data - see Introduction. Thus, we may write (3.22) without $`ϵ`$:

(3.23)
$$ff_0i\frac{\psi ^{(3)}(0)}{3!}f_1\left(\frac{\psi ^{(3)}(0)}{3!}\right)^2f_{21}\frac{\psi ^{(4)}(0)}{4!}f_{22}.$$

By using (3.23), the influence of the moments of order 3 and 4 on the price can be explicitly analyzed; and this influence is highly non-linear in $`(x,t)`$, since the functions in (3.23) are.

If $`P:=f`$ is the bond price, then we can derive from (3.23) similar approximate formulas for the yield and forward rate.

The final remark is: in order to find a current term of the asymptotics, we have to differentiate the previous terms, therefore an asymptotic solution with several terms may produce serious errors in the neighborhood of a point where the pay-off $`g`$ is not sufficiently smooth. Indeed, one can hardly expect that a formula which is polynomial in $`x`$, can give a high order approximation in this case. Hence, in a neighborhood of such a point, a different asymptotic formula should be written: see Section 6.

## 4. Bond pricing: Gaussian model, one-factor case

In this section, we recall the solution of problem (3.12)-(3.13) in the one-factor Gaussian case, when $`\psi (\xi )=i\mu \xi +\sigma ^2\xi ^2/2`$, and $`L=\mu _x+\frac{\sigma ^2}{2}_x^2`$, and indicate the properties of the solution which are crucial for the asymptotic method to work.

We assume that the interest rate is a quadratic function of the stochastic factor $`X(t)`$:

(4.1)
$$r(X(t))=R_0+2R_1X(t)+X(t)^2.$$

The dynamics of the stochastic factor, $`X`$, is governed by (3.1), where $`\stackrel{~}{\theta }`$ is a scalar function, and $`\kappa `$ is a positive scalar. The bond price is given by (3.4) with $`g(x)1`$, hence it is a bounded solution to problem (3.12)-(3.13) with $`g(x)=1`$ in the RHS of (3.13).

Set $`\tau =Tt`$, $`\theta (\tau )=\stackrel{~}{\theta }(Tt)+\mu `$, and with some abuse of notation, write $`f(x,\tau )`$ instead of $`f(x,T\tau )`$. We look for the bounded solution to the problem

(4.2) $`(_\tau +(\theta (\tau )\kappa )_x+{\displaystyle \frac{\sigma ^2}{2}}_x^2r(x))f(x,\tau )`$ $`=`$ $`0,\tau >0,`$
(4.3) $`f(x,0)`$ $`=`$ $`g(x)`$
in the form

(4.4)
$$f(x,\tau )=\mathrm{exp}\mathrm{\Phi }_0(x,\tau ),$$

where

(4.5)
$$\mathrm{\Phi }_0(x,\tau )=A(\tau )x^2+B(\tau )x+C(\tau ).$$

By substituting (4.4) into (4.2), we obtain

(4.6)
$$(\mathrm{exp}(\mathrm{\Phi }_0)\mathrm{exp}\mathrm{\Phi }_0)(x,\tau )r(x)=0,$$

where

(4.7)
$$=_\tau +(\theta (\tau )\kappa )_x+\frac{\sigma ^2}{2}_x^2,$$

subject to $`\mathrm{\Phi }_0(x,0)=0`$. Straightforward calculations yield the following system of ODE with zero initial data:

(4.8) $`A^{}(\tau )2\kappa A(\tau )+2\sigma ^2A(\tau )^21`$ $`=`$ $`0,`$
(4.9) $`B^{}(\tau )\kappa B(\tau )+2\sigma ^2A(\tau )B(\tau )+2\theta (\tau )A(\tau )2R_1`$ $`=`$ $`0,`$
(4.10) $`C^{}(\tau )+\sigma ^2A(\tau )+{\displaystyle \frac{\sigma ^2}{2}}B(\tau )^2+\theta (\tau )B(\tau )R_0`$ $`=`$ $`0.`$
Equation (4.8) is solved by separation of variables:

(4.11)
$$A(\tau )=A_1A_2\frac{1e^{\omega \tau }}{A_2A_1e^{\omega \tau }},$$

where $`A_1<0<A_2`$ are roots of the quadratic equation $`2\sigma ^2A^22\kappa A1=0`$, and $`\omega =2\sigma ^2(A_1A_2)<0`$. $`A(\tau )`$ having being found, we can calculate $`B(\tau )`$ from the linear equation (4.9):

(4.12)
$$B(\tau )=\frac{2e^{\omega _1\tau }(A_2I_1(\tau )A_1I_2(\tau ))}{(A_2A_1)(A_2A_1e^{\omega \tau })},$$

where $`\omega _1=2\sigma ^2A_1\kappa `$, and

$`I_1(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }(A_1\theta (s)R_1)e^{\omega _1s}𝑑s,`$
$`I_2(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }(A_2\theta (s)R_1)e^{(\omega \omega _1)s}𝑑s.`$
If $`\theta `$ is independent of $`\tau `$, then $`I_j`$ can be calculated explicitly:

$`I_1(\tau )`$ $`=`$ $`{\displaystyle \frac{A_1\theta R_1}{\omega _1}}(1e^{\omega _1\tau }),`$
$`I_2(\tau )`$ $`=`$ $`{\displaystyle \frac{A_2\theta R_1}{\omega _1\omega }}(1e^{(\omega \omega _1)\tau }),`$
and therefore

$`B(\tau )`$ $`=`$ $`{\displaystyle \frac{2}{(A_2A_1)(A_2A_1e^{\omega \tau })}}`$
$`\times \left[{\displaystyle \frac{A_2(A_1\theta R_1)}{\omega _1}}(e^{\omega _1\tau }1){\displaystyle \frac{A_1(A_2\theta R_1)}{\omega _1\omega }}(e^{\omega _1\tau }e^{\omega \tau })\right].`$
Finally, we find $`C(\tau )`$ from (4.10) by integration:

(4.14)
$$C(\tau )=_0^\tau \left[\sigma ^2A(s)+\frac{\sigma ^2}{2}B(s)^2+\theta (s)B(s)R_0\right]𝑑s.$$

If $`\theta `$ is constant, then $`C(\tau )`$ can be calculated explicitly. In order that $`B(\tau )`$ and $`C(\tau )`$ can be calculated explicitly, $`\theta `$ need not to be a constant; for instance, one can use exponential polynomials.

To end this section, we make the crucial remark on the properties of the solution. First, from (4.11),

(4.15)
$$A(\tau )<0,\tau >0,$$

and as $`\tau 0`$,

(4.16)
$$A(\tau )A_1A_2\frac{\omega }{A_2A_1}\tau =A_1A_22\sigma ^2\tau =\tau .$$

Hence, for any $`\tau (0,T]`$, $`f(x,\tau )`$ decays as $`\mathrm{exp}(\tau x^2)`$, as $`x\pm \mathrm{}`$, and $`\widehat{f}(\xi ,\tau )`$, the Fourier transform of $`f(x,\tau )`$ w.r.t. the first argument, decays as $`\tau ^{1/2}\mathrm{exp}(\xi ^2/(4\tau ))`$, as $`\xi \pm \mathrm{}`$. To be more specific,

$`\widehat{f}(\xi ,\tau )`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}e^{ix\xi +A(\tau )x^2+B(\tau )x+C(\tau )}𝑑x`$
$`=`$ $`{\displaystyle \frac{1}{\pi A(\tau )}}\mathrm{exp}[C(\tau )+(\xi +iB(\tau ))^2/(4A(\tau ))].`$
We conclude that for any $`N`$, in the region $`\tau (0,T]`$, $`|\xi |>ϵ^{1/2}`$,

(4.18)
$$|\xi ^N\widehat{f}(\xi ,\tau )|C_Ne^{\xi ^2/(8\tau )},$$

where $`C_N`$ is independent of $`\tau `$ and $`\xi `$. Notice that the RHS of (4.18) is negligible.

## 5. Bond pricing: Lévy model, one-factor case

### 5.1. The leading term of the asymptotics

We assume that (3.14)-(3.16) hold. Take $`\mu `$ and $`\sigma ^2`$ from (3.15), and denote by $`f_0`$ the solution to the Gaussian bond pricing problem (4.2)-(4.3); it is given by (4.4), (4.11), (4.12) and (4.14), and satisfies (4.15), (4.16), (4) and (4.18). Introduce $`f^1:=ff_0`$. Since $`f_0`$ and $`f`$ are solutions to problems (4.2)-(4.3) and (3.12)-(3.13), respectively, and (3.15) holds, we conclude that $`f^1`$ is the solution to the following problem: in the half-plane $`\tau >0`$,

(5.1)
$$(_\tau +(\theta (\tau )\kappa x)_x\psi (D_x)r(x))f^1(x,\tau )=𝒟_1(ϵ,D_x)f_0(x,\tau ),$$

where

$`𝒟_1(ϵ,\xi )`$ $`:=`$ $`\psi (ϵ,\xi )+{\displaystyle \frac{\sigma ^2}{2}}\xi ^2i\mu \xi `$
$`=`$ $`{\displaystyle \underset{j=3}{\overset{\mathrm{}}{}}}ϵ^{j2}k_j(i\xi )^j`$
$`=`$ $`ϵ{\displaystyle \underset{j=3}{\overset{\mathrm{}}{}}}ϵ^{j3}k_j(i\xi )^j.`$
We also have the initial condition

(5.2)
$$f^1(x,0)=0.$$

From (3.16) and (4.18), the following estimate for the RHS in (5.1) follows:

(5.3)
$$|𝒟_1(ϵ,\xi )\widehat{f}_0(\xi ,\tau )|C_0ϵ\tau ^{1/2}\mathrm{exp}(\xi ^2/(8\tau )),$$

where $`C_0`$ is independent of $`ϵ(0,1)`$ and $`\tau (0,T]`$. By making the inverse Fourier transform, we obtain

(5.4)
$$𝒟_1(ϵ,D_x)f_0(x,\tau )_{C(𝐑\times [0,T])}Cϵ,$$

where $`C`$ is independent of $`ϵ(0,1)`$. By applying the Feynman-Kac theorem to (5.1)-(5.2), the representation of $`f^1`$ in the form of the stochastic integral results:

$$f^1(x,\tau )=E_\tau \left[_\tau ^0\mathrm{exp}\left(_\tau ^sr(X(s^{}))𝑑s^{}\right)𝒟(ϵ,D_x)f_0(x,s)𝑑s\right],$$

and from (5.4), we derive an estimate

(5.5)
$$|f^1(x,\tau )|Cϵ,$$

where $`C`$ is independent of $`ϵ(0,1)`$, $`x𝐑`$ and $`\tau (0,T]`$.

### 5.2. First correction term

Estimate (5.5) shows that $`f_0`$ is indeed the leading term of the asymptotics of $`f`$ as $`ϵ0`$, and in view of (3.15), it is natural to look for the first correction term in the form $`ϵf_1`$, where $`f_1`$ is the solution to the following problem:

(5.6)
$$(_\tau +(\theta (\tau )\kappa x)_x+\frac{\sigma ^2}{2}_x^2r(x))f_1(x,\tau )=k_3_x^3f_0(x,\tau ),$$

in the half-plane $`\tau >0`$, subject to

(5.7)
$$f_1(x,0)=0.$$

We look for $`f_1`$ in the form

(5.8)
$$f_1(x,\tau )=k_3f_0(x,\tau )\stackrel{~}{f}_1(x,\tau )=k_3e^{\mathrm{\Phi }_0(x,\tau )}\stackrel{~}{f}_1(x,\tau ),$$

where $`\mathrm{\Phi }_0`$ is given by (4.5). By substituting into (5.6), we obtain that $`\stackrel{~}{f}_1`$ solves the problem

(5.9)
$$(+\sigma ^2(_x\mathrm{\Phi }_0)_x+e^{\mathrm{\Phi }_0}(e^{\mathrm{\Phi }_0})r(x))\stackrel{~}{f}_1=e^{\mathrm{\Phi }_0}_x^3e^{\mathrm{\Phi }_0},$$

in the half-plane $`\tau >0`$, subject to

(5.10)
$$\stackrel{~}{f}_1(x,\tau )=0,$$

where $``$ is defined by (4.7). Equation (4.6) allows us to simplify (5.9):

$$(+\sigma ^2(2A(\tau )x+B(\tau ))_x)\stackrel{~}{f}_1=e^{\mathrm{\Phi }_0}_x^3e^{\mathrm{\Phi }_0}.$$

We calculate the operator

$`e^{\mathrm{\Phi }_0}_x^3e^{\mathrm{\Phi }_0}`$ $`=`$ $`(e^{\mathrm{\Phi }_0}_xe^{\mathrm{\Phi }_0})^3`$
$`=`$ $`(_x+2A(\tau )x+B(\tau ))^3,`$
and by using $`_x1=0`$, rewrite (5.9) as

(5.11)
$$_1\stackrel{~}{f}_1(x,\tau )=\stackrel{~}{g}_0(x,\tau ),$$

where

$`\stackrel{~}{g}_0(x,\tau )`$ $`:=`$ $`8A^3x^3+12A^2Bx^2+(12A^2+6AB^2)x+6AB+B^3,`$
$`_1`$ $`:=`$ $`_\tau +(\theta _1(\tau )+\kappa _1(\tau )x)_x{\displaystyle \frac{\sigma ^2}{2}}_x^2,`$
$`\theta _1(\tau )`$ $`:=`$ $`\theta (\tau )+\sigma ^2B(\tau ),`$
$`\kappa _1(\tau )`$ $`:=`$ $`\kappa 2\sigma ^2A(\tau ),`$
and $`A=A(\tau )`$, $`B=B(\tau )`$. Clearly, we may look for the solution to (5.11) in the form of a polynomial in $`x`$, of degree 3, with coefficients vanishing at $`\tau =0`$:

(5.12)
$$\stackrel{~}{f}_1(x,\tau )=a_{13}(\tau )x^3+a_{12}(\tau )x^2+a_{11}(\tau )x+a_{10}(\tau ).$$

By substituting into (5.11), we obtain a system of linear ODE:

(5.13) $`a_{13}^{}+3\kappa _1a_{13}`$ $`=`$ $`8A^3,`$
(5.14) $`a_{12}^{}+2\kappa _1a_{12}3\theta _1a_{13}`$ $`=`$ $`12A^2B,`$
(5.15) $`a_{11}^{}+\kappa _1a_{11}2\theta _1a_{12}3\sigma ^2a_{13}`$ $`=`$ $`12A^2+6AB^2,`$
(5.16) $`a_{10}^{}\theta _1a_{11}\sigma ^2a_{12}`$ $`=`$ $`6AB+B^3,`$
which can easily be integrated step by step. Namely, let

(5.17)
$$\kappa _2(s)=_0^s\kappa _1(s^{})𝑑s^{};$$

then

(5.18) $`a_{13}(\tau )`$ $`=`$ $`8e^{3\kappa _2(\tau )}{\displaystyle _0^\tau }e^{3\kappa _2(s)}A(s)^3𝑑s,`$
(5.19) $`a_{12}(\tau )`$ $`=`$ $`e^{2\kappa _2(\tau )}{\displaystyle _0^\tau }e^{2\kappa _2(s)}(12A(s)^2B(s)+3\theta _1(s)a_{13}(s))𝑑s,`$
$`a_{11}(\tau )`$ $`=`$ $`e^{\kappa _2(\tau )}{\displaystyle _0^\tau }e^{\kappa _2(s)}(12A(s)^2+6A(s)B(s)^2`$
$`+2\theta _1(s)a_{12}(s)+3\sigma ^2a_{13}(s))ds,`$
(5.21) $`a_{10}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }(6A(s)B(s)+B(s)^3+\theta _1(s)a_{11}(s)+\sigma ^2a_{12}(s))𝑑s.`$
Formulas (5.8), (5.12) and (5.18)-(5.21) give the first order approximation

(5.22)
$$ff_0(1+ϵk_3\stackrel{~}{f}_1)$$

to the bond price. The proof similar to the one of estimate (5.5) albeit more involved shows that the error of approximation (5.22) is

(5.23)
$$|f(x,\tau )f_0(x,\tau )(1+ϵk_1\stackrel{~}{f}_1(x,\tau ))|Cϵ^2(1+|x|^2)^{3/2},$$

where $`C`$ is independent of $`ϵ(0,1)`$, $`x𝐑`$ and $`\tau (0,T]`$.

Unlike (5.5), we have a polynomially growing factor $`(1+|x|^2)^{3/2}`$ in the RHS of (5.23). Notice, however, that for practical purposes, one needs to know the bond price for small values of $`r(X(t))`$, hence for small values of $`X(t)`$, and therefore the polynomially growing factor $`(1+|x|^2)^{3/2}`$ does not matter much.

### 5.3. First correction term II: the derivation based on the change of variables

To simplify the calculation of the next terms of the asymptotics, it is advantageous to change the variables in equations similar to (5.11):

(5.24)
$$x=\theta _2(\tau )+e^{\kappa _2(\tau )}y,$$

where $`\kappa _2`$ is given by (5.17), and $`\theta _2`$ is the solution to the Cauchy problem

$`\theta _2^{}(\tau )\kappa _2(\tau )\theta _2(\tau )`$ $`=`$ $`\theta _1(\tau ),`$
$`\theta _2(0)`$ $`=`$ $`0,`$
that is,

$$\theta _2(\tau )=e^{\kappa _2(\tau )}_0^\tau e^{\kappa _2(s)}\theta _1(s)𝑑s.$$

The same change of the variables simplifies the calculation of $`\stackrel{~}{f}_1`$. Introduce an operator $`S`$ by $`S(f)(y,\tau )=f(x(y),\tau )`$. Under the change of variables (5.24), $`_\tau +(\theta _1(\tau )\kappa _1(\tau ))_x_\tau `$ and $`_xe^{\kappa _2(\tau )}_y`$, therefore

$$_2:=S^1_1S=_\tau \frac{\sigma ^2}{2}e^{2\kappa _2(\tau )}_y^2,$$

and we can rewrite (5.11) in the form

(5.25)
$$_2F_1(y,\tau )=G_0(y,\tau ),$$

where $`F_1=S\stackrel{~}{f}_1`$, $`G_0=S\stackrel{~}{g}_0`$. Clearly, $`G^0`$ is a polynomial in $`y`$ of the same order as $`\stackrel{~}{g}_0`$:

$$G_0(y,\tau )=G_{0,3}(\tau )y^3+G_{0,2}(\tau )y^2+G_{0,1}(\tau )y+G_{0,0}(\tau ),$$

and the coefficients $`G_{0,j}`$ can easily be calculated by using formulas for the coefficients of $`\stackrel{~}{g}_0`$ or, better, independently. Under the change of variables (5.24), $`_x+2A(\tau )x+B(\tau )`$ becomes

$$𝒟:=e^{\kappa _2(\tau )}(_y+A_1(\tau )y+B_1(\tau )),$$

where

$$A_1(\tau )=2e^{2\kappa _2(\tau )}A(\tau ),B_1(\tau )=e^{\kappa _2(\tau )}(B(\tau )2A(\tau )\theta _2(\tau )),$$

therefore

$`G_0`$ $`=`$ $`𝒟^31`$
$`=`$ $`e^{3\kappa _2}(A_1^3y^3+3A_1^2B_1y^2+(3A_1^2+3A_1B_1^2)y+3A_1B_1+B_1^3),`$
where $`\kappa _2=\kappa _2(\tau ),A_1=A_1(\tau )`$ and $`B_1=B_1(\tau )`$. The solution of (5.25) subject to $`F_1(y,0)=0`$ is a polynomial in $`y`$ of the same order as $`G_0`$:

$$F_1(y,\tau )=F_{1,3}(\tau )y^3+F_{1,2}(\tau )y^2+F_{1,1}(\tau )y+F_{1,0}(\tau ),$$

whose coefficients can easily be found by integration:

(5.27) $`F_{1,3}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }G_{0,3}(s)𝑑s,`$
(5.28) $`F_{1,2}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }G_{0,2}(s)𝑑s,`$
(5.29) $`F_{1,1}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }(G_{0,1}(s)+3\sigma ^2e^{2\kappa _2(\tau )}G_{0,3}(s))𝑑s,`$
(5.30) $`F_{1,0}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }(G_{0,0}(s)+\sigma ^2e^{2\kappa _2(\tau )}G_{0,2}(s))𝑑s.`$
After that we make the inverse change of variables $`y=e^{\kappa _2(\tau )}(x+\theta _2(\tau ))`$, and calculate $`\stackrel{~}{f}_1=S^1F_1`$.

### 5.4. Next terms of the asymptotics

Suppose that the approximation of order $`j1`$ has been found:

$$f=f_0\left(1+\underset{l=1}{\overset{j}{}}ϵ^lk_{l+2}\stackrel{~}{f}_l\right)=f_0\underset{l=0}{\overset{j}{}}ϵ^lk_{l+2}\stackrel{~}{f}_l,$$

where $`k_2=1,f_0=e^{\mathrm{\Phi }_0},\stackrel{~}{f}_01`$, and $`\stackrel{~}{f}_l`$, $`1lj`$, are polynomials in $`x`$ with coefficients depending on on $`\tau `$:

(5.31)
$$\stackrel{~}{f}_l(x,\tau )=\underset{s=0}{\overset{m_l}{}}a_{ls}(\tau )x^s.$$

We look for the next term of the asymptotics in the form $`ϵ^{j+1}k_{j+3}f_0\stackrel{~}{f}_{j+1}`$, where $`f_{j+1}:=f_0\stackrel{~}{f}_{j+1}`$ is the solution to the problem

(5.32) $`(r(x))f_{j+1}(x,\tau )`$ $`=`$ $`g_j(x,\tau ),\tau >0,`$
(5.33) $`f_{j+1}(x,0)`$ $`=`$ $`0,`$
where $`ϵ^{j+1}g_j`$ is the collection of terms of order $`ϵ^{j+1}`$ in the expression

$$\underset{p=3}{\overset{\mathrm{}}{}}ϵ^{p2}k_p_x^p\left(\underset{l=0}{\overset{j}{}}ϵ^lk_{l+2}f_0\stackrel{~}{f}_l\right),$$

that is,

$$g_j=\underset{p=3}{\overset{j+3}{}}k_pk_{j+5p}_x^p(f_0\stackrel{~}{f}_{j+3p}).$$

Set

$`\stackrel{~}{g}_j`$ $`=`$ $`e^{\mathrm{\Phi }_0}g_j`$
$`=`$ $`{\displaystyle \underset{p=3}{\overset{j+3}{}}}k_pk_{j+5p}(e^{\mathrm{\Phi }_0}_xe^{\mathrm{\Phi }_0})^p\stackrel{~}{f}_{j+3p}`$
$`=`$ $`{\displaystyle \underset{p=3}{\overset{j+3}{}}}k_pk_{j+5p}(_x+2A(\tau )x+B(\tau ))^p\stackrel{~}{f}_{j+3p},`$
and multiply (5.32) by $`e^{\mathrm{\Phi }_0}`$. We obtain

(5.35)
$$(+\sigma ^2(_x\mathrm{\Phi }_0)_x+e^{\mathrm{\Phi }_0}(e^{\mathrm{\Phi }_0})r(x))\stackrel{~}{f}_{j+1}=\stackrel{~}{g}_j.$$

Equation (4.6) allows us to simplify (5.35) and obtain a problem in the half-space $`\tau >0`$ with the unknown $`\stackrel{~}{f}_{j+1}`$:

(5.36) $`\left(_\tau +(\theta _1(\tau )\kappa _1(\tau )x)_x+{\displaystyle \frac{\sigma ^2}{2}}_x^2\right)\stackrel{~}{f}_{j+1}`$ $`=`$ $`\stackrel{~}{g}_j,\tau >0,`$
(5.37) $`\stackrel{~}{f}_{j+1}(x,0)`$ $`=`$ $`0.`$
From (5.4), $`\stackrel{~}{g}_j`$ is a polynomial in $`x`$ as well:

(5.38)
$$\stackrel{~}{g}_j(x,\tau )=\underset{s=0}{\overset{m_j^{}}{}}b_{j,s}(\tau )x^s,$$

where $`m_j^{}=\mathrm{max}_{0lj}(m_ll+j+3)`$, and any of the coefficients $`b_{j,s}`$ may be zero, that is, the order of $`\stackrel{~}{g}_j`$ may be less than $`m_j^{}`$. Denote by $`m_j`$ the order of $`\stackrel{~}{g}_j`$.

By making the change of variables (5.24), we simplify (5.36):

(5.39)
$$_2F_{j+1}=G_j,$$

where $`F_j=S\stackrel{~}{f}_j,G_j=S\stackrel{~}{g}_j`$. Since $`G_j`$ is a polynomial in $`y`$ of order $`m_j`$:

$$G_j(y,\tau )=\underset{l=0}{\overset{m_j}{}}G_{j,l}(\tau )y^l$$

(the coefficients $`G_{j,l}`$ will be calculated in Subsection 5.5), $`F_{j+1}`$ also is:

$$F_{j+1}(y,\tau )=\underset{l=0}{\overset{m_j}{}}F_{j+1,l}(\tau )y^l,$$

and the coefficients $`F_{j+1,l}`$ are easily found by integration:

(5.40) $`F_{j+1,l}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }G_{j,l}(s)𝑑s,l=m_j,m_j1;`$
(5.41) $`F_{j+1,l}(\tau )`$ $`=`$ $`{\displaystyle _0^\tau }\left({\displaystyle \frac{\sigma ^2}{2}}e^{2\kappa _2(s)}(l+2)(l+1)F_{j+1,l+2}(s)+G_{j,l}(s)\right)𝑑s,`$
for $`l=m_j2,m_j3,\mathrm{}0.`$

### 5.5. Calculation of $`G_{j,l}`$

We can rewrite (5.4) as

(5.42)
$$G_j=\underset{p=3}{\overset{j+3}{}}k_pk_{j+5p}𝒟^pF_{j+3p},j=0,1,\mathrm{}.$$

By using the initial data $`F_0=1`$, we find $`G_j`$ and $`F_{j+1}`$ from (5.42) and (5.40)-(5.41) step by step. In particular, $`G_0`$ is given by (5.3), and

(5.43)
$$G_1=k_3^2G_{11}+k_4G_{12},$$

where

(5.44)
$$G_{11}=𝒟^3F_1,\mathrm{and}G_{12}=𝒟^4F_0=𝒟^41$$

are polynomials of degree 6 and 4, respectively. Hence,

(5.45)
$$F_1=k_3^2F_{21}+k_4F_{22},$$

where $`F_{21}`$ and $`F_{22}`$ are polynomials of degree 6 and 4, respectively, which solve (5.39) with $`G_{11}`$ and $`G_{12}`$ in the RHS (and satisfy the initial condition $`F_{2j}(y,0)=0`$). The reader can use (5.40)-(5.41) to obtain explicit formulas for coefficients of $`F_{2j}`$ in terms of coefficients of $`G_{1j}`$. Notice that though the latter can be written explicitly, in practical implementation of the method, it is simpler to write a program which calculates the coefficients of a polynomial $`𝒟P`$, given coefficients of a polynomial $`P`$, and use this program to calculate $`G_{1j}`$ (and $`G_l`$ for $`l>1`$ should one wish it, though in applications, it seems not a reasonable thing to do).

### 5.6. Second order approximation for the bond price, yield and forward rate

We see that the second order approximation can be written in the form

(5.46)
$$fe^{\mathrm{\Phi }_0}(1+ϵk_3\stackrel{~}{f}_1+(ϵk_3)^2\stackrel{~}{f}_{21}+ϵ^2k_4\stackrel{~}{f}_{22}),$$

where $`\mathrm{\Phi }_0`$, $`\stackrel{~}{f}_1=S^1F_1`$, $`\stackrel{~}{f}_{21}=S^1F_{21}`$ and $`\stackrel{~}{f}_{22}=S^1F_{22}`$ are polynomials in $`x`$ with coefficients depending on $`\tau `$ \- and parameters $`\stackrel{~}{\theta }`$, $`\kappa `$, $`\mu :=i\psi ^{}(0)`$, $`\sigma ^2:=\psi ^{\prime \prime }(0)`$. For practical applications,(5.46) can be rewritten in the form

(5.47)
$$fe^{\mathrm{\Phi }_0}\left(1+K_3\stackrel{~}{f}_1+K_3^2\stackrel{~}{f}_{21}+K_4\stackrel{~}{f}_{22}\right),$$

where

$`K_3:`$ $`=`$ $`ϵk_3=i\psi ^{(3)}(0)/3!,`$
$`K_4:`$ $`=`$ $`ϵ^2k_4=\psi ^{(4)}(0)/4!`$
can be inferred from the data.

Denote by $`P:=f`$ is the price of the bond, and by $`P_0:=\mathrm{exp}(\mathrm{\Phi }_0)`$ the price in the Gaussian model; then (5.47) becomes

(5.48)
$$PP_0+K_3P_1+K_3^2P_{21}+K_4P_{22},$$

where

$$P_1=P_0\stackrel{~}{f}_1,P_{21}=P_0\stackrel{~}{f}_{21},P_{22}=P_0\stackrel{~}{f}_{22}.$$

By using the formulas for the yield

$$R(x,\tau )=\frac{\mathrm{ln}P(x,\tau )}{\tau }$$

and forward rate

$$F(x,\tau )=\frac{}{\tau }\mathrm{ln}P(x,\tau ),$$

we obtain approximate formulas

(5.49)
$$RR_0+K_3R_1+K_3^2R_{21}+K_4R_{22},$$

where

$`R_0(x,\tau )`$ $`=`$ $`\mathrm{\Phi }_0(x,\tau )/\tau ,`$
$`R_1(x,\tau )`$ $`=`$ $`\stackrel{~}{f}_1(x,\tau )/\tau ,`$
$`R_{21}(x,\tau )`$ $`=`$ $`(\stackrel{~}{f}_1(x,\tau )^2/2\stackrel{~}{f}_{21}(x,\tau ))/\tau ,`$
$`R_{22}(x,\tau )`$ $`=`$ $`\stackrel{~}{f}_{22}(x,\tau )/\tau ,`$
and

(5.50)
$$FF_0+K_3F_1+K_3^2F_{21}+K_4F_{22},$$

where

$`F_0(x,\tau )`$ $`=`$ $`{\displaystyle \frac{}{\tau }}\mathrm{\Phi }_0(x,\tau ),`$
$`F_1(x,\tau )`$ $`=`$ $`{\displaystyle \frac{}{\tau }}\stackrel{~}{f}_1(x,\tau ),`$
$`F_{21}(x,\tau )`$ $`=`$ $`{\displaystyle \frac{}{\tau }}(\stackrel{~}{f}_1(x,\tau )^2/2\stackrel{~}{f}_{21}(x,\tau )),`$
$`F_{22}(x,\tau )`$ $`=`$ $`{\displaystyle \frac{}{\tau }}\stackrel{~}{f}_{22}(x,\tau ).`$
### 5.7. Numerical examples

(The author thanks Nina Boyarchenko for writing the programs for numerical examples and checking the algebra in the previous two sections.) We take the simplest model for $`r`$: $`r=x^2`$, and constant $`\stackrel{~}{\theta }(t)=0.06`$, $`\kappa =0.3`$, $`\mu =0`$ and $`\sigma ^2=0.08`$. We also fix $`x=0.25`$, and study shapes of correction terms to the bond price, yield and forward rates in (5.48), (5.49) and (5.50) (see Figures 1–3).

From Figure 3, we clearly see that it is the first correction term, $`F_1`$, that can account for the hump of the forward rate curve - provided the skewness is negative and not too small in modulus. In the next three Figures 4–6, we plot the bond price, yield and forward rate; first, the leading term (dots), then the formula with the first correction term taken into account (solid line), and finally, the formula with the two correction terms (dotted line). We take the same parameters as above, and $`K_3=\sigma ^2/16=0.005`$, $`K_4=K_3/20=2.510^4`$. We see that fairly large skewness does produce a hump-shaped forward rate curve, when the Gaussian curve has no hump; the asymptotic formulas are applicable since $`K_3`$ is small w.r.t. $`\sigma ^2`$, and $`K_4`$ is small w.r.t. $`K_3`$.

In the last series of figures (Figures 7–9), we fix small $`K_4=\sigma ^2/200`$, and show how the shapes of the curves vary with the skewness.

Notice that if one fits the Gaussian QTSM and the non-Gaussian one to the same data set, then the leading Gaussian-like leading term in the Lévy model will be determined by slightly different drift and volatility parameters and the spot value of the factor, hence the pictures above do not describe quite accurately the difference between the Gaussian and Lévy models. However, the shapes of the curves in the Gaussian model and the leading Gaussian approximation in the non-Gaussian models are the same (and not differ much). Thus, one can use the pictures to get a feeling what a difference between the Gaussian and non-Gaussian QTSMs can be.

## 6. Extensions and ramifications

### 6.1. Multi-factor case

All the constructions in the previous two sections admit modifications for the multi-factor case. The only two differences are

* in the Gaussian approximation, the $`A(\tau )`$, $`B(\tau )`$ and $`C(\tau )`$ are matrix functions. The first two can be found by solving a system of linear ODE’s (for detailed realization in the Gaussian QTSM-model, see Kim (2003)), and then $`C(\tau )`$ is found by the integration;
* the correction terms are polynomials not in one factor but in several factors. This leads to systems of linear ODE’s whose unknowns are coefficients at factors $`x^\alpha =x_1^{\alpha _1}\mathrm{}x_n^{\alpha _n}`$ of order $`\alpha =\alpha _1+\mathrm{}+\alpha _nm_j`$, where $`m_j`$ depends on the step of the method.

Certainly, the systems of linear ODE to be solved become significantly larger but they remain linear systems; therefore, the numerical solution is fairly stable.

### 6.2. Pricing under historic and risk-neutral measures

Assume that under the historic measure, the dynamics of $`X_t`$ is given by

(6.1)
$$dX(t)=(\stackrel{~}{\theta }(t)\kappa X(t))dt+dZ^𝐏(t),$$

where $`\{Z^𝐏(t)\}`$ is an $`n`$-dimensional Lévy process with the characteristic exponent $`\psi ^𝐏`$. We assume that $`\psi ^𝐏`$ admits the analytic continuation into the tube domain $`𝐑^n+iU^𝐏`$, where $`U^𝐏`$ is an open set containing 0. To specify the interest rate dynamics under a risk-neutral measure, $`𝐐`$, we consider first the state price deflator in the form $`\pi _t=\mathrm{exp}(q_t)`$, where $`q_t`$ obeys the SDE

(6.2)
$$dq_t=r_tdt\mathrm{\Lambda },dZ^𝐏(t)\mathrm{\Lambda }_0t.$$

The vector $`\mathrm{\Lambda }𝐑^n`$ represent the market prices of risk of the factors. Notice that for processes with jumps, one cannot use an arbitrary $`\mathrm{\Lambda }`$. For the bond to be priced, it is necessary that $`\mathrm{\Lambda }\overline{U^𝐏}`$, and the condition $`\mathrm{\Lambda }U^𝐏`$ suffices for the bond and options on the yield to be priced (to be more specific, any payoff which admits a polynomial bound with respect to factors is admissible; if a payoff grows exponentially, then additional restrictions on $`\mathrm{\Lambda }`$ must be imposed).

Set $`\psi ^𝐐(\xi )=\psi ^𝐏(\xi +i\mathrm{\Lambda })\psi ^𝐏(i\mathrm{\Lambda })`$. It is easy to check that under the condition $`\mathrm{\Lambda }U^𝐏`$, $`\psi ^𝐐`$ is the characteristic exponent of a Lévy process, call it $`Z^𝐐`$. The new measure $`𝐐`$ is the Esscher transform of $`𝐏`$ popular in the literature on the pricing of options on stocks (see e.g. Eberlein et al. (1998) and Boyarchenko and Levendorskiǐ (2002b)). Notice that $`\psi ^𝐐`$ is analytic in the tube domain $`U^𝐐=U^𝐏\mathrm{\Lambda }\{0\}`$, and choose $`\mathrm{\Lambda }_0=\psi ^𝐏(i\mathrm{\Lambda })`$. The straightforward calculations show that the pricing formula

(6.3)
$$f(x,t)=E_t^𝐏\left[\frac{\pi _t}{\pi _T}g(X(T))|X(t)=x\right]$$

can be written as the pricing formula under the risk-neutral measure $`𝐐`$:

(6.4)
$$f(x,t)=E_t^𝐐\left[\mathrm{exp}\left(_t^Tr_s𝑑s\right)g(X(T))|X(t)=x\right].$$

It is possible that in some situations, the specification of risk (6.2) is not sufficiently flexible, and in the Gaussian QTSM, Ahn et al. (2002) consider a more general model for prices of risk. This more general specification is not applicable in Lévy models with exponentially decaying Lévy densities; however, there is an additional flexibility in model with jumps, which may provide any number of additional degrees of freedom: one can use the pricing formula (6.4) with any Lévy process $`Z^𝐐`$ provided the difference of Lévy densities of $`𝐏`$ and $`𝐐`$ is of finite variation (this can be shown as in the case of pricing of derivatives on a stock; see e.g. Carr et al. (2002)). An additional restriction should be taken into account if we want to use the approximate formulas of Sections 4–5: the Lévy densities under $`𝐏`$ and $`𝐐`$ should decay exponentially, and sufficiently fast.

Notice that contrary to the change of measure in the Gaussian model, the change of measure in the non-Gaussian model may lead to the changes (albeit small) in the instantaneous moments of order two.

### 6.3. Option pricing

Let $`f(x,t)`$ be the price of an European style derivative contract with the pay-off $`g(X(T_1))`$, where $`T_1<T`$. A typical example is a call option on the yield with the pay-off $`g(X(T_1))=\mathrm{max}\{R(X(T_1)K,0\}`$. For a fixed $`T_1<T`$, the formulas in the preceding section allows one to find an approximation to $`g(X(T_1))`$ as a function of the spot values of factors $`x=X(t)`$ (and of $`\tau =TT_1`$). In the Gaussian model, this is a quadratic function, and in the non-Gaussian model, it is a fourth-order polynomial with small coefficients of order 3 and 4. Hence, we can calculate the two roots of the equation $`R(x)=K`$, which are close to the two roots in the Gaussian approximation, by using simple perturbation technique. Denote these roots by $`x_ϵ^\pm (K)`$, and represent the pay-off function $`g(x)`$ in the form

$$g(x)=g^+(x)+g^{}(x),$$

where

$$g^+(x)=\mathrm{𝟏}_{[x_ϵ^+,+\mathrm{})}(x)(R(x)K),$$

$$g^{}(x)=\mathrm{𝟏}_{(\mathrm{},x_ϵ^{}]}(x)(R(x)K).$$

By using the perturbation technique once again, we can calculate the Fourier transforms of functions $`g^\pm `$, and reduce the problem of the calculation of $`f(x,t)`$ to the family of problems considered in Sections 4–5, with the pay-off $`e^{ix\xi }`$ instead of 1. The modification of the asymptotic calculations is straightforward albeit lengthy. In the end, we make the inverse Fourier transform, and obtain an asymptotic formula for the price $`f(x,t)`$, $`t<T_1`$. For a numerical realization of this formula, one needs to choose an appropriate grid in the frequency domain, and for each $`\xi `$ from the grid, solve a number of systems of ODE’s. The inverse FFT finishes the job. Notice that this is a variant of the standard transform method (see e.g. Duffie et al. (2000), Chacko and Das (2002) and the bibliography therein).

When using this scheme, one should remember that due to the non-smoothness of the pay-off at the money, the approximate formulas will not work well near expiry, especially near expiry and strike. Fortunately, near the strike, a different approximation – and much simpler one – can be derived. Fix $`x_0`$, and introduce

(6.5)
$$r(x_0;x)=R_0+R^{}(x_0)(xx_0).$$

Denote by $`f^0(x_0;x,t)`$ the solution to the affine model with the short rate modelled by (6.5), and the same dynamics of the factors under the risk-neutral measure and pay-off as in the initial QTSM.

###### Theorem 6.1.

$`f(x_0,t)=f^0(x_0;x_0,t)+o(\tau )+O(ϵ^{\mathrm{}})`$ as $`\tau =T_1t+0`$ and $`ϵ0`$.

The proof of Theorem 6.1 is relatively straightforward, and it relies on the exponential decay of the density of jumps. Since the solution in the affine model with jumps is well-known and fairly simple, we get an efficient approximation to the option price in the non-Gaussian QTSM. Should one wish it, the correction term to the formula $`f(x_0,t)f^0(x_0;x_0,t)`$ can be obtained as the price of a derivative security (in the same affine term structure model), which pays the stream of dividends at rate $`g(x)`$. For in-the-money options, one can use a constant function $`g(x_0)`$ instead of $`g(x)`$.

Finally, for out-of-the-money options, a simpler approximation can be derived, in the form

(6.6)
$$f(x,t)𝒞(g;x)\tau ,$$

where $`𝒞(g;x)`$ depends only on the Lévy density $`F(dx)`$ of the process but not on the Gaussian part and drift. For the case of the call option on the yield considered above, essentially the same calculations as in Levendorskiǐ (2003) give for $`x(x_ϵ^{},x_ϵ^+)`$

(6.7)
$$𝒞(x,t)=_{\mathrm{}}^+\mathrm{}g(x+y)F(dy),$$

and the explicit analytic expression in terms of parameters of model classes of RLPE processes can be derived (the modification of the calculations in Levendorskiǐ (2003) made for $`g(x)=(e^x1)^+`$ and $`g(x)=(1e^x)^+`$ is straightforward).

### 6.4. Parameters’ fitting

Far from expiry, the leading term of the yields depends on the instantaneous moments of order one and two, and for out-of-the-money options on yields, near expiry, the leading term depends only on the jump part of the process. This allows us to suggest the following scheme of the fitting the model to the data, under a risk-neutral measure.

1. Infer parameters of the Gaussian model (including the spot values of the factors) from the data on yields. Here one can use the efficient method of moments as in Ahn et al. (2002) or an estimation method integrated with the extraction of the state variables (the extended-Kalman-filter-based quasi-maximum-likelihood estimate) as in Kim (2003), say. We regard these parameters’ values as a zero-order approximation to the spot values of factors, drift, mean reverting and variance-covariance parameters of the non-Gaussian model.
2. Given the spot values of the factors, one can infer the conditional characteristic function of the process from empirical data as in Singleton (1999) and calculate the moments up to order 4, or infer these moments as in Collin-Dufresne et al. (2003). However, in order to identify the contribution of jumps more accurately, the following steps seems to be reasonable.
3. Choose a parametrized model for the Lévy density (or one of the standard Lévy models described in Section 2), and by using the prices on interest rate derivatives near expiry and approximate formulas near expiry, fit the parameters of the Lévy density, and calculate the instantaneous moments of order 3 and 4.
4. Calculate the correction terms to the yields in the Gaussian approximation by using the asymptotic formulas of Section 5 and the zero-order approximations for the parameters in the Gaussian approximation, and moments of order 3 and 4.
5. Subtract the correction terms from the data, and use the new data set in the Gaussian procedure to infer the corrected values of the spot factors and the first two moments.
6. Step 3 can be repeated in order to get a corrected specification of the Lévy density and moments of order 3 and 4.

Notice that in order that this procedure be consistent, the resulting moments of order 3 and 4 must be small relative to the smallest eigenvalue of the matrix of the instantaneous second moments.

## 7. Conclusion

We constructed a class of QTSM models with a regular Lévy process of exponential type in place of the Gaussian one in standard QTSM. By using the Feynman-Kac formula, we have reduced the pricing problem for an interest rate derivative to a boundary problem for a pseudo-differential operator. In the case of the bond, we have found an approximate solution to the boundary problem assuming that the tails of probability densities of a process decay sufficiently fast. The leading term of the approximate solution looks as in the Gaussian model (even when the underlying process has no Gaussian component), and the correction terms depend on skewness and kurtosis.

Numerical examples are produced to show that by changing skewness and kurtosis, various shapes of the forward rate curve can be obtained. In particular, negative skewness can produce a hump-shaped forward rate curve even when the Gaussian curve has no hump: the very non-Gaussianity of the process is (one of) causes of the hump of the forward rate curve. Bond prices and the yield curve also change but the types of the shape of the curves remain essentially the same.

We discussed possible choices of a risk-neutral measure, and indicated additional flexibility which the usage of jumps provides. A brief outline of asymptotic pricing in the multi-factor case, and the pricing of interest rate derivatives is given. We derived simple asymptotic formulas for option prices near expiry, and suggested a procedure of the fitting of the model to the data.

Notice that the use of pseudo-diffusions for option pricing instead of Gaussian models with approximately the same number of parameters is advantageous not only near expiry where the latter are expected to produce serious errors. If the number of parameters is approximately the same, then the number of factors in the former is smaller than in the latter, therefore the use of the inverse FFT simplifies significantly in the non-Gaussian model, and the numerical procedure becomes much more stable.

## Appendix A Proofs of technical results

### A.1. Proof of Theorem 3.1

By using the decomposition $`g=g_+g_{}`$, where $`g_+(x):=\mathrm{max}\{g(x),0\}`$ and $`g_{}=g_+g`$ are non-negative, we see that it suffices to prove Theorem 3.1 for continuous non-negative $`g`$. Fix $`\chi C^{\mathrm{}}(𝐑)`$ such that $`0\chi (x)1`$ for all $`x`$, $`\chi (x)=1,x1`$, $`\chi (x)=0,x2`$, and for any $`m>0`$, set $`g^m(x):=\chi (|x|/m)g(x)`$. Then $`g^m`$ is a continuous function with the compact support. Define $`f^m`$ by (3.4) with $`g^m`$ in the RHS. For any $`x`$, $`g^m(x)g(x)`$ as $`m\mathrm{}`$, and by the Monotone Convergence Theorem, $`f^m(x)f(x)`$. Notice that $`f^mf`$ in the sense of generalized functions: for any non-negative $`uC_0^{\mathrm{}}(𝐑^n\times (0,T))`$,

$$_{𝐑^n}f^m(x)u(x)𝑑x_{𝐑^n}f(x)u(x)𝑑x,m\mathrm{}.$$

Below we will show that

1. let $`g`$ be continuous and satisfy (3.5); then problem (3.12)-(3.13) has a unique continuous solution $`f(g;x,t)`$, which satisfies (3.11);
2. in the half-space $`t<T`$, $`f`$ is of the class $`C^{2,1}`$ w.r.t. $`(x,t)`$;
3. $`f(g^m;,)f(g;,)`$ as $`m\mathrm{}`$, in the sense of generalized functions;
4. if $`g`$ is a continuous function of the compact support, then $`f(x,t):=f(g;x,t)`$ is given by (3.4).

By (iv), $`f(g^m;,)=f^m(,)`$, by (iii) and (ii), the limit in (iii) is a continuous function, and since we already know that $`f^m(x,t)f(x,t)`$ pointwise, we conclude that $`f(,)=f(g;,)`$ solves the problem (3.12)-(3.13), and finish the proof of Theorem 3.1.

It remains to prove (i)-(iv). We start with the proof of (i). Assume that $`\kappa `$ is diagonalizable: there exists a matrix $`C`$ such that $`\kappa _C:=C^1\kappa C`$ is a diagonal matrix with the diagonal entries $`\kappa _j`$ (if $`\kappa `$ is not diagonalizable, an additional step is to be made - see the end of the proof of (i)). By making the change of variables $`x=Cy`$, we reduce to the case of the diagonal matrix $`\kappa (=\kappa _C)`$; the $`\psi (\xi )`$ becomes $`\psi _C(\xi ):=\psi ((C^{})^1\xi )`$. To simplify the notation below (and without loss of generality), we assume that $`\kappa `$ itself is diagonal. In (3.12)-(3.13), change the variables:

$$t=T\tau ,x_j=\stackrel{~}{\theta }_{2j}(\tau )+e^{\kappa _j\tau }y_j,j=1,\mathrm{},n,$$

where

$$\stackrel{~}{\theta }_{2j}(\tau )=e^{\tau \kappa _j}_0^\tau e^{s\kappa _j}\stackrel{~}{\theta }_j(Ts)𝑑s$$

is the solution to ODE

$$\stackrel{~}{\theta }_{2j}^{}(\tau )\kappa \stackrel{~}{\theta }_{2j}(\tau )=\stackrel{~}{\theta }_j(T\tau ),$$

subject to $`\stackrel{~}{\theta }_{2j}(0)=0`$, and set

$`v(y,\tau )`$ $`=`$ $`f(x,t),`$
$`e^{\tau \kappa }`$ $`=`$ $`\mathrm{diag}(e^{\kappa _j\tau }),`$
$`r_1(y,\tau )`$ $`=`$ $`r(\stackrel{~}{\theta }_2(\tau )+e^{\tau \kappa }y).`$
We obtain

(A.1) $`(_\tau +\psi (e^{\kappa \tau }D_y)+r_1(y,\tau ))v(y,\tau )`$ $`=`$ $`0,\tau >0,`$
(A.2) $`v(y,0)`$ $`=`$ $`g(y).`$
Notice that $`\stackrel{~}{\theta }_1C^2([0,T])`$ since $`\stackrel{~}{\theta }C([0,T])`$, and $`r_1`$ satisfies estimate

(A.3) $`|_y^\alpha _\tau ^sr_1(y,\tau )|`$ $``$ $`C_{\alpha ,s}(1+|y|)^{2|\alpha |},|\alpha |,s=0,1,2,`$
(A.4) $`c_0|y|^2C_0`$ $``$ $`r_1(y,\tau ),`$
where $`c_0>0`$ and $`C_0,C_{\alpha ,s}`$ depend on $`T`$ but not on $`x𝐑^n`$ and $`\tau [0,T]`$.

Estimates (3.8), (3.9), (3.10), and (A.3)-(A.4) allows one to apply the standard technique of construction of the inverse to the operator of a boundary problem for PDO to problem (A.1)-(A.2). This technique is based on the construction of an appropriate partition of unity, localization and patching of an approximate inverse from local inverses; for the realization of this general scheme for many classes of PDO see Levendorskiǐ (1993). In op. cit., boundary value problems in $`L_p`$-based spaces were considered whereas here we need corresponding results for $`C^s`$-based spaces. This modification is straightforward: see e.g. the modification in Barndorff-Nielsen and Levendorskiǐ (2001), for a different class of PDO.

In the result, we obtain that $`v`$, the continuous solution to problem (A.1)-(A.2), which admits estimate (3.11), exists and it is unique. Moreover, it is of the class $`C^{2,1}`$ in the half-plane $`\tau >0`$, and satisfies estimate

(A.5)
$$\underset{𝐑^n\times [0,T]}{sup}|e^{\omega |y|}v(y,\tau )|C\underset{𝐑^n}{sup}|e^{\omega |y|}g(y)|,$$

where $`C`$ depends on $`T`$, $`\kappa `$, $`\stackrel{~}{\theta }`$ and the constants in estimates for $`\psi `$ and $`r`$. By making the inverse changes of variables and unknowns, we obtain (i) and (ii).

If $`\kappa `$ is not diagonalizable, then prior to the change of variables $`x=Cy`$, an additional change of variables $`x_je^{\rho _j}x_j,j=1,\mathrm{},n`$, where $`\rho _j>0`$, is needed. The $`\kappa `$ will be replaced by $`\kappa \mathrm{diag}(\rho _j)`$, which generically has pairwise distinct eigenvalues and hence, diagonalizable; $`\stackrel{~}{\theta }`$ will change as well, and $`\psi (D_x)`$ becomes $`\psi (e^{\rho _1\tau }D_{x_1},\mathrm{},e^{\rho _n\tau }D_{x_n}).`$ After that we make the same changes of variables (using the new $`\kappa `$ and $`\stackrel{~}{\theta }`$).

To prove (iii), we take $`\omega _1(\omega ,\lambda )`$, and apply the argument above starting with $`gg^m`$ instead of $`g`$ and $`\omega _1`$ instead of $`\omega `$. Since

$$\underset{𝐑^n}{sup}|e^{\omega _1|x|}(g(x)g^m(x))|0\mathrm{as}m\mathrm{},$$

estimate (A.5) implies that

$$\underset{𝐑^n\times [0,T]}{sup}|e^{\omega _1|x|}(f(x,t)f^m(x,t))|0\mathrm{as}m\mathrm{},$$

which proves (iii).

It remains to prove (iv). We have seen that for a continuous $`g`$ with compact support, $`f`$ is continuous in the half-plane $`t<T`$, and of the class $`C^{2,1}`$ in the open half-plane. Moreover, $`f(x,\tau )`$ decays faster than $`e^{\omega |x|}`$ as $`x\mathrm{}`$, for any $`\omega >\lambda `$ (notice that a continuous $`g`$ of the compact support satisfies (3.5) for any $`\omega `$, and in order that the proof of estimate (3.11) remain valid, we may use any $`\omega >\lambda `$, negative ones in particular). Further, $`r`$ is continuous and semi-bounded from below. These conditions are more than sufficient for the Feynman-Kac theorem to be applicable (for instance, at this stage, we can repeat the proof on p.274 in Rogers and Williams (1994)), which gives (iv). Theorem 3.1 has been proved.

REFERENCES

Ahn, D.-H., R.F. Dittmar, and A.R. Gallant (2002): Quadratic term structure

models: theory and evidence, Review of Financial Studies 15:1, 243-288.

Ahn, D.-H., R.F. Dittmar, A.R. Gallant, and Bin Gao (2003): Pure bred or

hybrid?: Reproducing the volatility in term structure dynamics, Journal of Econo-

metrics 116 (2003), 147–180.

Barndorff-Nielsen, O.E. (1998): Processes of Normal Inverse Gaussian Type,

Finance and Stochastics 2, 41–68.

Barndorff-Nielsen, O.E. and W. Jiang (1998): An initial analysis of some Ger-

man stock price series, Working Paper Series 15. Aarhus: CAF Univ. of Aarhus/Aarhus

School of Business.

Barndorff-Nielsen, O.E., and S. Levendorskiǐ (2001): Feller Processes of Nor-

mal Inverse Gaussian Type, Quantitative Finance 1, 318-331.

Barndorff-Nielsen, O.E., E. Nicolato, and N. Shephard (2002): Some recent

developments in stochastic volatility modelling, Quantitative Finance 2, 11-23.

Barndorff-Nielsen, O.E. and N. Shephard (2001a): Normal modified stable

processes, Working paper MaPhySto: Aarhus University

Barndorff-Nielsen, O.E. and N. Shephard (2001b): Non-Gaussian Ornstein-

Uhlenbeck- based models and some of their uses in financial economics (with discus-

sion), J. Royal. Stat. Soc. 63, 167-241.

Bouchaud, J-P. and M. Potters, M. (1997): Theory of Financial Risk. Aléa-

Saclay, Eurolles: Paris

Boyarchenko, S.I. and S.Z. Levendorskiǐ (2000): Option pricing for truncated

Lévy processes, Intern. Journ. Theor. and Appl. Finance 3, 549-552.

Boyarchenko, S.I. and S.Z. Levendorskiǐ (2002a): Perpetual American options

under Lévy processes, SIAM J. Control and Optimization 40, 1663-1696.

Boyarchenko, S.I. and S.Z. Levendorskiǐ (2002b): Non-Gaussian Merton-Black-

Scholes theory. World Scientific: Singapore

Boyarchenko, S.I. and S.Z. Levendorskiǐ (2002c): Barrier options and touch and

out options under regular Lévy processes of exponential type, Annals of Applied Pro-

bability, 12:4, 1261–1298.

Boyarchenko, S.I. and S.Z. Levendorskiǐ (2002d): Pricing of perpetual Bermudan

options, Quantitative Finance, 2:6, 432–442.

Carr, P., H. Geman, D.B. Madan, M. Yor (2002): The fine structure of asset

returns: an empirical investigation, Journal of Business 75, 305-332.

Chacko, G., and S. Das (2002): Pricing interest rate derivatives: a general approach,

Review of Financial Studies 15:1, 195-241.

Chen, L. and H.V. Poor (2002): A general characterization of quadratic term struc-

ture models, Working paper

Chen, P. and O. Scaillet (2002): Linear-quadratic jump-diffusion modeling with

application to stochastic volatility, Working paper

Collin-Dufresne, P., R.S. Goldstein, and C. Johnes (2003): Identification

and estimation of “maximal” affine term structure models: an application to stochastic

volatility, Working paper

Cont, R., M. Potters, and J.-P. Bouchaud (1997): Scaling in stock market data:

stable laws and beyond, in Scale Invariance and beyond (Proceedings of the CNRS

Workshop on Scale Invariance, Les Houches, March 1997). B. Dubrulle, F. Graner,

and D. Sornette, eds., Springer: Berlin, 75–85.

Das, S.R. (2002): The surprise element: jumps in interest rates, Journal of Economet-

rics 106, 27-65.

Duffie, D., and R. Kan (1996): A yield-factor model of interest rates, Mathematical

Finance 6, 376-406.

Duffie, D., J. Pan, and K. Singleton (2000): Transform analysis and option

pricing for affine jump-diffusions, Econometrica 68, 1343-1376.

Duffie, D., D. Filipović, and W. Schachermayer (2002): Affine processes and

applications in Finance, Annals of Applied Probability

Eberlein, E., U. Keller, and K. Prause (1998): New insights into smile, mispric-

ing and value at risk: The hyperbolic model, Journ. of Business 71, 371–406.

Eberlein, E., and S. Raible (1999): Term structure models driven by general Lévy

process, Mathematical Finance 9 , 31-53.

Eberlein, E., and F. Ozkan (2001): The defaultable Lévy Term Structure: Ratings

and Restructuring, Preprint 71, University of Freiburg.

Fama, E.F. (1965): The behavior of stock market prices, Journ. of Business 38, 34–105.

Don H. Kim (2003): Time-varying risk and return in the quadratic-gaussian model of

the term structure, Working paper

Koponen, I. (1995): Analytic approach to the problem of convergence of truncated

Lévy flights towards the Gaussian stochastic process, Physics Review E 52, 1197–1199.

Kudryavzev, O.E., and S.Z. Levendorskiǐ (2002): Comparative study of first

touch digitals: normal inverse gaussian vs. gaussian modelling, Working paper

MaPhySto: Aarhus, October 2002.

Levendorskiǐ, S. (1993): Degenerate elliptic equations. Mathematics and its Applica-

tions, 258. Kluwer Academic Publishers Group: Dordrecht

Levendorskiǐ, S. (2003): The American put and European options near expiry, under

Lévy processes, Working paper

Madan, D.B., P. Carr, and E.C. Chang (1998): The variance Gamma process and

option pricing, European Finance Review 2, 79–105.

Mandelbrot, B.B. (1963): The variation of certain speculative prices, Journ. of

Business 36, 394-419.

Matacz, A. (2000): Financial modelling and option theory with the truncated Lévy

process, Intern. Journ. Theor. and Appl. Finance 3:1, 143-160.

Rogers, L.C.G., and D. Williams (1994): Diffusions, Markov processes, and mar-

tingales”, John Wiley and Sons: West Sussex UK

K. Sato (1999), Lévy processes and infinitely divisible distributions, Cambridge Univer-

sity Press: Cambridge

Singleton, K. (1999): Estimation of affine asset pricing models unsing the empirical

characteristic function, Journal of Econometrics 102, 111–141