key: cord-0148735-e5j2ay0g
authors: Hu, Yuan; Shirvani, Abootaleb; Lindquist, W. Brent; Fabozzi, Frank J.; Rachev, Svetlozar T.
title: Market Complete Option Valuation using a Jarrow-Rudd Pricing Tree with Skewness and Kurtosis
date: 2021-06-16
journal: nan
DOI: nan
sha: b29efdd719734792d624b4ebeb12f05c7c32b6dc
doc_id: 148735
cord_uid: e5j2ay0g

Applying the Cherny-Shiryaev-Yor invariance principle, we introduce a generalized Jarrow-Rudd (GJR) option pricing model with uncertainty driven by a skew random walk. The GJR pricing tree exhibits skewness and kurtosis in both the natural and risk-neutral world. We construct implied surfaces for the parameters determining the GJR tree. Motivated by Merton's pricing tree incorporating transaction costs, we extend the GJR pricing model to include a hedging cost. We demonstrate ways to fit the GJR pricing model to a market driver that influences the price dynamics of the underlying asset. We supplement our findings with numerical examples.

Pricing trees determined by a Markov chain have been studied in a number of academic papers. 1 Since the Markov chain pricing tree is constructed directly in the risk-neutral world in these papers, it is not clear what discrete pricing model in the natural world would evolve to a risk-neutral Markov chain pricing model in accordance with dynamic asset pricing theory. 2 In this work, we follow the classical binomial pricing model framework (see Cox et al. (1979) and Jarrow and Rudd (2012, Chapters V and VI) ). We begin with a Markov chain model in the natural world replicating a self-financing portfolio, and then transform to risk-neutral option pricing preserving market completeness. It is essential that our option pricing model is defined first in the natural world. If the option pricing model is placed directly in the risk-neutral world, and calibrated with market option data, no option market dislocation or option mispricing can be revealed. Option pricing models which do not start with modeling the underlying assets in the natural world and then, via risk-neutral valuation based on Black-Scholes-Merton dynamic pricing theory, pass to the risk-neutral world are meaningless, if not dangerous, in practical applications. Such pricing models are generally used to predict what the option traders believe the correct option prices are and not what option prices are actually fair. It is imperative that opinion prices be aligned with reliable spot market models. 3 To illustrate this claim, suppose that a trader would like to determine if potential mispricing in the option market exists due to a market bubble that the trader suspects will burst. 4 The trader's option pricing model, when calibrated to option market data, should recover, uniquely, the spot price dynamics (Ross, 2015) of the underlying asset (Kim et al. (2016) , Hu et al. (2020a) ). If that recovered spot price process does not conform with market data on the asset spot price, this could be a trading signal that there is potential option mispricing. 5 Our paper is close in spirit to Kijima and Yoshida (1993) , where an option pricing model with Markov chain stochastic volatility is introduced and the continuous-time limiting price process is determined as a subordinated geometric Brownian motion (GBM). 6 While the discrete-and continuous-time market models are incomplete in the paper by Kijima and Yoshida, in our paper we deal with complete market models, both in the discrete-and continuous-time settings. We extend the Jarrow and Rudd (JR) binomial pricing model (Jarrow and Rudd (2012) , Hull (2012, p. 442) , Kim et al. (2016 Kim et al. ( , 2019 ) with an additional parameter determining the skewness and excess kurtosis of the underlying asset return distribution. We refer to our extended model as the generalized Jarrow-Rudd (GJR) option pricing model.

