An Asymptotic Solution for Call Options on Zero-Coupon Bonds

: We present an asymptotic solution for call options on zero-coupon bonds, assuming a stochastic process for the price of the bond, rather than for interest rates in general. The stochastic process for the bond price incorporates dampening of the price return volatility based on the maturity of the bond. We derive the PDE in a similar way to Black and Scholes. Using a perturbation approach, we derive an asymptotic solution for the value of a call option. The result is interesting, as the leading order terms are equivalent to the Black–Scholes model and the additional next order terms provide an adjustment to Black–Scholes that results from the stochastic process for the price of the bond. In addition, based on the asymptotic solution, we derive delta, gamma, vega and theta solutions. We present some comparison values for the solution and the Greeks.


Introduction
Black and Scholes' [1] solution for European style options assumes a geometric Brownian motion process for the underlying stock price. Many additional asset classes were also assumed to follow this process, and so the Black-Scholes model was used to price options on a variety of asset classes. Since many asset classes have no maturity (stocks, commodity prices, foreign exchange rates), this process seems reasonable. However, it was clear from the outset that this process would not work for most bond prices. To price bonds, and other interest rate securities, researchers turned to modeling interest rates and using the resulting rates to price interest rate securities and their related derivative securities.
In this paper, we propose a process for a zero-coupon bond price directly, and using that process, we derive a partial differential equation (PDE) and an asymptotic solution for a call option on the bond. Our approach of modeling bond price dynamics rather than interest rate dynamics offers several important advantages. First, we are not restricted to assumptions regarding the relationship of interest rates in a term structure, that is, we do not need to assume rates fit a model, such as "Preferred Habitat." Second, modeling the bond dynamics is computationally simpler than many spot and term structure models with multiple factors. In order for spot and term structure models to provide for reasonable representation of interest rate dynamics, the models often require multiple uncertainty terms with non-constant volatility. This results in computationally complex formulations, which normally require numerical solutions. Finally, the resulting solution to our formulation is similar to the familiar Black-Scholes model for European style call options on stocks. While there are advantages and disadvantages to every approach, modeling bond dynamics directly is more closely aligned with the nature of trading in the asset itself; bonds, like stocks, trade in denominated currencies (models that focus on interest rate dynamics are one-step removed, in the sense that they are inputs into a pricing process for an underlying asset. If we were examining stocks, for instance, one possible way to price options on stocks would be to model dynamics for factors that affect stock price-perhaps a stochastic process for the market risk premium, as an example-then derive a stock price and then price the option. Black and Scholes decided to model the stock price directly. In that vein, we are returning to directly modeling the bond's price). In this fashion, our approach adds to the literature on option pricing where the underlying asset price is modeled directly.
In what follows, we provide a short review of extant interest rate models. We discuss a stochastic process for the bond price and derive the PDE. We then present the asymptotic solution for a European style call option. Finally, we discuss the properties of our model, including the commonly used Greeks.

Previous Literature: Interest Rate Models
Vasicek [2], Rendleman and Bartter [3], Black, Derman and Toy [4], Black and Karasinski [5], are examples of short rate models, where parameters are generally chosen to make the evolution of the spot rate fit the current term structure. In Vasicek [2], the bond's price follows from an evolution of a spot rate that follows a Markov process. Vasicek [2] uses the specific example of an Ornstein-Uhlenbeck process to derive the representation of the bond price. Rendleman and Bartter [3] discuss the pricing of bond options where interest rates follow a binomial process with parameters chosen so that the price relative (ratio of) interest rates approximate a lognormal distribution. In the study of Black, Derman, and Toy [4], all security prices and rates depend on only one factor, the short rate. The current structure of long rates and their estimated volatilities are used to construct a tree of possible future short rates, which may then be used to value interest-rate-sensitive securities. Black and Karasinski [5] present a short-term interest rate one-factor model for bond and option pricing that allows the target rate, mean reversion and local volatility to vary deterministically through time. The distribution of possible short rates is lognormal in their model, so the rate will not fall below zero. Jamshidian [6] derives a closed-form solution for an option on a zero-coupon bond assuming the process developed by Vasicek [2].
The studies of Ho and Lee [7], Heath, Jarrow and Morton [8,9], Duffie and Kan [10], Brace, Gatarek, and Muiela [11], and Musiela and Rutkowski [12] are examples of term structure models. These models take a given term structure and model its behavior. Ho and Lee [7] derive an arbitrage free discrete trading economy, where bond prices fluctuate according to a single factor path independent binomial process. Securities are priced relative to a term structure of interest rates (spot rate and forward rates). The path independence condition is equivalent to constant volatility in the forward rate process. Heath, Jarrow, and Morton [8,9] develop a framework for pricing interest rate claims based on the term structure of forward rates. The model (HJM for convenience) treats the entire forward rate curve at a point in time as an input, and evolves the curve according to a stochastic process. The stochastic process for the forward rate is the sum of a simple diffusion and a deterministic drift, and is more general than Ho and Lee [7]; the forward rate process may have non-constant volatilities for instance. Duffie and Kan [10] analyze a multifactor model of the term structure, where the yield factors form a Markov process. In Brace, Gatarek, and Muiela [11], a class of term structure models with a volatility of the lognormal type is analyzed in the general HJM framework. The corresponding market forward rates do not explode and are positive and mean reverting. Musiela and Rutkowski [12] describe the properties and problems associated with continuous-time presentations of term-structure models and specify the general properties for wide classes of models.
Two highly recommended, and excellent discussions of interest rate models, and their associated problems, are provided by Björk and Gaspar [13] and Rebonato [14]. Björk and Gaspar present a collection of important topics and conditions for arbitrage free pricing of interest rate claims, including a discussion of the market price of risk, short-rate models, and term structure models. For those just starting to study interest rate models, Sections 1 and 2 provide a host of useful theorems and definitions. In addition, there is an excellent discussion on the difficulties of operationalizing these models, and how to approach the realization of essentially infinite-dimensional models in finite dimensional space. Rebonato discusses the evolution of interest-rate models from the early short-rate models to more recent developments in term-structure models and provides a far more comprehensive discussion than what can be provided here. In fact, two quotes from Rebonato provide the inspiration for this work: Therefore, while today's models are indubitably more effective (at least in certain respects) than the early ones, and, therefore, one can certainly speak of an 'evolution' of interest-rate models, this does not necessary imply that certain choices, abandoned in the past, might not have ultimately led a more fruitful description of interest-rate derivatives products if they had been pursued more actively. [14] (p. 672) As well as the following: As for the log-normal assumption for the bond price, this was partly satisfactory, because it did not constrain the price to be smaller than face value, thereby assigning a non-zero probability to negative interest rates . . . What gave greater discomfort, however, was the so-called pull-to-par problem. [14] (p. 683) The pull-to-par issue is an important facet of bond returns, as noted in such works as that of Esquível, Gaspar, and Beleza Sousa [15]. Esquível et al. use pull-to-par returns to analyze implicit default propensities from observed bond prices. While their paper does not explicitly put forth a stochastic process to create bond prices, and instead is based on market data, they do note that one could imagine a stochastic process that might create price realizations: Although this is not necessary for what follows, we may suppose that it is a stochastic process-defined on a complete probability space (Ω, A, P)-from which some realization is observed. [15] (p. 82) The valuation of simple contingent claims under simple continuous time models is no longer novel; however, we believe that there is value in the return to the "pull-to-par" issue and the pursuit of a path "abandoned in the past." We propose an approach based on the bond's price, rather than a set of interest rates/forward rates, that better captures the pullto-par dynamics of a bond than the traditional geometric Brownian motion (we do note that the pull-to-par problem can be handled in the classical short-rate and term-structure models of interest rates, and so when we are referring to addressing the pull-to-par problem, it is with respect to a model of bond price dynamics. We thank an anonymous referee for pointing out our need to clarify this point).
By going back to the bond's price, we are not limited to processes for interest rates that, due to computational complexities, are typically limited to three factors (or less) to characterize a yield curve (partly to avoid overfitting problems and to allow for some tractability). There is no need to fit a curve or restrict the evolution of rates so that the resulting term structure fits with observable bond prices. By starting with a process for the bond, we can, in a fashion similar to Black and Scholes develop a solution that does not require assumptions about the term structure or the relation of rates across maturities. Our model would also be an example of a stochastic process that could characterize zero coupon bond prices as described by Esquível, Gaspar, and Beleza Sousa [15].

A Process for Zero-Coupon Bonds
To begin, we assume that the zero-coupon bond's price before maturity can be any positive value, including values higher than par (nominally negative interest rates). We make very few assumptions on the implied interest rates other than rates are made up of two components: a risk-free rate (time value of money) and a risk premium (if the bond is risky). We assume that interest rates in a term structure may or may not be related to such frameworks as those of the "expectations hypothesis," "preferred habitat theory," "liquidity preference," or "market segmentation," though we are not specifying any specific structure (an excellent discussion of these can be found in the study of Cox, Ingersoll, and Ross [16] in the introduction). As such, the rate implied by our zero-coupon bond price (the process of which is given later) can be expressed as an average of a current short rate and future short rates, if so chosen, that is, the bond's current price reflects the term structure to its maturity at the current time.

Model Definition
We start with dynamics for the zero-coupon bond price.
where B is the price of a zero-coupon bond that is worth F at maturity. T B is the maturity time for the bond, t is the current time, z is a Brownian motion, and σ is the volatility of the bond's return at full duration (initial duration, which for the zero-coupon bond equals the bonds maturity). The dz and σ are unique for each bond maturity, and are not to be seen as common terms across all maturities (as noted by an anonymous referee, there is empirical research indicating that only three or four independent stochastic components are needed to describe the market. Our assumption of different Brownian motions across the maturities is made to ensure no arbitrage). The first term (drift term) represents the pull-to-par effect that takes place if the implied interest rate on the bond does not change as time changes. The second term (shock term) represents the riskiness of the bond, which is a function of its duration (in this case, the remaining percent of initial duration). Therefore, the volatility constant, σ, is scaled down over time as the bond's remaining time to maturity becomes shorter. For simplicity of discussion, we refer to the model as the pull-to-par model (or PTP). Of course, there are many different ways in which the dynamics could be proposed. We chose to try and capture two important characteristics in the bond price dynamics (pull-to-par and duration), while adopting a form that would adjust for the problems of the Black and Scholes geometric Brownian motion described in Rebonato. Yet, the dynamics we chose are still suggestive of that in Black and Scholes, so the derivation and solution will be familiar to many. While our focus is on characterizing risky bonds, we do note that our process could conceptually represent a risk-free bond (a somewhat trivial result but important in the resulting model, which, like Black and Scholes, has a constant risk-free rate as part of the solution). In the case of a conceptually true risk-free bond, B r f , the volatility for it would be σ r f = 0, (there would be no uncertainty of any type-no risk premiums) and the term in (1) would result in the bond earning a constant continuously compounded risk-free In our model, in order for there to be no arbitrage, we need only assume that a risky bond, B(0, T), with volatility σ, has the following properties: The resulting initial (expected) rate of return on the risky bond, r B , could therefore be expressed as follows: where λ is the risk-premium(s) for bond B. The process in (1) therefore can represent both risk-free and risky bonds. Since our formulation is general, one could imagine the risk-premium representing such concepts as inflation risk, liquidity risk, default risk, etc. (or some combination of these and others). There are some limits to our model in conceptualizing the default risk since we do not specifically cover the recovery rates for bonds in default but, as noted later, our model would work in situations where the recovery rate is below the strike on the option.
The following exposition on the derivation of the governing PDE is for those interested in seeing the steps in comparison to the well-known Black-Scholes solution to a European option on stocks. Those familiar with the derivation might prefer to move directly to the solution in Section 4.
Suppose that there is a European style call option on this bond (a function of the bond price and time), denoted as follows: By Ito's Lemma, with Equation (1), we have the following: where C t and C B are partial derivatives of C with respect to t and B, respectively, and C BB is the second partial derivative of C with respect to B. We can form a hedge portfolio by combining N B shares of the bond and N C call options. The value of the hedge portfolio (V H ) is given by the following: Differentiating (4) and substituting for dB from (1) and dC from (3) and rearranging, we obtain the following: A perfectly hedged portfolio should be risk-free and earn the risk-free rate, r; therefore, the following holds: To hedge, we need to set the following condition: Substituting (4), (6) and (7) into (5), and dividing through by dt and N C , we obtain the following: In addition, we have the following initial and boundary conditions for (8): where E is the exercise price of the call option (Condition 8b could also be written as C(γ,t) = 0 for recovery value γ, where 0 ≤ γ < E < 1. Note, our model would not be appropriate in the case of default where the recovery value γ has the following property 0 ≤ E < γ < 1).

Solution
Under the assumption that T B is much larger than T C (i.e., T C T B 1), we use a regular perturbation method to find an asymptotic solution for (8). For the mathematical definition of T C T B 1 and the asymptotic solution procedure, see Kevorkian and Cole [17]. We further assume that T C = O(1) and look for an asymptotic solution of the following form: Which, for our PDE and its initial and boundary conditions, is equal to the following: The derivation of this solution is rather involved. We detail it in the Appendix A and omit it here for brevity. If we let τ = T C − t (time moves forward), this becomes the following: where M is defined as follows: Additionally, N(z) is defined as follows: and d 1 and d 2 are defined as follows: We also have α defined as: α = 2r σ 2 and erfc(z) as:

Discussion of Results
The first two terms in the solution, (Equation (10)) is equivalent to the Black-Scholes model. The additional terms provide an adjustment to Black-Scholes that results from the pull-to-par effect and dampening effect on volatility that are present in the original stochastic Equation (1). In particular, the third term represents the leading order approximation to that adjustment under the assumption that T B is much larger than T C . Mathematically, the solution in (10) decomposes the PTP model into a descending order of components with the leading order O(1), a Black-Scholes model. This is given by the first two terms in (10) and it is governed by the homogeneous Black-Scholes Equation (A1a), subject to the initial and boundary conditions (8a, b). The next order O 1 T B in (10), which is the leading adjustment to the Black-Scholes model, is governed by a forced Black-Scholes Equation (A1b), subject to zero initial and boundary conditions. The forcing term here is given by the gamma of the Black-Scholes model, multiplied by σ 2 B 2 t. Given these effects, the PTP model results in a solution that is less than that of a pure Black-Scholes geometric Brownian motion solution.
To illustrate the differences, we have included a few computations, which are not meant to be exhaustive (see Table 1). Using the daily 20-year constant maturity treasury yields for the year 2020 (Treasury.Gov, 1 April 2021), we created daily zero-coupon bond prices, from which we then derived the estimated daily price returns. The choice of treasury data was two-fold: (1) the data are easily accessed by all interested, and (2) since our paper is not focused on any particular set of risk premiums, treasury bonds, while generally considered default free, are nonetheless risky (it is generally commonly held that treasury debt has at least inflation risk). The data were largely used to give us a set of inputs for illustrative purposes. The annualized standard deviation of the returns serves as a basis for the volatility input in the tables that follow. The annualized standard deviation of the hypothetical 20-year zero-coupon bond (assuming 256 trading days) was 20.20%. Given this, table computations have volatilities of 15%, 20%, and 25% for the option values and the Greeks. For the year 2020, short-term treasury rates (3 months, 6 months, and 1 year) ranged from near zero to just over 1.5%, with an average of slightly below 0.5%. Based on this, an input of 0.5% for the short-term risk-free rate was used as an input to the model.  13.1937 In addition to the hypothetical values, an analysis of the contribution of the terms O(1) and O 1 T B are provided in Table 2. As noted earlier, our approximation solution assumes that the maturity of the bond is considerably larger than that of the option.
We illustrate the contribution for a range of bond maturities against option maturities of 3 months, 6 months and 1 year. As the bond maturity increases, the contribution of the O 1 T B generally decreases.

Properties of Greeks
Based on the solution in Equation (10), we derived the standard "Greeks" associated with the call option. The solutions are presented below: The call option's delta is equal to the following: As with the options value, we also calculated values for the call's delta. Table 3 illustrates these for the same strike, short rate, and set of volatilities and option maturities.
As with the value, the delta is generally a bit lower than that of an equivalent Black-Scholes model. Graphically, Figures 1 and 2 present the deltas for differing option maturities and volatilities, respectively.     The call option's gamma is equal to the following:  The call option's gamma is equal to the following: Table 4 presents the gamma values with the same underlying data used for the values and deltas. As with the delta, the gamma is also somewhat lower than that of an equivalent Black-Scholes model. Graphically, Figures 3 and 4 present the gammas for differing option maturities and volatilities, respectively. The call option's theta is equal to the following:   Table 5 presents the theta values for our set of inputs. Similar to the delta and gamma results, the theta results are lower in magnitude than their Black-Scholes counterparts.  The gamma results present a "corner" for the at-the-money options. This is likely due to the method of derivation, i.e., the approximate solution. In solving the differential equation at order O 1 T B , the Laplace transform method was used, and the solutions for intervals x > 0 and x < 0 were obtained first, then matched at x = 0 (i.e., B = 80).
In reality, there probably should be no "corner"; this will likely disappear with the addition of higher order terms (though computationally, these terms are quite time-consuming to derive).
The call option's theta is equal to the following: Table 5 presents the theta values for our set of inputs. Similar to the delta and gamma results, the theta results are lower in magnitude than their Black-Scholes counterparts. Figures 5 and 6 illustrate the thetas for different option maturities and volatilities, respectively.    The call option's vega is: Finally, Table 6 presents the vega values for our set of inputs. Not surprisingly, the results mirror those of the other Greeks, with the PTP model values being lower than their Black-Scholes counterparts.  0  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99  100 Theta Bond Price Thetas for 0.5 year 80 strike call option on 20 year bond Volatility 15% Volatility 20% Volatility 25% The call option's vega is: Finally, Table 6 presents the vega values for our set of inputs. Not surprisingly, the results mirror those of the other Greeks, with the PTP model values being lower than their Black-Scholes counterparts. Figures 7 and 8 illustrate the thetas for different option maturities and volatilities, respectively.

Conclusions
We derive an asymptotic solution for call options on zero-coupon bonds, assuming a stochastic process for the price of the bond, rather than for interest rates in general. Our approach of modeling bond price dynamics has several important advantages: (1) we are not restricted to assumptions regarding the relationship of interest rates in a term structure; (2) modeling the bond dynamics is computationally simpler than many spot and term structure models with multiple factors; and (3) the resulting solution to our formulation is similar to the familiar Black-Scholes model for European-styled call options on stocks. In our model, we constrain the stochastic process to have a pull-to-par effect and price volatility scales by the remaining original duration of the bond. In this manner, our process more realistically represents the price path of a bond than otherwise given by the geometric Brownian motion proposed in Black and Scholes. The pull-to-par and dampening effect on the volatility results in prices and Greeks that are below those of the Black-Scholes model, in general.

Conclusions
We derive an asymptotic solution for call options on zero-coupon bonds, assuming a stochastic process for the price of the bond, rather than for interest rates in general. Our approach of modeling bond price dynamics has several important advantages: (1) we are not restricted to assumptions regarding the relationship of interest rates in a term structure; (2) modeling the bond dynamics is computationally simpler than many spot and term structure models with multiple factors; and (3) the resulting solution to our formulation is similar to the familiar Black-Scholes model for European-styled call options on stocks. In our model, we constrain the stochastic process to have a pull-to-par effect and price volatility scales by the remaining original duration of the bond. In this manner, our process more realistically represents the price path of a bond than otherwise given by the geometric Brownian motion proposed in Black and Scholes. The pull-to-par and dampening effect on the volatility results in prices and Greeks that are below those of the Black-Scholes model, in general.

Data Availability Statement:
The data that support the findings of this study are available from the United States Treasury at Treasury.Gov. These data were obtained from the following resources available in the public domain: https://www.treasury.gov/resource-center/data-chart-center/interestrates/pages/textview.aspx?data=yield. The data was accessed on 15 February 2021.

Conflicts of Interest:
The authors have no relevant financial or non-financial competing interests to report.

Appendix A. Derivation of Asymptotic Solution
Here, we derive the asymptotic solution to the PDE (8), subject to the initial and boundary conditions (8a) and (8b).
The assumption is as follows: T B is much larger than T C (for example T B is 60 months and T C is one month). We consider the solution as C(B, t; T B ), and we look for a solution in the following form (regular perturbation): Substituting in Equation (8), we have the following: All of the terms independent of T B can be collected as follows: All of the terms multiplied by 1 T B can be collected as follows: All of the terms multiplied by 1 can be collected as follows: The two conditions (8a) and (8b) therefore become the following: The solution to Equation (A1a) with the conditions C 0 (B, T C ) = Max(B − E, 0) and C 0 (0, t) = 0 is the same as that for Black-Scholes, as follows: where N(z) is defined as follows: d 1 is defined as follows: and d 2 is defined as follows: Next, we turn to Equation (A1b) with the conditions C 1 (B, T C ) = 0 and C 1 (0, t) = 0: First, we make the following change of variables in Equation (A1b): The PDE for v 1 (x, τ) is (after letting α = 2r σ 2 ) as follows: In order to eliminate the v 1 and v 1 x terms in Equation (A7), we let the following hold: The PDE for u 1 (x, τ) is: where u 0 (x, τ) is the Black-Scholes solution presented earlier after all the substitutions above, i.e., the following: Taking the derivatives of u 0 (x, τ) with respect to x in the term u 0 xx − αu 0 x + α 2 −1 4 u 0 in Equation (A9) and simplifying, we find the following: Therefore, Equation (A9) becomes the following: From the initial condition, C 1 (B, T C ) = 0, in (A2) we find the initial condition for u 1 (x, τ) as u 1 (x, 0) = 0. Now, we take the Laplace transform U 1 (x, s) = L{u 1 (x, τ)}.
We note that the following formulas and properties will be useful: Taking the Laplace transform of Equation (A10), using (A11b)-(A11d) and simplifying, we obtain the following: For x > 0, we have the following solution form for (A12): where U 1p is a particular solution to (A12). In this case U 1p is the following: The boundedness of U 1 (x, s) as x → ∞ requires A 1 = 0. So, we have the following: Similarly, for x < 0 we have the following solution form for (A12): x e x √ s The continuity of U 1 (x, s) at x = 0 gives B 1 = A 2 . The continuity of U 1 x at x = 0 gives B 1 = 3 4σ 2 s (5/2) − T C 2s (3/2) . Therefore, we have the following: we use (A11a), (A11d) and (A14), and then simplify to obtain the following: Applying (A11b), (A11d), (A14), and (A15) to the inverse Laplace transform of Equation (A13) and then simplifying we find the following: Now, the function values of C 1 (B, t) can be computed as follows: Given B and t, we compute τ and x using (A6), i.e., the following: , and x = ln B E then u 1 (x, τ) can be calculated using (A16), and v 1 (x, τ) from (A8a). Finally, C 1 (B, t) is available using (A6), i.e., C 1 (B, t) = Ev 1 (x, τ). Therefore, our final perturbation solution, C(B, t; T B ) = C 0 (B, t) + 1 T B C 1 (B, t) + O 1 T B 2 becomes the following: where M is defined as follows: