]> High-resolution path-integral development of financial options

High-resolution path-integral development of financial options

1. INTRODUCTION
1.1. Background
1.2. Our Approach
1.3. Outline of Paper
2. DATA
2.1. Eurodollars
2.2. Futures
2.3. Options on Futures
2.4. Front/Back Month Contracts
2.5. Stochastic Volatility
3. MODELS
3.1. Random walk model
3.2. Black-Scholes (BS) Theory
3.3. Some Key Issues in Derivation of BS
3.4. Models
3.4.1. Various Diffusions
4. STATISTICAL MECHANICS OF FINANCIAL MARKETS (SMFM)
4.1. Statistical Mechanics of Large Systems
4.2. Correlations
5. ADAPTIVE SIMULATED ANNEALING (ASA) FITS
5.1. ASA Outline
5.1.1. General description
5.1.2. Mathematical outline
5.1.3. Reannealing
5.1.4. Quenching
5.2. -Indicator of Market Contexts
6. PATH-INTEGRAL (PATHINT) DEVELOPMENT
6.1. PATHINT Outline
6.2. Development of Long-Time Probabilities
6.3. Dependence of Probabilities on and
6.4. Two-Factor Volatility and PATHINT Modifications
7. CALCULATION OF DERIVATIVES
7.1. Primary Use of Probabilities For European Options
7.2. PATHINT Baselined to CRR and BS
7.3. Two-Factor PATHINT Baselined to One-Factor PATHINT
8. CONCLUSIONS
ACKNOWLEDGMENTS
REFERENCES
FIGURE CAPTIONS
TABLE CAPTIONS
Figure 1
Figure 2
Figure 3
Table 1

Lester Ingber

Lester Ingber Research
POB 06440 Sears Tower, Chicago, IL 60606
and
DRW Investments LLC
311 S Wacker Dr, Ste 900, Chicago, IL 60606

ingber@ingber.com, ingber@alumni.caltech.edu

+1.312.542.1006 Voice
+1.801.881.1088 Fax

ABSTRACT

The Black-Scholes theory of option pricing has been considered for many years as an important but very approximate zeroth-order description of actual market behavior. We generalize the functional form of the diffusion of these systems and also consider multi-factor models including stochastic volatility. Daily Eurodollar futures prices and implied volatilities are fit to determine exponents of functional behavior of diffusions using methods of global optimization, Adaptive Simulated Annealing (ASA), to generate tight fits across moving time windows of Eurodollar contracts. These short-time fitted distributions are then developed into long-time distributions using a robust non-Monte Carlo path-integral algorithm, PATHINT, to generate prices and derivatives commonly used by option traders.

Keywords: options; eurodollar; volatility; path integral; optimization; statistical mechanics

PACS Nos.: 05.10.-a, 02.70.-c, 82.20.wt, 89.90.+n

1. INTRODUCTION

1.1. Background

There always is much interest in developing more sophisticated pricing models for financial instruments. In particular, there currently is much interest in improving option pricing models, particularly with respect to stochastic variables [1-4].
The standard Black-Scholes (BS) theory assumes a lognormal distribution of market prices, i.e., a diffusion linearly proportional to the market price. However, many texts include outlines of more general diffusions proportional to an arbitrary power of the market price [5].
The above aspects of stochastic volatility and of more general functional dependencies of diffusions are most often “swept under the rug” of a simple lognormal form. Experienced traders often use their own intuition to put volatility “smiles” into the BS theoretical constant coefficient in the BS lognormal distribution to compensate for these aspects.
It is generally acknowledged that since the market crash of 1987, markets have been increasingly difficult to describe using the BS model, and so better modelling and computational techniques should be used traders [6], although in practice simple BS models are the rule rather than the exception simply because they are easy to use [7]. To a large extent, previous modelling that has included stochastic volatility and multiple factors has been driven more by the desire to either delve into mathematics tangential to these issues, or to deal only with models that can accommodate closed-form algebraic expressions. We do not see much of the philosophy in the literature that has long driven the natural sciences: to respect first raw data, secondly models of raw data, and finally the use of numerical techniques that do not excessively distort models for the sake of ease of analysis and speed of computation. Indeed, very often the reverse set of priorities is seen in mathematical finance.

1.2. Our Approach

We have addressed the above issues in detail within the framework of a previously developed statistical mechanics of financial markets (SMFM) [8-13].
Our approach requires three sensible parts. Part one is the formulation of the model, which to some extent also involves specification of the specific market(s) data to be addressed. Part two is the fitting of the model to specific market data. Part three is the use of the resulting model to calculate option prices and their Greeks (partial derivatives of the prices with respect to their independent variables), which are used as risk parameters by traders. Each part requires some specific numerical tuning to the market under consideration.
The first part was to develop the algebraic model to replace/generalize BS, including the possibility of also addressing how to handle data regions not previously observed in trading. This is not absurd; current BS models perform integrals that must include a much influence from fat tails that include data regions never seen or likely to be seen in real-world markets. There are some issues as to whether we should take seriously the notion that the market is strongly driven by some element of a “self-fulfilling prophesy” by the BS model [14], but in any case our models have parameters to handle a wide range of possible cases that might arise.
We have developed two parallel tracks starting with part one, a one-factor and a two-factor model. The two-factor model includes stochastic volatility. At first we sensed the need to develop this two-factor model, and we now see that this is at the least an important benchmark against which to judge the worth of the one-factor model.
The second part was to fit the actual raw data so we can come up with real distributions. Some tests illustrated that standard quasi-linear fitting routines, could not get the proper fits, and so we used a more powerful global optimization, Adaptive Simulated Annealing (ASA) [15]. Tuning and selection of the time periods to perform the fits to the data were not trivial aspects of this research. Practical decisions had to be made on the time span of data to be fit and how to aggregate the fits to get sensible “fair values” for reasonable standard deviations of the exponents in the diffusions.
The third part was to develop Greeks and risk parameters from these distributions without making premature approximations just to ease the analysis. Perhaps someday, simple approximations and intuitions similar to what traders now use for BS models will be available for these models, but we do not think the best approach is to start out with such approximations until we first see proper calculations, especially in this uncharted territory. When it seemed that Cox-Ross-Rubenstein (CRR) standard tree codes (discretized approximations to partial differential equations) [16] were not stable for general exponents, i.e., for other than the lognormal case, we turned to a PATHINT code developed a decade ago for some hard nonlinear multifactor problems [17], e.g., combat analyses [18], neuroscience [19,20], and potentially chaotic systems [21,22]. In 1990 and 1991 papers on financial applications, it was mentioned how these techniques could be used for stochastic interest rates and bonds [9,10]. The modifications required here for one-factor European and then American cases went surprisingly smoothly; we still had to tune the meshes, etc. The two-factor model presented a technical problem to the algorithm, which we have reasonably handled using a combination of selection of the model in part one and a reasonable approach to developing the meshes.
More detailed numerical calculations than contained in this paper will be presented in a subsequent paper [23].

1.3. Outline of Paper

Section 1 is this introduction. Section 2 describes the nature of Eurodollar (ED) futures data and the evidence for stochastic volatility. Section 3 outlines the algebra of modelling options, including the standard BS theory and our generalizations. Section 4 outlines the three equivalent mathematical representations used by SMFM; this is required to understand the development of the short-time distribution that defines the cost function we derive for global optimization, as well as the numerical methods we have developed to calculate the long-time evolution of these short-time distributions. Section 5 outlines ASA and explains its use to fit short-time probability distributions defined by our models to the Eurodollar data; we offer the fitted exponent in the diffusion as a new important technical indicator of market behavior. Section 6 outlines PATHINT and explains its use to develop long-time probability distributions from the fitted short-time probability distributions, for both the one-factor and two-factor tracks. Section 7 describes how we use these long-time probability distributions to calculate European and American option prices and Greeks; here we give numerical tests of our approach to BS CRR algorithms. Section 8 is our conclusion.

2. DATA

2.1. Eurodollars

Eurodollars are fixed-rate time deposits held primarily by overseas banks, but denominated in US dollars. They are not subject to US banking regulations and therefore tend to have a tighter bid-ask spread than deposits held in the United States [24].

2.2. Futures

The three-month Eurodollar futures contract is one of the most actively traded futures markets in the world. The contract is quoted as an index where the yield is equal to the Eurodollar price subtracted from 100. This yield is equal to the fixed rate of interest paid by Eurodollar time deposits upon maturity and is expressed as an annualized interest rate based on a 360-day year. The Eurodollar futures are cash settled based on the 90-day London Interbank Offer Rate (LIBOR). A “notional” principal amount of $1 million, is used to determine the change in the total interest payable on a hypothetical underlying time deposit, but is never actually paid or received [24].
Currently a total of 40 quarterly Eurodollar futures contracts (or ten years worth) are listed, with expirations annually in March, June, September and December.

2.3. Options on Futures

The options traded on the Eurodollar futures include not only 18 months of options expiring at the same time as the underlying future, but also various short dated options which themselves expire up to one year prior to the expiration of the underlying futures contract.

2.4. Front/Back Month Contracts

For purposes of risk minimization, as discussed in a previous paper [4], traders put on spreads across a variety of option contracts. One common example is to trade the spread on contracts expiring one year apart, where the future closer to expiration is referred to as the front month contract, and the future expiring one year later is called the back month. The availability of short dated or “mid-curve” options which are based on an underlying back month futures contract, but expire at the same time as the front month, allow one to trade the volatility ratios of the front and back month futures contracts without having to take the time differences in option expirations into consideration. We studied the volatilities of these types of front and back month contracts. Here, we give analyses with respect only to quarterly data longer than six months from expiration.

2.5. Stochastic Volatility

Below we develop two-factor models to address stochastic volatility. In a previous paper, we have performed empirical studies of Eurodollar futures to support the necessity of dealing with these issues [4].

3. MODELS

3.1. Random walk model

The use of Brownian motion as a model for financial systems is generally attributed to Bachelier [25], though he incorrectly intuited that the noise scaled linearly instead of as the square root relative to the random log-price variable. Einstein is generally credited with using the correct mathematical description in a larger physical context of statistical systems. However, several studies imply that changing prices of many markets do not follow a random walk, that they may have long-term dependences in price correlations, and that they may not be efficient in quickly arbitraging new information [26-28]. A random walk for returns, rate of change of prices over prices, is described by a Langevin equation with simple additive noise $η$ , typically representing the continual random influx of information into the market.

\dot{Γ} = - γ_{1} + γ_{2} η,

\dot{Γ} = d Γ / d t,

< η (t) >_{η} = 0, < η (t), η (t') >_{η} = δ (t - t'),

(1)

where $γ_{1}$ and $γ_{2}$ are constants, and $Γ$ is the logarithm of (scaled) price. Price, although the most dramatic observable, may not be the only appropriate dependent variable or order parameter for the system of markets [29]. This possibility has also been called the “semistrong form of the efficient market hypothesis” [26].
The generalization of this approach to include multivariate nonlinear nonequilibrium markets led to a model of statistical mechanics of financial markets (SMFM) [8].

3.2. Black-Scholes (BS) Theory

The standard partial-differential equation used to formulate most variants of Black-Scholes (BS) models describing the market value of an option, $V$ , is

\frac{\partial V}{\partial t} + \frac{1}{2} σ^{2} S^{2} \frac{\partial^{2} V}{\partial S^{2}} + r S \frac{\partial V}{\partial S} - r V = 0,

(2)

where $S$ is the asset price, and $σ$ is the standard deviation, or volatility of $S$ , and $r$ is the short-term interest rate. The solution depends on boundary conditions, subject to a number of interpretations, some requiring minor transformations of the basic BS equation or its solution. For example, the basic equation can apply to a number of one-dimensional models of interpretations of prices given to $V$ , e.g., puts or calls, and to $S$ , e.g., stocks or futures, dividends, etc.
For instance, if $V$ is set to $C$ , a call on an European option with exercise price $X$ with maturity at $T$ , the solution is

C (S, t) = S N (d_{1}) - X e^{- r (T - t)} N (d_{2}),

d_{1} = \frac{ln (S / X) + (r + \frac{1}{2} σ^{2}) (T - t)}{σ (T - t)^{1 / 2}},

d_{2} = \frac{ln (S / X) + (r - \frac{1}{2} σ^{2}) (T - t)}{σ (T - t)^{1 / 2}} .

(3)

In practice, the volatility $σ$ is the least known parameter in this equation, and its estimation is generally the most important part of pricing options. Usually the volatility is given in a yearly basis, baselined to some standard, e.g., 252 trading days per year, or 360 or 365 calendar days. Therefore, all values of volatility given in the graphs in this paper, based on daily data, would be annualized by multiplying the standard deviations of the yields by $\sqrt{252} = 15.87$ We have used this factor to present our implied volatilities as daily movements.

3.3. Some Key Issues in Derivation of BS

The basic BS model considers a portfolio in terms of delta ( $Δ$ ),

Π = V - Δ S,

(4)

in a market with Gaussian-Markovian (“white”) noise $X$ and drift $μ$ ,

\frac{d S}{S} = σ d X + μ d t,

(5)

where $V (S, t)$ inherits a random process from $S$ ,

d V = σ S \frac{\partial V}{\partial S} d X + (μ S \frac{\partial V}{\partial S} + \frac{1}{2} σ^{2} S^{2} \frac{\partial^{2} V}{\partial S^{2}} + \frac{\partial V}{\partial t}) d t .

(6)

This yields

d Π = σ (\frac{\partial V}{\partial S} - Δ) d X + (μ S \frac{\partial V}{\partial S} + \frac{1}{2} σ^{2} S^{2} \frac{\partial^{2} V}{\partial S^{2}} + \frac{\partial V}{\partial t} - μ Δ S) d t .

(7)

The expected risk-neutral return of $Π$ is

d Π = r Π d t = r (V - Δ S) d t .

(8)

Options $V$ on futures $F$ can be derived, e.g., using simple transformations to take cost of carry into consideration, such as

F = S e^{r (T - t)},

(9)

and setting

d Π = r V d t .

(10)

The corresponding BS equation for futures $F$ is

\frac{\partial V}{\partial t} + \frac{1}{2} σ^{2} F^{2} \frac{\partial^{2} V}{\partial S^{2}} - r V = 0 .

(11)

At least two advantages are present if $Δ$ is chosen such that

Δ = \frac{\partial V}{\partial S} .

(12)

Then, the portfolio can be instantaneously “risk-neutral,” in terms of zeroing the coefficient of $X$ , as well as independent of the direction of market, in terms of zeroing the coefficient of $μ$ For the above example of $V$ = $C$ ,

Δ = N (d_{1}) .

(13)

Other trading strategies based on this simple model use similar constructs as risk parameters, e.g., gamma ( $Γ$ ), theta ( $Θ$ ), vega ( $Υ$ ), rho ( $ρ$ ) [5],

Γ = \frac{\partial^{2} Π}{\partial S^{2}},

Θ = \frac{\partial Π}{\partial t},

Υ = \frac{\partial Π}{\partial σ},

ρ = \frac{\partial Π}{\partial r} .

(14)

The BS equation, Eq. (2), may be written as

Θ + r S Δ + \frac{1}{2} (σ S)^{2} Γ = r f .

(15)

3.4. Models

Our two-factor model includes stochastic volatility $σ$ of the underlying $S$ ,

d S = μ d t + σ F (S, S_{0}, S_{\infty}, x, y) d z_{S},

d σ = ν d t + ε d z_{σ},

< d z_{i} > = 0, i = {S, σ},

< d z_{i} (t) d z_{j} (t') > = d t δ (t - t'), i = j,

< d z_{i} (t) d z_{j} (t') > = ρ d t δ (t - t'), i \neq j,

F (S, S_{0}, S_{\infty}, x, y) = {\begin{array}{l} S, \\ S^{x} S_{0}^{1 - x}, \\ S^{y} S_{0}^{1 - x} S_{\infty}^{x - y}, \end{array} \begin{array}{l} S < S_{0} \\ S_{0} \leq S \leq S_{\infty} \\ S > S_{\infty} \end{array}, (null)

(16)

where $S_{0}$ and $S_{\infty}$ are selected to lie outside the data region used to fit the other parameters, e.g., $S_{0} = 1$ and $S_{\infty} = 20$ for fits to Eurodollar futures which historically have a very tight range relative to other markets. We have used the Black-Scholes form $F = S$ inside $S < S_{0}$ to obtain the usual benefits, e.g., no negative prices as the distribution is naturally excluded from $S < 0$ and preservation of put-call parity. Put-call parity for European options is derived quite independent of any mathematical model of options [5]. In its simplest form, it is given by

c + X e^{- r (T - t)} = p + S,

(17)

where $c$ ( $p$ ) is the fair price of a call (put), $X$ is the strike price, $r$ is the risk-free interest rate, $t$ is the present time, $T$ is the time of expiration, and $S$ is the underlying market. We have taken $y = 0$ , a normal distribution, to reflect total ignorance of markets outside the range of $S > S_{\infty}$ The one-factor model just assumes a constant $σ$ It is often noted that BS models incorrectly include untenable contributions from large $S$ regions because of their fat tails [30]. (If we wished to handle negative interest rates, ED prices > 100, we would move shift the $S = 0$ axis to some $S < 0$ value.)
We found that the abrupt, albeit continuous, changes across $S_{0}$ especially for $x \leq 0$ did not cause any similar effects in the distributions evolved using these diffusions, as reported below.
The formula for pricing an option $P$ , derived in a Black-Scholes generalized framework after factoring out interest-rate discounting, is equivalent to using the form

d S = μ S d t + σ F (S, S_{0}, S_{\infty}, x, y) d z_{S},

d σ = ν d t + ε d z_{σ} .

(18)

We experimented with some alternative functional forms, primarily to apply some smooth cutoffs across the above three regions of $S$ For example, we used $F'$ , a function $F$ designed to revert to the lognormal Black-Scholes model in several limits,

F' (S, S_{0}, S_{\infty}, x) = S C_{0} + (1 - C_{0}) (S^{x} S_{0}^{1 - x} C_{\infty} + S_{0} (1 - C_{\infty})),

C_{0} = exp [- {(\frac{S}{S_{0}} \frac{| 1 - x |}{1 + | 1 - x |})}^{| 2 - x | + 1}],

C_{\infty} = exp [- {(\frac{S}{S_{\infty}})}^{2}],

lim_{S \to \infty, x \neq 1} F' (S, S_{0}, S_{\infty}, x) = S_{0} = c o n s t a n t,

lim_{S \to 0^{+}} F' (S, S_{0}, S_{\infty}, x) = lim_{x \to 1} F' (S, S_{0}, S_{\infty}, x) = S .

(19)

However, our fits were most sensitive to the data when we permitted the central region to be pure $S^{x}$ using $F$ above.

3.4.1. Various Diffusions

Fig. 1. gives examples of $F (S, S_{0}, S_{\infty}, x, y) d z_{S}$ for $x$ in {-1, 0, 1, 2}. The other parameters are $S = 5$ , $S_{0} = 0.5$ , $S_{\infty} = 20$ , $y = 0$

4. STATISTICAL MECHANICS OF FINANCIAL MARKETS (SMFM)

4.1. Statistical Mechanics of Large Systems

Aggregation problems in nonlinear nonequilibrium systems typically are “solved” (accommodated) by having new entities/languages developed at these disparate scales in order to efficiently pass information back and forth. This is quite different from the nature of quasi-equilibrium quasi-linear systems, where thermodynamic or cybernetic approaches are possible. These approaches typically fail for nonequilibrium nonlinear systems.
Many systems are aptly modeled in terms of multivariate differential rate-equations, known as Langevin equations,

{\dot{M}}^{G} = f^{G} + {\hat{g}}_{j}^{G} η^{j}, (G = 1, . . ., Λ), (j = 1, . . ., N),

{\dot{M}}^{G} = d M^{G} / d t,

< η^{j} (t) >_{η} = 0, < η^{j} (t), η^{j'} (t') >_{η} = δ^{j j'} δ (t - t'),

(20)

where $f^{G}$ and ${\hat{g}}_{j}^{G}$ are generally nonlinear functions of mesoscopic order parameters $M^{G}$ , $j$ is a microscopic index indicating the source of fluctuations, and $N \geq Λ$ The Einstein convention of summing over repeated indices is used. Vertical bars on an index, e.g., |j|, imply no sum is to be taken on repeated indices.
Via a somewhat lengthy, albeit instructive calculation, outlined in several other papers [8,10,31], involving an intermediate derivation of a corresponding Fokker-Planck or Sch¨odinger-type equation for the conditional probability distribution $P [M (t) | M (t_{0})]$ , the Langevin rate Eq. (20) is developed into the more useful probability distribution for $M^{G}$ at long-time macroscopic time event $t = (u + 1) θ + t_{0}$ , in terms of a Stratonovich path-integral over mesoscopic Gaussian conditional probabilities [32-36]. Here, macroscopic variables are defined as the long-time limit of the evolving mesoscopic system.
The corresponding Sch¨odinger-type equation is [34,35]

\partial P / \partial t = \frac{1}{2} (g^{G G'} P)_{, G G'} - (g^{G} P)_{, G} + V,

g^{G G'} = k_{T} δ^{j k} {\hat{g}}_{j}^{G} {\hat{g}}_{k}^{G'},

g^{G} = f^{G} + \frac{1}{2} δ^{j k} {\hat{g}}_{j}^{G'} {\hat{g}}_{k, G'}^{G},

[. . .]_{, G} = \partial [. . .] / \partial M^{G} .

(21)

This is properly referred to as a Fokker-Planck equation when $V \equiv 0$ Note that although the partial differential Eq. (21) contains information regarding $M^{G}$ as in the stochastic differential Eq. (20), all references to $j$ have been properly averaged over. I.e., ${\hat{g}}_{j}^{G}$ in Eq. (20) is an entity with parameters in both microscopic and mesoscopic spaces, but $M$ is a purely mesoscopic variable, and this is more clearly reflected in Eq. (21).
The path integral representation is given in terms of the “Feynman” Lagrangian

L

P [M_{t} | M_{t_{0}}] d M (t) = \int . . . \int \underline{D} M exp (- S) δ [M (t_{0}) = M_{0}] δ [M (t) = M_{t}],

S = k_{T}^{- 1} min \int_{t_{0}}^{t} d t' L,

\underline{D} M = lim_{u \to \infty} Π_{ρ = 1}^{u + 1} g^{1 / 2} \underset{G}{Π} (2 π θ)^{- 1 / 2} d M_{ρ}^{G},

L ({\dot{M}}^{G}, M^{G}, t) = \frac{1}{2} ({\dot{M}}^{G} - h^{G}) g_{G G'} ({\dot{M}}^{G'} - h^{G'}) + \frac{1}{2} {h^{G}}_{; G} + R / 6 - V,

h^{G} = g^{G} - \frac{1}{2} g^{- 1 / 2} (g^{1 / 2} g^{G G'})_{, G'},

g_{G G'} = (g^{G G'})^{- 1},

g = det (g_{G G'}),

{h^{G}}_{; G} = h_{, G}^{G} + Γ_{G F}^{F} h^{G} = g^{- 1 / 2} (g^{1 / 2} h^{G})_{, G},

Γ_{J K}^{F} \equiv g^{L F} [J K, L] = g^{L F} (g_{J L, K} + g_{K L, J} - g_{J K, L}),

R = g^{J L} R_{J L} = g^{J L} g^{J K} R_{F J K L},

R_{F J K L} = \frac{1}{2} (g_{F K, J L} - g_{J K, F L} - g_{F L, J K} + g_{J L, F K}) + g_{M N} (Γ_{F K}^{M} Γ_{J L}^{N} - Γ_{F L}^{M} Γ_{J K}^{N}) .

(22)

Mesoscopic variables have been defined as $M^{G}$ in the Langevin and Fokker-Planck representations, in terms of their development from the microscopic system labeled by $j$ The Riemannian curvature term $R$ arises from nonlinear $g_{G G'}$ , which is a bona fide metric of this space [34]. Even if a stationary solution, i.e., ${\dot{M}}^{G} = 0$ , is ultimately sought, a necessarily prior stochastic treatment of ${\dot{M}}^{G}$ terms gives rise to these Riemannian “corrections.” Even for a constant metric, the term ${h^{G}}_{; G}$ contributes to $L$ for a nonlinear mean $h^{G}$ $V$ may include terms such as $\underset{T'}{Σ} J_{T' G} M^{G}$ , where the Lagrange multipliers $J_{T' G}$ are constraints on $M^{G}$ , which are advantageously modeled as extrinsic sources in this representation; they too may be time-dependent.
For our purposes, the above Feynman Lagrangian defines a kernel of the short-time conditional probability distribution, in the curved space defined by the metric, in the limit of continuous time, whose iteration yields the solution of the previous partial differential equation Eq. (21). This differs from the Lagrangian which satisfies the requirement that the action is stationary to the first order in $d t$ — the WKBJ approximation, but which does not include the first-order correction to the WKBJ approximation as does the Feynman Lagrangian. This latter Lagrangian differs from the Feynman Lagrangian, essentially by replacing $R / 6$ above by $R / 12$ [37]. In this sense, the WKBJ Lagrangian is more useful for some theoretical discussions [38]. However, the use of the Feynman Lagrangian coincides with the numerical method we present here using the PATHINT code.
Using the variational principle, $J_{T G}$ may also be used to constrain $M^{G}$ to regions where they are empirically bound. More complicated constraints may be affixed to $L$ using methods of optimal control theory [39]. With respect to a steady state $\bar{P}$ , when it exists, the information gain in state $P$ is defined by

Υ [P] = \int . . . \int \underline{D} M' P ln (P / \bar{P}),

\underline{D} M' = \underline{D} M / d M_{u + 1} .

(23)

In the economics literature, there appears to be sentiment to define Eq. (20) by the Ito, rather than the Stratonovich prescription. It is true that Ito integrals have Martingale properties not possessed by Stratonovich integrals [40] which leads to risk-neural theorems for markets [41,42], but the nature of the proper mathematics should eventually be determined by proper aggregation of relatively microscopic models of markets. It should be noted that virtually all investigations of other physical systems, which are also continuous time models of discrete processes, conclude that the Stratonovich interpretation coincides with reality, when multiplicative noise with zero correlation time, modeled in terms of white noise $η^{j}$ , is properly considered as the limit of real noise with finite correlation time [43]. The path integral succinctly demonstrates the difference between the two: The Ito prescription corresponds to the prepoint discretization of $L$ , wherein $θ \dot{M} (t) \to M_{ρ + 1} - M_{ρ}$ and $M (t) \to M_{ρ}$ The Stratonovich prescription corresponds to the midpoint discretization of $L$ , wherein $θ \dot{M} (t) \to M_{ρ + 1} - M_{ρ}$ and $M (t) \to \frac{1}{2} (M_{ρ + 1} + M_{ρ})$ In terms of the functions appearing in the Fokker-Planck Eq. (21), the Ito prescription of the prepoint discretized Lagrangian, $L_{I}$ , is relatively simple, albeit deceptively so because of its nonstandard calculus.

L_{I} ({\dot{M}}^{G}, M^{G}, t) = \frac{1}{2} ({\dot{M}}^{G} - g^{G}) g_{G G'} ({\dot{M}}^{G'} - g^{G'}) - V .

(24)

In the absence of a nonphenomenological microscopic theory, the difference between a Ito prescription and a Stratonovich prescription is simply a transformed drift [37].
There are several other advantages to Eq. (22) over Eq. (20). Extrema and most probable states of $M^{G}$ , $≪ M^{G} ≫$ , are simply derived by a variational principle, similar to conditions sought in previous studies [44]. In the Stratonovich prescription, necessary, albeit not sufficient, conditions are given by

δ_{G} L = L_{, G} - L_{, \dot{G} : t} = 0,

L_{, \dot{G} : t} = L_{, \dot{G} G'} {\dot{M}}^{G'} + L_{, \dot{G} \dot{G}'} {\overset{..}{M}}^{G'} .

(25)

For stationary states, ${\dot{M}}^{G} = 0$ , and $\partial \bar{L} / \partial {\bar{M}}^{G} = 0$ defines $≪ {\bar{M}}^{G} ≫$ , where the bars identify stationary variables; in this case, the macroscopic variables are equal to their mesoscopic counterparts. [Note that $\bar{L}$ is not the stationary solution of the system, e.g., to Eq. (21) with $\partial P / \partial t = 0$ However, in some cases [45], $\bar{L}$ is a definite aid to finding such stationary states.] Many times only properties of stationary states are examined, but here a temporal dependence is included. E.g., the ${\dot{M}}^{G}$ terms in $L$ permit steady states and their fluctuations to be investigated in a nonequilibrium context. Note that Eq. (25) must be derived from the path integral, Eq. (22), which is at least one reason to justify its development.

4.2. Correlations

Correlations between variables are modeled explicitly in the Lagrangian as a parameter usually designated $ρ$ (not to be confused with the Rho Greek calculated for options). This section uses a simple two-factor model to develop the correspondence between the correlation $ρ$ in the Lagrangian and that among the commonly written Weiner distributions $d z$
Consider coupled stochastic differential equations

d r = f^{r} (r, l) d t + {\hat{g}}^{r} (r, l) σ_{1} d z_{1},

d l = f^{l} (r, l) d t + {\hat{g}}^{l} (r, l) σ_{2} d z_{2},

< d z_{i} > = 0, i = {1, 2},

< d z_{i} (t) d z_{j} (t') > = d t δ (t - t'), i = j,

< d z_{i} (t) d z_{j} (t') > = ρ d t δ (t - t'), i \neq j,

δ (t - t') = {\begin{array}{c} 0,, \\ 1, \end{array} \begin{array}{c} t \neq t', \\ t = t', \end{array} (null)

(26)

where $< . >$ denotes expectations.
These can be rewritten as Langevin equations (in the It^ prepoint discretization)

d r / d t = f^{r} + {\hat{g}}^{r} σ_{1} (γ^{+} n_{1} + s g n ρ γ^{-} n_{2}),

d l / d t = g^{l} + {\hat{g}}^{l} σ_{2} (s g n ρ γ^{-} n_{1} + γ^{+} n_{2}),

γ^{\pm} = \frac{1}{\sqrt{2}} [1 \pm (1 - ρ^{2})^{1 / 2}]^{1 / 2},

n_{i} = (d t)^{1 / 2} p_{i},

(27)

where $p_{1}$ and $p_{2}$ are independent [0,1] Gaussian distributions.
The equivalent short-time probability distribution, $P$ , for the above set of equations is

P = g^{1 / 2} (2 π d t)^{- 1 / 2} exp (- L d t),

L = \frac{1}{2} F^{†} \underline{g} F,

F = (\begin{matrix} d r / d t - f^{r}) \\ d l / d t - f^{l}) \end{matrix}),

g = det (\underline{g}),

k = 1 - ρ^{2} .

(28)

\underline{g}

, the metric in ${r, l}$ -space, is the inverse of the covariance matrix,

{\underline{g}}^{- 1} = (\begin{matrix} ({\hat{g}}^{r} σ_{1})^{2} & ρ {\hat{g}}^{r} {\hat{g}}^{l} σ_{1} σ_{2} \\ ρ {\hat{g}}^{r} {\hat{g}}^{l} σ_{1} σ_{2} & ({\hat{g}}^{l} σ_{2})^{2} \end{matrix}) .

(29)

The above also corrects previous papers which inadvertently dropped the $s g n$ factors in the above [9,10,17].

5. ADAPTIVE SIMULATED ANNEALING (ASA) FITS

5.1. ASA Outline

The algorithm developed which is now called Adaptive Simulated Annealing (ASA) [46] fits short-time probability distributions to observed data, using a maximum likelihood technique on the Lagrangian. This algorithm has been developed to fit observed data to a theoretical cost function over a $D$ -dimensional parameter space [46], adapting for varying sensitivities of parameters during the fit. The ASA code can be obtained at no charge, via WWW from http://www.ingber.com/ or via FTP from ftp.ingber.com [15].

5.1.1. General description

Simulated annealing (SA) was developed in 1983 to deal with highly nonlinear problems [47], as an extension of a Monte-Carlo importance-sampling technique developed in 1953 for chemical physics problems. It helps to visualize the problems presented by such complex systems as a geographical terrain. For example, consider a mountain range, with two “parameters,” e.g., along the North−South and East−West directions. We wish to find the lowest valley in this terrain. SA approaches this problem similar to using a bouncing ball that can bounce over mountains from valley to valley. We start at a high “temperature,” where the temperature is an SA parameter that mimics the effect of a fast moving particle in a hot object like a hot molten metal, thereby permitting the ball to make very high bounces and being able to bounce over any mountain to access any valley, given enough bounces. As the temperature is made relatively colder, the ball cannot bounce so high, and it also can settle to become trapped in relatively smaller ranges of valleys.
We imagine that our mountain range is aptly described by a “cost function.” We define probability distributions of the two directional parameters, called generating distributions since they generate possible valleys or states we are to explore. We define another distribution, called the acceptance distribution, which depends on the difference of cost functions of the present generated valley we are to explore and the last saved lowest valley. The acceptance distribution decides probabilistically whether to stay in a new lower valley or to bounce out of it. All the generating and acceptance distributions depend on temperatures.
In 1984 [48], it was established that SA possessed a proof that, by carefully controlling the rates of cooling of temperatures, it could statistically find the best minimum, e.g., the lowest valley of our example above. This was good news for people trying to solve hard problems which could not be solved by other algorithms. The bad news was that the guarantee was only good if they were willing to run SA forever. In 1987, a method of fast annealing (FA) was developed [49], which permitted lowering the temperature exponentially faster, thereby statistically guaranteeing that the minimum could be found in some finite time. However, that time still could be quite long. Shortly thereafter, Very Fast Simulated Reannealing (VFSR) was developed in 1987 [46], now called Adaptive Simulated Annealing (ASA), which is exponentially faster than FA.
ASA has been applied to many problems by many people in many disciplines [50-52]. The feedback of many users regularly scrutinizing the source code ensures its soundness as it becomes more flexible and powerful.

5.1.2. Mathematical outline

ASA considers a parameter $α_{k}^{i}$ in dimension $i$ generated at annealing-time $k$ with the range

α_{k}^{i} \in [A_{i}, B_{i}],

(30)

calculated with the random variable $y^{i}$ ,

α_{k + 1}^{i} = α_{k}^{i} + y^{i} (B_{i} - A_{i}),

y^{i} \in [- 1, 1] .

(31)

The generating function $g_{T} (y)$ is defined,

g_{T} (y) = Π_{i = 1}^{D} \frac{1}{2 (| y^{i} | + T_{i}) ln (1 + 1 / T_{i})} \equiv Π_{i = 1}^{D} g_{T}^{i} (y^{i}),

(32)

where the subscript $i$ on $T_{i}$ specifies the parameter index, and the $k$ -dependence in $T_{i} (k)$ for the annealing schedule has been dropped for brevity. Its cumulative probability distribution is

G_{T} (y) = \int_{- 1}^{y^{1}} . . . \int_{- 1}^{y^{D}} d y'^{1} . . . d y'^{D} g_{T} (y') \equiv Π_{i = 1}^{D} G_{T}^{i} (y^{i}),

G_{T}^{i} (y^{i}) = \frac{1}{2} + \frac{s g n (y^{i})}{2} \frac{ln (1 + | y^{i} | / T_{i})}{ln (1 + 1 / T_{i})} .

(33)

y^{i}

is generated from a $u^{i}$ from the uniform distribution

u^{i} \in U [0, 1],

y^{i} = s g n (u^{i} - \frac{1}{2}) T_{i} [(1 + 1 / T_{i})^{| 2 u^{i} - 1 |} - 1] .

(34)

It is straightforward to calculate that for an annealing schedule for

T_{i}

T_{i} (k) = T_{0 i} exp (- c_{i} k^{1 / D}),

(35)

a global minima statistically can be obtained. I.e.,

Σ_{k_{0}}^{\infty} g_{k} ≈ Σ_{k_{0}}^{\infty} [Π_{i = 1}^{D} \frac{1}{2 | y^{i} | c_{i}}] \frac{1}{k} = \infty .

(36)

Control can be taken over $c_{i}$ , such that

T_{f i} = T_{0 i} exp (- m_{i}) w h e n k_{f} = exp n_{i},

c_{i} = m_{i} exp (- n_{i} / D),

(37)

where $m_{i}$ and $n_{i}$ can be considered “free” parameters to help tune ASA for specific problems.
ASA has over 100 OPTIONS available for tuning. A few important ones were used in this project.

5.1.3. Reannealing

Whenever doing a multi-dimensional search in the course of a complex nonlinear physical problem, inevitably one must deal with different changing sensitivities of the $α^{i}$ in the search. At any given annealing-time, the range over which the relatively insensitive parameters are being searched can be “stretched out” relative to the ranges of the more sensitive parameters. This can be accomplished by periodically rescaling the annealing-time $k$ , essentially reannealing, every hundred or so acceptance-events (or at some user-defined modulus of the number of accepted or generated states), in terms of the sensitivities $s_{i}$ calculated at the most current minimum value of the cost function, $C$ ,

s_{i} = \partial C / \partial α^{i} .

(38)

In terms of the largest $s_{i}$ = $s_{max}$ , a default rescaling is performed for each $k_{i}$ of each parameter dimension, whereby a new index $k'_{i}$ is calculated from each $k_{i}$ ,

k_{i} \to k'_{i},

T'_{i k'} = T_{i k} (s_{max} / s_{i}),

k'_{i} = (ln (T_{i 0} / T_{i k'}) / c_{i})^{D} .

(39)

T_{i 0}

is set to unity to begin the search, which is ample to span each parameter dimension.

5.1.4. Quenching

Another adaptive feature of ASA is its ability to perform quenching in a methodical fashion. This is applied by noting that the temperature schedule above can be redefined as

T_{i} (k_{i}) = T_{0 i} exp (- c_{i} k_{i}^{Q_{i} / D}),

c_{i} = m_{i} exp (- n_{i} Q_{i} / D),

(40)

in terms of the “quenching factor” $Q_{i}$ The sampling proof fails if $Q_{i} > 1$ as

\underset{k}{Σ} \overset{D}{Π} 1 / k^{Q_{i} / D} = \underset{k}{Σ} 1 / k^{Q_{i}} < \infty .

(41)

This simple calculation shows how the “curse of dimensionality” arises, and also gives a possible way of living with this disease. In ASA, the influence of large dimensions becomes clearly focussed on the exponential of the power of $k$ being $1 / D$ , as the annealing required to properly sample the space becomes prohibitively slow. So, if resources cannot be committed to properly sample the space, then for some systems perhaps the next best procedure may be to turn on quenching, whereby $Q_{i}$ can become on the order of the size of number of dimensions.
The scale of the power of $1 / D$ temperature schedule used for the acceptance function can be altered in a similar fashion. However, this does not affect the annealing proof of ASA, and so this may used without damaging the sampling property.

5.2. -Indicator of Market Contexts

Our studies of contexts of markets well recognized by option traders to have significantly different volatility behavior show that the exponents $x$ are reasonably faithful indicators defining these different contexts.
We feel the two-factor model is more accurate because the data indeed demonstrate stochastic volatility [4]. We also note that the two-factor $x$ fit by ASA across the last few years. This is not true of the one-factor ASA fitted $x$ the Black-Scholes $σ$ as a parameter, but rather calculate as historical volatility during all runs. Some results of two-factor studies and one-factor studies using a Black-Scholes $σ$ have been reported elsewhere [13].
Since $σ$ is not widely traded and arbitraged, to fit the two-factor model, we calculate this quantity as an historical volatility for both its prepoint and postpoint values. Some previous studies used a scaled implied volatility (which is calculated from a BS model). We use a standard deviation $σ'$ ,

σ' = S t d D e v (d S / F (S, S_{0}, S_{\infty}, x, y)) .

(42)

In the one-factor model, it does not make good numerical sense to have two free parameters in one term, i.e., $σ$ and $x$ , as these cannot be fit very well within the variance the data. Instead, one method is to take guidance from the two-factor results, to set a scale for an effective $σ$ , and then fit the parameter $x$ Another method it apply the above StdDev as a proxy for $σ$ Some motivation for this approach is given by considering collapsing a two-factor stochastic volatility model in one-factor model: The one-factor model now has an integral over the stochastic process in its diffusion term. The is integral is what we are approximating by using the standard deviation of a moving window of the data.

6. PATH-INTEGRAL (PATHINT) DEVELOPMENT

6.1. PATHINT Outline

The fits described above clearly demonstrate the need to incorporate stochastic volatility in option pricing models. If one-factor fits are desired, e.g., for efficiency of calculation, then at the least the exponent of price $x$ should be permitted to freely adapt to the data. In either case, it is required to develop a full set of Greeks for trading. To meet these needs, we have used a path-integral code, PATHINT, described below, with great success. At this time, the two-factor code takes too long to run for daily use, but it proves to be a good weekly baseline for the one-factor code.
The PATHINT algorithm develops the long-time probability distribution from the Lagrangian fit by the first optimization code. A robust and accurate histogram-based (non-Monte Carlo) path-integral algorithm to calculate the long-time probability distribution has been developed to handle nonlinear Lagrangians [18-20,22,53-55],
The histogram procedure recognizes that the distribution can be numerically approximated to a high degree of accuracy as sum of rectangles at points $M_{i}$ of height $P_{i}$ and width $Δ M_{i}$ For convenience, just consider a one-dimensional system. The above path-integral representation can be rewritten, for each of its intermediate integrals, as

P (M; t + Δ t) = \int d M' [g_{s}^{1 / 2} (2 π Δ t)^{- 1 / 2} exp (- L_{s} Δ t)] P (M'; t) = \int d M' G (M, M'; Δ t) P (M'; t),

P (M; t) = Σ_{i = 1}^{N} π (M - M_{i}) P_{i} (t),

π (M - M_{i}) = {\begin{array}{l} 0, (M_{i} - \frac{1}{2} Δ M_{i - 1}) \leq M \leq (M_{i} + \frac{1}{2} Δ M_{i}), \\ 1, o t h e r w i s e, \end{array}, (null)

(43)

which yields

P_{i} (t + Δ t) = T_{i j} (Δ t) P_{j} (t),

T_{i j} (Δ t) = \frac{2}{Δ M_{i - 1} + Δ M_{i}} \int_{M_{i} - Δ M_{i - 1} / 2}^{M_{i} + Δ M_{i} / 2} d M \int_{M_{j} - Δ M_{j - 1} / 2}^{M_{j} + Δ M_{j} / 2} d M' G (M, M'; Δ t) .

(44)

T_{i j}

is a banded matrix representing the Gaussian nature of the short-time probability centered about the (varying) drift.
Care must be used in developing the mesh in $Δ M^{G}$ , which is strongly dependent on the diagonal elements of the diffusion matrix, e.g.,

Δ M^{G} ≈ (Δ t g^{| G | | G |})^{1 / 2} .

(45)

Presently, this constrains the dependence of the covariance of each variable to be a nonlinear function of that variable, albeit arbitrarily nonlinear, in order to present a straightforward rectangular underlying mesh. Below we address how we have handled this problem in our two-factor stochastic-volatility model.
Fitting data with the short-time probability distribution, effectively using an integral over this epoch, permits the use of coarser meshes than the corresponding stochastic differential equation. The coarser resolution is appropriate, typically required, for numerical solution of the time-dependent path-integral: By considering the contributions to the first and second moments of $Δ M^{G}$ for small time slices $θ$ , conditions on the time and variable meshes can be derived [53]. The time slice essentially is determined by $θ \leq {\bar{L}}^{- 1}$ , where $\bar{L}$ is the “static” Lagrangian with $d M^{G} / d t = 0$ , throughout the ranges of $M^{G}$ giving the most important contributions to the probability distribution $P$ The variable mesh, a function of $M^{G}$ , is optimally chosen such that $Δ M^{G}$ is measured by the covariance $g^{G G'}$ , or $Δ M^{G} ≈ (g^{G G} θ)^{1 / 2}$
If the histogram method is further developed into a trapezoidal integration, then more accuracy can be expected at more complex boundaries [54]. Such problems does not arise here, and 6−7 significant figure accuracy is easily achieved provided great care is taken to develop the mesh as described above. For example, after setting the initial-condition discretized delta function $M S U P G$ at the prepoint of an interval, the mesh going forward in $M^{G}$ is simply calculated stepwise using

Δ M^{G} = (g^{G G} θ)^{1 / 2} .

(46)

However, going backwards in $M^{G}$ , an iterative procedure was used at each step, starting with an estimate from the prepoint and going forward again, until there was no mismatch. That much care is required for the mesh was observed in the original Wehner-Wolfer paper [53].
It is important to stress that very good numerical accuracy is required to get very good Greeks required for real-world trading. Many authors develop very efficient numerical schemes to get reasonable prices to 2 or 3 significant figures, but these methods often are not very good to enough significant figures to get good precision for the Greeks. Typical Monte Carlo methods are notorious for giving such poor results after very long computer runs. In particular, we do not believe that good Greeks required for trading can be obtained by using meshes obtained by other simpler numerical algorithms [56].
The PATHINT algorithm in its present form can “theoretically” handle any n-factor model subject to its diffusion-mesh constraints. In practice, the calculation of 3-factor and 4-factor models likely will wait until giga-hertz speeds and giga-byte RAM are commonplace.

6.2. Development of Long-Time Probabilities

The noise determined empirically as the diffusion of the data is the same, independent of $x$ within our approach. Therefore, we scale different exponents such that the the diffusions, the square of the “basis-point volatilities” (BPV), are scaled to be equivalent. Then, there is not a very drastic change in option prices for different exponents $x$ for the strike $X$ set to the $S$ underlying, the at-the-money (ATM) strike. This is not the case for out of the money (OTM) or in the money (ITM) strikes, e.g., when exercising the strike would generate loss or profit, resp. This implies that current pricing models are not radically mispricing the markets, but there still are significant changes in Greeks using more sophisticated models.

6.3. Dependence of Probabilities on and

Fig. 2. gives examples of the short-time distribution evolved out to $T = 0.5$ year for $x$ in {-1, 0, 1, 2}, with 500 intermediate epochs/foldings, and BS $σ = 0.0075$ Each calculation scales $σ$ by multiplying by $S / F (S, S_{0}, S_{\infty}, x, y)$
Fig. 3. gives an example of a two-factor distribution evolved out to $T = 0.5$ year for $x = 0.7$

6.4. Two-Factor Volatility and PATHINT Modifications

In our two-factor model, the mesh of $S$ would depend on $σ$ and cause some problems in any PATHINT grid to be developed in $S$ - $σ$
For some time we have considered how to handle this generic problem for $n$ -factor multivariate systems with truly multivariate diffusions within the framework of PATHINT. In one case, we have taken advantage of the Riemannian invariance of the probability distribution as discussed above, to transform to a system where the diffusions have only “diagonal” multiplicative dependence [19]. However, this leads to cumbersome numerical problems with the transformed boundary conditions [20]. Another method, not yet fully tested, is to develop a tiling of diagonal meshes for each factor $i$ that often are suitable for off-diagonal regions in an $n$ -factor system, e.g.,

Δ M_{k}^{i} = 2^{m_{k}^{i}} Δ M_{0}^{i},

Δ M_{0}^{i} ≈ \sqrt{g_{k_{0}}^{| i | | i |} Δ t},

(47)

where the mesh of variable $i$ at a given point labeled by $k$ is an exponentiation of 2, labeled by $m_{k}^{i}$ ; the integral power $m_{k}^{i}$ is determined such that it gives a good approximation to the diagonal mesh given by the one-factor PATHINT mesh conditions, in terms of some minimal mesh $Δ M_{0}^{i}$ , throughout regions of the Lagrangian giving most important contributions to the distribution as predetermined by a scan of the system. This tiling of the kernel is to be used together with interpolation of intermediate distributions.
The results of our study here are that, after the at-the-money BPV are scaled to be equivalent, there is not a very drastic change in the one-factor ATM Greeks developed below. Therefore, while we have not at all changed the functional dependence of the Lagrangian on $S$ and $σ$ , we have determined our meshes using a diffusion for the $S$ equation as $σ_{0} F (S, S_{0}, S_{\infty}, x, y)$ , where $σ_{0}$ is determined by the same BPV-equivalent condition as imposed on the one-factor models. This seems to work very well, especially since we have taken our $σ$ equation to be normal with a limited range of influence in the calculations. Future work yet has to establish a more definitive distribution for $σ$

7. CALCULATION OF DERIVATIVES

7.1. Primary Use of Probabilities For European Options

We can use PATHINT to develop the distribution of the option value back in time from expiration. This is the standard approach used by CRR, explicit and implicit Crank-Nicolson models, etc [57].
For European options, we also take advantage of the accuracy of PATHINT enhanced by normalizing the distribution as well as the kernel at each iteration (though in these studies this was not required after normalizing the kernel). Therefore, we have calculated our European option prices and Greeks using the most elementary and intuitive definition of the option’s price $V$ [58], which is the expected value

V = < max [z (S - X), 0] >, {\begin{matrix} z = 1, & c a l l \\ z = - 1, & p u t \end{matrix}, (null)

(48)

where $X$ is the strike price, and the expected value $< . >$ is taken with respect to the risk-neutral distribution of the underlying market $S$ It should be noted that, while the standard approach of developing the option price delivers at the present time a range of underlying values for a given strike, our approach delivers a more practical complete range of strikes (as many as 50−60 for Eurodollar options) for a given underlying at the present time, resulting in a greatly enhanced numerical efficiency. The risk-neutral distribution is effectively calculated taking the drift as the cost-of-carry $b$ times $S$ , using the above arguments leading to the BS formula. We have designed our codes to use parameters risk-free-rate $r$ and cost-of-carry $b$ such that

b = r, c o s t o f c a r r y o n n o n d i v i d e n d s t o c k,

b = r - q, c o s t o f c a r r y o n s t o c k p a y i n g d i v i d e n d y i e l d q,

b = 0, c o s t o f c a r r y o n f u t u r e c o n t r a c t,

b = r - r_{f}, c o s t o f c a r r y o n c u r r e n c y w i t h r a t e r_{f},

(49)

which is similar to how generalized European BS codes use $b$ and $r$ [59].
Using this approach, the European price $V_{E}$ is calculated as

V_{E} = < max [z (e^{(b - r) T} S_{T} - e^{- r T} X), 0] > .

(50)

The American price $V_{A}$ must be calculated using a different kernel going back in time from expiration, using as “initial conditions” the option values used in the above average. This kernel is the transposed matrix used for the European case, and includes additional drift and “potential” terms due to the need to develop this back in time. This can be understood as requiring the adjoint partial differential equation or a postpoint Lagrangian in real time.
That is, a forward equation for the conditional probability distribution $P [M (t + d t), t + d t | M (t), t]$ is

\partial P / \partial t = \frac{1}{2} (g^{G G'} P)_{, G G'} - (g^{G} P)_{, G} + V,

(51)

where the partial derivatives with respect to $M^{G}$ act on the postpoint $M^{G} (t + d t)$ A backward equation for the conditional probability distribution $P [M (t + d t), t + d t | M (t), t]$ is

\partial P / \partial t = \frac{1}{2} g^{G G'} P_{, G G'} - g^{G} P_{, G} + V,

(52)

where the partial derivatives with respect to $M^{G}$ act on the prepoint $M^{G} (t)$ The forward equation has a particularly simple form in the mathematically equivalent prepoint path integral representation.
Above, we have described how the forward distribution at present time $T_{0}$ is evolved using PATHINT to the time of expiration, $P (T)$ , e.g., using a path-integral kernel $K$ folded over $n$ epochs, where it is folded with the function $O$ ,

V = < O (T) P (T) > = < O (K^{n} P (t_{0})) >,

O (T) = max [z (e^{(b - r) T} S_{T} - e^{- r T} X), 0],

(53)

to determine the European value at the present time of the calls and puts at different strike values $X$ An equivalent calculation can be performed by using the backward equation, expressed in terms of the “equivalent” kernel $K^{†}$ acting on $O$ ,

V = < O (T) P (T) > = < (K^{n} O (T)) P (t_{0}) > .

(54)

It is convenient to use the simple prepoint representation for the Lagrangian, so the backward kernel is first re-expressed as a forward kernel by bringing the diffusions and drifts inside the partial derivatives, giving a transformed adjoint kernel $K^{†}$
The above mathematics is easily tested by calculating European options going forwards and backwards. For American options, while performing the backwards calculation, at each point of mesh, the options price evolving backwards from $T$ is tested and calculated as

O_{n e w} = max [(S - X), (e^{- r d t} O_{o l d})] .

(55)

The Greeks ${Δ, Γ, Θ}$ are directly taken off the final developed option at the present time, since the $M^{G}$ mesh is available for all derivatives. We get excellent results for all Greeks. Note that for CRR trees, only one point of mesh is at the present time, so $Δ$ requires moving back one epoch and $Γ$ requires moving back two epochs, unless the present time is pushed back additional epochs, etc.

7.2. PATHINT Baselined to CRR and BS

The CRR method is a simple binomial tree which in a specific limit approaches the BS partial differential equation. It has the virtues of being fast and readily accommodates European and American calculations. However, it suffers a number of well-known numerical problems, e.g., a systematic bias in the tree approximation and an oscillatory error as a function of the number of intermediate epochs/iterations in its time mesh. Some Greeks like ${Δ, Γ, Θ}$ can be directly taken off the tree used for pricing with reasonable approximations (at epochs just before the actual current time). The first problem for American options can be alleviated somewhat by using the variant method [5],

C R R_{v a r i a n t} = C R R_{A m e r i c a n} - C R R_{E u r o p e a n} + B S .

(56)

The second problem can be alleviated somewhat by averaging runs of $n$ iterations with runs of $n + 1$ iterations [60]. This four-fold increase of runs is rarely used, though perhaps it should be more often. Furthermore, if increased accuracy in price is needed in order to take numerical derivatives, typically 200−300 iterations should be used for expirations some months away, not 30−70 as too often used in practice.
When taking numerical derivatives there can arise a need to tune increments taken for the differentials. For some Greeks like $Δ$ and $Γ$ the size of the best differentials to use may vary with strikes that have different sensitivities to parts of the underlying distribution. One method of building in some adaptive flexibility across many such strikes is to increase the order of terms used in taking numerical derivatives. (This was not required for results presented here.) For example, it is straightforward to verify that, while the central difference

\frac{d f}{d x} = \frac{f (x + d x) - f (x - d x)}{2 d x}

(57)

is typically good to $o ((d x)^{3})$ ,

\frac{d f}{d x} = \frac{- f (x + 2 d x) + 8 f (x + d x) - 8 f (x - d x) + f (x - 2 d x)}{12 d x}

(58)

is typically good to $o ((d x)^{5})$ Similarly, while

\frac{d^{2} f}{d x^{2}} = \frac{f (x + d x) - 2 f (x) + f (x - d x)}{d x d x}

(59)

is typically good to $((d x)^{4})$ ,

\frac{d^{2} f}{d x^{2}} = \frac{- d (x + 2 d x) + 16 f (x + d x) - 30 f (x) + 16 f (x - d x) - f (x - 2 d x)}{d x d x}

(60)

is typically good to $((d x)^{6})$
Table 1 gives an example of baselining our one-factor PATHINT code to the CRR and BS codes using the above safeguards for an option whose American price is the same as its European counterpart, a typical circumstance [59]. In the literature, the CRR code is most often taken as the benchmark for American calculations. We took the number of intermediate epochs/points to be 300 for each calculation. Parameters used for this particular ATM call are $T = 0.5$ years, $r = 0.05$ , $b = 0$ , $σ = 0.10$

Table 1.

Tests with American CRR and American PATHINT lead to results with the same degrees of accuracy.

7.3. Two-Factor PATHINT Baselined to One-Factor PATHINT

Previous papers and tests have demonstrated that the two-factor PATHINT code performs as expected. The code was developed with only a few lines to be changed for running any $n$ -factor problem, i.e., of course after coding the Lagrangian specific to a given system. Tests were performed by combining two one-factor problems, and there is no loss of accuracy. However, here we are making some additional mesh approximations as discussed above to accommodate $σ$ in the $S$ diffusion. This seems quite reasonable, but there is no sure test of the accuracy. We indeed see that the ATM results are very close across $x$ similar to our ATM comparisons between BS and our one-factor PATHINT results for various $x$ scaling is performed to have all models used the same BPV (using the $σ_{0}$ procedure for the mesh as described above for the two-factor model).
The logical extension of Greeks for the two-factor model is to develop derivatives of price with respect to $ρ$ and $ε$ in $σ$ volatility equation. However, we did not find a bona fide two-factor proxy for the one-factor $Υ$ , the derivative of price with respect to the one-factor $σ$ constant. We get very good ATM $Υ$ comparisons between BS and our one-factor models with various $x$ We tried simply multiplying the noise in the two-factor stochastic volatility in the price equation by a parameter with deviations from 1 to get numerical derivatives of PATHINT solutions, and this gave somewhat good agreement with the ATM BPV-scaled BS $Υ$ within a couple of significant figures. Perhaps this is not too surprising, especially given the correlation substantial $ρ$ between the price and volatility equations which we do not neglect.

8. CONCLUSIONS

The results of our study are that, after the at-the-money basis-point volatilities are scaled to be equivalent, there is only a very small change in option prices for different exponents $x$ , both for the one-factor and two-factor models. There still are significant differences in Greeks using more sophisticated models, especially for out-of-the-money options. This implies that current pricing models are not radically mispricing the markets.
Our studies point to contexts of markets well recognized by option traders to have significantly different volatility behavior. Suppression of stochastic volatility in the one-factor model just leaks out into stochasticity of parameters in the model, e.g., especially in $x$ , unless additional numerical methods are employed, e.g., using an adaptive standard deviation. Our studies show that the two-factor exponents $x$ are reasonably faithful indicators defining these different contexts. The $x$ -exponents in the two-factor fits are quite stable. As such, especially the two-factor $x$ can be considered as a “context indicator” over a longer time scale than other indicators typically used by traders. The two-factor fits also exhibit differences due to the $σ$ parameters, including the $ρ$ correlations, in accord with the sense traders have about the nature of changing volatilities across this time period.
Modern methods of developing multivariate nonlinear multiplicative Gaussian-Markovian systems are quite important, as there are many such systems and the mathematics must be diligently exercised if such models are to faithfully represent the data they describe. Similarly, sophisticated numerical techniques, e.g., global optimization and path integration are important tools to carry out the modeling and fitting to data without compromising the model, e.g., by unwarranted quasi-linear approximations. Three quite different systems have benefited from this approach:
The large-scale modeling of neocortical interactions has benefited from the use of intuitive constructs that yet are faithful to the complex algebra describing this multiple-scaled complex system. For example, canonical-momenta indicators have been successfully applied to multivariate financial markets.
It is clear that ASA optimization and PATHINT path-integral tools are very useful to develop the algebra of statistical mechanics for a large class of nonlinear stochastic systems encountered in finance. However, it also is clear that each system typically presents its own non-typical unique character and this must be included in any such analysis. A virtue of this statistical mechanics approach and these associated tools is they appear to be flexible and robust to handle quite different systems.

ACKNOWLEDGMENTS

I thank Donald Wilson for his support and Jennifer Wilson for collaboration running ASA and PATHINT calculations. Implied volatility and yield data was extracted from the MIM database of Logical Information Machines (LIM).

REFERENCES

1.		L. Ederington and W. Guan, Is implied volatility an informationally efficient and effective predictor of future volatility?, U Oklahoma, Norman, OK, (1998).
2.		G. Bakshi, C. Cao, and Z. Chen, Pricing and hedging long-term options, Pennsylvania State U, University Park, PA, (1998).
3.		L. Ingber, Some Applications of Statistical Mechanics of Financial Markets, LIR-98-1-SASMFM, Lester Ingber Research, Chicago, IL, (1998).
4.		L. Ingber and J.K. Wilson, Volatility of volatility of financial markets, Mathl. Computer Modelling 29 (5), 39-57 (1999).
5.		J.C. Hull, Options, Futures, and Other Derivatives, Third Edition, Prentice Hall, Upper Saddle River, NJ, (1997).
6.		J.C. Jackwerth, Recovering risk aversion from option prices and realized returns, London Business School, London, UK, (1998).
7.		F.A. Longstaff and E.S. Schwartz, Valuing American options by simulation: A simple least-squares approach, Capital Management Sciences, Los Angeles, CA, (1998).
8.		L. Ingber, Statistical mechanics of nonlinear nonequilibrium financial markets, Math. Modelling 5 (6), 343-361 (1984).
9.		L. Ingber, Statistical mechanical aids to calculating term structure models, Phys. Rev. A 42 (12), 7057-7064 (1990).
10.		L. Ingber, M.F. Wehner, G.M. Jabbour, and T.M. Barnhill, Application of statistical mechanics methodology to term-structure bond-pricing models, Mathl. Comput. Modelling 15 (11), 77-98 (1991).
11.		L. Ingber, Statistical mechanics of nonlinear nonequilibrium financial markets: Applications to optimized trading, Mathl. Computer Modelling 23 (7), 101-121 (1996).
12.		L. Ingber, Canonical momenta indicators of financial markets and neocortical EEG, in Progress in Neural Information Processing, (Edited by S.-I. Amari, L. Xu, I. King, and K.-S. Leung), pp. 777-784, Springer, New York, (1996).
13.		L. Ingber and J.K. Wilson, Statistical mechanics of financial markets: Exponential modifications to Black-Scholes, Mathl. Computer Modelling 31 (8/9), 167-192 (2000).
14.		C. Azariadis, Self-fulfilling prophecies, J. Econ. Theory 25, 380-396 (1981).
15.		L. Ingber, Adaptive Simulated Annealing (ASA), Global optimization C-code, Caltech Alumni Association, Pasadena, CA, (1993).
16.		J. C. Cox, S. A. Ross, and M. Rubenstein, Option pricing: A simplified approach, J. Fin. Econ. 7, 229-263 (1979).
17.		L. Ingber, Data mining and knowledge discovery via statistical mechanics in nonlinear stochastic systems, Mathl. Computer Modelling 27 (3), 9-31 (1998).
18.		L. Ingber, H. Fujio, and M.F. Wehner, Mathematical comparison of combat computer models to exercise data, Mathl. Comput. Modelling 15 (1), 65-90 (1991).
19.		L. Ingber, Statistical mechanics of neocortical interactions: Path-integral evolution of short-term memory, Phys. Rev. E 49 (5B), 4652-4664 (1994).
20.		L. Ingber and P.L. Nunez, Statistical mechanics of neocortical interactions: High resolution path-integral calculation of short-term memory, Phys. Rev. E 51 (5), 5074-5083 (1995).
21.		L. Ingber, R. Srinivasan, and P.L. Nunez, Path-integral evolution of chaos embedded in noise: Duffing neocortical analog, Mathl. Computer Modelling 23 (3), 43-53 (1996).
22.		L. Ingber, Path-integral evolution of multivariate systems with moderate noise, Phys. Rev. E 51 (2), 1616-1619 (1995).
23.		L. Ingber and R. Mondescu, (in preparation), (2000).
24.		Federal Reserve Bank, Instruments of the Money Markets, Seventh Edition, Federal Reserve Bank of Richmond, Richmond, VA, (1993).
25.		L. Bachelier, Th´orie de la Sp´culation, Annales de l’Ecole Normale Sup*´rieure* 17, 21-86 (1900).
26.		M. C. Jensen, Some anomalous evidence regarding market efficiency, an editorial introduction, J. Finan. Econ. 6, 95-101 (1978).
27.		B. B. Mandelbrot, When can price be arbitraged efficiently? A limit to the validity of the random walk and martingale models, Rev. Econ. Statist. 53, 225-236 (1971).
28.		S. J. Taylor, Tests of the random walk hypothesis against a price-trend hypothesis, J. Finan. Quant. Anal. 17, 37-61 (1982).
29.		P. Brown, A. W. Kleidon, and T. A. Marsh, New evidence on the nature of size-related anomalies in stock prices, J. Fin. Econ. 12, 33-56 (1983).
30.		F.J. Fabozzi, Treasury Securities and Derivatives, Fabozzi Associates, New Hope, PA, (1998).
31.		L. Ingber, Statistical mechanics of neocortical interactions: A scaling paradigm applied to electroencephalography, Phys. Rev. A 44 (6), 4017-4060 (1991).
32.		K.S. Cheng, Quantization of a general dynamical system by Feynman’s path integration formulation, J. Math. Phys. 13, 1723-1726 (1972).
33.		H. Dekker, Functional integration and the Onsager-Machlup Lagrangian for continuous Markov processes in Riemannian geometries, Phys. Rev. A 19, 2102-2111 (1979).
34.		R. Graham, Path-integral methods in nonequilibrium thermodynamics and statistics, in Stochastic Processes in Nonequilibrium Systems, (Edited by L. Garrido, P. Seglar, and P.J. Shepherd), pp. 82-138, Springer, New York, NY, (1978).
35.		F. Langouche, D. Roekaerts, and E. Tirapegui, Discretization problems of functional integrals in phase space, Phys. Rev. D 20, 419-432 (1979).
36.		F. Langouche, D. Roekaerts, and E. Tirapegui, Short derivation of Feynman Lagrangian for general diffusion process, J. Phys. A 113, 449-452 (1980).
37.		F. Langouche, D. Roekaerts, and E. Tirapegui, Functional Integration and Semiclassical Expansions, Reidel, Dordrecht, The Netherlands, (1982).
38.		M. Rosa-Clot and S. Taddei, A path integral approach to derivative security pricing: I. Formalism and analytic results, INFN, Firenze, Italy, (1999).
39.		P. Hagedorn, Non-Linear Oscillations, Oxford Univ., New York, NY, (1981).
40.		B. Oksendal, Stochastic Differential Equations, Springer, New York, NY, (1998).
41.		J.M. Harrison and D. Kreps, Martingales and arbitrage in multiperiod securities markets, J. Econ. Theory 20, 381-408 (1979).
42.		S.R. Pliska, Introduction to Mathematical Finance, Blackwell, Oxford, UK, (1997).
43.		C.W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Springer-Verlag, Berlin, Germany, (1983).
44.		R. C. Merton, An intertemporal capital asset pricing model, Econometrica 41, 867-887 (1973).
45.		L. Ingber, Statistical mechanics of neocortical interactions: Stability and duration of the 7+−2 rule of short-term-memory capacity, Phys. Rev. A 31, 1183-1186 (1985).
46.		L. Ingber, Very fast simulated re-annealing, Mathl. Comput. Modelling 12 (8), 967-973 (1989).
47.		S. Kirkpatrick, C.D. Gelatt, Jr., and M.P. Vecchi, Optimization by simulated annealing, Science 220 (4598), 671-680 (1983).
48.		S. Geman and D. Geman, Stochastic relaxation, Gibbs distribution and the Bayesian restoration in images, IEEE Trans. Patt. Anal. Mac. Int. 6 (6), 721-741 (1984).
49.		H. Szu and R. Hartley, Fast simulated annealing, Phys. Lett. A 122 (3-4), 157-162 (1987).
50.		M. Wofsey, Technology: Shortcut tests validity of complicated formulas, The Wall Street Journal 222 (60), B1 (1993).
51.		L. Ingber, Simulated annealing: Practice versus theory, Mathl. Comput. Modelling 18 (11), 29-57 (1993).
52.		L. Ingber, Adaptive simulated annealing (ASA): Lessons learned, Control and Cybernetics 25 (1), 33-54 (1996).
53.		M.F. Wehner and W.G. Wolfer, Numerical evaluation of path-integral solutions to Fokker-Planck equations. I., Phys. Rev. A 27, 2663-2670 (1983).
54.		M.F. Wehner and W.G. Wolfer, Numerical evaluation of path-integral solutions to Fokker-Planck equations. II. Restricted stochastic processes, Phys. Rev. A 28, 3003-3011 (1983).
55.		M.F. Wehner and W.G. Wolfer, Numerical evaluation of path integral solutions to Fokker-Planck equations. III. Time and functionally dependent coefficients, Phys. Rev. A 35, 1795-1801 (1987).
56.		M. Rosa-Clot and S. Taddei, A path integral approach to derivative security pricing: II. Numerical methods, INFN, Firenze, Italy, (1999).
57.		P. Wilmott, S. Howison, and J. Dewynne, The Mathematics of Financial Derivatives, Cambridge U Press, Cambridge, (1995).
58.		L. Ingber, A simple options training model, Mathl. Computer Modelling 30 (5-6), 167-182 (1999).
59.		E.G. Haug, The Complete Guide to Option Pricing Formulas, McGraw-Hill, New York, NY, (1997).
60.		M. Broadie and J. Detemple, Recent advances in numerical methods for pricing derivative securities, in Numerical Methods in Finance, (Edited by L.C.G Rogers and D. Talay), pp. 43-66, Cambridge University Press, Cambridge, UK, (1997).

FIGURE CAPTIONS

FIG. 1. (a) $F (S, S_{0}, S_{\infty}, x, y)$ for $x = 1$ , the Black-Scholes case. The other parameters are $S = 5$ , $S_{0} = 0.5$ , $S_{\infty} = 20$ , $y = 0$ (b) $F (S, S_{0}, S_{\infty}, x, y)$ for $x = 0$ , the normal distribution. (c) $F (S, S_{0}, S_{\infty}, x, y)$ for $x = - 1$ (d) $F (S, S_{0}, S_{\infty}, x, y)$ for $x = 2$
FIG. 2. The short-time probability distribution at time $T = 0.5$ years for $x = 1$ , the (truncated) Black-Scholes distribution. The short-time probability distribution at time $T = 0.5$ years for $x = 0$ , the normal distribution. The short-time probability distribution at time $T = 0.5$ years for $x = - 1$ The short-time probability distribution at time $T = 0.5$ years for $x = 2$
FIG. 3. A two-factor distribution evolved out to $T = 0.5$ year for $x = 0.7$

TABLE CAPTIONS

Table 1. Calculation of prices and Greeks are given for closed form BS (only valid for European options), binomial tree $C R R_{E u r o p e a n}$ , $C R R_{A m e r i c a n}$ , $C R R_{v a r i a n t}$ , and PATHINT. As verified by calculation, the American option would not be exercised early, so the PATHINT results are identical to the European option. The $C R R_{A m e r i c a n}$ differs somewhat from the $C R R_{E u r o p e a n}$ due to the discrete nature of the calculation. All CRR calculations include averaging over 300 and 301 iterations to minimize oscillatory errors.