- Email: [email protected]
Quantifying risks with exact analytical solutions of derivative pricing distribution
Physica A 471 (2017) 757–766
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Quantifying risks with exact analytical solutions of derivative pricing distribution Kun Zhang a,1 , Jing Liu b,1 , Erkang Wang a , Jin Wang a,b,c,d,∗ a
State Key Laboratory of Electroanalytical Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun, Jilin, 130022, PR China b
Quantitative Finance Center, Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY, 11794, United States c Department of Chemistry and Physics, Stony Brook University, Stony Brook, NY, 11794, United States d
Harriman School of Business Management, Stony Brook University, Stony Brook, NY, 11794, United States
highlights • Path integral approach to the statistical distributions of option pricing. • Analytical expressions of statistical distribution of bond and bond option pricing. • Statistical fluctuations and fatty tails of bond and bond option pricing.
article
info
Article history: Received 15 May 2016 Received in revised form 30 September 2016 Available online 28 December 2016 Keywords: Option pricing Distribution Risk
abstract Derivative (i.e. option) pricing is essential for modern financial instrumentations. Despite of the previous efforts, the exact analytical forms of the derivative pricing distributions are still challenging to obtain. In this study, we established a quantitative framework using path integrals to obtain the exact analytical solutions of the statistical distribution for bond and bond option pricing for the Vasicek model. We discuss the importance of statistical fluctuations away from the expected option pricing characterized by the distribution tail and their associations to value at risk (VaR). The framework established here is general and can be applied to other financial derivatives for quantifying the underlying statistical distributions. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Quantitative finance started from the study of stochastic processes. The foundation for understanding the stochastic processes was laid down in the last century thorough Brownian motion [1] by Bachelier, Einstein and Wiener [2–4]. In 1951, Ito established the stochastic calculus for Brownian motion [5]. The stock price fluctuations in time on the financial market can be described by a stochastic process. Later in 1950s, Samuelson laid down the foundation of option pricing theories via expectations by quantifying both macro- and microeconomics through mathematics [6]. In 1973, Black, Scholes and Merton derived the Black–Scholes partial differential equations for option pricing [7,8].
∗
Corresponding author at: Department of Chemistry and Physics, Stony Brook University, Stony Brook, NY, 11794, United States. E-mail address: [email protected] (J. Wang).
1 Kun Zhang and Jing Liu contributed equally to this work. http://dx.doi.org/10.1016/j.physa.2016.12.044 0378-4371/© 2016 Elsevier B.V. All rights reserved.
758
K. Zhang et al. / Physica A 471 (2017) 757–766
In 1977, Boyle related the pricing of options to the simulations of random asset paths, assessing the important roles of Monte Carlo simulations in finance [9]. More recently, in 1994, Derman and Kani [10], Dupire [11] and Rubinstein [12] considered the stochastic nature of the volatility and generalized the Black–Scholes equation in this situation. In 1977, Vasicek investigated the pricing interest rate products such as bonds. He modeled a short-term interest rate as a random walk and found that the interest rate derivatives can be valued using equations similar to the Black–Scholes equation [13]. In 1986, Ho and Lee found a way of improving Vasicek model with the idea of yield curve fitting or calibration [14]. In 1992, Heath, Jarrow and Morton modeled the random motion of the whole yield curve and generalized Ho and Lee’s approach [15]. In addition to the expected pricing discussed, quantifying the statistical fluctuations of equity and bond products as well as their derivatives is critical for addressing the associated risks. Basel committee on banking supervision published on international convergence of capital measurement and capital standards in 2006 [16]. Therefore, banks are now required to assess the risk quantitatively. Moreover, as one can see clearly from the recent financial crisis started from Wall Street and propagated through the whole world, the quantitative risk assessment is crucial for financial stability. Progresses have been made on the equity (stock) pricing models. Furthermore, analytical models of option pricing for both equity and credit so far have often been focused on the expected pricing. The distributions of option pricing revealing high order statistical fluctuations critical for risk assessments, are still challenging to obtain analytically [17–23] despite of numerical efforts such as Monte Carlo approach [9]. In this study, we established a quantitative framework using path integrals to obtain an exact analytical formula for the statistical distribution of the bond and bond option pricing for the Vasicek model. We discuss the importance of statistical fluctuations away from the expected option pricing characterized by the tail of the distribution and the associations to value at risk (VaR) critical in the practice of risk control. Quantifying statistical distributions is also very important in characterizing complex physical and biological systems as well as single molecules dynamics [24–26]. 2. Model In order to explore the statistical distribution of the option pricing, it is important to start with a suitable method. In this study, we use a path integral method [24–33]. This method has its advantages that the expected pricing and associated distribution can all be quantitatively addressed in a relatively straight forward way. We start with a bond model. Bond is a contract, paid for up-front, that yields a known amount on a known date in the future, the maturity date T . Bonds maybe issued by either governments or companies. The main purpose of a bond issue is the raising of capital, and the up-front premium can be thought of as a loan to the government or the company. The issue related to valuing bond pricing lies on how much one should pay to get a guaranteed amount in certain times. In regarding to the bond options, one aims to find fair values of the contract. In this section, we will discuss the path integral formulation and the expected pricing of bond and bond options. In the results and discussion section, we will concentrate on the distribution of the derivative pricing of bond and bond options. 2.1. The Vasicek model of bond with fluctuating interest rate The interest rate for bonds is not a constant and stochastically fluctuates. The Vasicek model for bonds [13] is usually studied from a stochastic differential equation for the short term interest rate r [13,34,35]: dr dt
= a(b − r (t )) + σ f (t ).
(1)
The rate dynamics (change) is controlled by two terms: the drift term (first term) and the stochastic term (second term). The parameters a and b model the mean reversion, σ is the volatility of the short term interest rate. a, b and σ are constants in time. The whole dynamics is parameterized in terms of the continuous time variable t. Looking at the drift term, the interest rate will tend to relax (decrease) towards the mean for large r and move up on average for small r. The ‘‘force’’ f (t ) represents the stochastic fluctuations:
⟨f (t )f (t )⟩ = δ(t − t ),
⟨f (t )⟩ = 0.
(2)
The brackets ⟨⟩ indicate a mean value with respect to a Gaussian distribution of zero mean and variance of one. The interest rate term structure models (including Vasicek model) are for pricing risk-free government bonds. For corporate bonds, additional credit spread and credit risk model are needed. 2.2. The functional for the short term interest rate From Gaussian noise properties of the stochastic force, we can write down the distribution of ‘‘force’’ f as: 1 p({f (s)}) = e− 2
dsf 2 (s)
.
(3)
K. Zhang et al. / Physica A 471 (2017) 757–766
759
What we are interested in is the distribution of the interest rate r. One can substitute stochastic force f from Eq. (1) and taking care of the variable transformation (Jacobian) to obtain the probability density functional for the short rate r (s) [32]: p({r (s)}) = Ne
−
1 2σ 2
2 ds dr −a(b−r (s)) ds
.
(4)
We will be only interested in the family of short rates within a specific time window, e.g. between t and T . Therefore, the integration on the right hand side. of Eq. (4) has to be constrained to the interval[t , T ]. 3. Statistical methods 3.1. The expected bond pricing as a path integral Now we discuss the bond price as a specific example for a claim. The bond price today, i.e. at time t, for a monetary unit promised at some later time T can be stated as an expectation value of the discount function (from the maturity time T to the current time t). Let P (t , T ) = e−
T
r (s)ds
t
: When the maturity time T is reached, then the bond price is the one promised (here it is one for
convenience). But at current time t less than T , the promised amount will be discounted by a factor of e− words, the bond is worth less at present before the maturity time.
T ⟨P (t , T )⟩ = EQ e− t
dsr (s)
|r (t ) = r .
T t
r (s)ds
. In other
(5)
EQ is the expectation value. The classical Vasicek approach (as described in [33,34]) has been to write down the corresponding PDE for P (t , T ) rather than evaluating the expectation value in Eq. (5) directly. Using the probability density functional in Eq. (4), we are now in the position to evaluate the expectation value in Eq. (4) for P (t , T ). The result is well known. It is explicitly given as follows, see details in Appendix A, [28–32]:
⟨P ⟩ = A(t , T ) exp(−B(t , T )r )
(6)
with
A(t , T ) = exp (−B(t , T ) + T − t )
σ2 2a2
σ2 2 −b − B (t , T ) 4a
(7)
and B(t , T ) =
1 − e−a(T −t ) a
.
(8)
3.2. The expected bond option pricing For an example of bond option, we look at the European call. A European call option is the call option at time 0 with a strike price k, expiring at time t on a bond, with maturity time T with T > t. In other words, for a European call option of the bound, if the bond price P is less than the striking price k at time t, then one gets nothing. On the other hand, if the bond price P is larger than the striking price k at time t, then one gets the difference between the bond price P and the striking price k. The expected price of a European call at time 0 with strike k expiring at time t on zero bond with maturity time T , with T > t is given by,
t ⟨c ⟩ = E e− 0
r (s)ds
max(P (t , T ) − k; 0)|r (0) = r
where P (t , T ) = A(t , T )e−B(t ,T )r . The constraint on P-k is transformed to r: r (t ) <
1 B(t , T )
ln
A(t , T ) k
= K.
Thus
t ⟨c ⟩ = E e− 0
r (s)ds+ln(A(t ,T )e−B(t ,T )r (t ) −k)
|r (0) = r .
760
K. Zhang et al. / Physica A 471 (2017) 757–766
We can substitute x(s) = b − r (s), and the constraint translates into x(t ) > b − K . The final analytical result is known and given [28–32] as follows:
⟨c ⟩ = c (0, T )N (h) − kc (0, t )N (h − σp )
(9)
where N (x) is the cumulative normal distribution and h=
1
σp
c (0, T )
+
kc (0, t )
σp
(10)
2
and
σp =
σ a
(1 − e
−a(T −t )
1 − e−2at
)
(11)
2a
where
c (0, T ) = EQ e−
T t
dsr (s)
c (t , T )|r (0) = r
(12)
c (t , T ) = A(t , T )e−B(t ,T )r
(13)
and
c (0, t ) = EQ e−
T t
dsr (s)
|r (0) = r .
(14)
See details in Appendix A [28–32]. Here using a path integral formulation, the formula derived earlier using partial differential equation for expected bond pricing with Vasicek model can be recovered [32,34–36]. The analytical tractability of the bond option formulas for Vasicek dynamics is dependent on the simple form of the equation for the bond price itself. 3.3. The derivative price distribution from path integral method The advantage of the path integral method is that we can obtain not only the expectation value but also directly the statistical distribution of the option pricing. We will show below the examples of bond pricing and bond option pricing to illustrate the power of path integral method in quantifying and obtaining the analytical formula for the statistical distribution of option pricing. 3.4. The bond price distribution from path integral method Since from previous section we know how to compute the expectation value of the bond price P, we now turn to the question how we will obtain the statistical distribution of the bond price P. T
As discussed, P (t , T ) = e− t r (s)ds is a discount function as discussed and its expectation values gives the bond price. The interest rate, however, follows stochastic fluctuations characterized by the probability density functional p({r (s)}) given by T Eq. (4). If we know the statistical distribution of z = − t r (s)ds since P (t , T ) = exp[z ], we will know the distribution of
T
P (t , T ). The distribution of z can be calculated by f (z ) = ⟨δ[z − (− t r (s)ds)]⟩. This gives the density or probability of the z. The expectation is over the realizations of r, p({r (s)}). The physical meaning is clear. The density or probability of z is obtained by counting the events or peaks of z characterized by the delta function at locations with different realizations T of r. Fortunately, the delta function can be transformed into the exponential in Fourier space k, exp[ik(z + t r (s)ds)]. This allows us to map the current problem to essentially almost the same problem of expected pricing done before with path T
integrals. Notice that the expected bond price is obtained by expectation of the discount function P (t , T ) = e− t r (s)ds over T p({r (s)}) whereas now for distribution we are performing the expectation of exp[ik(z + t r (s)ds)] over p({r (s)}) with an additional Fourier integral k to carry out. Therefore, the probability density function of z and therefore P can be obtained. We give some derivation steps to follow in the result section (for those who care about the detailed procedures, see details in Appendix A, for those who do not care about intermediate steps, jump directly to Eq. (18) for the results). 3.5. The bond option price distribution from path integral method We now turn to the question how we will obtain the statistical distribution of the bond option price P. For a European call option of the bond with maturity time T , if bond price P is larger than the striking price k at time t, then one gets the difference between the bond price P and the striking price k. Otherwise, one gets nothing. This provides a constraint in the expectation over all the possible realizations or interest rate r. This is the main difference between the bond pricing and bond option pricing. In other words, the bond option price c (before the expectation) follows from a discount
K. Zhang et al. / Physica A 471 (2017) 757–766
function as discussed e−
t
0 r (s)ds
761
relaxing from striking time t to current time (t = 0) with a constraint P larger than k at t
striking time (option price (before expectation) c = e− 0 r (s)ds max[P (t , T ) − k; 0]). While the interest rate follows a stochastic fluctuations characterized by the probability density functional p({r (s)}) given t by Eq. (4). If we know the distribution of z = ln(P (t , T ) − k) − 0 r (s)ds, we will know the distribution of c since c = exp[z ].
t
The distribution of z can be calculated by f (z ) = δ[z − (ln(P (t , T ) − k) − 0 r (s)ds)], This gives the density or probability of the z. The expectation is over the realizations of r, p({r (s)}). The physical meaning is clear. The density or probability of z is obtained by counting the events or peaks of z characterized by the delta function at locations with different realizations of r. Furthermore, the delta function can be transformed into the exponential in Fourier space g, exp[ig {z − (ln(P (t , T ) − t k) − 0 r (s)ds)}]. This allows us to map the current problem to essentially almost the same problem of expected option pricing done before with path integrals. Notice that the expected option price is obtained by performing expectations t
over the discount function c = e− 0 r (s)ds max[P (t , T ) − k; 0] (with the constraint P larger than k) over the realizations of interest rate with probability density p({r (s)}) whereas now for distribution we are performing the expectation values t of exp[ig {z − (ln(P (t , T ) − k) − 0 r (s)ds)}] over p({r (s)}) (with a constraint P larger than k) with an additional Fourier integral g to carry out. Therefore, the probability density function of z and therefore option price c can be obtained. We give some derivation steps to follow in the result section (for those who care about the detailed procedures, see details in Appendix A. For those who do not care about intermediate steps, jump directly to Eq. (26) for the results). 4. Results and discussion 4.1. The results and derivations for bond price distribution from path integral method
We are now trying to find out the distribution of P (t , T ), that is f (P ). First we define z as, z = − and z is related by P = ez ⇔ z = ln P ⇔ dz =
1
T t
r (s)ds we see that P
dP .
P
We can also see the relationship between the distribution of z, f (z ) and the distribution of P, f (P ) : f (P )dP = f (z )dz f (P ) = f (z ) f (P ) =
1 P
dz dP
f (z ).
Thus, to find out f (P ), we only need to find out f (z ).
f (z ) = δ
T
r (s)ds
z+
(15)
t
and we can write down the path integral representation for the distribution of z
f (z ) =
1 2π
∞ −∞
∞
dk −∞ dr (T )
r (T ) r (t )
∞
−∞ dr (T )
Dr (s) exp − 2σ 2
r (T )
1
T t
ds
Dr (s) exp − 2σ1 2 r (t )
dr (s) ds
T t
ds
− a(b − r (s))
dr (s) ds
2
+ ik z + 2 − a(b − r (s))
T t
dsr (s) .
(16)
+∞
Here, we use δ(x) = 21π −∞ eikx dk to calculate f (z ). We see that it is similar to the formula used for expected value of bond pricing except with one addition integration over k. −a(T −t )
Define, B(t , T ) = 1−e a we perform the integrals with respect to r (T ) and path integrals with respect to r (s) by applying the generating function for the harmonic oscillator [27] (details in Appendix A). We obtain: f (z ) =
X Y
=
1 2π
∞
dk exp −
−∞
σ2
σ2
2a
4a
(−B(t , T ) + (T − t ))k2 + 2
k2 B2 (t , T ) − ikB(t , T )x(t ) + ikb(T − t ) + ikz
where x(t ) = b − r (t ). Now, we only have the Gaussian integration with respect to k for f (z ) to carry out:
f (z ) =
1
π W (t , T )
exp
− ((−B(t , T ) + T − t )b + B(t , T )r + z )2 W (t , T )
2 2 where W (t , T ) = 2aσ2 (−B(t , T ) + T − t ) − σa B2 (t , T ).
Since f (P ) =
1 f P
(z ) and z = ln P.
(17)
762
K. Zhang et al. / Physica A 471 (2017) 757–766
Fig. 1. Statistical distribution of bond price with respect to (a) interest reversion rate a, (b = 0.01, r = 0.01, σ = 0.03, t = 10, T = 30); (b) interest reversion b, (a = 0.3, r = 0.01, σ = 0.03, t = 10, T = 30); (c) interest rate r, (a = 0.3, b = 0.01, σ = 0.03, t = 10, T = 30); (d) volatility, (a = 0.3, b = 0.01, r = 0.01, t = 10, T = 30); (e) current time t, (a = 0.3, b = 0.01, r = 0.01, σ = 0.05, T = 30); (f) maturity time T . (a = 0.3, b = 0.01, r = 0.01, σ = 0.05, t = 10).
Therefore, we finally obtain the distribution of bond pricing. f (P ) =
1 P
− ((−B(t , T ) + T − t )b + B(t , T )r + ln P )2 exp π W (t , T ) W (t , T ) 1
.
(18)
Eq. (18) gives the analytical solution of the distribution for option pricing. A(t , T ) and B(t , T ) are given by Eqs. (7) and (8). As we can see, the bond price follows a log normal distribution, with the average and width W (t , T ) controlled by the interest rate r at time t and associated mean reversion a and b, interest rate fluctuation strength or spread, σ current time t, and maturity time T . We plot the distribution of f (P ) with respect to mean reversion a, b, the current interest rate r, volatility, σ current time t and maturity time T in Fig. 1. We found that as a, b, r and t increase, the tails of the distribution become less prominent. On the other hand, as σ and T increase, the tails of the distribution become more prominent. This implies the more volatility and longer maturity time, the more statistical fluctuations of the bond pricing is. On the other hand, the closer of current time t to the maturity time T and larger the interest rate and associated mean reversion, the bond pricing has less fluctuations and is more determined. This provides the quantitative origin of the bond pricing fluctuations and important for quantitative risk control from the statistical distribution perspective. 4.2. The results and derivations for bond option price distribution from path integral method Let us consider the option price distribution of a European call at time 0 with strike k expiring at time t on zero bond with maturity time of T , with T > t. t
Let option pricing (before expectation) be c (0) = e− 0 r (s)ds max(P (t , T ) − k; 0) with r (0) = r, so that the option pricing A(t ,T ) is given by the discount function with the constraint written as: r (t ) < B(t1,T ) ln( k ) = K . We need to find out the distribution of option pricing c (0), that is f (c (0)). For convenience, we use f (c ) instead of f (c (0)) in the text. t First, we define variable z, so that, z = ln(P (t , T ) − k) − 0 r (s)ds then c and z are related, c = ez ⇔ z = ln c ⇔ dz =
1 c
dc .
K. Zhang et al. / Physica A 471 (2017) 757–766
763
In general, we can find the relationship between distribution of option pricing c, f (c ) and distribution of z, f (z ): f (c )dc = f (z )dz dz f (c ) = f (z ) dc 1 f (c ) = f (z ). c
+∞
To find out f (c ), we only need to find out f (z ). Since δ(x) = 21π −∞ eigx dg, we have
t r (s)ds δ z − ln(P (t , T ) − k) − 0 t r (s)ds δ z − ln(P (t , T ) − k) + 0 t ∞ 1 r (s)ds dg exp ig z − ln(P (t , T ) − k) + 2π −∞ 0 ∞ t 1 dg exp ig z − ln(P (t , T ) − k) + r (s)ds . 2π −∞ 0
f (z ) =
= = =
(19)
Therefore, the distribution of z, f (z ) is given in the path integral formulation as: 1 2π
f (z ) =
∞
K
−∞
dg −∞ dr (t )
r (t )
Dr (s) exp − 2σ1 2
r (0)
∞
−∞ dr (t )
t 0
ds
dr (s) ds
− a(b − r (s))
r (t )
Dr (s) exp − 2σ1 2 r (0)
t 0
ds
dr (s) ds
t + ig z − ln(P (t , T ) − k) + 0 r (s)ds . 2 − a(b − r (s))
2
(20) We can see that this formula is similar to the one used in the expectation of option pricing except the additional integration with respect to g. Substitute x(s) = b − r (s) and by applying the generating function for the harmonic oscillator [27] (see details in Appendix A), we obtain, f (z ) = XY , where X =
1
a
∞
∞
σ2 dx(t ) exp − 3 (e−at − 1 + at )g 2
dg 2π σ 2 sinh(at ) −∞ 2a b−K a a 2 − 2 (x2 (t ) cosh(at ) − 2x(t )x(0) + 2(eat − 1)(igC (x(t ) + x(0)) − g 2 C 2 )) − x (t ) 2σ sinh(at ) 2σ 2 2π
+ igz + igbt − ig ln(A(t , T )e−B(t ,T )(b−x(t )) − k)
(21)
and
Y =
a
∞
2π σ 2 sinh(at ) −∞
dx(t ) exp −
a 2σ
x2 (t ) − 2
a 2σ 2 sinh(at )
(x2 (t ) cosh(at ) − 2x(t )x(0))
(22)
−a(T −t )
with C = σ2a (1 − e−at ) and B(t , T ) = 1−e a . Now we perform the Gaussian integration with respect to g, we obtain: 2
X =
a 2π σ 2 sinh(at )
−
× exp
1
πV
∞
dx(t ) b−K
1−eat x(t )+x(0) a 1+eat
2 + bt + z − ln(A(t , T )e−B(t ,T )(b−x(t )) − k) V
−
a 2σ 2 sinh(at )
where V = − 2aσ3 ( 2
2(eat −1) eat +1
(x2 (t ) cosh(at ) − 2x(t )x(0)) −
− at ) =
2σ 2 a2
(2ah + t ), ah =
1 a
a 2σ 2
e−at −1 . e−at +1
x2 (t )
(23)
764
K. Zhang et al. / Physica A 471 (2017) 757–766
Also, we can perform the Gaussian integration with respect to x(t ). We obtain
Y =
a
π σ 2 (1 − e−2at )
2π σ 2 sinh(at )
a
ax2 (0)
exp −
σ 2 (1 − e2at )
.
(24)
Substitute back, r (t ) = b − x(t ), r (0) = r = b − x(0) f (z ) =
1
π
ap V
K
dr (t ) exp
−(−ah (r (t ) + r − 2b) + bt + z − ln(A(t , T )e−B(t ,T )r (t ) − k))2 V
−∞
− ap (b − r (t ) − e−at (b − r ))2
(25)
−at
1 where, ap = σ 2 (1−ae−2at ) , ah = 1a ee−at − . +1 Therefore we finally obtain the analytical form of the distribution of bond option pricing c as follows,
f (c ) =
1
πc
ap V
K
dr (t ) exp
−(−ah (r (t ) + r − 2b) + bt + ln c − ln(A(t , T )e−B(t ,T )r (t ) − k))2 V
−∞
− ap (b − r (t ) − e−at (b − r ))2 .
(26)
Eq. (26) gives the exact analytical solution of the distribution for option pricing (up to a one dimensional integral with respect to r (t )). A(t , T ) and B(t , T ) are given in Eqs. (7) and (8). We plot the distribution of option pricing f (c ) with respect to mean reversion a, b, current interest rate r (r (0)), volatility σ , striking time t and maturity time T in Fig. 2. We found that as a, b, r and t increase, the tails of the distribution become less prominent. On the other hand, as σ and T increase, the tails of the distribution become more prominent. This implies the more volatility and longer maturity time, the more statistical fluctuations of the option pricing is. On the other hand, the closer the striking time to the maturity time and the larger the interest rate and associated mean reversion, the option pricing has less statistical fluctuations and more determined values. This provides the quantitative origin of the bond option pricing fluctuations and important for quantitative risk control from the statistical distributionperspective. Once the statistical fluctuations are known, the underlying fluctuations are known and can be quantified. The value at risk (VaR) is often used to quantify the risk and is simply the variation or the second order moments which can be directly calculated from the distribution. The skew (3rd order moment) and kurtosis (4th order moment) are sometimes used for characterizing the asymmetry of the underlying observables. They are nothing but the third order and fourth order statistical fluctuations which can also be easily obtained from the distribution. For Gaussian distribution, the variable itself should not have skewness and kurtosis. Here the stochastic variable is the interest rate. However, the quantity we are interested in is the option pricing which is a nonlinear function of interest rate and no symmetry on option pricing is guaranteed. Therefore, the associated distribution of the option pricing is not expected to be Gaussian. The skewness and kurtosis of the option pricing are not expected to be zero for this reason. As one can see, statistical fluctuations can be completely characterized and quantified by the underlying distribution function. The statistical distribution can help to quantify the high order fluctuations. The high order fluctuations in option pricing can be obtained by integrating the same order of power for the option pricing over the corresponding statistical distribution. In principle once the distribution is known, the information on the higher order fluctuations such as VaR can be obtained. In practice, however, option prices are often used to calibrate the term structure model parameters (such as a and b in the Vasicek model), Fitting the mean option pricing is already a challenging. It may be even more challenging to obtain VaR directly. In strong fluctuation case, the expected or average pricing may not faithfully characterize the system since there are large statistical fluctuations around it. Instead, high order moments become more important. The distribution often has fatty tails. This implies the ‘‘intermittency’’ where rare event (low probability event at the tail) [24–26] can have significant effects on the nature of the system. It is exactly this rare event poses the high risk which we need to capture in order to prevent it to happen in finance. As we can see one cannot find such rare risk event by just looking at the expected values or the average. High order statistics through distribution is the only way to quantify the whole event. Such strongly fluctuating systems can be found not only in financial world as we have witnessed recently, but also every where else in physical and biological systems such as biomolecules, networks, weather, earth quakes, and evolution etc. [24–26,37–39].
K. Zhang et al. / Physica A 471 (2017) 757–766
765
Fig. 2. Statistical distribution of bond option price with respect to (a) interest reversion rate a, (b = 0.01, r = 0.01, σ = 0.2, t = 10, T = 30); (b) interest reversion b, (a = 0.35, r = 0.01, σ = 0.2, t = 10, T = 30); (c) interest rate r, (a = 0.35, b = 0.01, σ = 0.2, t = 10, T = 30); (d) volatility, (a = 0.35, b = 0.01, r = 0.01, t = 10, T = 30); (e) striking time t, (a = 0.35, b = 0.01, r = 0.01, σ = 0.2, T = 30); (f) maturity time T . (a = 0.35, b = 0.01, r = 0.01, σ = 0.2, t = 10).
5. Conclusions In this study, we obtained exact analytical solutions of the statistical distribution of the bond and bond option pricing for the Vasicek model with a path integral framework. This approach provides a quantitative way to obtain the underlying distribution of the derivative pricing which is critical for the financial risk control. We discuss the importance of statistical fluctuations away from the expected option pricing characterized by the tail of the distribution and the associations to value at risk (VaR). The framework established here is general and can be applied not only to financial derivatives but also to wide variety of problems in complex physical and biological systems, and single molecules for quantifying the underlying intrinsic statistical distributions. Acknowledgments We would like to thank Dr. Xidi Wang from Citibank, Dr. Jun Zang from Morgan-Chase Bank, Dr. Joy Zhang from Credit Suisse Bank, Peter G. Wolynes from Rice University and Prof. Raphael Douady from Stony Brook University/CNRS for interesting discussions. We would like to thank NSF (Grant No. NSF-PHY 76066) and NSFC (China) (Grant No. 91430217) for financial support. Appendix A. Supplementary data Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.physa.2016.12.044. References [1] R. Brown, A Brief Account of Microscopical Observations, London, (1827). [2] L. Bachelier, Theorie de la Speculation, Jacques Gabay Pub., Paris, 1995. [3] A. Einstein, On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat, Ann. Phys. 17 (1905) 549–560. [4] N. Wiener, The average of an analytic functional, Proc. Natl. Acad. Sci. 7 (1921) 253–260.
766 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
K. Zhang et al. / Physica A 471 (2017) 757–766 K. Ito, On stochastic differential equations, Mem. Amer. Math. Soc. 4 (1951) 1–51. P. Samuelson, Rational theory of warrant pricing, Ind. Manag. Rev. 6 (1965) 13–39. F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Polit. Econ. 81 (1973) 637–659. R.C. Merton, Theory of rational option pricing, Bell J. Econ. Manage. Sci. 4 (1973) 141–183. P. Boyle, Options: a Monte Carlo approach, J. Financ. Econ. 4 (1977) 323–338. E. Derman, I. Kani, Riding on a smile, Risk 7 (1994) 32–39. B. Dupire, Pricing with a smile, Risk 7 (1994) 18–20. M. Rubinstein, Implied binomial trees, J. Finance 69 (1994) 771–818. O.A. Vasicek, An equilibrium characterization of the term structure, J. Financ. Econ. 5 (1977) 177–188. T. Ho, S. Lee, Term structure movements and pricing interest rate contingent claims, J. Finance 42 (1986) 1129–1142. D. Heath, R. Jarrow, A. Morton, Bond Pricing and the term structure of interest rates: a new methodology, Econometrica 60 (1992) 77–105. Basel Committee on Banking Supervision (2006) International Convergence of Capital Measurement and Capital Standards: A Revised Framework Comprehensive Version Bank for International Settlements 1–347. P. Carr, M. Schroder, On the valuation of arithmetic-average Asian options: the Geman-Yor Laplace transform revisited. arXiv preprint math/0102080, 2001. M. Yor, On some exponential functionals of Brownian motion, Adv. Appl. Probab. 24 (1992) 509–531. R. Douady, Closed form formulas for exotic options and their lifetime distribution, Int. J. Theor. Appl. Finance 2 (01) (1999) 17–42. H. Geman, N. El Karoui, J.C. Rochet, Changes of numeraire, changes of probability measure and option pricing, J. Appl. Probab. (1995) 443–458. M. Jeanblanc, M. Yor, M. Chesney, Mathematical Methods for Financial Markets, Springer Science & Business Media, 2009. D. Brigo, F. Mercurio, Interest Rate Models: Theory and Practice, Springer Finance, Berlin, Heidelberg, New York, 2001. T. Bielecki, M. Rutkowski, Credit Risk: Modeling, Valuation and Hedging, Spinger Finance, Berlin, Heidelberg, New York, 2002. J. Wang, P.G. Wolynes, Intermittency of single molecule reaction dynamics in fluctuating environments, Phys. Rev. Lett. 74 (1995) 4317. S.S. Plotkin, P.G. Wolynes, Buffed energy landscapes: Another solution to the kinetic paradoxes of protein folding, Proc. Natl. Acad. Sci. 100 (8) (2003) 4417–4422. L. Mahadevani, J.M. Deutch, Influence of feedback on the stochastic evolution of simple climate systems, Proc. R. Soc. A 466 (2010) 993–1003. R.P. Feynman, A.R. Hibbs, Quantum Mechanics and Path Integral, McGraw-Hill Book Company, New York, 1965. B.E. Baaquie, Quantum Finance, Cambridge University Press, 2004. V. Linetsky, The path integral approach to financial modeling and options pricing, Comput. Econ. 11 (1998) 129. J. Dash, Path Integrals and Options. Part I, CNRS Preprint CPT88/PE.2206, 1988. J. Dash, Path Integrals and Options. Part II, CNRS Preprint CPT89/PE.2333, 1989. M. Otto, Using path integrals to price interest rate derivatives, 1999. arXiv:cond-mat/9812318v2. H.D. Feng, K. Zhang, J. Wang, Non-equilibrium transition state rate theory, Chem. Sci. (2014) (5) 3761–3769. P. Wilmott, Paul Wilmott on Quantitative Finance. I and II, John Wiley & Sons, New York, 2000. R. Rebonato, Interest-Rate Option Models, John Wiley & Sons, Chichester, 1996. F. Jamshidian, An exact bond option pricing formula, J. Finance 44 (1989) 205. M. Sasai, P.G. Wolynes, Stochastic gene expression as a many-body problem, Proc. Natl. Acad. Sci. 100 (2003) 2374–2379. M. Bottiglieri, C. Godano, On-offintermittency in earthquake occurrence, Phys. Rev. E 75 (2007) 026101–5. S.J. Gould, N. Eldredge, Punctuated equilibria: the tempo and mode of evolution reconsidered, Paleobiology 3 (1977) 115–151.
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Quantifying risks with exact analytical solutions of derivative pricing distribution Kun Zhang a,1 , Jing Liu b,1 , Erkang Wang a , Jin Wang a,b,c,d,∗ a
State Key Laboratory of Electroanalytical Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun, Jilin, 130022, PR China b
Quantitative Finance Center, Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY, 11794, United States c Department of Chemistry and Physics, Stony Brook University, Stony Brook, NY, 11794, United States d
Harriman School of Business Management, Stony Brook University, Stony Brook, NY, 11794, United States
highlights • Path integral approach to the statistical distributions of option pricing. • Analytical expressions of statistical distribution of bond and bond option pricing. • Statistical fluctuations and fatty tails of bond and bond option pricing.
article
info
Article history: Received 15 May 2016 Received in revised form 30 September 2016 Available online 28 December 2016 Keywords: Option pricing Distribution Risk
abstract Derivative (i.e. option) pricing is essential for modern financial instrumentations. Despite of the previous efforts, the exact analytical forms of the derivative pricing distributions are still challenging to obtain. In this study, we established a quantitative framework using path integrals to obtain the exact analytical solutions of the statistical distribution for bond and bond option pricing for the Vasicek model. We discuss the importance of statistical fluctuations away from the expected option pricing characterized by the distribution tail and their associations to value at risk (VaR). The framework established here is general and can be applied to other financial derivatives for quantifying the underlying statistical distributions. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Quantitative finance started from the study of stochastic processes. The foundation for understanding the stochastic processes was laid down in the last century thorough Brownian motion [1] by Bachelier, Einstein and Wiener [2–4]. In 1951, Ito established the stochastic calculus for Brownian motion [5]. The stock price fluctuations in time on the financial market can be described by a stochastic process. Later in 1950s, Samuelson laid down the foundation of option pricing theories via expectations by quantifying both macro- and microeconomics through mathematics [6]. In 1973, Black, Scholes and Merton derived the Black–Scholes partial differential equations for option pricing [7,8].
∗
Corresponding author at: Department of Chemistry and Physics, Stony Brook University, Stony Brook, NY, 11794, United States. E-mail address: [email protected] (J. Wang).
1 Kun Zhang and Jing Liu contributed equally to this work. http://dx.doi.org/10.1016/j.physa.2016.12.044 0378-4371/© 2016 Elsevier B.V. All rights reserved.
758
K. Zhang et al. / Physica A 471 (2017) 757–766
In 1977, Boyle related the pricing of options to the simulations of random asset paths, assessing the important roles of Monte Carlo simulations in finance [9]. More recently, in 1994, Derman and Kani [10], Dupire [11] and Rubinstein [12] considered the stochastic nature of the volatility and generalized the Black–Scholes equation in this situation. In 1977, Vasicek investigated the pricing interest rate products such as bonds. He modeled a short-term interest rate as a random walk and found that the interest rate derivatives can be valued using equations similar to the Black–Scholes equation [13]. In 1986, Ho and Lee found a way of improving Vasicek model with the idea of yield curve fitting or calibration [14]. In 1992, Heath, Jarrow and Morton modeled the random motion of the whole yield curve and generalized Ho and Lee’s approach [15]. In addition to the expected pricing discussed, quantifying the statistical fluctuations of equity and bond products as well as their derivatives is critical for addressing the associated risks. Basel committee on banking supervision published on international convergence of capital measurement and capital standards in 2006 [16]. Therefore, banks are now required to assess the risk quantitatively. Moreover, as one can see clearly from the recent financial crisis started from Wall Street and propagated through the whole world, the quantitative risk assessment is crucial for financial stability. Progresses have been made on the equity (stock) pricing models. Furthermore, analytical models of option pricing for both equity and credit so far have often been focused on the expected pricing. The distributions of option pricing revealing high order statistical fluctuations critical for risk assessments, are still challenging to obtain analytically [17–23] despite of numerical efforts such as Monte Carlo approach [9]. In this study, we established a quantitative framework using path integrals to obtain an exact analytical formula for the statistical distribution of the bond and bond option pricing for the Vasicek model. We discuss the importance of statistical fluctuations away from the expected option pricing characterized by the tail of the distribution and the associations to value at risk (VaR) critical in the practice of risk control. Quantifying statistical distributions is also very important in characterizing complex physical and biological systems as well as single molecules dynamics [24–26]. 2. Model In order to explore the statistical distribution of the option pricing, it is important to start with a suitable method. In this study, we use a path integral method [24–33]. This method has its advantages that the expected pricing and associated distribution can all be quantitatively addressed in a relatively straight forward way. We start with a bond model. Bond is a contract, paid for up-front, that yields a known amount on a known date in the future, the maturity date T . Bonds maybe issued by either governments or companies. The main purpose of a bond issue is the raising of capital, and the up-front premium can be thought of as a loan to the government or the company. The issue related to valuing bond pricing lies on how much one should pay to get a guaranteed amount in certain times. In regarding to the bond options, one aims to find fair values of the contract. In this section, we will discuss the path integral formulation and the expected pricing of bond and bond options. In the results and discussion section, we will concentrate on the distribution of the derivative pricing of bond and bond options. 2.1. The Vasicek model of bond with fluctuating interest rate The interest rate for bonds is not a constant and stochastically fluctuates. The Vasicek model for bonds [13] is usually studied from a stochastic differential equation for the short term interest rate r [13,34,35]: dr dt
= a(b − r (t )) + σ f (t ).
(1)
The rate dynamics (change) is controlled by two terms: the drift term (first term) and the stochastic term (second term). The parameters a and b model the mean reversion, σ is the volatility of the short term interest rate. a, b and σ are constants in time. The whole dynamics is parameterized in terms of the continuous time variable t. Looking at the drift term, the interest rate will tend to relax (decrease) towards the mean for large r and move up on average for small r. The ‘‘force’’ f (t ) represents the stochastic fluctuations:
⟨f (t )f (t )⟩ = δ(t − t ),
⟨f (t )⟩ = 0.
(2)
The brackets ⟨⟩ indicate a mean value with respect to a Gaussian distribution of zero mean and variance of one. The interest rate term structure models (including Vasicek model) are for pricing risk-free government bonds. For corporate bonds, additional credit spread and credit risk model are needed. 2.2. The functional for the short term interest rate From Gaussian noise properties of the stochastic force, we can write down the distribution of ‘‘force’’ f as: 1 p({f (s)}) = e− 2
dsf 2 (s)
.
(3)
K. Zhang et al. / Physica A 471 (2017) 757–766
759
What we are interested in is the distribution of the interest rate r. One can substitute stochastic force f from Eq. (1) and taking care of the variable transformation (Jacobian) to obtain the probability density functional for the short rate r (s) [32]: p({r (s)}) = Ne
−
1 2σ 2
2 ds dr −a(b−r (s)) ds
.
(4)
We will be only interested in the family of short rates within a specific time window, e.g. between t and T . Therefore, the integration on the right hand side. of Eq. (4) has to be constrained to the interval[t , T ]. 3. Statistical methods 3.1. The expected bond pricing as a path integral Now we discuss the bond price as a specific example for a claim. The bond price today, i.e. at time t, for a monetary unit promised at some later time T can be stated as an expectation value of the discount function (from the maturity time T to the current time t). Let P (t , T ) = e−
T
r (s)ds
t
: When the maturity time T is reached, then the bond price is the one promised (here it is one for
convenience). But at current time t less than T , the promised amount will be discounted by a factor of e− words, the bond is worth less at present before the maturity time.
T ⟨P (t , T )⟩ = EQ e− t
dsr (s)
|r (t ) = r .
T t
r (s)ds
. In other
(5)
EQ is the expectation value. The classical Vasicek approach (as described in [33,34]) has been to write down the corresponding PDE for P (t , T ) rather than evaluating the expectation value in Eq. (5) directly. Using the probability density functional in Eq. (4), we are now in the position to evaluate the expectation value in Eq. (4) for P (t , T ). The result is well known. It is explicitly given as follows, see details in Appendix A, [28–32]:
⟨P ⟩ = A(t , T ) exp(−B(t , T )r )
(6)
with
A(t , T ) = exp (−B(t , T ) + T − t )
σ2 2a2
σ2 2 −b − B (t , T ) 4a
(7)
and B(t , T ) =
1 − e−a(T −t ) a
.
(8)
3.2. The expected bond option pricing For an example of bond option, we look at the European call. A European call option is the call option at time 0 with a strike price k, expiring at time t on a bond, with maturity time T with T > t. In other words, for a European call option of the bound, if the bond price P is less than the striking price k at time t, then one gets nothing. On the other hand, if the bond price P is larger than the striking price k at time t, then one gets the difference between the bond price P and the striking price k. The expected price of a European call at time 0 with strike k expiring at time t on zero bond with maturity time T , with T > t is given by,
t ⟨c ⟩ = E e− 0
r (s)ds
max(P (t , T ) − k; 0)|r (0) = r
where P (t , T ) = A(t , T )e−B(t ,T )r . The constraint on P-k is transformed to r: r (t ) <
1 B(t , T )
ln
A(t , T ) k
= K.
Thus
t ⟨c ⟩ = E e− 0
r (s)ds+ln(A(t ,T )e−B(t ,T )r (t ) −k)
|r (0) = r .
760
K. Zhang et al. / Physica A 471 (2017) 757–766
We can substitute x(s) = b − r (s), and the constraint translates into x(t ) > b − K . The final analytical result is known and given [28–32] as follows:
⟨c ⟩ = c (0, T )N (h) − kc (0, t )N (h − σp )
(9)
where N (x) is the cumulative normal distribution and h=
1
σp
c (0, T )
+
kc (0, t )
σp
(10)
2
and
σp =
σ a
(1 − e
−a(T −t )
1 − e−2at
)
(11)
2a
where
c (0, T ) = EQ e−
T t
dsr (s)
c (t , T )|r (0) = r
(12)
c (t , T ) = A(t , T )e−B(t ,T )r
(13)
and
c (0, t ) = EQ e−
T t
dsr (s)
|r (0) = r .
(14)
See details in Appendix A [28–32]. Here using a path integral formulation, the formula derived earlier using partial differential equation for expected bond pricing with Vasicek model can be recovered [32,34–36]. The analytical tractability of the bond option formulas for Vasicek dynamics is dependent on the simple form of the equation for the bond price itself. 3.3. The derivative price distribution from path integral method The advantage of the path integral method is that we can obtain not only the expectation value but also directly the statistical distribution of the option pricing. We will show below the examples of bond pricing and bond option pricing to illustrate the power of path integral method in quantifying and obtaining the analytical formula for the statistical distribution of option pricing. 3.4. The bond price distribution from path integral method Since from previous section we know how to compute the expectation value of the bond price P, we now turn to the question how we will obtain the statistical distribution of the bond price P. T
As discussed, P (t , T ) = e− t r (s)ds is a discount function as discussed and its expectation values gives the bond price. The interest rate, however, follows stochastic fluctuations characterized by the probability density functional p({r (s)}) given by T Eq. (4). If we know the statistical distribution of z = − t r (s)ds since P (t , T ) = exp[z ], we will know the distribution of
T
P (t , T ). The distribution of z can be calculated by f (z ) = ⟨δ[z − (− t r (s)ds)]⟩. This gives the density or probability of the z. The expectation is over the realizations of r, p({r (s)}). The physical meaning is clear. The density or probability of z is obtained by counting the events or peaks of z characterized by the delta function at locations with different realizations T of r. Fortunately, the delta function can be transformed into the exponential in Fourier space k, exp[ik(z + t r (s)ds)]. This allows us to map the current problem to essentially almost the same problem of expected pricing done before with path T
integrals. Notice that the expected bond price is obtained by expectation of the discount function P (t , T ) = e− t r (s)ds over T p({r (s)}) whereas now for distribution we are performing the expectation of exp[ik(z + t r (s)ds)] over p({r (s)}) with an additional Fourier integral k to carry out. Therefore, the probability density function of z and therefore P can be obtained. We give some derivation steps to follow in the result section (for those who care about the detailed procedures, see details in Appendix A, for those who do not care about intermediate steps, jump directly to Eq. (18) for the results). 3.5. The bond option price distribution from path integral method We now turn to the question how we will obtain the statistical distribution of the bond option price P. For a European call option of the bond with maturity time T , if bond price P is larger than the striking price k at time t, then one gets the difference between the bond price P and the striking price k. Otherwise, one gets nothing. This provides a constraint in the expectation over all the possible realizations or interest rate r. This is the main difference between the bond pricing and bond option pricing. In other words, the bond option price c (before the expectation) follows from a discount
K. Zhang et al. / Physica A 471 (2017) 757–766
function as discussed e−
t
0 r (s)ds
761
relaxing from striking time t to current time (t = 0) with a constraint P larger than k at t
striking time (option price (before expectation) c = e− 0 r (s)ds max[P (t , T ) − k; 0]). While the interest rate follows a stochastic fluctuations characterized by the probability density functional p({r (s)}) given t by Eq. (4). If we know the distribution of z = ln(P (t , T ) − k) − 0 r (s)ds, we will know the distribution of c since c = exp[z ].
t
The distribution of z can be calculated by f (z ) = δ[z − (ln(P (t , T ) − k) − 0 r (s)ds)], This gives the density or probability of the z. The expectation is over the realizations of r, p({r (s)}). The physical meaning is clear. The density or probability of z is obtained by counting the events or peaks of z characterized by the delta function at locations with different realizations of r. Furthermore, the delta function can be transformed into the exponential in Fourier space g, exp[ig {z − (ln(P (t , T ) − t k) − 0 r (s)ds)}]. This allows us to map the current problem to essentially almost the same problem of expected option pricing done before with path integrals. Notice that the expected option price is obtained by performing expectations t
over the discount function c = e− 0 r (s)ds max[P (t , T ) − k; 0] (with the constraint P larger than k) over the realizations of interest rate with probability density p({r (s)}) whereas now for distribution we are performing the expectation values t of exp[ig {z − (ln(P (t , T ) − k) − 0 r (s)ds)}] over p({r (s)}) (with a constraint P larger than k) with an additional Fourier integral g to carry out. Therefore, the probability density function of z and therefore option price c can be obtained. We give some derivation steps to follow in the result section (for those who care about the detailed procedures, see details in Appendix A. For those who do not care about intermediate steps, jump directly to Eq. (26) for the results). 4. Results and discussion 4.1. The results and derivations for bond price distribution from path integral method
We are now trying to find out the distribution of P (t , T ), that is f (P ). First we define z as, z = − and z is related by P = ez ⇔ z = ln P ⇔ dz =
1
T t
r (s)ds we see that P
dP .
P
We can also see the relationship between the distribution of z, f (z ) and the distribution of P, f (P ) : f (P )dP = f (z )dz f (P ) = f (z ) f (P ) =
1 P
dz dP
f (z ).
Thus, to find out f (P ), we only need to find out f (z ).
f (z ) = δ
T
r (s)ds
z+
(15)
t
and we can write down the path integral representation for the distribution of z
f (z ) =
1 2π
∞ −∞
∞
dk −∞ dr (T )
r (T ) r (t )
∞
−∞ dr (T )
Dr (s) exp − 2σ 2
r (T )
1
T t
ds
Dr (s) exp − 2σ1 2 r (t )
dr (s) ds
T t
ds
− a(b − r (s))
dr (s) ds
2
+ ik z + 2 − a(b − r (s))
T t
dsr (s) .
(16)
+∞
Here, we use δ(x) = 21π −∞ eikx dk to calculate f (z ). We see that it is similar to the formula used for expected value of bond pricing except with one addition integration over k. −a(T −t )
Define, B(t , T ) = 1−e a we perform the integrals with respect to r (T ) and path integrals with respect to r (s) by applying the generating function for the harmonic oscillator [27] (details in Appendix A). We obtain: f (z ) =
X Y
=
1 2π
∞
dk exp −
−∞
σ2
σ2
2a
4a
(−B(t , T ) + (T − t ))k2 + 2
k2 B2 (t , T ) − ikB(t , T )x(t ) + ikb(T − t ) + ikz
where x(t ) = b − r (t ). Now, we only have the Gaussian integration with respect to k for f (z ) to carry out:
f (z ) =
1
π W (t , T )
exp
− ((−B(t , T ) + T − t )b + B(t , T )r + z )2 W (t , T )
2 2 where W (t , T ) = 2aσ2 (−B(t , T ) + T − t ) − σa B2 (t , T ).
Since f (P ) =
1 f P
(z ) and z = ln P.
(17)
762
K. Zhang et al. / Physica A 471 (2017) 757–766
Fig. 1. Statistical distribution of bond price with respect to (a) interest reversion rate a, (b = 0.01, r = 0.01, σ = 0.03, t = 10, T = 30); (b) interest reversion b, (a = 0.3, r = 0.01, σ = 0.03, t = 10, T = 30); (c) interest rate r, (a = 0.3, b = 0.01, σ = 0.03, t = 10, T = 30); (d) volatility, (a = 0.3, b = 0.01, r = 0.01, t = 10, T = 30); (e) current time t, (a = 0.3, b = 0.01, r = 0.01, σ = 0.05, T = 30); (f) maturity time T . (a = 0.3, b = 0.01, r = 0.01, σ = 0.05, t = 10).
Therefore, we finally obtain the distribution of bond pricing. f (P ) =
1 P
− ((−B(t , T ) + T − t )b + B(t , T )r + ln P )2 exp π W (t , T ) W (t , T ) 1
.
(18)
Eq. (18) gives the analytical solution of the distribution for option pricing. A(t , T ) and B(t , T ) are given by Eqs. (7) and (8). As we can see, the bond price follows a log normal distribution, with the average and width W (t , T ) controlled by the interest rate r at time t and associated mean reversion a and b, interest rate fluctuation strength or spread, σ current time t, and maturity time T . We plot the distribution of f (P ) with respect to mean reversion a, b, the current interest rate r, volatility, σ current time t and maturity time T in Fig. 1. We found that as a, b, r and t increase, the tails of the distribution become less prominent. On the other hand, as σ and T increase, the tails of the distribution become more prominent. This implies the more volatility and longer maturity time, the more statistical fluctuations of the bond pricing is. On the other hand, the closer of current time t to the maturity time T and larger the interest rate and associated mean reversion, the bond pricing has less fluctuations and is more determined. This provides the quantitative origin of the bond pricing fluctuations and important for quantitative risk control from the statistical distribution perspective. 4.2. The results and derivations for bond option price distribution from path integral method Let us consider the option price distribution of a European call at time 0 with strike k expiring at time t on zero bond with maturity time of T , with T > t. t
Let option pricing (before expectation) be c (0) = e− 0 r (s)ds max(P (t , T ) − k; 0) with r (0) = r, so that the option pricing A(t ,T ) is given by the discount function with the constraint written as: r (t ) < B(t1,T ) ln( k ) = K . We need to find out the distribution of option pricing c (0), that is f (c (0)). For convenience, we use f (c ) instead of f (c (0)) in the text. t First, we define variable z, so that, z = ln(P (t , T ) − k) − 0 r (s)ds then c and z are related, c = ez ⇔ z = ln c ⇔ dz =
1 c
dc .
K. Zhang et al. / Physica A 471 (2017) 757–766
763
In general, we can find the relationship between distribution of option pricing c, f (c ) and distribution of z, f (z ): f (c )dc = f (z )dz dz f (c ) = f (z ) dc 1 f (c ) = f (z ). c
+∞
To find out f (c ), we only need to find out f (z ). Since δ(x) = 21π −∞ eigx dg, we have
t r (s)ds δ z − ln(P (t , T ) − k) − 0 t r (s)ds δ z − ln(P (t , T ) − k) + 0 t ∞ 1 r (s)ds dg exp ig z − ln(P (t , T ) − k) + 2π −∞ 0 ∞ t 1 dg exp ig z − ln(P (t , T ) − k) + r (s)ds . 2π −∞ 0
f (z ) =
= = =
(19)
Therefore, the distribution of z, f (z ) is given in the path integral formulation as: 1 2π
f (z ) =
∞
K
−∞
dg −∞ dr (t )
r (t )
Dr (s) exp − 2σ1 2
r (0)
∞
−∞ dr (t )
t 0
ds
dr (s) ds
− a(b − r (s))
r (t )
Dr (s) exp − 2σ1 2 r (0)
t 0
ds
dr (s) ds
t + ig z − ln(P (t , T ) − k) + 0 r (s)ds . 2 − a(b − r (s))
2
(20) We can see that this formula is similar to the one used in the expectation of option pricing except the additional integration with respect to g. Substitute x(s) = b − r (s) and by applying the generating function for the harmonic oscillator [27] (see details in Appendix A), we obtain, f (z ) = XY , where X =
1
a
∞
∞
σ2 dx(t ) exp − 3 (e−at − 1 + at )g 2
dg 2π σ 2 sinh(at ) −∞ 2a b−K a a 2 − 2 (x2 (t ) cosh(at ) − 2x(t )x(0) + 2(eat − 1)(igC (x(t ) + x(0)) − g 2 C 2 )) − x (t ) 2σ sinh(at ) 2σ 2 2π
+ igz + igbt − ig ln(A(t , T )e−B(t ,T )(b−x(t )) − k)
(21)
and
Y =
a
∞
2π σ 2 sinh(at ) −∞
dx(t ) exp −
a 2σ
x2 (t ) − 2
a 2σ 2 sinh(at )
(x2 (t ) cosh(at ) − 2x(t )x(0))
(22)
−a(T −t )
with C = σ2a (1 − e−at ) and B(t , T ) = 1−e a . Now we perform the Gaussian integration with respect to g, we obtain: 2
X =
a 2π σ 2 sinh(at )
−
× exp
1
πV
∞
dx(t ) b−K
1−eat x(t )+x(0) a 1+eat
2 + bt + z − ln(A(t , T )e−B(t ,T )(b−x(t )) − k) V
−
a 2σ 2 sinh(at )
where V = − 2aσ3 ( 2
2(eat −1) eat +1
(x2 (t ) cosh(at ) − 2x(t )x(0)) −
− at ) =
2σ 2 a2
(2ah + t ), ah =
1 a
a 2σ 2
e−at −1 . e−at +1
x2 (t )
(23)
764
K. Zhang et al. / Physica A 471 (2017) 757–766
Also, we can perform the Gaussian integration with respect to x(t ). We obtain
Y =
a
π σ 2 (1 − e−2at )
2π σ 2 sinh(at )
a
ax2 (0)
exp −
σ 2 (1 − e2at )
.
(24)
Substitute back, r (t ) = b − x(t ), r (0) = r = b − x(0) f (z ) =
1
π
ap V
K
dr (t ) exp
−(−ah (r (t ) + r − 2b) + bt + z − ln(A(t , T )e−B(t ,T )r (t ) − k))2 V
−∞
− ap (b − r (t ) − e−at (b − r ))2
(25)
−at
1 where, ap = σ 2 (1−ae−2at ) , ah = 1a ee−at − . +1 Therefore we finally obtain the analytical form of the distribution of bond option pricing c as follows,
f (c ) =
1
πc
ap V
K
dr (t ) exp
−(−ah (r (t ) + r − 2b) + bt + ln c − ln(A(t , T )e−B(t ,T )r (t ) − k))2 V
−∞
− ap (b − r (t ) − e−at (b − r ))2 .
(26)
Eq. (26) gives the exact analytical solution of the distribution for option pricing (up to a one dimensional integral with respect to r (t )). A(t , T ) and B(t , T ) are given in Eqs. (7) and (8). We plot the distribution of option pricing f (c ) with respect to mean reversion a, b, current interest rate r (r (0)), volatility σ , striking time t and maturity time T in Fig. 2. We found that as a, b, r and t increase, the tails of the distribution become less prominent. On the other hand, as σ and T increase, the tails of the distribution become more prominent. This implies the more volatility and longer maturity time, the more statistical fluctuations of the option pricing is. On the other hand, the closer the striking time to the maturity time and the larger the interest rate and associated mean reversion, the option pricing has less statistical fluctuations and more determined values. This provides the quantitative origin of the bond option pricing fluctuations and important for quantitative risk control from the statistical distributionperspective. Once the statistical fluctuations are known, the underlying fluctuations are known and can be quantified. The value at risk (VaR) is often used to quantify the risk and is simply the variation or the second order moments which can be directly calculated from the distribution. The skew (3rd order moment) and kurtosis (4th order moment) are sometimes used for characterizing the asymmetry of the underlying observables. They are nothing but the third order and fourth order statistical fluctuations which can also be easily obtained from the distribution. For Gaussian distribution, the variable itself should not have skewness and kurtosis. Here the stochastic variable is the interest rate. However, the quantity we are interested in is the option pricing which is a nonlinear function of interest rate and no symmetry on option pricing is guaranteed. Therefore, the associated distribution of the option pricing is not expected to be Gaussian. The skewness and kurtosis of the option pricing are not expected to be zero for this reason. As one can see, statistical fluctuations can be completely characterized and quantified by the underlying distribution function. The statistical distribution can help to quantify the high order fluctuations. The high order fluctuations in option pricing can be obtained by integrating the same order of power for the option pricing over the corresponding statistical distribution. In principle once the distribution is known, the information on the higher order fluctuations such as VaR can be obtained. In practice, however, option prices are often used to calibrate the term structure model parameters (such as a and b in the Vasicek model), Fitting the mean option pricing is already a challenging. It may be even more challenging to obtain VaR directly. In strong fluctuation case, the expected or average pricing may not faithfully characterize the system since there are large statistical fluctuations around it. Instead, high order moments become more important. The distribution often has fatty tails. This implies the ‘‘intermittency’’ where rare event (low probability event at the tail) [24–26] can have significant effects on the nature of the system. It is exactly this rare event poses the high risk which we need to capture in order to prevent it to happen in finance. As we can see one cannot find such rare risk event by just looking at the expected values or the average. High order statistics through distribution is the only way to quantify the whole event. Such strongly fluctuating systems can be found not only in financial world as we have witnessed recently, but also every where else in physical and biological systems such as biomolecules, networks, weather, earth quakes, and evolution etc. [24–26,37–39].
K. Zhang et al. / Physica A 471 (2017) 757–766
765
Fig. 2. Statistical distribution of bond option price with respect to (a) interest reversion rate a, (b = 0.01, r = 0.01, σ = 0.2, t = 10, T = 30); (b) interest reversion b, (a = 0.35, r = 0.01, σ = 0.2, t = 10, T = 30); (c) interest rate r, (a = 0.35, b = 0.01, σ = 0.2, t = 10, T = 30); (d) volatility, (a = 0.35, b = 0.01, r = 0.01, t = 10, T = 30); (e) striking time t, (a = 0.35, b = 0.01, r = 0.01, σ = 0.2, T = 30); (f) maturity time T . (a = 0.35, b = 0.01, r = 0.01, σ = 0.2, t = 10).
5. Conclusions In this study, we obtained exact analytical solutions of the statistical distribution of the bond and bond option pricing for the Vasicek model with a path integral framework. This approach provides a quantitative way to obtain the underlying distribution of the derivative pricing which is critical for the financial risk control. We discuss the importance of statistical fluctuations away from the expected option pricing characterized by the tail of the distribution and the associations to value at risk (VaR). The framework established here is general and can be applied not only to financial derivatives but also to wide variety of problems in complex physical and biological systems, and single molecules for quantifying the underlying intrinsic statistical distributions. Acknowledgments We would like to thank Dr. Xidi Wang from Citibank, Dr. Jun Zang from Morgan-Chase Bank, Dr. Joy Zhang from Credit Suisse Bank, Peter G. Wolynes from Rice University and Prof. Raphael Douady from Stony Brook University/CNRS for interesting discussions. We would like to thank NSF (Grant No. NSF-PHY 76066) and NSFC (China) (Grant No. 91430217) for financial support. Appendix A. Supplementary data Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.physa.2016.12.044. References [1] R. Brown, A Brief Account of Microscopical Observations, London, (1827). [2] L. Bachelier, Theorie de la Speculation, Jacques Gabay Pub., Paris, 1995. [3] A. Einstein, On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat, Ann. Phys. 17 (1905) 549–560. [4] N. Wiener, The average of an analytic functional, Proc. Natl. Acad. Sci. 7 (1921) 253–260.
766 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
K. Zhang et al. / Physica A 471 (2017) 757–766 K. Ito, On stochastic differential equations, Mem. Amer. Math. Soc. 4 (1951) 1–51. P. Samuelson, Rational theory of warrant pricing, Ind. Manag. Rev. 6 (1965) 13–39. F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Polit. Econ. 81 (1973) 637–659. R.C. Merton, Theory of rational option pricing, Bell J. Econ. Manage. Sci. 4 (1973) 141–183. P. Boyle, Options: a Monte Carlo approach, J. Financ. Econ. 4 (1977) 323–338. E. Derman, I. Kani, Riding on a smile, Risk 7 (1994) 32–39. B. Dupire, Pricing with a smile, Risk 7 (1994) 18–20. M. Rubinstein, Implied binomial trees, J. Finance 69 (1994) 771–818. O.A. Vasicek, An equilibrium characterization of the term structure, J. Financ. Econ. 5 (1977) 177–188. T. Ho, S. Lee, Term structure movements and pricing interest rate contingent claims, J. Finance 42 (1986) 1129–1142. D. Heath, R. Jarrow, A. Morton, Bond Pricing and the term structure of interest rates: a new methodology, Econometrica 60 (1992) 77–105. Basel Committee on Banking Supervision (2006) International Convergence of Capital Measurement and Capital Standards: A Revised Framework Comprehensive Version Bank for International Settlements 1–347. P. Carr, M. Schroder, On the valuation of arithmetic-average Asian options: the Geman-Yor Laplace transform revisited. arXiv preprint math/0102080, 2001. M. Yor, On some exponential functionals of Brownian motion, Adv. Appl. Probab. 24 (1992) 509–531. R. Douady, Closed form formulas for exotic options and their lifetime distribution, Int. J. Theor. Appl. Finance 2 (01) (1999) 17–42. H. Geman, N. El Karoui, J.C. Rochet, Changes of numeraire, changes of probability measure and option pricing, J. Appl. Probab. (1995) 443–458. M. Jeanblanc, M. Yor, M. Chesney, Mathematical Methods for Financial Markets, Springer Science & Business Media, 2009. D. Brigo, F. Mercurio, Interest Rate Models: Theory and Practice, Springer Finance, Berlin, Heidelberg, New York, 2001. T. Bielecki, M. Rutkowski, Credit Risk: Modeling, Valuation and Hedging, Spinger Finance, Berlin, Heidelberg, New York, 2002. J. Wang, P.G. Wolynes, Intermittency of single molecule reaction dynamics in fluctuating environments, Phys. Rev. Lett. 74 (1995) 4317. S.S. Plotkin, P.G. Wolynes, Buffed energy landscapes: Another solution to the kinetic paradoxes of protein folding, Proc. Natl. Acad. Sci. 100 (8) (2003) 4417–4422. L. Mahadevani, J.M. Deutch, Influence of feedback on the stochastic evolution of simple climate systems, Proc. R. Soc. A 466 (2010) 993–1003. R.P. Feynman, A.R. Hibbs, Quantum Mechanics and Path Integral, McGraw-Hill Book Company, New York, 1965. B.E. Baaquie, Quantum Finance, Cambridge University Press, 2004. V. Linetsky, The path integral approach to financial modeling and options pricing, Comput. Econ. 11 (1998) 129. J. Dash, Path Integrals and Options. Part I, CNRS Preprint CPT88/PE.2206, 1988. J. Dash, Path Integrals and Options. Part II, CNRS Preprint CPT89/PE.2333, 1989. M. Otto, Using path integrals to price interest rate derivatives, 1999. arXiv:cond-mat/9812318v2. H.D. Feng, K. Zhang, J. Wang, Non-equilibrium transition state rate theory, Chem. Sci. (2014) (5) 3761–3769. P. Wilmott, Paul Wilmott on Quantitative Finance. I and II, John Wiley & Sons, New York, 2000. R. Rebonato, Interest-Rate Option Models, John Wiley & Sons, Chichester, 1996. F. Jamshidian, An exact bond option pricing formula, J. Finance 44 (1989) 205. M. Sasai, P.G. Wolynes, Stochastic gene expression as a many-body problem, Proc. Natl. Acad. Sci. 100 (2003) 2374–2379. M. Bottiglieri, C. Godano, On-offintermittency in earthquake occurrence, Phys. Rev. E 75 (2007) 026101–5. S.J. Gould, N. Eldredge, Punctuated equilibria: the tempo and mode of evolution reconsidered, Paleobiology 3 (1977) 115–151.