{'fraction_non_alphanumeric': 0.04787778838642479, 'fraction_numerical': 0.018106067181668257, 'mean_word_length': 4.227404718693285, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaluating Callable and Putable Bonds: An Eigenfunction Expansion Approach

We propose an efficient method to evaluate callable and putable bonds under a wide class of interest rate models, including the popular short rate diffusion models, as well as their time changed versions with jumps. The method is based on the eigenfunction expansion of the pricing operator. Given the set of call and put dates, the callable and putable bond pricing function is the value function of a stochastic game with stopping times. Under some technical conditions, it is shown to have an eigenfunction expansion in eigenfunctions of the pricing operator with the expansion coefficients determined through a backward recursion. For popular short rate diffusion models, such as CIR, Vasicek, 3/2, the method is orders of magnitude faster than the alternative approaches in the literature. In contrast to the alternative approaches in the literature that have so far been limited to diffusions, the method is equally applicable to short rate jump-diffusion and pure jump models constructed from diffusion models by Bochner's subordination with a L\'{e}vy subordinator.


Introduction
A large fraction of all corporate and sovereign bond issues in the global financial markets have embedded options. The call option allows the bond issuer, such as a corporation or a government, to buy the bond back from the bond holder (call the bond) for pre-specified call prices at some pre-specified times prior to bond's maturity. This allows the bond issuer to refinance the bond if interest rates decline. The put option allows the bond holder to sell (put) the bond back to the bond issuer for pre-specified put prices at some pre-specified times prior to maturity. This allows the bond holder to re-invest the proceeds into a bond with a higher coupon if interest rates rise. The bond with both a call and a put option can be analyzed as an instance of a stochastic game with stopping times (also known as Dynkin games as they have been introduce by Dynkin (1969) as a generalization of optimal stopping problems) driven by the underlying stochastic interest rate model. The bond issuer and the bond holder are opposing players whose opposing optimal strategies are to minimize and to maximize the bond value, respectively. When the call and put decisions can be made at discrete times (typically an advance notice has to be given to the other party ahead of a coupon payment date, when the option is exercised), this sets up a stochastic game with stopping times in discrete time. The value function and the optimal call and put policies can then be determined by solving Bellman's dynamic programming backward induction, starting from maturity and rolling back recursively through the decision times. Developing solution methods for this problem is of significant practical importance.
The pricing of bonds with embedded options has attracted considerable interest in the literature over the years. It goes back to Brennan and Schwartz (1977) who used a finite-difference approach in timehomogeneous diffusion models. Later Büttler (1995) showed that when finite-difference methods were used to price callable bonds under the Vasicek model, the presence of slowly decaying oscillations in the solution after each coupon/call date resulted in poor numerical accuracy. This led Büttler and Waldvogel (1996) (BW) to develop an alternative method for pricing callable bonds under the Vasicek and CIR models utilizing the explicit form of the Green's function in these models. The method relies on the interpolation of the value function and on the numerical quadrature procedure for the integration involving the value function and the Green's function. More recently d 'Halluin et al. (2001) (DFVL) showed that finite-difference methods for these problems can be stabilized via van Leer flux limiter, appropriate non-uniform time stepping schemes, and careful consideration of boundary conditions and presented numerical experiments demonstrating that properly formulated finite-difference methods, in fact, outperformed alternative approaches. Other recent works on applying numerical PDE methods to callable and putable bonds in diffusion interest rate models include Farto and Vázquez (2005) and de Frutos (2008). In Farto and Vázquez (2005), the convection dominated diffusion equation is solved numerically by combining the characteristics method with piecewise-linear Lagrange finite elements. de Frutos (2008) (F) proposes a spectral numerical method for pricing callable bonds. The holding value function is approximated as a finite summation involving the Laguerre polynomials, and the problem is converted to a stiff system of ordinary differential equations (ODEs) for the time-dependent coefficients of the Laguerre expansion. Ben-Ameur et al. (2007) (BBKL) propose an alternative dynamic programming approach not based on PDEs in which a piecewise linear approximation is used for the value function and the exact transition probability is used to compute the expectation of the discounted piecewise linear approximation of the value function in the cases of CIR and Vasicek short rate models where the exact analytical expressions are available. All of these papers provide numerical experiments illustrating computational performance of their respective methods on the same example of a Swiss callable bond. This provides a natural comparison benchmark.
The present paper proposes an efficient method to evaluate callable and putable bonds under a wider class of interest rate models than any of the previous approaches, including the popular short rate diffusion models, as well as their time changed versions with jumps, including both jump-diffusion and pure jump models with state-dependent jumps. The method is based on the eigenfunction expansion of the pricing operator. We show that, under some technical conditions, the callable and putable bond pricing function has an eigenfunction expansion in eigenfunctions of the pricing operator with the expansion coefficients determined through a backward recursion. For popular short rate diffusion models, such as Cox-Ingersoll-Ross (CIR) (Cox et al., 1985), Vasicek (Vasicek, 1977), 3/2 (Ahn and Gao, 1999), we demonstrate on the test case of the Swiss bond used in the previous studies that the method is orders of magnitude faster than the alternative approaches in the literature. In addition, in contrast to the alternative approaches in the literature that have so far been limited to diffusions, the method is equally applicable to short rate jumpdiffusion and pure jump models constructed from diffusion models by Bochner's subordination with a Lévy subordinator.
The strength of the eigenfunction expansion method is that the value function is constructed globally in state and time with no need for discretization of either state or time variables. The only approximations involved in the computation scheme based on the method are the truncation of the infinite eigenfunction expansion (that, under some technical conditions, is uniformly convergent with uniformly controlled truncation error) and the numerical solution of a non-linear equation to determine the stopping boundary at each step of the backward recursion solved by the fast-converging bisection algorithm. Another strength of the method is that it can be seamlessly applied to both jump-diffusion and pure jump interest rate models obtained from diffusion models by subordination. Semi-analytical methods of BW and BBKL are not suitable to handling jump-diffusion and pure jump models since no analytical solutions are available for transition probabilities and Green's functions in these models. Numerical PDE methods, such as finitedifference and finite element methods, can, in principle, be extended to handle partial integro-differential equations (PIDE) arising in jump-diffusion and pure jump models, but at substantial costs both in the theoretical complexity and structure of the schemes and in their computational implementation and computational performance. Moreover, most of the existing instances of applying numerical PIDE methods in computational finance have been limited to Lévy processes with state-independent jumps. In contrast, the eigenfunction expansion method is capable of handling models with state-dependent jumps, such as meanreverting jumps in the interest rate. Both at the theoretical and computational level, moving from a pure diffusion model to a jump-diffusion or a pure jump model obtained from the diffusion model by the time change with respect to a Lévy subordinator amounts to no more than replacing the eigenvalues e −λnt of the pricing operator in the original diffusion model with the eigenvalues e −φ(λn)t of the pricing operator in the time changed model, where φ(λ) is the Laplace exponent of the Lévy subordinator. Remarkably, this insight goes back to the seminal work of Bochner (1949) introducing the idea of subordination (see page 370). It has been applied in probability theory (Albeverio and Rüdiger (2003); Chen and Song (2005)) and in finance (Albanese and Kuznetsov (2004); Boyarchenko and Levendorskiǐ (2006); Mendoza-Arriaga et al. (2010); Mendoza-Arriaga and Linetsky (2011); Li and Linetsky (2011)). The idea to use subordinated diffusions to build financial models goes back to Barndorff-Nielsen and Levendorskiǐ (2001), who considered NIG-like Feller processes for option pricing in the pseudo-differential operator framework.
Surveys on the application of eigenfunction expansions to the valuation of European-style derivatives can be found in Linetsky (2004Linetsky ( , 2008, where extensive bibliographies are given. Applications to interest rate models and the valuation of bonds without embedded options in particular can be found in Lewis (1998); Davydov and Linetsky (2003); Linetsky (2004, 2007); Boyarchenko and Levendorskiǐ (2006). Applications of eigenfunction expansions to European-style derivatives in models with jumps constructed by time changing diffusions with Lévy subordinators can be found in Albanese and Kuznetsov (2004), Boyarchenko and Levendorskiǐ (2006), Mendoza-Arriaga et al. (2010), Mendoza-Arriaga and Linetsky (2011), and Li and Linetsky (2011). The reference Boyarchenko and Levendorskiǐ (2006) is particularly relevant to our paper, as they also consider interest rate models based on subordinated diffusions.
The rest of the paper is organized as follows. Section 2.1 describes the general framework for the application of eigenfunction expansions to short rate diffusion models. Section 2.2 describes short rate models with jumps constructed by time-changing diffusion models with a Lévy subordinator. Section 3 presents our eigenfunction expansion method for solving the dynamic programming backward induction for callable and putable bonds. Section 4 presents examples of eigenfunction expansions for CIR, subordinate CIR, Vasicek, subordinate Vasicek, the 3/2 model, and the subordinate 3/2 model. Section 5 presents numerical experiments demonstrating computational performance of the method on the test case considered in the previous callable bond literature. Appendix contains selected proofs.

Short Rate Diffusion Models
Let {X t , t ≥ 0} be a one-dimensional, time-homogeneous regular (i.e. it reaches every point in (l, r) with positive probability) diffusion process on the interval I ⊆ R with (finite or infinite) endpoints l and r, −∞ ≤ l < r ≤ ∞. An endpoint is unattainable if it is a natural or an entrance boundary and is attainable if it is an exit or a regular boundary (see pp.14-15 of Borodin and Salminen (2002) for Feller's classification of boundaries for one-dimensional diffusions). In this paper we assume that the endpoints are either unattainable (and, thus, not included in the state space, so the interval I is open at such an endpoint) or regular and specified as instantaneously reflecting (the endpoint is included in the state space in that case, so the interval is closed at such an endpoint). We assume that the diffusion is conservative, that is P t (x, I) = 1 for each t ≥ 0 and x ∈ I, where P t (x, A) is the transition function from the initial state x to the Borel set A ⊆ I in time t. Thus the process X has infinite lifetime. We assume that the volatility σ(x) of X is positive and continuous on the open interval (l, r) and the drift µ(x) is continuous on (l, r).
We further assume that the instantaneous interest rate (the short rate) r t at time t is a function of the state variable X t and is given by r t = r(X t ). We assume that r(x) is continuous on (l, r). The continuity assumptions for σ, µ and r are not necessary, but simplify exposition in what follows. Consider the family of pricing operators (mathematically, Feynman-Kac (FK) operators) {P r t , t ≥ 0} defined by where E x is the expectation operator with respect to the probability measure P x of the process X starting at x ∈ I at time zero. Since in this paper we are interested in pricing, we always work with risk-neutral probabilities chosen by the market. The pricing operator P r t maps future payoff functions of the future state at time t into present value functions of the present state at time zero by discounting from time t to time zero and taking the expectation conditional on the current state at time zero. Under our assumptions, these operators form a strongly continuous semigroup on the Banach space C b (I) of bounded continuous functions on I with the supremum norm (see, e.g., Borodin and Salminen (2002)). The infinitesimal generator G r of this semigroup acts on twice continuously differentiable functions on I with compact supports by the second-order differential operator (the Sturm-Liouville operator): Define s(x) and m(x) to be the scale and speed densities of the diffusion process X: where x 0 ∈ (l, r) is an arbitrary point (see Karlin and Taylor (1981), Borodin and Salminen (2002) for discussions of the scale function and the speed measure of the one-dimensional diffusion). The infinitesimal generator can then be re-written in the formally self-adjoint form: Under our assumptions, the generator G r and the FK semigroup {P r t , t ≥ 0} can be extended to a self-adjoint operator and the symmetric strongly-continuous semigroup in the Hilbert space L 2 (I, m) of functions on I square-integrable with the speed measure m(dx) = m(x)dx and endowed with the inner product (f, g) = I f (x)g(x)m(x) dx and norm f = (f, f ). Thus, the Spectral Theorem for self-adjoint operators in Hilbert spaces can now be applied to write down the spectral decomposition of the generator and the semigroup. The spectral representation for one-dimensional diffusions goes back to the classical work of McKean (1956). We refer the reader to Linetsky (2004Linetsky ( , 2008 for surveys of applications in finance.
In this paper we further assume that the negative of the infinitesimal generator −G r has a purely discrete spectrum in L 2 (I, m) bounded from below. Sufficient conditions for the purely discrete spectrum in terms of the behavior of the functions σ, µ and r near the boundaries l and r are given in Linetsky (2004Linetsky ( , 2008. When the spectrum of −G r is purely discrete and bounded from below, the FK semigroup has the eigenfunction expansion of the form: for any f ∈ L 2 (I, m), where {λ n } ∞ n=0 such that λ 0 < λ 1 < · · · , lim n↑∞ λ n = ∞, are the eigenvalues of −G r and ϕ n are the corresponding eigenfunctions normalized so that ϕ n 2 = 1 (for future convenience we index the eigenvalues and eigenfunctions starting from n = 0 rather than n = 1). The eigenfunctions form a complete orthonormal basis in L 2 (I, m). We also assume that the eigenvalues satisfy the condition: for all t > 0, so that the FK semigroup is trace class (see Section 7.2 in Davies (2007)). We recall that the semigroup of a one-dimensional diffusion has a symmetric density p t (x, y) with respect to the speed measure m(dx) = m(x)dx that is a continuous function in t, x and y (see p.149 in Ito and McKean (1974) or p.13 of Borodin and Salminen (2002)). Hence, according to Theorem 7.2.5 in Davies (2007), the eigenfunctions ϕ n (x) are continuous functions with the global estimate |ϕ n (x)| ≤ e λnt/2 p t (x, x) for all t > 0, and the density p t (x, y) has an eigenfunction expansion for all t > 0 that converges uniformly in x and y on compacts. This ensures that, in addition to the L 2 convergence, the eigenfunction expansion (1) converges uniformly in x on compacts for all f ∈ L 2 (I, m) and t > 0. This follows from the Cauchy-Schwartz bound for the expansion coefficients |f n | ≤ f , the eigenfunction estimate, and the trace class condition (2). Since we are interested in bond pricing in this paper, we assume that the constant payoffs are in the Hilbert space L 2 (I, m), i.e. 1 ∈ L 2 (I, m). This is equivalent to assuming that the speed measure m is a finite measure on I, m(I) < ∞. In that case, the speed density can be normalized to one to be a probability density and, thus, serves as the steady state density of the underlying process X. Then the present value at time zero of a zero-coupon bond with unit face value and maturity t ≥ 0 when the underlying process is in state x, X 0 = x, has the eigenfunction expansion given by: with the expansion coefficients p n = (1, ϕ n ). Under our assumptions, the expansion converges uniformly in x on compacts for all t > 0. Virtually all popular short rate models in the financial economics literature fit into the general framework described above, including Cox-Ingersoll-Ross (CIR) model (Cox et al., 1985), Vasicek model (Vasicek, 1977), the 3/2 model (Ahn and Gao, 1999), Black's model of interest rates as options (Gorovoi and Linetsky, 2004), and the quadratic model (Beaglehole and Tenney, 1992;Leippold and Wu, 2002). Note that we have not made the assumption that r(x) is non-negative to accommodate the Vasicek model. If we make that assumption, then G r is positive semi-definite (so that −G r is negative semi-definite), the FK semigroup (P r t ) t≥0 is a contraction semigroup on L 2 (I, m), and λ 0 ≥ 0. To accommodate the Vasicek model, we made a weaker assumption that the spectrum of −G r is bounded from below, rather than non-negative.

Short Rate Models with Jumps Constructed by Subordination
A subordinator {T t , t ≥ 0} is a nondecreasing Lévy process with the Laplace transform given by the Lévy-Kchintchine formula with the Laplace exponent φ(λ), with nonnegative drift γ ≥ 0, and Lévy measure ν(ds) satisfying the integrability condition ∞ 0 (s ∧ 1)ν(ds) < ∞. For any set A ⊂ R bounded away from zero, jumps of sizes in A arrive according to a Poisson process with the arrival rate ν(A). If ν is a finite measure on (0, ∞), the subordinator is a compound Poisson process plus drift at the rate γ. If ν(0, ∞) = ∞, the subordinator is a jump process with infinite activity and drift γ. If γ = 0, it is a pure jump process. Examples of subordinators important in applications include compound Poisson processes with exponential or gamma distributed jump sizes, inverse Gaussian (IG) subordinators (Barndorff-Nielsen, 1998), and gamma subordinators (Madan et al., 1998), with the latter two having infinite activity. These examples are special cases of subordinators with Lévy measures of the form ν(ds) = Cs −p−1 e −ηs ds with C > 0, η > 0, and p < 1. The case with p ∈ (0, 1) are the tempered stable subordinators (the limiting cases with η = 0 are stable subordinators). The special case with p = 1/2 is the IG subordinator. The limiting case with p = 0 is the Gamma subordinator. The Laplace exponent is given by: As an example, the Lévy measure and Laplace exponent for an IG subordinator parameterized with µ and ν, the mean and variance of an IG process at time one, t = 1, are given by: Further mathematical details on subordinators can be found in Schilling et al. (2010) and on financial applications in Cont and Tankov (2004).
Since subordinators are non-negative, non-decreasing processes, they can be used as stochastic time changes to time change other processes. This procedure is known as Bochner's subordination and goes back to Bochner (1949Bochner ( , 1955. In particular, time changing a Markov process with a subordinator yields another Markov process whose semigroup and infinitesimal generator are given by Phillips' theorem (Phillips (1952) (2011), and Li and Linetsky (2011).
In particular, we can construct new short rate models with jumps from diffusion short rate models as follows. Let X be an ergodic one-dimensional diffusion on I with volatility σ(x) and drift µ(x) and with the stationary density given by the (normalized) speed density m as described in the previous section. Let r(x) be the function defining the diffusion short rate model as in section 2.1 and {P r t , t ≥ 0} the corresponding FK semigroup. Then a new short rate model with jumps is obtained by subordinating the original diffusion short rate model with respect to a given subordinator T with the Laplace exponent φ(λ) by defining a new semigroup (the superscript φ signifies that the subordination is performed with respect to the subordinator with the Laplace exponent φ): where Y t is a new Markov jump-diffusion process on I with the infinitesimal generator G φ acting on twicedifferentiable functions with compact supports as an integro-differential operator (see Mendoza-Arriaga et al. (2010) and Mendoza-Arriaga and Linetsky (2011)): where the drift with respect to the truncation function 1 {|y−x|≤1} is and π φ (x, y) is the symmetric state-dependent Lévy density: where p s (x, y) is the density of the original FK semigroup (P r t ) t≥0 of the pure diffusion short rate model, and γ and ν are the drift and the Lévy measure of the subordinator. When γ > 0, Y is a jump-diffusion. When γ = 0, Y is a pure jump process. Y has the same steady state density m(x) as the original diffusion X. The short rate in (7) is the function of the state variable Y t so that r t = r φ (Y t ) with the function r φ (x) given by: where P (s, x) is the price (4) of the s-maturity zero-coupon bond at time zero when the state of the underlying diffusion X 0 = x. The form of the generator G φ of Y and the function r φ follow from Phillips' theorem that characterizes the subordinate semigroup and its infinitesimal generator. Here we subordinate the pricing (FK) semigroup of the diffusion process X with the discount rate r t = r(X t ) with respect to a given subordinator to obtain a new semigroup interpreted as the pricing (FK) semigroup of a new short rate model r t = r φ (Y t ) driven by the state variable Y t following a Markov process with jumps.
Remark 2.1. Mathematically, the subordination of the FK semigroup can be interpreted as follows. First formulate the FK semigroup P r of the original conservative diffusion X with the discount rate r(x) as the transition semigroup of the diffusionX with killing at the rate r(x) (cf. Section II.4 on pp.27-28 in Borodin and Salminen (2002) for the connection between discounting and killing). Then construct a new processX φ t :=X Tt by time changingX with the subordinator T . Use Phillips' theorem to write down its infinitesimal generator and, thus, its local characteristics (diffusion, drift, state-dependent Lévy measure, and state-dependent killing rate r φ (x)). Then formulate the transition semigroup of the processX φ as the FK semigroup of a conservative process Y with the generator given above and with the discount rate r φ (x) given above. The formulation of the application of Phillips' theorem to this situation is given in Mendoza-Arriaga et al. (2010) and we do not repeat it here to save space.
The subordinate FK semigroup {P r,φ t , t ≥ 0} is also a symmetric strongly continuous semigroup of operators on L 2 (I, m) (Chen, 2005), and, under the assumptions we have made about the semigroup {P r t , t ≥ 0} in the previous section, it possesses an eigenfunction expansion in the same eigenfunctions ϕ n (x) with λ n in Eq.(1) replaced with λ φ n := φ(λ n ), where φ(λ) is the Laplace exponent of the subordinator: for all t > 0, then the semigroup {P r,φ t , t ≥ 0} is also trace class. If we further assume that the eigenfunctions of the original pure diffusion FK semigroup have a bound independent of n on each compact interval K = [a, b] ⊂ (l, r), i.e, |ϕ n (x)| ≤ C K for all n, where the constants C K may depend on the interval K but are independent of n, then these assumptions ensures that, in addition to the L 2 convergence, the eigenfunction expansion of the subordinate semigroup (8) converges uniformly in x on compacts for all f ∈ L 2 (I, m) and t > 0 and that the subordinate semigroup {P r,φ t , t ≥ 0} also has a continuous density with respect to m(x)dx with the eigenfunction expansion for all t > 0 uniformly convergent on compacts in x and y. For the present value at time zero of a zero-coupon bond with unit face value and maturity t ≥ 0 when the underlying state variable has initial value of Y 0 = x, we then obtain the eigenfunction expansion given by: with the expansion coefficients p n = (1, ϕ n ). Under the assumption (9) on the growth φ(λ) and the bound on eigenfunctions independent of n, the expansion converges uniformly in x on compacts for all t > 0.
Using the subordination approach, we can extend all the diffusion short rate models popular in financial economics to jump-diffusion and pure jump models, in particular constructing subordinate CIR (SubCIR), subordinate Vasicek (SubVasicek), etc. Subordinate models allow for jumps in the interest rate dynamics. Moreover, if the diffusion process is mean-reverting, the subordinate process will have jumps that are also mean-reverting (see Li and Linetsky (2011) for the proof in the subordinate Ornstein-Uhlenbeck context). While adding jumps improves the model's realism and flexibility, remarkably, the analytical and computational framework remains entirely unchanged, as the only modification required in the eigenfunction expansion is the replacement of λ n in Eq.(1) with λ φ n = φ(λ n ) in Eq. (8). Remark 2.2. Matching the initial yield curve. The time-homogeneous short rate models discussed in sections 2.1 and 2.2 can be extended to match any initial term structure of interest rates as proposed by Brigo and Mercurio (2001) and commonly done in fixed income market practice by adding a deterministic function of time to the short rate process in the extended diffusion and the subordinated diffusion models, respectively. The function can then be chosen so that the initial zero-coupon bond prices of all maturities in the extended model match the zero-coupon prices consistent with the given initial term structure of interest rates (given yield curve). The callable and putable bond pricing developed in this paper can then be immediately extended to this class of models. To simplify notation, we do not explicitly consider this extension in what follows and assume the short rate model is time homogeneous.

The Eigenfunction Expansion Method for Callable and Putable Bonds
The call option allows the bond issuer to buy the bond back from the bond holder (call the bond) for prespecified call prices at some pre-specified times prior to maturity. This allows the bond issuer to refinance the bond if interest rates decline. The put option allows the bond holder to sell (put) the bond back to the bond issuer for pre-specified put prices at some pre-specified times prior to maturity. This allows the bond holder to re-invest the proceeds into a bond with higher coupon if interest rates rise. We assume that the bond principal is equal to one dollar and the bond pays coupons of C dollars on the coupon dates t i , i = 1, ..., k. Let the bond issue date and maturity date be t 0 = 0 and t k = T , respectively. After some initial protection period from t 0 until t k * , the call and put options can be exercised at the subsequent coupon dates t i , i = k * , ..., k − 1, prior to maturity t k . We also assume that there are notice periods of lengths δ so that the option exercise decision to exercise at the coupon date t i has to be made at an earlier time τ i = t i − δ so that the adequate advance notice of duration δ can be given to the other party. Denote the call and put prices at times t i as K c i and K p i , respectively. It is assumed that K c i > K p i . This sets up an optimal stopping game in discrete time with finite horizon where one player (the bond holder) chooses a stopping time to maximize the bond value and the other player (the bond issuer) chooses a stopping time to minimize the bond value.
We assume that the short rate is where Y t is a subordinate diffusion, as discussed in section 2. We note that the pure diffusion model can be viewed a special case with the trivial time change T t = t with φ(λ) = λ, γ = 1 and ν ≡ 0. We assume that all the assumptions made in section 2 are in force. Then the pricing operator P r,φ t has the eigenfunction expansion (11). To simplify notation, we drop the superscripts r and φ and simply write P t for the pricing operator P r,φ t in what follows. Let V (t, x) be the value of the bond at time t ∈ [0, T ] when the underlying state is x. Since the decisions to exercise the call and put options are made at times τ i = t i − δ prior to the coupon dates, the present values at τ i of the call and put prices at time t i must be compared to the holding (continuation) values of the bond at time τ i . The discounted value of the call price at time τ i is K c i P (δ, x), where P (δ, x) is the value of the zero-coupon bond with time to maturity equal to the notice period δ and unit face value when the state variable is in the state x at time τ i . The expected discounted value of the put price at time τ i is K p i P (δ, x). At maturity the bond value is equal to its principal plus the last coupon, denote the bond's value at time τ i and C i (x) denote the bond's holding (continuation) value at time τ i ex-coupon at time t i , i = k * , ..., k − 1 (the assumption is that the next coupon C is always paid at time t i , whether or not decisions are made to exercise any of the two options at time τ i ). Let V 0 (x) := V (0, x) denote the value of the bond at the time of issue.
Assuming both players behave rationally and maximize the value of their assets and minimize the value of their liabilities, in state x the bond issuer will exercise the call option at time , the bond holder will exercise the put option at time , and there will be no exercise of either option at time . By the assumption that K c i > K p i the simultaneous exercise of both options is never optimal. Further set Then the value of the bond with call and put options satisfies the following Bellman's dynamic programming backward induction (see also Ben-Ameur et al. (2007)): Assuming that for each i = k * , ..., k − 1 each of the two equations has at most one solution in I denoted by x c i and x p i , respectively, setting x c i := l if the first equation has no solution in I, setting x p i := r if the second equation has no solution in I, and noting that x c i < x p i , the backward induction (14) can be re-written in the form: for all k * ≤ i ≤ k − 1 and x ∈ I, where 1 {·} denotes the indicator function, and 1 {x<x c i } ≡ 0 if x c i = l and 1 {x>x p i } ≡ 0 if x p i = r by convention. This backward induction can be solved by a variety of computational methods in the literature, as discussed in the introduction and in section 5, based on the different methods to approximate the pricing operator P t appearing in Eqs. (13) and (15). In this paper we follow the approach based on representing the pricing operator P t by its eigenfunction expansion. Our main result is the following theorem that summarizes our eigenfunction expansion method for the valuation of callable and putable bonds.
Theorem 3.1. Suppose that m is a finite measure on I, m(I) < ∞.
(i) The value function V 0 (x) and the value functions V i (x) and the continuation value functions C i (x) are in L 2 (I, m) for all i = k * , ..., k − 1.
(ii) The continuation value functions have the following eigenfunction expansions: The value function at the time of the bond issue has the following eigenfunction expansion: (iii) The eigenfunction expansions (18) and (19) converge uniformly in x on compacts in I.
(iv) If each of the two equations in (16) has at most one solution in I, then the eigenfunction expansion coefficients in (18) and (19) satisfy the following backward recursion: c k n = (1 + C)p n , n = 0, 1, 2, ..., where p n = (1, ϕ n ) are the expansion coefficients of the unit payoff appearing in the eigenfunction expansion of the zero-coupon bond (11), and for each i = k * , ..., k − 1, for n = 0, 1, 2, ..., where x c i and x p i are as previously defined, and for l ≤ x < y ≤ r we introduced the following notation where 1 (x,y) = 1 (x,y) (z) is the indicator function of the interval (x, y) and P (t) = P (t, x) is the value function of the zero-coupon bond with time to maturity t ≥ 0 and unit face value when the underlying state variable is x.
Proof. See the Appendix.
Theorem 3.1 reduces the solution of the backward induction for callable and putable bonds to the recursion for the expansion coefficients (20), together with finding the roots of equations (16) by a numerical root finding algorithm, such as bisection. The continuation value C i (x) on the right hand side of the equations (16) is given by the eigenfunction expansion (18) with the coefficients determined on the previous step of the recursion. The eigenfunction expansion can be truncated at a finite level, and the truncation error is uniformly controlled due to the uniform convergence of the expansions when condition (9) is satisfied.
Remark 3.1. When the state process is an affine diffusion, such as CIR or Vasicek, the zero-coupon bond value function has the exponential-affine form in the state variable and the integral in Eq. (22) can be written as and in some cases can be calculated in closed form. Generally, it can be calculated using the eigenfunction expansion for the zero-coupon bond value function: where π m,n (x, y) are defined in Eq.(21).
Remark 3.2. If the bond has only the call option and no put option, then the game reduces to the optimal stopping problem for the bond issuer. In that case in Eq.(20) the term with K p i p n (x i p , r) is absent and π m,n (x c i , x p i ) is replaced with π m,n (x c i , r) for all i. Similarly, if the bond has only the put option and no call option, the game reduces to the optimal stopping problem for the bond holder and in Eq.(20) the term with K c i p n (l, x i c ) is absent and π m,n (x c i , x p i ) is replaced with π m,n (l, x p i ) for all i.
The condition that each of the two equations in (16) has at most one solution in I generally needs to be checked case by case. For CIR, Vasicek, and 3/2 models considered in section 4 the condition can be verified by the following proposition. For CIR and Vasicek, (26) below follows from Theorem 1.1 in Ikeda and Watanabe (1977) while (26) holds for the 3/2 model since the diffusion process X(t) in the 3/2 model can be written as X(t) = 1/Y (t), where Y (t) is a CIR process satisfying the Feller condition.
Proposition 3.2. Suppose that the discount rate r(x) is a non-decreasing function. Let (Ω, F , P) be a complete probability space with right continuous increasing family (F t ) t≥0 of sub σ-fields of F each containing P -null sets and let B t be a one-dimensional F t -Brownian motion. Let σ(x) and µ(x) be continuous. Let X 1 (t) and X 2 (t) be processes started at different initial states such that X 1 (0) < X 2 (0) and then each of the two equations in (16) has at most one solution in I.
Proof. See the Appendix.
Remark 3.3. We emphasize that uniqueness of roots of (16) is not a requirement for our method to work. In fact, it is one of the strengths of our approach that it can handle just as easily more general cases with multiple break-even points and, hence, early exercise regions that are not necessarily one-sided and, in general, can be unions of multiple intervals. Proving the uniqueness of the break-even point provides a convenience for numerical implementation, as we can stop after finding a single root. In general, without the proof of uniqueness, a more thorough numerical investigation of the functions is required in each case to either establish uniqueness or determine multiple roots. While we have been able to prove in Proposition 3.2 uniqueness for pure diffusion short rate models, we have been unable to extend the proof to the case of subordinate diffusions, as it is based on classical SDE comparison results that to the best of our knowledge are not available in general for subordinate diffusions. Nevertheless, in our extensive numerical experimentation in all cases of subordinate diffusions we have considered, we have observed similar behavior of functions in (16) that lead to unique solutions. We thus conjecture that uniqueness also holds for subordinate diffusions, perhaps subject to some condition on the subordinator.

CIR and SubCIR
In the Cox-Ingersoll-Ross (CIR) model (Cox et al., 1985), the short rate follows the CIR diffusion with drift µ(x) = κ(θ − x) and volatility σ(x) = σ √ x where κ > 0, θ > 0, and σ > 0 are the rate of mean reversion, the long run mean and volatility, respectively. In this case r(x) = x. When Feller's condition 2κθ/σ 2 ≥ 1 is satisfied, the origin is an unattainable entrance boundary and infinity is an unattainable natural boundary. In this case I = (0, ∞). When Feller's condition is not satisfied, the origin is an attainable regular boundary and is specified as instantaneously reflecting. In this case I = [0, ∞). The CIR speed density reads m(x) = 2 σ 2 x b−1 e − 2κx σ 2 . The celebrated CIR zero-coupon bond pricing formula is: where A(t) = 2γe (κ+γ)t/2 (γ + κ)(e γt − 1) + 2γ b , B(t) = 2(e γt − 1) (γ + κ)(e γt − 1) + 2γ , The eigenfunction expansion (4) of the CIR zero-coupon bond pricing function is given in Davydov and Linetsky (2003). In this case the eigenfunctions, eigenvalues, and the expansion coefficients for the unit payoff are: n (x) are the generalized Laguerre polynomials and b and γ as defined in (28). The CIR eigenfunctions are continuous and have a bound independent of n on each compact interval K = [a, b] ⊂ (l, r), i.e. |ϕ n (x)| ≤ C K for all n, where the constant C K is independent of n, since by inequality (27a) on p.53 of Nikiforov and Uvarov (1988), for any compact interval in I, the CIR eigenfunctions satisfy the bound |ϕ n (x)| ≤ Cn −1/4 , where the constant C is independent of n (but depends on the interval).
The quantities (21) and (22) in the CIR model can be calculated as follows: where we introduced the following notation: y α e −sy L (α) n (y) dy.
In the calculation of (30) we used the explicit expression for the CIR zero-coupon bond pricing function (27) as in Eq.(24), rather than its eigenfunction expansion. The quantities a For n ≥ 1, For m = n, a (α) where γ(α + 1, x) = x 0 e −y y α dy is the lower incomplete gamma function. The coefficients b Proof. See the Appendix.
The Laguerre polynomials of degree α satisfy the following classical recursion (Lebedev, 1965, Eq. 4.18.1) 1 (x) = −x + α + 1. Then for any N the quantities {a (α) n,n (x), 0 ≤ n ≤ N } can be efficiently computed recursively in the following order: The computation of a (α) n,m (x) with n = m can be done directly using (31) and (32)   n (x), 0 ≤ n ≤ N } can be efficiently computed recursively in the order: N (x) For the SubCIR model the explicit bond pricing formula similar to (27) is not available, and we use the eigenfunction expansion (11) of the SubCIR zero-coupon bond price as in Eq.(25). The expression (30) is then replaced with: In Theorem 3.1 λ n are now the eigenvalues of the SubCIR model related by φ(λ n ) to the eigenvalues of the CIR diffusion model.

Vasicek and SubVasicek
In Vasicek model (Vasicek, 1977), the short rate follows the OU diffusion with drift µ(x) = κ(θ − x) with κ > 0 and θ > 0 and constant volatility σ > 0. In this case I = R and r(x) = x, both boundaries at plus and minus infinity are unattainable natural boundaries, and the process can get negative. However, when θ and the initial state x 0 are sufficiently above zero and κ > 0 is sufficiently large, the probability of the rate falling below zero is relatively small due to mean reversion pulling the process back towards the positive long run mean as it approaches zero from above. The Vasicek speed density is a Gaussian density . The celebrated Vasicek zero-coupon bond pricing formula has the same exponential affine form as the CIR (27) with The eigenfunction expansion (4) of the Vasicek zero-coupon bond pricing function is given in Gorovoi and Linetsky (2004). In this case the eigenfunctions, eigenvalues, and the expansion coefficients for the unit payoff are: where H n (x) are Hermite polynomials. The Vasicek eigenfunctions are continuous and have a bound independent of n on each compact interval K = [a, b] ⊂ (l, r), i.e. |ϕ n (x)| ≤ C K for all n, where the constant C K is independent of n, since by inequality (28a) on p.53 of Nikiforov and Uvarov (1988), for any compact interval in I, the Vasicek eigenfunctions satisfy the bound |ϕ n (x)| ≤ Cn −1/4 , where the constant C is independent of n.
The quantities (21) and (22) in the Vasicek model can be calculated as follows: where we introduced the following notation: In the calculation of (37) we used the explicit expression for the Vasicek zero-coupon bond pricing function as in (24), rather than its eigenfunction expansion. The quantities a n,m and b n can be computed efficiently. The coefficients of the Hermite polynomial H n (x) can be computed from the recursive equation (Lebedev, 1965, Eq. 4.10.1) Proposition 4.2. The quantities a m,n (x) can be computed as follows. For m = n a n,m ( (38) a n,n (x) can be computed recursively as follows: where Φ(x) is the standard normal cumulative distribution function. The quantities b n (s, x) can be computed as follows.
where Erf(x) is the error function, and Proof. See the Appendix.
For the SubVasicek model the explicit bond pricing formula similar to (23) is not available, and we use the eigenfunction expansion (11) of the SubVasicek zero-coupon bond price instead. The expression (37) is then replaced with: In the recursion (20) λ n are the eigenvalues of the SubVasicek model related by φ(λ n ) to the eigenvalues of the Vasicek diffusion model.  (Gradshteyn and Ryzhik, 2007, Eq. 7.374.6, p.803 These coefficients can be used to evaluate callable bonds without the put option, similar to remark 4.1 for the CIR.

The 3/2 and Sub-3/2 Model
In this model, the short rate process is a diffusion on (0, ∞) with infinitesimal parameters σ( where κ, θ, and σ are positive constant parameters. This process was proposed by Cox et al. (1985) as a model for the inflation rate. Ahn and Gao (1999) propose this process as a model for the short rate and show that this model is empirically more plausible than the square-root model.
The eigenfunction expansion (4) of the zero-coupon bond pricing bond is given in Linetsky (2004). In this case the eigenfunctions, eigenvalues, and the expansion coefficients for the unit payoff are: The 3/2 eigenfunctions are continuous on I and have a bound independent of n on each compact interval K = [a, b] ⊂ (l, r), i.e. |ϕ n (x)| ≤ C K for all n, where the constant C K is independent of n, since by inequality (27a) on p.53 of Nikiforov and Uvarov (1988), for any compact interval in I, the 3/2 eigenfunctions satisfy the bound |ϕ n (x)| ≤ Cn −1/4 , where the constant C is independent of n.
The quantities (21) and (22) in the 3/2 and Sub-3/2 model can be calculated as follows: where we introduced the following notation: For the Sub-3/2 model, λ n are the eigenvalues of the Sub-3/2 model related by φ(λ n ) to the eigenvalues of the 3/2 diffusion model. The expression for a n,k (x) is same as the one for the CIR model. Hence, the quantities a n,k (x) can be computed using the method given in the CIR section.

Computational Results
This section shows our computational results for CIR, Vasicek, SubCIR, and SubVasicek models. We consider the callable bond example that has been extensively used in the literature starting from Büttler and Waldvogel (1996) and including d' Halluin et al. (2001), Ben-Ameur et al. (2007), andde Frutos (2008), as the test case to compare computational performance of a number of computational approaches to the callable bond valuation. The callable bond was issued by Swiss Confederation in 1987 with maturity in 2012. At the time of valuation considered in Büttler and Waldvogel (1996) and the subsequent papers, the remaining time to maturity of the bond was t k = 20.172 years with k = 21 remaining annual coupons of 4.25% per annum. The notice period is 2 months, δ = 0.1666. The protection period is t k * = 10.172 with k * = 11. There are ten early exercise dates, t 11 = 10.172, t 12 = 11.172,..., t 20 = 19.172. The call prices corresponding to these dates are given in Table 1. The bond did not include a put option. We use the values of the parameters κ, σ, θ for the CIR and Vasicek models estimated in Büttler and Waldvogel (1996) and used in the subsequent papers in the literature. They are given in Table 2. For the SubCIR and SubVasicek models, we used the same parameter values for the underlying CIR and Vasicek diffusions, while specifying the subordinator to be the inverse Gaussian (IG) subordinator with drift (IG Lévy measure given in (6)). The IG parameter ν was set to 1. The subordinator drift γ and the IG parameter µ were chosen so that E[T t ] = t to normalize the time change. For the jump-diffusion case, we used γ = 0.5 and µ = 0.5. For the pure jump case, we used γ = 0 and µ = 1.
In the process of finding the break-even point x c i at each step of the recursion, the infinite series for the continuation value given in (18) needs to be truncated at some finite level. At time t 0 , the series in (19) also needs to be truncated. In the recursion formula for the expansion coefficients in (20), only the previous expansion coefficients that were calculated at the previous step are used in the eigenfunction expansion. For the subordinated models, we also compute the zero-coupon bonds by the eigenfunction expansions that are also truncated at some finite level. We used an adaptive truncation strategy that truncated the expansion after a user-specified relative error tolerance ǫ was reached in each instance of the series evaluation by comparing with ǫ the ratio of the next term and the sum of the next two terms relative to the sum of all the previous terms.
In order to find the break-even point at each decision point in time, the bisection method was used. The break-even short rates are shown in table 3 (JD stands for jump-diffusion, and PJ for pure jump). For Vasicek and CIR models, r(x) = x, so the short rate is equal to the state value. For subordinated models, the short rate is given by the function r φ (x) of the state variable. For the CIR model, we start by checking the boundary at zero to see if the expected discounted value of the strike is greater than the continuation value, in which case there is no non-negative break-even point at that decision time instance. Otherwise, there is a unique non-negative break-even point for the CIR model. For the Vasicek model, there always is a unique break-even point. The break-even point was found by the bisection method until the length of the search interval became less than 10 −7 . Table 4 shows computational results for all the models considered in this section with the initial short rate r = 0.05. The first column indicates the absolute pricing error in pricing the callable bond. During the process of finding the break-even point, the truncation level is determined for each evaluation of the continuation value. The second column gives the average truncation level N at each decision point at time τ 20 , τ 19 , ..., τ 11 , and at t 0 (the average of truncations levels as determined by our adaptive truncation algorithm in evaluating the expansion of the continuation value needed for each step of the bisection algorithm). The third column shows the maximum truncation level N at those times. The fourth column shows the CPU time of our algorithm implemented in C using the GNU Scientific Library (GSL) and compiled with gcc and executed on a 2.4 GHz Intel Core i3 370M processor. The CPU time includes the time taken for any precomputations required. For CIR and Vasicek models, the CPU time to price the callable bond to approximately five correct decimal points (convergence of 10 −5 in the tables) was about one millisecond. For comparison, de Frutos (2008) reported CPU times of 0.75 seconds using Matlab on a 3GHz processor, while Ben-Ameur et al. (2007) reported CPU times of 2 to 3 seconds using C on a 2.0 GHz Pentium 4 processor. Thus, the eigenfunction expansion approach is approximately three orders of magnitude faster in this instance of pricing the callable bond.
For the subordinated models with jumps, jump-diffusion models required slightly longer CPU times than pure diffusion models, while pure jump models required slightly longer times than jump-diffusions. This is due to the replacement of the diffusion eigenvalues λ n with eigenvalues φ(λ n ) of subordinated processes that slows down the eigenvalue growth and, hence, requires more terms in the expansions, as evidenced in Table 4. Still, the algorithm reached the pricing error of under 10 −5 in 2.2 and 2.5 milliseconds under pure jump CIR and Vasicek models. Tables 5 and 6 show the computed values of the callable bond in comparison to other methods. The columns BW, DFVL, BBKL, and F refer to Büttler and Waldvogel (1996), d 'Halluin et al. (2001), Ben-Ameur et al. (2007), andde Frutos (2008). Table 7 shows the results for the subordinated models in the present paper. We stress that neither of the alternative approaches in the literature is capable of handling jump-diffusion and pure jump models with state dependent jumps. The remarkable advantage of the eigenfunction expansion method is that it is entirely straightforward to move from pure diffusion models to jump-diffusion and pure jump models obtained by subordination by simply replacing the diffusion eigenvalues λ n with subordinate eigenvalues φ(λ n ).
While Tables 1-7 provide results for the bond with the call option only to facilitate comparisons with the literature, Tables 8-10 provide the corresponding results for the bond that is both callable and putable. To generate this example of a bond with both options, we added the put option to the callable bond considered previously in this section. Table 8 gives the put prices we have assumed. Table 9 presents results for breakeven short rates for call and put options under the range of short rate models considered in this paper. Table  10 presents the corresponding prices of the bond with both call and put options.

Conclusion
This paper proposed an efficient method to evaluate bonds with embedded options under a wide class of interest rate models, including the popular short rate diffusion models, as well as their time changed versions with jumps. The method is based on the eigenfunction expansion of the pricing operator. Given the set of call and put dates, the callable and putable bond pricing function is the value function of a stochastic game with stopping times. Under some technical conditions, it is shown to have an eigenfunction expansion in eigenfunctions of the pricing operator with the expansion coefficients determined through a backward recursion. For CIR and Vasicek the method is orders of magnitude faster than the alternative approaches in the literature. In contrast to the alternative approaches in the literature that have so far been limited to diffusions, the method is equally applicable to short rate jump-diffusion and pure jump models constructed from diffusion models by Bochner's subordination with a Lévy subordinator. In future work we plan to apply the eigenfunction expansion method of this paper to convertible bonds, where the stock price process is the stochastic variable driving the conversion and call decisions.

Appendix
Proof of Theorem 3.1. (i) For any f ∈ L 2 (I, m), P t f ∈ L 2 (I, m) for any t ≥ 0. Since m is a finite measure on I, 1 ∈ L 2 (I, m), so P (t, x) ∈ L 2 (I, m) for any t ≥ 0. Then by (12) to (14), C k−1 (x) and V k−1 (x) are in L 2 (I, m). By induction using (13), (14), and (15), it can be shown that the value function V 0 (x) and the value functions V i (x) and the continuation value function C i (x) are in L 2 (I, m) for all i = k * , ..., k − 1.
(ii) The expressions for the continuation value function and the value function at the time of the bond issue are given in (13) and (15). By part (i), the eigenfunction expansion can be obtained from (1) or (8), where c i+1 n = (V i+1 , ϕ n ), i = k * , ..., k − 1, and c k * n = (V k * , ϕ n ) for n = 0, 1, 2, .... (iii) Under the conditions given in section 2.1 or 2.2, the eigenfunction expansion for the density given in (3) or (10) holds. Then for f ∈ L 2 (I, m) and any x ∈ I, The first inequality follows from the Cauchy-Schwartz inequality, and the last inequality follows from the bounds on eigenfunctions described in section 2.1 or 2.2 (e.g. |ϕ n (x)| ≤ e λnt/2 p t (x, x) or |ϕ n (x)| ≤ C K ) and the trace class condition (2) or (9). Hence, the eigenfunction expansion converges pointwise to P t f (x) for any x ∈ I. The eigenfunction expansions converge uniformly in x on compacts in I by the following: Let K be any compact subset of I.
The last expression goes to 0 as M goes to ∞ by the bounds on |ϕ n (x)| and the trace class condition. Therefore, the eigenfunction expansions converge uniformly in x on compacts in I.
(iv) By (12), the values of c k n are given by the coefficients of the eigenfunction expansion of the zero-coupon bond: c k n = (1 + C)p n , n = 0, 1, 2, ... Suppose that we know the values of the coefficients c i+1 n , n = 0, 1, 2, .... Denote 1 (x,y) as an indicator function that is 1 on the interval (x, y) and 0 otherwise.
The last equality follows from (f, g) = ∞ n=0 f n g n , where f n = (f, ϕ n ) and g n = (g, ϕ n ) for f, g ∈ L 2 (I, m).
Proof of Proposition 3.2. Suppose that (26) holds. From this we get that, with probability one, Hence, the zero-coupon bond function P (t, x) is a positive decreasing function. We can also show that the continuation value function C i (x) and bond value function V i (x) are positive decreasing functions of x for all i = k * , ..., k − 1. For i = k − 1, C k−1 (x) is a constant multiple of a zerocoupon bond function by (12) and (13), so C k−1 (x) is a positive decreasing function. By (14), the bond value function V i (x) is a positive decreasing function if C i (x) is a positive decreasing function since we already showed that zero-coupon bond functions are positive decreasing functions. It remains to show that for i = k * , ..., k − 2 the continuation function C i (x) is a positive decreasing function given that V i+1 (x) is a positive decreasing function. This is shown by the following. By (26), for X 1 (0) < X 2 (0). We show that there can be at most one solution x such that where K is a constant (either K c i or K p i ), for each i = k * , ..., k − 1. It is first shown that there is at most one solution to the equation given by Suppose that x k−1 is a solution to the equation. By (13), where g(x) = E e − t k t k−1 r(X(s)) ds X(t k−1 ) = x is a decreasing function and E t k−1 denotes the expectation under the t k−1 -forward adjusted measure. Denote the probability measure under the t k−1 -forward adjusted measure as Q k−1 . Q k−1 and P are equivalent probability measures, so P[X 1 (t) ≤ X 2 (t) for all t ≥ 0] = 1 implies that Q k−1 [X 1 (t) ≤ X 2 (t) for all t ≥ 0] = 1. Then E t k−1 [g(X(t k−1 ))|X(τ k−1 ) = x] is a decreasing function of x. At x = x k−1 , It is shown next that there is at most one solution to the equation given by (13), x is a decreasing function and E ti denotes the expectation under the t i -forward adjusted measure. Denote the probability measure under the t i -forward adjusted measure as Q i . Q i and P are equivalent probability measures, so Proof of Proposition 4.1. For Laguerre polynomials, the forward shift property is d dx n−1 (x) and the backward shift property is To derive the recursion for a (α) n,m (x), we first apply the backward shift, integrate by parts and then apply the forward shift.
For n ≥ 1, m ≥ 1, m = n, by solving the above equation for a α+1 n,m (x) and equating the expressions for a α+1 n,m (x) and a α+1 m,n (x), we have For n ≥ 1, using the backward shift property, we have For n = m, we get the following from (45): where γ(α + 1, x) = x 0 e −y y α dy is the lower incomplete gamma function. The coefficients b (α) n (x) are computed from equations (46) and (47) below. For n ≥ 1, Proof of Proposition 4.2. For Hermite polynomials, the forward shift property is d dx H n (x) = 2nH n−1 (x) and the backward shift property is To derive the recursion for a n,m (x), we first apply the backward shift, integrate by parts and then apply the forward shift.
If m = n, then subtracting (48) from (49) and rearranging, we have a n,m ( a n,n (x) can be computed recursively as follows: where Φ(x) is the cumulative distribution function of a standard normal distribution.