As in the original JR model (Kim et al., 2016 (Kim et al., , 2019 Hu et al., 2020a,b) , our GJR pricing model preserves market completeness. 7 In the GJR pricing model, the embedded Markov chain driving the discrete underlying price process is a skew random walk which, in the limit, becomes a skew Brownian motion (SBM) after the necessary scalar normalization. 8 In discrete time, the distributional mapping between the GJR price dynamics of the underlying asset and its risk-neutral dynamics is one-to-one. 9 The GJR pricing model 1 See, for example, Duan and Simonato (2001) , Simonato (2011) , Bhat and Kumar (2012) , Fuh et al. (2012) , Van and Elliot (2012a,b) , and Fan et al. (2016) . In Bhat and Kumar (2012) , the authors claim that the limiting return distribution is a mixture of normal distributions. A mixture of different normal distributions is not infinitely divisible (Steutel and van Harn, 2004, Chapter VI) , and it is therefore not clear that their limiting continuous-time model is arbitrage-free. 2 Black and Scholes (1973) , Merton (1973) , Delbaen and Schachermayer (1994) , Delbaen and Schachermayer (1998) , and Duffie (2001, Chapter 6) .

3 Black (1975) , Brenner and Galai (1984) , Melick and Thomas (1997) , Hilber et al. (2009) , Pasquariello (2014) , and Ross (2015) . 4 The trader can decide to ride the bubble as long as possible, as many other option traders will do, with the hope of unwinding option positions before the bubble bursts. In this case, the trader does not need to worry about what the true spot market model is that corresponds to the risk-neutral option pricing model. For an extensive study of market bubbles and related option markets, we refer to Heston et al. (2007) , Jarrow et al. (2010) , and Vogel (2018) .

5 In this setting model risk does exist. That is why in real trading a suite of option pricing models is used. 6 We employ the following abbreviations throughout this paper. Each abbreviation is also defined the first time it is referenced. BM: Brownian motion; CPM: continuous-time option pricing model; CSYIP: Cherny-Shiryaev-Yor Invariance Principle; DPM: discrete-time option pricing model; ECC: European contingency claim; GBM: geometric Brownian motion; GJR: generalized Jarrow-Rudd; HTC: hedging transaction cost; JR: Jarrow-Rudd; RMSE: relative mean-square error; SBM: skew Brownian motion; w.p.: with probability.

7 As far the authors of this paper are aware, all discrete market models in the literature exhibiting skewness and kurtosis lead to incomplete market models. Thus, the problem of determining the unique spot price tree corresponding to the risk-neutral tree chosen in any of those papers will be unsolved.

8 See Harrison and Shepp (1981) , Cherny et al. (2003) , Revuz and Yor (1994, Chapters VII and X) , and Corns and Satchell (2007) . 9 That is, for every fixed trading frequency, the probability law of the risk-neutral tree for the underlying asset uniquely allows us to study option pricing in a Markov chain market model with transaction costs. Our market model has a relatively simple parametric form and is easily calibrated to option data, as illustrated by numerical examples. Our paper proceeds as follows. In section 2 we introduce the Cherny-Shiryaev-Yor invariance principle (CSYIP) (Cherny et al., 2003) for SBM. As with the use of the Donsker-Prokhorov invariance principle 10 in the Cox-Ross-Rubinstein pricing tree (Cox et al., 1979) and in the JR pricing tree, we apply CSYIP to introduce the GJR pricing tree and obtain the limiting continuous-time price dynamics. The mean return µ is retained as a parameter in the GJR model. The GJR tree includes an additional parameter β governing the skewness and kurtosis of the GJR tree in the natural world. While the GJR tree ultimately converges to a GBM as in the classical JR pricing tree, the GJR pre-limiting behavior is that of geometric SBM. We investigate numerically the pre-limiting behavior of the GJR tree to obtain estimates of the smallest number, n, of trading intervals to maturity T = n∆t, n ↑ ∞, at which the skewness and excess kurtosis vanish from the GJR tree distribution. In section 3, we determine the risk-neutral probabilities in the GJR pricing model which, together with the volatility σ and the risk-free rate r f , depend on µ and β. Using daily closing price data for the SPDR S&P 500 ETF Trust fund (SPY) and call option data for the same underlying asset, we estimate the implied µ, β and σ surfaces. Motivated by Merton's binomial tree model with transaction costs (Merton, 1990, Chapter 14) , in section 4 we extend the GJR pricing tree model to include a hedging transaction cost term. Using the SPY option data, we estimate option transaction costs, the implied transaction cost surface, and the impact transaction costs have on the implied µ, β and σ surfaces. In section 5, we explore possibilities for fitting the GJR model to a market driver that affects the price dynamics of the underlying asset. We estimate the parameters for the new pricing tree model in a numerical example with the underlying asset being the stock of Microsoft Corporation (MSFT). We explore both an endogenous and exogenous approach. The exogenous approach allows for greater generalization of the market driver; in a numerical example we assume the asset returns are dependent on Fama-French five-factor loading values (Fama and French, 2015) . In section 6, we derive the risk-neutral probabilities for the extended pricing tree introduced in section 5 and estimate the corresponding implied volatility surface. Chapter 7 concludes of our work.

In this section, we introduce the GJR pricing tree model. To do so, we describe SBM and then apply the CSYIP to formulate a new path-dependent pricing model defined by a recombined tree that generalizes the JR binomial option pricing model. 11 In contrast to the classical Cox-Ross-Rubinstein and JR binomial pricing models, in the GJR pricing model the pre-limiting pricing process is determined by SBM. That is, for a moderately small trading frequency ∆t, the GJR dynamics exhibit the properties of geometric SBM. This feature leads to a more flexible discrete pricing model for the underlying asset price behavior for any realistically small trading frequency. However, as ∆t ↓ 0 our pricing tree converges to a GBM, as in the traditional JR pricing model. As shown in Hu et al. (2020a,b) , discrete-time option pricing models contain considerably more information than continuous-time pricing models.

Due to the option trader's presumed ability to trade continuously with no transaction costs, in a continuous-time option pricing model the information about the underlying stock mean return and stockprice direction at a given trading frequency is lost. That is the main reason that we emphasize a discrete-time option pricing model, rather than paying more attention to the limiting continuous-time option price dynamics. In the GJR option pricing model, new features, such as skewness and excess kurtosis 12 , of the underlying stock return distribution will be present.

determines the probability law of the asset's pricing tree in the natural world. 10 Donsker (1951) , Prokhorov (1956) , Billingsley (1999, section 14) , Gikhman and Skorokhod (1969, Chapter IX) , Skorokhod (2005, section 5.3 .3), Davydov and Rotar (2008) 11 See Jarrow and Rudd (2012); Hull (2012, p. 442); Kim et al. (2016 Kim et al. ( , 2019 . 12 Discrete-time option pricing models with underlying stock return distributions exhibiting skewness and excess kurtosis are known; see, for example, Yamada and Primbs (2004) . But the Yamada-Primbs' pricing model is based on mean-variance hedging in incomplete markets. Our GJR pricing model is based on no-arbitrage asset valuation arguments leading to a complete market model.

We start with the definition of a SBM. 13 Let B = {B t , t ≥ 0} be a standard Brownian motion (BM) generating a stochastic basis (Ω, F = {F t = σ(B u , u ≤ t) t≥0 }, P). Let α ∈ [0, 1] and set 14

= 0 is a SBM with parameter α having the following properties. 15

(v) B (α) is a semimartingale satisfying the strong Markov property.

(vi) When t ≥ 0, sample paths of B (α) can be generated using the representation 17

is the conditional density f

(viii) Corns and Satchell (2007) developed a continuous-time option pricing model based on geometric Az-

where B 1,t and B 2,t are two independent BMs. Corns and Satchell (2007) showed that A (δ) d = B (α) with α = (1 + δ)/2.

13 See Itô and McKean (1996, section 4.2, problem 1, p. 115 ); Harrison and Shepp (1981) ; Revuz and Yor (1994, Chapters VII and X ) ; Lang (1995) ; Lejay (2006) ; Corns and Satchell (2007) ; Cherny et al. (2003) ; Ramirez (2011) ; Atar and Budhiraja (2015) ; Trutnau et al. (2015) ; and Li (2019) .

14 See Cherny et al. (2003) . 15 See Cherny et al. (2003, Chapter 4) and Corns and Satchell (2007) . 16 Here, and in what follows, d = stands for "equal in distribution" or "equal in probability law". 17 See Itô and McKean (1996, section 4.2, problem 1, p. 115) , Lejay (2006) and Corns and Satchell (2007) . We denote "with probability" as "w.p." (ix) The moment-generating function of B (α) t has the form

Note that odd moments depend on α through the term 2α − 1, while even moments are independent of α.

From (2), the mean, variance, skewness and excess kurtosis of B (α) are given by

Note that the skewness and excess kurtosis are time-independent quantities. The behavior of these four quantities is demonstrated in Fig. 1 . 

We start with the formulation of the CSYIP for SBM with the definition of a piecewise continuous function. PSC(i) ∪ ∞ n=1 J n = R; PSC(ii) for every compact interval J there exists n ∈ N such that ∪ n k=1 J k ⊇ J; and PSC(iii) on each J n , n ∈ N , h : J n → R is continuous and has finite limits at those endpoints of J n which do not belong to J n .

We call M (α) a skew random walk with parameter α. 

, t ≥ 0, to be the random process with piecewise linear trajectories having vertexes

, t ≥ 0, to be a random process with the piecewise linear trajectories having vertexes k/n, C

19 See Harrison and Shepp (1981) and Cherny et al. (2003) . 20 See Theorem 4.1 in Cherny et al. (2003) .

We begin the construction of a binomial pricing tree by applying CSYIP using only the lower moment SBM process B (α,n) t . We will add the higher moment C (α,n) t in section 5.2. We fix n ∈ N , the time interval ∆t = T /n, and β ∈ −1/ √ ∆t, 1/ √ ∆t . Set

If we consider the process

then D (α∆t) weakly converges in the Skorokhod space D[0, T ] to the BM B as ∆t ↓ 0. However, for fixed ∆t, D (α∆t) exhibits the properties of a SBM with parameter α ∆t . From (8), (7) and Fig. 2 , we see that β will affect the moments of the process B 

with mean E B

We next define the GJR pricing tree, determined by the skew random walk M (α∆t) . 22 GJR Pricing Tree: Let µ > 0 and σ > 0. For k = 1, . . . , n, n∆t = T , α ∆t = (1 + β √ ∆t)/2 and skew random walk M (α∆t) , define the GJR pricing recombined tree by

where v k = kµ + σβ 2k/π − 1 . We study the limiting behavior of this tree as ∆t ↓ 0. Note that, to leading order in ∆t,

Consider the cumulative log-return R

k∆t /S 0 , k = 1, . . . , n. From (10) and (11) we obtain

From (12) and the CSYIP it follows that, for a fixed but relatively small time-increment ∆t, the pricing tree (10) approximates

However, if we let ∆t ↓ 0, then the D[0, T ]-process generated by the tree (10) will ultimately converge to a GBM, that is, S t = S 0 exp(µt + σB t ), t ∈ [0, T ].

21 As ∆t = 1/n, we can assume that the value of E B

22 If α ∆t = 1/2 (i.e. β = 0), we obtain the JR binomial tree. Jarrow and Rudd (2012) and Hull (2012, p. 442 ) defined the JR binomial tree directly in the risk-neutral world, that is, when µ = r f , where r f is the risk-free rate. Kim et al. (2016) showed that JR binomial tree option pricing model provides the fastest rate of convergence to the corresponding GBM in the risk-neutral world.

As a numerical example to investigate the pre-limiting behavior of the GJR pricing tree, we use SPY daily closing prices. Price data covering N = 7136 trading days over the period 1/29/1993 to 6/1/2021 was obtained from Bloomberg Professional Services. We denote SPY daily closing prices as S (SPY) k∆t and cumulative log-returns as R

, where k = 1, ..., N, and ∆t = 1/252. We first determine a length, L, for a moving-window estimator that provides relatively strong explanatory power and acceptable precision for computing the parameters µ, σ and β, by using the following procedure based on equations (12). For each tested length L and for each moving window w i , i = 1, . . . , N − L + 1, repeat the following steps.

Step 1: Fit a robust linear regression using a logistic weight function 23 to the SPY return data in w i using the model 24 ln R

to produce estimates for the value ofσ (α,wi) and the error terms e

(1,wi) k∆t , k = 1, ..., L.

Step 2: Withσ (α,wi) estimated from Step 1, apply the conditional least squares optimization,

to determineβ (α,wi) ,μ (α,wi) and the second error sequence e (2,wi) k∆t , k = 1, ..., L.

Step 3: Define to be the random variable having the sample 23 See Holland and Welsch (2007) and Pregibon (1981) . 24 Since the cumulative return is not stationary, we employ the model Var R The results clearly show a sharp drop inμ (α) and increased volatility (σ (α) ) resulting from the 2008 global financial crisis. The most recent peak inσ (α) in 2020 corresponds to the market reaction to the Covid-19 pandemic. In Fig. 4c , we note that optimum values forβ (α) change rapidly and frequently hit the limits ±1/ √ ∆t. These results suggest that the least squares optimization in Step 2 may be improved with a smaller value of ∆t, which would require intra-day data. We therefore smooth the time-series data in Figs. 4a to 4c using one-year moving averages. The time-series for the smoothed parameters, denotedσ (α) ,μ (α) andβ (α) , are presented in Figs. 4d to 4f. The impact of the global financial crises is retained in the smoothed series. As the market disruption due to the Covid-19 pandemic was of shorter duration, the pandemic impact is lessened in the averaged data. Most significantly, the averaged valuesβ (α) are better behaved.

Using the last one-year estimation window (from 6/2/2020 to 6/1/2021), we obtain the estimatesσ (α) = 0.151,μ (α) = 0.119,β (α) = −0.978, andᾱ = (1 +β (α) √ ∆t)/2 = 0.469 for the date 6/1/2021. Based on these estimates, Fig. 5 shows a constructed, recombined, GJR price tree comprised of 30 price trajectories (10). 

We now consider the following discrete version of Corns' and Satchell's extension 25 of the Black-Scholes-Merton market model. 26 Suppose the dynamics 27 of the risky asset follows the GJR pricing tree (10). The dynamics 28 of the riskless asset is given by

where S T . Following the construction of the Cox-Ross-Rubinstein and JR binomial pricing models, our next goal is to use the pricing tree (10) to derive the price dynamics of f (n) k∆t , k = 0, 1, . . . , n. We start by forming the replicating risk-neutral portfolio P conditionally on F k∆t , the replicating portfolio should be riskless, that is P (7) and (10), and defining 30

conditionally on F k∆t , by (12) it follows that

The risk-neutrality assumption implies that D 

Furthermore, given F k∆t , the portfolio P (n) (k+1)∆t = P (n,u) (k+1)∆t is riskless, leading to f (n,u)

where the risk-neural probability q n,k+1 for the time period [k∆t, (k + 1)∆t) is given by 

. Constructing a binomial model using A (δ) instead of B (δ) would require two independent pricing trees to model the discrete price dynamics of the underlying asset, which is not desirable for derivative hedging. 26 For the Black-Scholes-Merton market model, see Black and Scholes (1973) and Duffie (2001, Chapter 6) . 27 We call the risky asset a stock, and it will be denoted by S. 28 We call the riskless asset a bond, and it will be denoted by B. 29 We use the terms ECC, option, and derivative interchangeably. 30 Here, "u" stands for upward movement and "d" for downward movement in the binomial tree.

To leading term in ∆t, q n,k+1 has the form,

where θ = µ − r f + σ 2 /2 /σ is the market price of risk. For β = 0, we obtain the risk-neutral probabilities under the JR binomial model (see Kim et al. (2016 Kim et al. ( , 2019 ). Now we consider the risk-neutral dynamics of the stock S under the GJR pricing tree model (10) and the ECC dynamics (17). Given F k∆t , to leading order in ∆t,

where k = 1, . . . , n−1. Conditionally on F k∆t , the discrete risk-neutral return R (n;q) (k+1)∆t = ln S (n;q) (k+1)∆t /S (n;q) k∆t , has mean E R (n;q) (k+1)∆t = r f − σ 2 /2 ∆t and variance Var R (n;q) (k+1)∆t = σ 2 ∆t, and for γ > 2, E R 

Following this framework, we use the data from section 2.4 and the market option prices 31 for the underlying SPY asset to compute implied µ, β, and σ surfaces. Let C is the spot price of the underlying SPY ETF; K is the strike price; T is the terminal time; and r f is the risk-free rate. 32 For brevity, we simply refer to C (SPY,Market) i . Applying equations (18), (19) and (20), we construct the GJR pricing tree's theoretical call option price for day t, denoted as C (SPY,GJR) i S (SPY) 0 , K, T, t, r f , ρ For brevity, we refer to this as C (SPY,GJR) i (ρ), where ρ = µ, β or σ is the parameter indicating the surface to be computed. For example, to obtain the implied µ surface, we use the estimated values forσ andβ from section 2.4 and designate ρ = ρ (µ) = (µ |σ (α) = 0.151,β (α) = −0.978). We estimate the value for ρ (µ) i for the i th call option contract bŷ

We proceed analogously using parameters ρ (β) and ρ (σ) to compute the β and σ surfaces. Fig. 6 shows the results for the implied µ and β surfaces, plotted in terms of moneyness K/S and time to maturity T (in days), where K is the strike price and S is the current spot price of the underlying asset. Of the two, the µ surface has the more complex behavior with different time-development of the surface as K/S ranges from "in the money" to "out of the money" values. In contrast, values ofβ decrease on both sides of a ridge of values that generally aligns with K/S ∼ 1.1. Values of β, ranging from −0.978 to 0.147, are mostly negative indicating a negatively skewed pricing tree. Values of β increase with time at any constant value of K/S. 31 The SPY call option data was collected from Bloomberg Professional Services on 6/1/2021, 19:32 EST. The data set includes call options, with all strike values, having expiration date no later than 12/31/2021. In total the data involves M = 1, 913 contracts with valid bid and ask quotes. The SPY spot price for this date and time was $419.67. 32 We use 10-year Treasury yield curve rates for risk-free rate values; the annual rate was r f = 1.62% on 6/1/2021. To evaluate the implied σ surface, we also calculated the Black-Scholes implied volatility surface 33 σ (BLS) for the same option prices, and considered their difference using the percent deviation, Dev (α,BLS) = 100 σ − σ (BLS) /σ (BLS) . These surfaces are shown in Fig. 7 . Both surfaces, σ and σ (BLS) , show similar volatility smiles and roughly similar values when K/S > 1. However, the σ surface increases more rapidly than σ (BLS) as K/S moves deeper into the money. 

There is a vast literature on option pricing incorporating transaction costs. 34 Discrete-time option pricing models have been studied under various assumptions regarding the types of the transaction costs incurred in trading the replicating self-financing portfolios. 35 Here, we extend the GJR option pricing model to the case when the hedging is subject to transaction costs. Our approach in extending the GJR pricing tree model to include transaction costs is based on Merton's binomial option pricing model with transaction costs. 36 33 See Hull (2012, Chapter 15) . 34 Some basic references are: Leland (1985) ; Hodges and Neuberger (1989) ; Boyle and Vorst (1992) ; Davis et al. (1993) ; Edirisinghe et al. (1993) ; Kabanov and Safarian (1997) ; Broadie et al. (1998) ; Kabanov and Stricker (2001) ; Palmer (2001) ; Lai and Lim (2009); and Guasoni et al. (2012) . 35 See Merton (1990, Chapter 14) , Stettner (1997) , Palmer (2001) , Delbaen et al. (2002) , Melnikov and Petrachenko (2005) , and Chen et al. (2008) . 36 See Merton (1990, Chapter 14) .

Consider a market with three securities (S, B, C) consisting of: (a) a stock S with pricing tree {S (n) k∆t , k = 1, . . . , n, n∆t = T } given by (10), (b) a bond B with price dynamics given by (16), and (c) an ECC C with price dynamics given by (17) . In contrast to the previous section, we assume that the hedger (the ECC-contract seller) trades the replicating risk-neutral portfolio P (n)

k∆t at a cost. Namely, given F k∆t = σ M (α∆t) 1 , . . . , M (α∆t) k , the hedged portfolio at (k + 1)∆t is determined by

In (22)

is the hedging transaction cost (HTC). From (22), to leading order in ∆t, it follows that

As in (18) and (19), we obtain the risk-neutral valuation of the ECC,

where the risk-neutral probability q (λ)

n,k+1 has the form

with θ (λ (0) ) = µ + σ 2 /2 − r f / 1 + λ (0) /σ. 37 In the special case, λ (0) = λ (1) = 0, (25) coincides with (19). Now consider the risk-neutral dynamics of the stock S in the presence of HTC. Conditional on F k∆t , k = 1, . . . , n − 1, the risk-neutral value of the stock price at (k + 1)∆t is determined by S (n;q,λ)

where q (λ) n,k+1 is given by (25). Conditionally on F k∆t , the discrete risk-neutral log-return R (n;q,λ) (k+1)∆t = ln S (n;q,λ) (k+1)∆t /S (n;q,λ) k∆t has mean and variance E R (n;q,λ)

and S (q,λ)

where B k∆t is similar to formula (14.2a) in Merton (1990, Chapter 14) .

However, for every fixed ∆t, the risk-neutral tree (20) will exhibit the behavior of a geometric SBM and the value of λ (1) becomes relevant, as evident from (24).

Remark: It appears natural to assume that HTC could be asymmetric; that is, instead of (22), we could assume that P (n,u,λ)

with λ

A close inspection shows that, as ∆t ↓ 0, (29) leads to a pricing model which is not arbitrage free. However, for fixed ∆t, the model (29) is arbitrage free and an expression for the risk-neutral probability q (λ) n,k+1 similar to formula (24) can be readily obtained. 38

As in section 3.1, we numerically illustrate the HTC and its limiting behavior using the same data set, including closing and call option prices, for the underlying asset SPY . Following the notation in section 3.1, we set the parameter ρ = ρ (λ) = λ (0) , λ (1) |μ (α) ,β (α) ,σ (α) . From (24), (25) and (26), we construct the GJR tree option price C (SPY,GJR) i (ρ (λ) ). We find the optimal solution for ρ (λ) by minimizating the relative mean-square error (relMSE) min λ (0) >0, λ (1) ∈R relMSE = min

The resulting minimizing values are relMSE = 0.433 for λ (0) = 28.8 and λ (1) = 0.297. 39 Using the parameters ρ (λ (0) ) = (λ (0) |μ (α) ,β (α) ,σ (α) , λ (1) = 0.297) and ρ (λ (1) ) = (λ (1) |μ (α) ,β (α) ,σ (α) , λ (0) = 28.8), we reran the minimization problem (21) and computed the λ (0) and λ (1) surfaces shown in Fig. 8 . Both surfaces have a similar ridge shape, but of vastly different scales. The value of λ (1) is constant through five computed digits, whereas λ (0) varies significantly in the first computed digit. Thus, for these values ofμ (α) ,β (α) andσ (α) , the constant value model, λ ∆t = λ (0) , would suffice in (22). To ascertain the effects of adding HTC, we computed the implied µ, β, and σ surfaces using the minimizing values for λ (0) and λ (1) . (For example, for the µ surface we estimated values for ρ (µ) = (µ |σ (α) = 0.151,β (α) = −0.978, λ (0) = 28.8, λ (1) = 0.297) using (21).) To compare the results with those reported in section 3.1, we also computed the percent deviation surfaces, ρ (dev) = 100 ρ (HTC) − ρ (GJR) /ρ (GJR) for ρ = {µ, σ}, where ρ (GJR) refers to the surface computed for the GJR model without HTC in section 3.1 and ρ (HTC) refers to the surface computed for the GJR model with HTC. To avoid division by zero, for the β surface comparison we plot the difference surface β (diff) = β (HTC) − β (GJR) . The results are shown in Fig. 9 . Inclusion of HTC has the effect of differentially lowering µ values (as expected since the trader is expending money on fees), most noticeably along the coordinate K/S ∼ 1.2. Similarly expected, inclusion of HTC differentially raises the volatility surface. HTC differentially raises β values, with greater increases arising as K/S values move more "into the money", resulting in larger regions of the phase space where the pricing tree is positively skewed. 38 Discrete-time market models, which are free of arbitrages, but for which the corresponding continuous-time market model is not arbitrage-free, are known; see, for example, Karandikar and Rachev (1995) and Hurst et al. (1999) . In those models, the discrete equivalent martingale measure explodes in the limit. Similarly, if we assume that (29) holds, the limiting martingale measure does not exist. 39 The three terms on the right-hand side of (22) have the form DS + λD∆S − f , where the transaction cost term is λD∆S. The estimated values of λ (0) and λ (1) yield λ ≈ λ (0) = 28.8. If daily stock price changes for the SPY fund are around 0.1% of the fund price (i.e. ∆S ∼ 10 −3 S), the transaction cost term would have the value 0.0288DS, indicating transactions costs contribute ∼3% to the total price of the hedged portfolio. 

In this section we explore possibilities for fitting the pricing tree model, S (n) k∆t , k = 1, . . . , n, to a market driver that affects the price dynamics, S (n,S) k∆t , of stock S. We consider both endogenous and exogenous approaches. In the endogenous approach, the discrete dynamics for the price of S is modeled by (j−1)∆t is the market's daily log-return for stock S. In (31), the parameters σ (0,S) , µ (0,S) and β (0,S) are determined for the initial time k∆t = 0 from equations (14) and (15) using the cumulative return data for S applicable at the chosen start date k∆t = 0. The values obtained from (14) and (15) are also smoothed using the moving-window averaging procedure introduced in section 2.4. This will be demonstrated explicitly by numerical example in section 5.1.

The endogenous approach is limited while the exogenous approach assumes that the M (α∆t) k are determined by the upward and downward movements of a more general market driver. For concreteness, we initially assume the market driver is a factor model for the price dynamics of S; in particular we consider the Fama-French five-factor asset pricing model. 

In (32), r f,j∆t is the risk-free rate in period [j∆t, (j + 1)∆t) and the remaining Fama-French factors are:

41 See Fama and French (2012 , 2015 , 2017 . Figure 9 : The implied (a) µ, (b) β and (c) σ surfaces generated by the GJR pricing tree with HTC added. The percent deviation surfaces µ (dev) and σ (dev) are shown in (d) and (f) respectively, while (e) displays the difference surface β (diff) . Each surface is plotted as a function of moneyness K/S and time to maturity T (in days). Let j∆t, j = 1, . . . , m denote a time period τ 1 and j = m + 1, . . . , m + n denote the following time period τ 2 . Letr (n,F ) j∆t , j = 1, . . . , m, denote sample return values (factor returns) for (32) obtained from the regression (33) over the period τ 1 , andR (n,F ) j∆t denote the resultant cumulative factor returns. Using the cumulative factor returns on the left-hand side of (14) and (15), and the moving-window averaging procedure introduced in section 2.4, we obtain estimatesσ (F ) ,μ (F ) andβ (F ) for the values of the parameters σ, µ and β at timestep m∆t. Let k∆t, k = 1, . . . , n, where k = j − m, label the trading days in time period τ 2 . Then σ (F ) ,μ (F ) andβ (F ) correspond to parameter values for k∆t = 0.

We model the discrete dynamics for the factor-driven price over the period τ 2 by S (n,exo) k∆t is then chosen as that path that minimizes the relMSE between the two sides of (34) when the stock prices S (n,S) k∆t are used on the left-hand side of (34). This exogenous method is also demonstrated in section 5.1.

We demonstrate the endogenous and exogenous approaches using the daily return series for the stock S = MSFT. The price data set 42 covers the period 4/30/2015 tthrough 4/30/2021. Since the estimation (equations (14) and (15)) and averaging procedures of section 2.4 for σ, µ and β require one-year moving windows, we divide the data set into two time periods, τ 1 = 4/31/2015 through 4/28/2017 and τ 2 = 5/1/2017 through 4/30/2021. MSFT data on the close of 4/28/2017 in τ 1 serves as t = 0 data for τ 2 ; the closing price of MSFT on 4/28/2017 serves as S 0 for the period τ 2 and cumulative returns over τ 2 are then computed relative to this S 0 . And applying the estimation and averaging procedures of section 2.4 to the cumulative log-return data for MSFT over τ 1 , we obtain estimatesσ (0,MSFT) ,μ (0,MSFT) andβ (0,MSFT) determined for the close of trading on 4/28/2017. These values, along withᾱ (0,MSFT) and the value of the relMSE, are given in Table 1 For the exogenous approach, we use the Fama-French five-factor model outlined in the previous section as the market driver. The Fama-French factor returns, r (n,F ) k∆t in (32), were computed over the period τ 1 using a robust form for the regression (33). 43 The values of the coefficients, as well as the R-square and root mean square-error (RMSE) values, are given in Table 2 . The sample returnr (n,F ) k∆t and priceŜ (n,F ) k∆t series computed from these coefficients are presented in Fig. 10 where they are compared with the market return r are indicative of issues with applying market driver analyses to stock returns. Since the Fama-French model is based upon a regression fit to return values, errors produced in return values have large-time decay structure and lead ultimately to significant cumulative differences in price performance.

The parameter estimatesμ (F ) ,σ (F ) ,β (F ) , andᾱ (F ) obtained for 4/28/2017 using the cumulative Fama-French returns over τ 1 are given in Table 1 . These values generally differ by a factor of ∼ 3 compared to the corresponding endogenous variables, with the exception of a factor of ∼ 13 difference betweenσ (F ) and σ (0,MSFT) . Using the estimated valueᾱ (F ) , we generate 10 6 scenarios of the sample path M (ᾱ (F ) ) k covering 42 MSFT market closing prices from Bloomberg Professional Services. Return values for the Fama-French factors from the U.S. Research Returns Data Center (https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data Library.html).

43 Motivated by results in Knez and Ready (1997) , we choose robust regression in estimating the parameters of the Fama-French five-factor model. × Ö R-square RMSE −6.5 · 10 −4 1.2 · 10 −2 −3.9 · 10 −3 −4.0 · 10 −4 3.0 · 10 −3 −5.8 · 10 −3 0.58 9.1 · 10 −3 Fig. 11 . Both estimates produce reasonable agreement to the actual price data; however both significantly smooth and decrease the impact of the Covid-19 pandemic; the exogenous method showing less impact from the pandemic than the endogenous method.

To demonstrate the application of the model (35), (36), we use (36) to estimate the MSFT log-return r (n,MSFT) k∆t over the period τ 2 . To capture likely heavy-tailed behavior, we identify h(·) as the Student's t probability density function with κ degrees of freedom. For convenience, we utilize the optimal path M (ᾱ (F ) ,exo) k determined from the exogenous procedure in section 5.1 for the required path in (36). To establish the model parameters, v, σ, γ, and κ, we construct the conditional least squares minimization problem, min ν∈R,σ∈R,γ∈R,κ∈ [5, 30] r (n,MSFT) k∆t − r (α∆t,n;h) k∆t 2 2 .

Table 3 provides the resulting parameter estimates. The model (35), (36) with these parameter estimates will be used further in section 6.1.

v σ γ κ RMSE 0.272 3.44 · 10 −5 4.22 · 10 −2 6.24 1.86 · 10 −2 

Given the filtration F (n;M (α ∆t ) ) = {F where the risk-neutral probability, q (α∆t,n;h) 

In the option pricing literature since 1990, academic work on continuous-time option pricing models (CPMs) has greatly overshadowed the work on discrete-time option pricing models (DPMs). Moreover, with very few exceptions, DPMs are used only to approximate the CPMs (presumably to decrease the numerical complexity in generating Monte Carlo price trajectories). Those DPMs are placed directly into the riskneutral world and do not address an issue that is very important for every option trader: "Which DPM in the natural world is uniquely determined by the DPM introduced directly in the risk-neutral world?" The reason for the "subordinated" role of the DPM's is the "obsessive love" of academics for CPMs, arising from the alluring mathematical beauty of the semimartingale theory 48 of the Strasbourg school and the truly paramount Fundamental Theorem of Asset Pricing. 49 However, crucial information found in DPMs in the natural world, such as the mean-return parameter, the probability for stock-upturn, and the distributional skewness and kurtosis, is lost in CPMs. Most notably, in CPMs the mean-return parameter is lost due to the assumed ability of the hedger to trade continuously in time. DPMs with defined dynamics in the natural world do not have this issue. 50 That observation was the main motivation for this paper. Specifically, we have extended the classical Jarrow-Rudd pricing tree to include skewness and kurtosis in the underlying asset's return distribution in both the natural and in the risk-neutral world. We have extended Merton's option pricing tree model with hedging-transaction costs to our new generalized Jarrow Rudd (GJR) option pricing tree model. We have applied the Cherny-Shiryaev-Yor invariance principle and the Fama-French fivefactor model to further extend the GJR pricing tree model to cover path-dependent options. Our numerical examples include estimation of the implied surfaces of all parameters in the GJR pricing trees.

On the multi-dimensional skew Brownian motion

Option pricing under a normal mixture distribution derived from the Narkov tree model

Convergence of Probability Measures

Fact and fantasy in the use of options

The pricing of options and corporate liabilities

Optimal replication in discrete time with transaction costs

On measuring the risk of common stocks implied by options prices: A note

Optimal replication of contingent claims under portfolio constraints

The least cost super replicating portfolio in the Boyle-Vorst model with transaction costs

Limit behavior of the "horizontal-vertical" random walk and some extensions of the Donsker-Prokhorov invariance principle

Jacod and Shiryaev (1987), and Protter

Skew Brownian motion and pricing European options

Options pricing: a simplified approach

European option pricing with transaction costs

On a non-classical invariance principle

Hedging under transaction costs in currency markets: A discrete-time model

A general version of the fundamental theorem of asset pricing

The fundamental theorem of asset pricing for unbounded stochastic processes

An invariant principle for certain probability limit theorems. Memoirs of the

American option pricing under GARCH by a Markov chain approximation

Dynamic Asset Pricing Theory, Third Edition

Optimal replication of options with transactions costs and trading restrictions

Séminaire de Probabilités

A five-factor asset pricing model

On a Markov chain approximation method for option pricing with regime switching

Option pricing with Markov switching

Introduction to the Theory of Random Processes

The fundamental theorem of asset pricing under transaction costs

On skew Brownian motion

Options and bubbles

Numerical methods for Lévy processes

Optimal replication of contingent claims under transactions costs

Robust regression using iteratively reweighted least-squares

Option pricing in markets with informed traders

Option pricing incorporating factor dynamics in complete markets

Options, Futures, and Other Derivatives, Eighth Edition

Option pricing for a logstable asset price model

Diffusion Processes and Their Sample Paths

Limit Theorems for Stochastic Processes

Asset price bubbles in incomplete markets

Option Pricing

On Icland's strategy of option pricing with transaction costs

The Harrison-Pliska arbitrage pricing theorem under transaction costs

A generalized binomial model and option pricing formulae for subordinated stock-price processes

A simple option pricing model with Markovian volatilities

Multi-purpose binomial model: fitting all moments to the underlying Brownian motion

Enhancing binomial and trinomial option pricing models

On the robustness of size and book-to-market in cross-sectional

Option hedging theory under transaction costs

Effective conductivity and skew Brownian motion

On the constructions of the skew Brownian motion

Option pricing and replication with transactions costs

On general skew Brownian motions

Recovering an asset's implied pdf from option prices: an application to crude oil during the Gulf Crisis

On option pricing in binomial market with transaction costs

Theory of rational option pricing

Continuous-Time Finance

A note on the Boyle-Vorst discrete-time option pricing model with transactions costs

Financial market dislocations. The Review of Financial Studies

Logistic regression diagnostics

Convergence of random processes and limit theorems in probability theory. Theory of Probability and its Applications

Stochastic Integration and Differential Equations

Multi-skewed Brownian motion and diffusion in layered media

Continuous Martingales and Brownian Motion

The recovery theorem

Computing American option prices in the lognormal jump-diffusion framework with a Markov chain

Basic Principles and Applications of Probability Theory

Option pricing in the CRR model with proportional transaction costs: A cone transformation approach

Infinitely Divisibility of Probability Distributions on the Real Line

On countably skewed Brownian motion with accumulation point

American option prices in a Markov chain market model

Asset pricing using finite state Markov chain stochastic discount functions. Stochastic Analysis and Applications

Financial Market Bubbles and Crashes

Properties of multinomial lattices with cumulants for option pricing and hedging

(α,n) t For a general market driver, we extend to the full functionality of the CYSIP by adding the higher moment SBM process C (α,n) t into our model. 44 According to (34), S (n,exo) k∆t is adapted to the filtration F (n;M (α ∆t ) ) = F (n;M (α ∆t ) ) k = σ M (α∆t) j ; j = 1, . . . , k , k = 1, . . . , n, F , t ≥ 0 to be a random process with the piecewise linear trajectories having vertexes k/n, C (α∆t,n;h) k/n , k ∈ N , n ∈ N , where C (α∆t,n;h) k/n = Y (α∆t,n;h) k/n . Then, for a fixed, relatively small value of ∆t, the bivariate process BNext, consider the following two processes in D[0, T ]:Then, for large n, the bivariate process BUltimately, as n ↑ ∞, Bwhere v ∈ R, σ ∈ R\{0} and γ ∈ R are parameters determining the dynamics of the stock price as a function of the index dynamics. The discrete dynamics of the stock log-return is given byThus, if γ = 0, the stock log-return r = S (α∆t,n;h)where k = 1, . . . , n, n∆t = T . Let h(x) ≥ 0, x ∈ R, be a chosen probability density function. 46 DefineNote that η For example, h is the Student's t density function with κ-degrees of freedom.(k+1)∆t , for the time period [k∆t, (k + 1)∆t) is given byTo leading order in ∆t we have q (α∆t,n;h)Note that, conditionally on F 

Using the daily closing prices and the corresponding log-returns for MSFT, we apply the full GJR model to generate a binomial tree to construct call option prices using (40) Table 3 .(a) (b) (c) Figure 12 : (a) The implied σ (ν,γ,κ,ᾱ∆t) surface generated by the path-dependent GJR pricing tree. (b) The implied σ (BLS) surface using the Black-Scholes formula for implied volatility. (c) The percent deviation surface Dev (ν,γ,κ,ᾱ∆t,BLS) . Each surface is plotted as functions of moneyness K/S and time to maturity T (in days).As in the numerical examples in sections 3.1 and 4.1, we compute the implied volatility, σ (ν,γ,κ,ᾱ∆t) , based on the path-dependent GJR pricing tree. We also estimate the Black-Scholes implied volatility σ (BLS) surface for the same option prices, and define the percent deviation between the two surfaces, Dev (ν,γ,κ,ᾱ∆t,BLS) = 100 σ (ν,γ,κ,ᾱ∆t) − σ (BLS) /σ (BLS) . Fig. 12 shows the implied volatility surfaces plotted as functions of time to maturity T (in days) and moneyness K/S. For maturity times 50 days, both the GJR and the Black-Scholes surfaces contain volatility smiles and the deviation between the two surfaces is relatively small. Beyond 50 days, the Black-Scholes implied volatility surface flattens, while the GJR surface flattens only in the "out-of-the-money" region, leading to large differences between the two over the "in-the-money" region.