

$\stackrel{\text{.}}{\Gamma}=d\Gamma /dt\text{\hspace{0.5em}},$ 


$<\eta \left(t\right){>}_{\eta}=0\text{\hspace{0.5em}},\text{\hspace{0.5em}}<\eta \left(t\right),\eta (t\prime ){>}_{\eta}=\delta (tt\prime )\text{\hspace{0.5em}},$ 

(1)
where${\gamma}_{1}$
and${\gamma}_{2}$
are constants,
and$\Gamma $
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. BlackScholes (BS) Theory
The standard
partialdifferential equation used to formulate most
variants of BlackScholes (BS) models describing the market
value of an
option,$V$
, is 

$\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}+rS\frac{\partial V}{\partial S}rV=0\text{\hspace{0.5em}},$ 

(2)
where$S$
is the asset price,
and$\sigma $
is the standard deviation, or volatility
of$S$
,
and$r$
is the shortterm 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 onedimensional 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)=SN\left({d}_{1}\right)X{e}^{r(Tt)}N\left({d}_{2}\right)\text{\hspace{0.5em}},$ 


${d}_{1}=\frac{\text{ln}\left(S/X\right)+(r+\frac{1}{2}{\sigma}^{2})(Tt)}{\sigma (Tt{)}^{1/2}}\text{\hspace{0.5em}},$ 


${d}_{2}=\frac{\text{ln}\left(S/X\right)+(r\frac{1}{2}{\sigma}^{2})(Tt)}{\sigma (Tt{)}^{1/2}}\text{\hspace{0.5em}}.$ 

(3)
In practice, the
volatility$\sigma $
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}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}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
($\Delta $
), 

$\Pi =V\Delta S\text{\hspace{0.5em}},$ 

(4)
in a market with GaussianMarkovian
(“white”)
noise$X$
and
drift$\mu $
, 

$\frac{dS}{S}=\sigma dX+\mu dt\text{\hspace{0.5em}},$ 

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

$dV=\sigma S\frac{\partial V}{\partial S}dX+\left(\mu S\frac{\partial V}{\partial S}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}+\frac{\partial V}{\partial t}\right)dt\text{\hspace{0.5em}}.$ 

(6)
This yields 

$d\Pi =\sigma \left(\frac{\partial V}{\partial S}\Delta \right)dX+\left(\mu S\frac{\partial V}{\partial S}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}+\frac{\partial V}{\partial t}\mu \Delta S\right)dt\text{\hspace{0.5em}}.$ 

(7)
The expected riskneutral return
of$\Pi $
is 

$d\Pi =r\Pi dt=r(V\Delta S)dt\text{\hspace{0.5em}}.$ 

(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(Tt)}\text{\hspace{0.5em}},$ 

(9)
and setting 

$d\Pi =rV\text{\hspace{0.17em}}dt\text{\hspace{0.5em}}.$ 

(10)
The corresponding BS equation for
futures$F$
is 

$\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma}^{2}{F}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}rV=0\text{\hspace{0.5em}}.$ 

(11)
At least two advantages are present
if$\Delta $
is chosen such that 

$\Delta =\frac{\partial V}{\partial S}\text{\hspace{0.5em}}.$ 

(12)
Then, the portfolio can be instantaneously
“riskneutral,” in terms of zeroing the
coefficient
of$X$
, as well as independent of the direction of market, in
terms of zeroing the coefficient
of$\mu $
For the above example
of$V$
=$C$
, 

$\Delta =N\left({d}_{1}\right)\text{\hspace{0.5em}}.$ 

(13)
Other trading strategies based on this simple model use
similar constructs as risk parameters, e.g., gamma
($\Gamma $
), theta
($\Theta $
), vega
(${\rm Y}$
), rho
($\rho $
) [5], 

$\Gamma =\frac{{\partial}^{2}\Pi}{\partial {S}^{2}}\text{\hspace{0.5em}},$ 


$\Theta =\frac{\partial \Pi}{\partial t}\text{\hspace{0.5em}},$ 


${\rm Y}=\frac{\partial \Pi}{\partial \sigma}\text{\hspace{0.5em}},$ 


$\rho =\frac{\partial \Pi}{\partial r}\text{\hspace{0.5em}}.$ 

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

$\Theta +rS\Delta +\frac{1}{2}(\sigma S{)}^{2}\Gamma =rf\text{\hspace{0.5em}}.$ 

(15)
3.4. Models
Our twofactor model includes
stochastic
volatility$\sigma $
of the
underlying$S$
, 

$dS=\mu \text{\hspace{0.17em}}dt+\sigma \text{\hspace{0.17em}}F(S,{S}_{0},{S}_{\infty},x,y)\text{\hspace{0.17em}}d{z}_{S}\text{\hspace{0.5em}},$ 


$d\sigma =\nu \text{\hspace{0.17em}}dt+\epsilon \text{\hspace{0.17em}}d{z}_{\sigma}\text{\hspace{0.5em}},$ 


$<d{z}_{i}>\text{\hspace{0.17em}}=0\text{\hspace{0.5em}},\text{\hspace{0.5em}}i=\text{{}S,\sigma \text{}}\text{\hspace{0.5em}},$ 


$<d{z}_{i}\left(t\right)\text{\hspace{0.17em}}d{z}_{j}(t\prime )>\text{\hspace{0.17em}}=\text{\hspace{0.17em}}dt\text{\hspace{0.17em}}\delta (tt\prime )\text{\hspace{0.5em}},\text{\hspace{0.5em}}i=j\text{\hspace{0.5em}},$ 


$<d{z}_{i}\left(t\right)\text{\hspace{0.17em}}d{z}_{j}(t\prime )>\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\rho \text{\hspace{0.17em}}dt\text{\hspace{0.17em}}\delta (tt\prime )\text{\hspace{0.5em}},\text{\hspace{0.5em}}i\ne j\text{\hspace{0.5em}},$ 


$F(S,{S}_{0},{S}_{\infty},x,y)=\{\begin{array}{l}S,\\ {S}^{x}{S}_{0}^{1x},\\ {S}^{y}{S}_{0}^{1x}{S}_{\infty}^{xy},\end{array}\begin{array}{l}\text{\hspace{0.5em}}\text{\hspace{0.5em}}S<{S}_{0}\\ \text{\hspace{0.5em}}\text{\hspace{0.5em}}{S}_{0}\le S\le {S}_{\infty}\\ \text{\hspace{0.5em}}\text{\hspace{0.5em}}S>{S}_{\infty}\end{array}\text{\hspace{0.5em}},(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
BlackScholes
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 putcall parity. Putcall 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(Tt)}=p+S\text{\hspace{0.5em}},$ 

(17)
where$c$
($p$
) is the fair price of a call
(put),$X$
is the strike
price,$r$
is the riskfree 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 onefactor model just assumes a
constant$\sigma $
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\le 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 BlackScholes generalized framework after
factoring out interestrate discounting, is equivalent to
using the form 

$dS=\mu \text{\hspace{0.17em}}S\text{\hspace{0.17em}}dt+\sigma \text{\hspace{0.17em}}F(S,{S}_{0},{S}_{\infty},x,y)\text{\hspace{0.17em}}d{z}_{S}\text{\hspace{0.5em}},$ 


$d\sigma =\nu \text{\hspace{0.17em}}dt+\epsilon \text{\hspace{0.17em}}d{z}_{\sigma}\text{\hspace{0.5em}}.$ 

(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\prime $
, a
function$F$
designed to revert to the lognormal BlackScholes model in
several limits, 

$F\prime (S,{S}_{0},{S}_{\infty},x)=S\text{\hspace{0.17em}}{C}_{0}+(1{C}_{0})\text{\hspace{0.17em}}({S}^{x}\text{\hspace{0.17em}}{S}_{0}^{1x}\text{\hspace{0.17em}}{C}_{\infty}+{S}_{0}(1{C}_{\infty}))\text{\hspace{0.5em}},$ 


${C}_{0}=\text{exp}\left[{\left(\frac{S}{{S}_{0}}\text{\hspace{0.17em}}\frac{1x}{1+1x}\right)}^{2x+1}\right]\text{\hspace{0.5em}},$ 


${C}_{\infty}=\text{exp}\left[{\left(\frac{S}{{S}_{\infty}}\right)}^{2}\right]\text{\hspace{0.5em}},$ 


$\underset{S\to \infty ,\text{\hspace{0.17em}}x\ne 1}{\text{lim}}F\prime (S,{S}_{0},{S}_{\infty},x)={S}_{0}=constant\text{\hspace{0.5em}},$ 


$\underset{S\to {0}^{+}}{\text{lim}}F\prime (S,{S}_{0},{S}_{\infty},x)=\underset{x\to 1}{\text{lim}}F\prime (S,{S}_{0},{S}_{\infty},x)=S\text{\hspace{0.5em}}.$ 

(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)\text{\hspace{0.17em}}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 quasiequilibrium
quasilinear 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 rateequations, known as Langevin
equations, 

${\stackrel{\text{.}}{M}}^{G}={f}^{G}+{\widehat{g}}_{j}^{G}{\eta}^{j}\text{\hspace{0.5em}},(G=1,...,\Lambda )\text{\hspace{0.5em}},\text{\hspace{0.5em}}(j=1,...,N)\text{\hspace{0.5em}},\text{\hspace{0.5em}}$ 


${\stackrel{\text{.}}{M}}^{G}=d{M}^{G}/dt\text{\hspace{0.5em}},$ 


$<{\eta}^{j}\left(t\right){>}_{\eta}=0\text{\hspace{0.5em}},\text{\hspace{0.5em}}<{\eta}^{j}\left(t\right),{\eta}^{j\prime}(t\prime ){>}_{\eta}={\delta}^{jj\prime}\delta (tt\prime )\text{\hspace{0.5em}},$ 

(20)
where${f}^{G}$
and${\widehat{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\ge \Lambda $
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 FokkerPlanck
or Sch¨odingertype equation for the conditional
probability
distribution$P\left[M\right(t\left)M\right({t}_{0}\left)\right]$
, the Langevin rate Eq. (20) is developed into the more
useful probability distribution
for${M}^{G}$
at longtime macroscopic time
event$t=(u+1)\theta +{t}_{0}$
, in terms of a Stratonovich pathintegral over mesoscopic
Gaussian conditional probabilities [3236]. Here,
macroscopic variables are defined as the longtime limit of
the evolving mesoscopic system.
The corresponding Sch¨odingertype equation
is [34,35] 

$\partial P/\partial t=\frac{1}{2}({g}^{GG\prime}P{)}_{,GG\prime}({g}^{G}P{)}_{,G}+V\text{\hspace{0.5em}},$ 


${g}^{GG\prime}={k}_{T}{\delta}^{jk}{\widehat{g}}_{j}^{G}{\widehat{g}}_{k}^{G\prime}\text{\hspace{0.5em}},$ 


${g}^{G}={f}^{G}+\frac{1}{2}{\delta}^{jk}{\widehat{g}}_{j}^{G\prime}{\widehat{g}}_{k,G\prime}^{G}\text{\hspace{0.5em}},$ 


$[...{]}_{,G}=\partial [...]/\partial {M}^{G}\text{\hspace{0.5em}}.$ 

(21)
This is properly referred to as a FokkerPlanck 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.,${\widehat{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\left[{M}_{t}{M}_{{t}_{0}}\right]dM\left(t\right)=\int \text{\hspace{0.17em}}...\text{\hspace{0.17em}}\int \text{\hspace{0.17em}}\underset{\_}{D}M\text{exp}(S)\delta \left[M\right({t}_{0})={M}_{0}]\delta \left[M\right(t)={M}_{t}]\text{\hspace{0.5em}},$ 


$S={k}_{T}^{1}\text{min}\text{\hspace{0.17em}}\underset{{t}_{0}}{\overset{t}{\int}}\text{\hspace{0.17em}}dt\prime L\text{\hspace{0.5em}},$ 


$\underset{\_}{D}M=\underset{u\to \infty}{\text{lim}}\text{\hspace{0.17em}}\underset{\rho =1}{\overset{u+1}{\Pi}}\text{\hspace{0.17em}}{g}^{1/2}\text{\hspace{0.17em}}\underset{G}{\Pi}\text{\hspace{0.17em}}(2\pi \theta {)}^{1/2}d{M}_{\rho}^{G}\text{\hspace{0.5em}},$ 


$L({\stackrel{\text{.}}{M}}^{G},{M}^{G},t)=\frac{1}{2}({\stackrel{\text{.}}{M}}^{G}{h}^{G}){g}_{GG\prime}({\stackrel{\text{.}}{M}}^{G\prime}{h}^{G\prime})+\frac{1}{2}{{h}^{G}}_{;G}+R/6V\text{\hspace{0.5em}},$ 


${h}^{G}={g}^{G}\frac{1}{2}{g}^{1/2}({g}^{1/2}{g}^{GG\prime}{)}_{,G\prime}\text{\hspace{0.5em}},$ 


${g}_{GG\prime}=({g}^{GG\prime}{)}^{1}\text{\hspace{0.5em}},$ 


$g=\text{det}\left({g}_{GG\prime}\right)\text{\hspace{0.5em}},$ 


${{h}^{G}}_{;G}={h}_{,G}^{G}+{\Gamma}_{GF}^{F}{h}^{G}={g}^{1/2}({g}^{1/2}{h}^{G}{)}_{,G}\text{\hspace{0.5em}},$ 


${\Gamma}_{JK}^{F}\equiv {g}^{LF}[JK,L]={g}^{LF}({g}_{JL,K}+{g}_{KL,J}{g}_{JK,L})\text{\hspace{0.5em}},$ 


$R={g}^{JL}{R}_{JL}={g}^{JL}{g}^{JK}{R}_{FJKL}\text{\hspace{0.5em}},$ 


${R}_{FJKL}=\frac{1}{2}({g}_{FK,JL}{g}_{JK,FL}{g}_{FL,JK}+{g}_{JL,FK})+{g}_{MN}({\Gamma}_{FK}^{M}{\Gamma}_{JL}^{N}{\Gamma}_{FL}^{M}{\Gamma}_{JK}^{N})\text{\hspace{0.5em}}.$ 

(22)
Mesoscopic variables have been defined
as${M}^{G}$
in the Langevin and FokkerPlanck representations, in terms
of their development from the microscopic system labeled
by$j$
The Riemannian curvature
term$R$
arises from
nonlinear${g}_{GG\prime}$
, which is a bona fide metric of this space [34]. Even
if a stationary solution,
i.e.,${\stackrel{\text{.}}{M}}^{G}=0$
, is ultimately sought, a necessarily prior stochastic
treatment
of${\stackrel{\text{.}}{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\prime}{\Sigma}\text{\hspace{0.17em}}{J}_{T\prime G}{M}^{G}$
, where the Lagrange
multipliers${J}_{T\prime G}$
are constraints
on${M}^{G}$
, which are advantageously modeled as extrinsic sources in
this representation; they too may be timedependent.
For our purposes, the above Feynman Lagrangian defines a
kernel of the shorttime 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$dt$
— the WKBJ approximation, but which does not include
the firstorder 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}_{TG}$
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$\overline{P}$
, when it exists, the information gain in
state$P$
is defined by 

${\rm Y}\left[P\right]=\text{\hspace{0.17em}}\int \text{\hspace{0.17em}}...\text{\hspace{0.17em}}\int \text{\hspace{0.17em}}\underset{\_}{D}M\prime \text{\hspace{0.17em}}P\text{\hspace{0.17em}}\text{ln}\text{\hspace{0.17em}}\left(P/\overline{P}\right)\text{\hspace{0.5em}},$ 


$\underset{\_}{D}M\prime =\underset{\_}{D}M/d{M}_{u+1}\text{\hspace{0.5em}}.$ 

(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 riskneural 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${\eta}^{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$\theta \stackrel{\text{.}}{M}\left(t\right)\to {M}_{\rho +1}{M}_{\rho}$
and$M\left(t\right)\to {M}_{\rho}$
The Stratonovich prescription corresponds to the midpoint
discretization
of$L$
,
wherein$\theta \stackrel{\text{.}}{M}\left(t\right)\to {M}_{\rho +1}{M}_{\rho}$
and$M\left(t\right)\to \frac{1}{2}({M}_{\rho +1}+{M}_{\rho})$
In terms of the functions appearing in the FokkerPlanck 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}({\stackrel{\text{.}}{M}}^{G},{M}^{G},t)=\frac{1}{2}({\stackrel{\text{.}}{M}}^{G}{g}^{G}){g}_{GG\prime}({\stackrel{\text{.}}{M}}^{G\prime}{g}^{G\prime})V\text{\hspace{0.5em}}.$ 

(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}$
,$\ll {M}^{G}\gg $
, 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 

${\delta}_{G}L={L}_{,G}{L}_{,\stackrel{\text{.}}{G}:t}=0\text{\hspace{0.5em}},$ 


${L}_{,\stackrel{\text{.}}{G}:t}={L}_{,\stackrel{\text{.}}{G}G\prime}{\stackrel{\text{.}}{M}}^{G\prime}+{L}_{,\stackrel{\text{.}}{G}\stackrel{\text{.}}{G}\prime}{\stackrel{\text{..}}{M}}^{G\prime}\text{\hspace{0.5em}}.$ 

(25)
For stationary
states,${\stackrel{\text{.}}{M}}^{G}=0$
,
and$\partial \overline{L}/\partial {\overline{M}}^{G}=0$
defines$\ll {\overline{M}}^{G}\gg $
, where the bars identify stationary variables; in this
case, the macroscopic variables are equal to their
mesoscopic counterparts. [Note
that$\overline{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],$\overline{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${\stackrel{\text{.}}{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$\rho $
(not to be confused with the Rho Greek calculated for
options). This section uses a simple twofactor model to
develop the correspondence between the
correlation$\rho $
in the Lagrangian and that among the commonly written Weiner
distributions$dz$
Consider coupled stochastic differential equations 

$dr={f}^{r}(r,l)dt+{\widehat{g}}^{r}(r,l){\sigma}_{1}d{z}_{1}\text{\hspace{0.5em}},$ 


$dl={f}^{l}(r,l)dt+{\widehat{g}}^{l}(r,l){\sigma}_{2}d{z}_{2}\text{\hspace{0.5em}},$ 


$<d{z}_{i}>=0\text{\hspace{0.5em}},\text{\hspace{0.5em}}i=\text{{}1,2\text{}}\text{\hspace{0.5em}},$ 


$<d{z}_{i}\left(t\right)d{z}_{j}(t\prime )>=dt\delta (tt\prime )\text{\hspace{0.5em}},\text{\hspace{0.5em}}i=j\text{\hspace{0.5em}},$ 


$<d{z}_{i}\left(t\right)d{z}_{j}(t\prime )>=\rho dt\delta (tt\prime )\text{\hspace{0.5em}},\text{\hspace{0.5em}}i\ne j\text{\hspace{0.5em}},$ 


$\delta (tt\prime )=\{\begin{array}{c}0\text{\hspace{0.17em}},,\\ 1\text{\hspace{0.17em}},\end{array}\text{\hspace{0.5em}}\text{\hspace{0.5em}}\begin{array}{c}t\ne t\prime \text{\hspace{0.5em}},\\ t=t\prime \text{\hspace{0.5em}},\end{array}(null)$ 

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

$dr/dt={f}^{r}+{\widehat{g}}^{r}{\sigma}_{1}({\gamma}^{+}{n}_{1}+sgn\rho \text{\hspace{0.17em}}{\gamma}^{}{n}_{2})\text{\hspace{0.5em}},$ 


$dl/dt={g}^{l}+{\widehat{g}}^{l}{\sigma}_{2}(sgn\rho \text{\hspace{0.17em}}{\gamma}^{}{n}_{1}+{\gamma}^{+}{n}_{2})\text{\hspace{0.5em}},$ 


${\gamma}^{\pm}=\frac{1}{\sqrt{2}}[1\pm (1{\rho}^{2}{)}^{1/2}{]}^{1/2}\text{\hspace{0.5em}},$ 


${n}_{i}=(dt{)}^{1/2}{p}_{i}\text{\hspace{0.5em}},$ 

(27)
where${p}_{1}$
and${p}_{2}$
are independent [0,1] Gaussian distributions.
The equivalent shorttime probability
distribution,$P$
, for the above set of equations is 

$P={g}^{1/2}(2\pi dt{)}^{1/2}\text{exp}(Ldt)\text{\hspace{0.5em}},$ 


$L=\frac{1}{2}{F}^{\u2020}\underset{\_}{g}F\text{\hspace{0.5em}},$ 


$F=\left(\begin{array}{c}\hfill dr/dt{f}^{r})\hfill \\ \hfill dl/dt{f}^{l})\hfill \end{array}\right)\text{\hspace{0.5em}},$ 


$g=\text{det}\left(\underset{\_}{g}\right)\text{\hspace{0.5em}},$ 


$k=1{\rho}^{2}\text{\hspace{0.5em}}.$ 

(28) $\underset{\_}{g}$
, the metric
in$\text{{}r,l\text{}}$
space, is the inverse of the covariance matrix, 

${\underset{\_}{g}}^{1}=\left(\begin{array}{cc}({\widehat{g}}^{r}{\sigma}_{1}{)}^{2}\hfill & \hfill \rho {\widehat{g}}^{r}{\widehat{g}}^{l}{\sigma}_{1}{\sigma}_{2}\\ \rho {\widehat{g}}^{r}{\widehat{g}}^{l}{\sigma}_{1}{\sigma}_{2}\hfill & \hfill ({\widehat{g}}^{l}{\sigma}_{2}{)}^{2}\end{array}\right)\text{\hspace{0.5em}}.$ 

(29)
The above also corrects previous papers which
inadvertently dropped
the$sgn$
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
shorttime 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 MonteCarlo
importancesampling 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 [5052]. 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${\alpha}_{k}^{i}$
in
dimension$i$
generated at
annealingtime$k$
with the range 

${\alpha}_{k}^{i}\in [{A}_{i},{B}_{i}]\text{\hspace{0.5em}},$ 

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

${\alpha}_{k+1}^{i}={\alpha}_{k}^{i}+{y}^{i}({B}_{i}{A}_{i})\text{\hspace{0.5em}},$ 


${y}^{i}\in [1,1]\text{\hspace{0.5em}}.$ 

(31)
The generating
function${g}_{T}\left(y\right)$
is defined, 

${g}_{T}\left(y\right)=\text{\hspace{0.17em}}\underset{i=1}{\overset{D}{\Pi}}\text{\hspace{0.17em}}\frac{1}{2({y}^{i}+{T}_{i})\text{\hspace{0.17em}}\text{ln}(1+1/{T}_{i})}\equiv \text{\hspace{0.17em}}\underset{i=1}{\overset{D}{\Pi}}\text{\hspace{0.17em}}{g}_{T}^{i}\left({y}^{i}\right)\text{\hspace{0.5em}},$ 

(32)
where the
subscript$i$
on${T}_{i}$
specifies the parameter index, and
the$k$
dependence
in${T}_{i}\left(k\right)$
for the annealing schedule has been dropped for brevity. Its
cumulative probability distribution is 

${G}_{T}\left(y\right)=\text{\hspace{0.17em}}\underset{1}{\overset{{y}^{1}}{\int}}...\underset{1}{\overset{{y}^{D}}{\int}}\text{\hspace{0.17em}}dy{\prime}^{1}...dy{\prime}^{D}\text{\hspace{0.17em}}{g}_{T}(y\prime )\equiv \text{\hspace{0.17em}}\underset{i=1}{\overset{D}{\Pi}}\text{\hspace{0.17em}}{G}_{T}^{i}\left({y}^{i}\right)\text{\hspace{0.5em}},$ 


${G}_{T}^{i}\left({y}^{i}\right)=\frac{1}{2}+\frac{\text{\hspace{0.17em}}sgn\text{\hspace{0.17em}}\left({y}^{i}\right)}{2}\frac{\text{ln}(1+{y}^{i}/{T}_{i})}{\text{ln}(1+1/{T}_{i})}\text{\hspace{0.5em}}.$ 

(33) ${y}^{i}$
is generated from
a${u}^{i}$
from the uniform distribution 

${u}^{i}\in U[0,1]\text{\hspace{0.5em}},$ 


${y}^{i}=\text{\hspace{0.17em}}sgn\text{\hspace{0.17em}}({u}^{i}\frac{1}{2}){T}_{i}\left[\right(1+1/{T}_{i}{)}^{2{u}^{i}1}1]\text{\hspace{0.5em}}.$ 

(34)
It is straightforward to calculate that for an annealing
schedule for 
${T}_{i}$

${T}_{i}\left(k\right)={T}_{0i}\text{exp}({c}_{i}{k}^{1/D})\text{\hspace{0.5em}},$ 

(35)
a global minima statistically can be obtained. I.e., 

$\underset{{k}_{0}}{\overset{\infty}{\Sigma}}\text{\hspace{0.17em}}{g}_{k}\text{\u2248}\text{\hspace{0.17em}}\underset{{k}_{0}}{\overset{\infty}{\Sigma}}\text{\hspace{0.17em}}\left[\text{\hspace{0.17em}}\underset{i=1}{\overset{D}{\Pi}}\text{\hspace{0.17em}}\frac{1}{2{y}^{i}{c}_{i}}\right]\frac{1}{k}=\infty \text{\hspace{0.5em}}.$ 

(36)
Control can be taken
over${c}_{i}$
, such that 

${T}_{fi}={T}_{0i}\text{exp}({m}_{i})\text{\hspace{0.5em}}when\text{\hspace{0.5em}}{k}_{f}=\text{exp}{n}_{i}\text{\hspace{0.5em}},$ 


${c}_{i}={m}_{i}\text{exp}({n}_{i}/D)\text{\hspace{0.5em}},$ 

(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
multidimensional search in the course of a complex
nonlinear physical problem, inevitably one must deal with
different changing sensitivities of
the${\alpha}^{i}$
in the search. At any given annealingtime, 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
annealingtime$k$
, essentially reannealing, every hundred or so
acceptanceevents (or at some userdefined 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 {\alpha}^{i}\text{\hspace{0.5em}}.$ 

(38)
In terms of the
largest${s}_{i}$
=${s}_{\text{max}}$
, a default rescaling is performed for
each${k}_{i}$
of each parameter dimension, whereby a new
index$k{\prime}_{i}$
is calculated from
each${k}_{i}$
, 

${k}_{i}\to k{\prime}_{i}\text{\hspace{0.5em}},$ 


$T{\prime}_{ik\prime}={T}_{ik}\left({s}_{\text{max}}/{s}_{i}\right)\text{\hspace{0.5em}},$ 


$k{\prime}_{i}=(\text{ln}\left({T}_{i0}/{T}_{ik\prime}\right)/{c}_{i}{)}^{D}\text{\hspace{0.5em}}.$ 

(39) ${T}_{i0}$
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}\left({k}_{i}\right)={T}_{0i}\text{exp}({c}_{i}{k}_{i}^{{Q}_{i}/D})\text{\hspace{0.5em}},$ 


${c}_{i}={m}_{i}\text{exp}({n}_{i}{Q}_{i}/D)\text{\hspace{0.5em}},$ 

(40)
in terms of the “quenching
factor”${Q}_{i}$
The sampling proof fails
if${Q}_{i}>1$
as 

$\underset{k}{\Sigma}\stackrel{D}{\Pi}1/{k}^{{Q}_{i}/D}=\underset{k}{\Sigma}1/{k}^{{Q}_{i}}<\infty \text{\hspace{0.5em}}.$ 

(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 twofactor model is more accurate because the
data indeed demonstrate stochastic volatility [4]. We
also note that the
twofactor$x$
fit by ASA across the last few years. This is not true of
the onefactor ASA
fitted$x$
the
BlackScholes$\sigma $
as a parameter, but rather calculate as historical
volatility during all runs. Some results of twofactor
studies and onefactor studies using a
BlackScholes$\sigma $
have been reported elsewhere [13].
Since$\sigma $
is not widely traded and arbitraged, to fit the twofactor
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$\sigma \prime $
, 

$\sigma \prime =StdDev(dS\text{\hspace{0.17em}}/\text{\hspace{0.17em}}F(S,{S}_{0},{S}_{\infty},x,y))\text{\hspace{0.5em}}.$ 

(42)
In the onefactor model, it does not make good numerical
sense to have two free parameters in one term,
i.e.,$\sigma $
and$x$
, as these cannot be fit very well within the variance the
data. Instead, one method is to take guidance from the
twofactor results, to set a scale for an
effective$\sigma $
, and then fit the
parameter$x$
Another method it apply the above StdDev as a proxy
for$\sigma $
Some motivation for this approach is given by considering
collapsing a twofactor stochastic volatility model in
onefactor model: The onefactor 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. PATHINTEGRAL (PATHINT) DEVELOPMENT
6.1. PATHINT Outline
The fits described above clearly
demonstrate the need to incorporate stochastic volatility in
option pricing models. If onefactor 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 pathintegral
code, PATHINT, described below, with great success. At this
time, the twofactor code takes too long to run for daily
use, but it proves to be a good weekly baseline for the
onefactor code.
The PATHINT algorithm develops the longtime probability
distribution from the Lagrangian fit by the first
optimization code. A robust and accurate histogrambased
(nonMonte Carlo) pathintegral algorithm to calculate the
longtime probability distribution has been developed to
handle nonlinear Lagrangians [1820,22,5355],
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$\Delta {M}_{i}$
For convenience, just consider a onedimensional system. The
above pathintegral representation can be rewritten, for
each of its intermediate integrals, as 

$P(M;t+\Delta t)=\text{\hspace{0.17em}}\int \text{\hspace{0.17em}}dM\prime \left[{g}_{s}^{1/2}\right(2\pi \Delta t{)}^{1/2}\text{exp}({L}_{s}\Delta t)\left]P\right(M\prime ;t)=\text{\hspace{0.17em}}\int \text{\hspace{0.17em}}dM\prime G(M,M\prime ;\Delta t\left)P\right(M\prime ;t)\text{\hspace{0.5em}},$ 


$P(M;t)=\underset{i=1}{\overset{N}{\Sigma}}\pi (M{M}_{i}){P}_{i}\left(t\right)\text{\hspace{0.5em}},$ 


$\pi (M{M}_{i})=\{\begin{array}{l}0\text{\hspace{0.17em}},\text{\hspace{0.5em}}({M}_{i}\frac{1}{2}\Delta {M}_{i1})\le M\le ({M}_{i}+\frac{1}{2}\Delta {M}_{i})\text{\hspace{0.5em}},\\ 1\text{\hspace{0.17em}},\text{\hspace{0.5em}}otherwise\text{\hspace{0.17em}},\end{array}\text{\hspace{0.5em}},(null)$ 

(43)
which yields 

${P}_{i}(t+\Delta t)={T}_{ij}(\Delta t){P}_{j}\left(t\right)\text{\hspace{0.5em}},$ 


${T}_{ij}(\Delta t)=\frac{2}{\Delta {M}_{i1}+\Delta {M}_{i}}\text{\hspace{0.17em}}\int \underset{{M}_{i}\Delta {M}_{i1}/2}{\overset{{M}_{i}+\Delta {M}_{i}/2}{\text{\hspace{0.17em}}}}dM\text{\hspace{0.17em}}\int \underset{{M}_{j}\Delta {M}_{j1}/2}{\overset{{M}_{j}+\Delta {M}_{j}/2}{\text{\hspace{0.17em}}}}dM\prime G(M,M\prime ;\Delta t)\text{\hspace{0.5em}}.$ 

(44) ${T}_{ij}$
is a banded matrix representing the Gaussian nature of
the shorttime probability centered about the (varying)
drift.
Care must be used in developing the mesh
in$\Delta {M}^{G}$
, which is strongly dependent on the diagonal elements of
the diffusion matrix, e.g., 

$\Delta {M}^{G}\text{\u2248}(\Delta t{g}^{GG}{)}^{1/2}\text{\hspace{0.5em}}.$ 

(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
twofactor stochasticvolatility model.
Fitting data with the shorttime 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 timedependent pathintegral: By considering the
contributions to the first and second moments
of$\Delta {M}^{G}$
for small time
slices$\theta $
, conditions on the time and variable meshes can be
derived [53]. The time slice essentially is determined
by$\theta \text{\hspace{0.5em}}\le \text{\hspace{0.5em}}{\overline{L}}^{1}$
,
where$\overline{L}$
is the “static” Lagrangian
with$d{M}^{G}/dt\text{\hspace{0.5em}}=\text{\hspace{0.5em}}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$\Delta {M}^{G}$
is measured by the
covariance${g}^{GG\prime}$
,
or$\Delta {M}^{G}\text{\u2248}({g}^{GG}\theta {)}^{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
initialcondition discretized delta
function$MSUPG$
at the prepoint of an interval, the mesh going forward
in${M}^{G}$
is simply calculated stepwise using 

$\Delta {M}^{G}=({g}^{GG}\theta {)}^{1/2}\text{\hspace{0.5em}}.$ 

(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 WehnerWolfer
paper [53].
It is important to stress that very good numerical accuracy
is required to get very good Greeks required for realworld
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 nfactor model
subject to its diffusionmesh constraints. In practice, the
calculation of 3factor and 4factor models likely will wait
until gigahertz speeds and gigabyte RAM are
commonplace.
6.2. Development of LongTime 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
“basispoint 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 atthemoney (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 shorttime distribution evolved out
to$T=0.5$
year
for$x$
in {1, 0, 1, 2}, with 500 intermediate epochs/foldings, and
BS$\sigma =0.0075$
Each calculation
scales$\sigma $
by multiplying
by$S/F(S,{S}_{0},{S}_{\infty},x,y)$
Fig. 3. gives an example of a
twofactor distribution evolved out
to$T=0.5$
year
for$x=0.7$
6.4. TwoFactor Volatility and PATHINT Modifications
In our twofactor model, the
mesh
of$S$
would depend
on$\sigma $
and cause some problems in any PATHINT grid to be developed
in$S$
$\sigma $
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 offdiagonal regions in
an$n$
factor system, e.g., 

$\Delta {M}_{k}^{i}={2}^{{m}_{k}^{i}}\Delta {M}_{0}^{i}\text{\hspace{0.5em}},$ 


$\Delta {M}_{0}^{i}\text{\u2248}\sqrt{{g}_{{k}_{0}}^{ii}\Delta t}\text{\hspace{0.5em}},$ 

(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 onefactor PATHINT mesh
conditions, in terms of some minimal
mesh$\Delta {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
atthemoney BPV are scaled to be equivalent, there is not a
very drastic change in the onefactor ATM Greeks developed
below. Therefore, while we have not at all changed the
functional dependence of the Lagrangian
on$S$
and$\sigma $
, we have determined our meshes using a diffusion for
the$S$
equation
as${\sigma}_{0}\text{\hspace{0.17em}}F(S,{S}_{0},{S}_{\infty},x,y)$
,
where${\sigma}_{0}$
is determined by the same BPVequivalent condition as
imposed on the onefactor models. This seems to work very
well, especially since we have taken
our$\sigma $
equation to be normal with a limited range of influence in
the calculations. Future work yet has to establish a more
definitive distribution
for$\sigma $
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 CrankNicolson 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=\text{\hspace{0.17em}}<\text{max}\left[z\right(SX),0]>\text{\hspace{0.5em}},\text{\hspace{0.5em}}\{\begin{array}{cc}z=1\text{\hspace{0.5em}},\hfill & call\hfill \\ z=1\text{\hspace{0.5em}},\hfill & put\hfill \end{array}\text{\hspace{0.5em}},(null)$ 

(48)
where$X$
is the strike price, and the expected
value$<.>$
is taken with respect to the riskneutral
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 riskneutral distribution
is effectively calculated taking the drift as the
costofcarry$b$
times$S$
, using the above arguments leading to the BS formula. We
have designed our codes to use parameters
riskfreerate$r$
and
costofcarry$b$
such that 

$b=r,\text{\hspace{0.5em}}\text{\hspace{0.5em}}cost\text{\hspace{0.5em}}of\text{\hspace{0.5em}}carry\text{\hspace{0.5em}}on\text{\hspace{0.5em}}nondividend\text{\hspace{0.5em}}stock\text{\hspace{0.5em}},$ 


$b=rq,\text{\hspace{0.5em}}cost\text{\hspace{0.5em}}of\text{\hspace{0.5em}}carry\text{\hspace{0.5em}}on\text{\hspace{0.5em}}stock\text{\hspace{0.5em}}paying\text{\hspace{0.5em}}dividend\text{\hspace{0.5em}}yield\text{\hspace{0.5em}}q\text{\hspace{0.5em}},$ 


$b=0,\text{\hspace{0.5em}}cost\text{\hspace{0.5em}}of\text{\hspace{0.5em}}carry\text{\hspace{0.5em}}on\text{\hspace{0.5em}}future\text{\hspace{0.5em}}contract\text{\hspace{0.5em}},$ 


$b=r{r}_{f},\text{\hspace{0.5em}}cost\text{\hspace{0.5em}}of\text{\hspace{0.5em}}carry\text{\hspace{0.5em}}on\text{\hspace{0.5em}}currency\text{\hspace{0.5em}}with\text{\hspace{0.5em}}rate\text{\hspace{0.5em}}{r}_{f}\text{\hspace{0.5em}},$ 

(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}=\text{\hspace{0.17em}}<\text{max}[z({e}^{(br)T}{S}_{T}{e}^{rT}X),0]>\text{\hspace{0.5em}}.$ 

(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\left[M\right(t+dt),t+dtM(t),t]$
is 

$\partial P/\partial t=\frac{1}{2}({g}^{GG\prime}P{)}_{,GG\prime}({g}^{G}P{)}_{,G}+V\text{\hspace{0.5em}},$ 

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

$\partial P/\partial t=\frac{1}{2}{g}^{GG\prime}{P}_{,GG\prime}{g}^{G}{P}_{,G}+V\text{\hspace{0.5em}},$ 

(52)
where the partial derivatives with respect
to${M}^{G}$
act on the
prepoint${M}^{G}\left(t\right)$
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\left(T\right)$
, e.g., using a pathintegral
kernel$K$
folded
over$n$
epochs, where it is folded with the
function$O$
, 

$V=\text{\hspace{0.17em}}<O\left(T\right)\text{\hspace{0.17em}}P\left(T\right)>\text{\hspace{0.17em}}=\text{\hspace{0.17em}}<O({K}^{n}P\left({t}_{0}\right))>\text{\hspace{0.5em}},$ 


$O\left(T\right)=\text{max}[z({e}^{(br)T}{S}_{T}{e}^{rT}X),0]\text{\hspace{0.5em}},$ 

(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}^{\u2020}$
acting
on$O$
, 

$V=\text{\hspace{0.17em}}<O\left(T\right)\text{\hspace{0.17em}}P\left(T\right)>\text{\hspace{0.17em}}=\text{\hspace{0.17em}}<({K}^{n}O\left(T\right))P\left({t}_{0}\right)>\text{\hspace{0.5em}}.$ 

(54)
It is convenient to use the simple prepoint
representation for the Lagrangian, so the backward kernel is
first reexpressed as a forward kernel by bringing the
diffusions and drifts inside the partial derivatives, giving
a transformed adjoint
kernel${K}^{\u2020}$
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}_{new}=\text{max}\left[\right(SX),({e}^{rdt}{O}_{old}\left)\right]\text{\hspace{0.5em}}.$ 

(55)
The
Greeks$\text{{}\Delta ,\Gamma ,\Theta \text{}}$
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$\Delta $
requires moving back one epoch
and$\Gamma $
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 wellknown
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$\text{{}\Delta ,\Gamma ,\Theta \text{}}$
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], 

$CR{R}_{variant}=CR{R}_{American}CR{R}_{European}+BS\text{\hspace{0.5em}}.$ 

(56)
The second problem can be alleviated somewhat by
averaging runs
of$n$
iterations with runs
of$n+1$
iterations [60]. This fourfold 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$\Delta $
and$\Gamma $
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{df}{dx}=\frac{f(x+dx)f(xdx)}{2dx}$ 

(57)
is typically good
to$o((dx{)}^{3})$
, 

$\frac{df}{dx}=\frac{f(x+2dx)+8f(x+dx)8f(xdx)+f(x2dx)}{12dx}$ 

(58)
is typically good
to$o((dx{)}^{5})$
Similarly, while 

$\frac{{d}^{2}f}{d{x}^{2}}=\frac{f(x+dx)2f\left(x\right)+f(xdx)}{dx\text{\hspace{0.17em}}dx}$ 

(59)
is typically good
to$((dx{)}^{4})$
, 

$\frac{{d}^{2}f}{d{x}^{2}}=\frac{d(x+2dx)+16f(x+dx)30f\left(x\right)+16f(xdx)f(x2dx)}{dx\text{\hspace{0.17em}}dx}$ 

(60)
is typically good
to$((dx{)}^{6})$
Table 1 gives an example of baselining our onefactor
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$
,$\sigma =0.10$
Table 1.
Tests with American CRR and American PATHINT lead to
results with the same degrees of accuracy.
7.3. TwoFactor PATHINT Baselined to OneFactor PATHINT
Previous papers and tests have
demonstrated that the twofactor 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 onefactor problems, and there is no loss of
accuracy. However, here we are making some additional mesh
approximations as discussed above to
accommodate$\sigma $
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 onefactor
PATHINT results for
various$x$
scaling is performed to have all models used the same BPV
(using
the${\sigma}_{0}$
procedure for the mesh as described above for the twofactor
model).
The logical extension of Greeks for the twofactor model is
to develop derivatives of price with respect
to$\rho $
and$\epsilon $
in$\sigma $
volatility equation. However, we did not find a bona
fide twofactor proxy for the
onefactor${\rm Y}$
, the derivative of price with respect to the
onefactor$\sigma $
constant. We get very good
ATM${\rm Y}$
comparisons between BS and our onefactor models with
various$x$
We tried simply multiplying the noise in the twofactor
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 BPVscaled
BS${\rm Y}$
within a couple of significant figures. Perhaps this is not
too surprising, especially given the correlation
substantial$\rho $
between the price and volatility equations which we do not
neglect.
8. CONCLUSIONS
The results of our study are
that, after the atthemoney basispoint volatilities are
scaled to be equivalent, there is only a very small change
in option prices for different
exponents$x$
, both for the onefactor and twofactor models. There still
are significant differences in Greeks using more
sophisticated models, especially for outofthemoney
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
onefactor 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 twofactor
exponents$x$
are reasonably faithful indicators defining these different
contexts.
The$x$
exponents in the twofactor fits are quite stable. As such,
especially the
twofactor$x$
can be considered as a “context indicator” over
a longer time scale than other indicators typically used by
traders. The twofactor fits also exhibit differences due to
the$\sigma $
parameters, including
the$\rho $
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 GaussianMarkovian 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 quasilinear
approximations. Three quite different systems have benefited
from this approach:
The largescale modeling of neocortical interactions has
benefited from the use of intuitive constructs that yet are
faithful to the complex algebra describing this
multiplescaled complex system. For example,
canonicalmomenta indicators have been successfully applied
to multivariate financial markets.
It is clear that ASA optimization and PATHINT pathintegral
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 nontypical 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
longterm options, Pennsylvania State U, University Park,
PA, (1998). 
3. 

L. Ingber, Some Applications of Statistical Mechanics of
Financial Markets, LIR981SASMFM, Lester Ingber Research,
Chicago, IL, (1998). 
4. 

L. Ingber and J.K. Wilson, Volatility of volatility of
financial markets, Mathl. Computer Modelling
29 (5), 3957 (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 leastsquares approach,
Capital Management Sciences, Los Angeles, CA, (1998). 
8. 

L. Ingber, Statistical mechanics of nonlinear
nonequilibrium financial markets, Math. Modelling
5 (6), 343361 (1984). 
9. 

L. Ingber, Statistical mechanical aids to calculating
term structure models, Phys. Rev. A
42 (12), 70577064 (1990). 
10. 

L. Ingber, M.F. Wehner, G.M. Jabbour, and T.M. Barnhill,
Application of statistical mechanics methodology to
termstructure bondpricing models, Mathl. Comput.
Modelling
15 (11), 7798 (1991). 
11. 

L. Ingber, Statistical mechanics of nonlinear
nonequilibrium financial markets: Applications to optimized
trading, Mathl. Computer Modelling
23 (7), 101121 (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. 777784, Springer, New York,
(1996). 
13. 

L. Ingber and J.K. Wilson, Statistical mechanics of
financial markets: Exponential modifications to
BlackScholes, Mathl. Computer Modelling
31 (8/9), 167192 (2000). 
14. 

C. Azariadis, Selffulfilling prophecies, J. Econ.
Theory 25, 380396 (1981). 
15. 

L. Ingber, Adaptive Simulated Annealing (ASA), Global
optimization Ccode, 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, 229263 (1979). 
17. 

L. Ingber, Data mining and knowledge discovery via
statistical mechanics in nonlinear stochastic systems,
Mathl. Computer Modelling
27 (3), 931 (1998). 
18. 

L. Ingber, H. Fujio, and M.F. Wehner, Mathematical
comparison of combat computer models to exercise data,
Mathl. Comput. Modelling
15 (1), 6590 (1991). 
19. 

L. Ingber, Statistical mechanics of neocortical
interactions: Pathintegral evolution of shortterm memory,
Phys. Rev. E
49 (5B), 46524664 (1994). 
20. 

L. Ingber and P.L. Nunez, Statistical mechanics of
neocortical interactions: High resolution pathintegral
calculation of shortterm memory, Phys. Rev. E
51 (5), 50745083 (1995). 
21. 

L. Ingber, R. Srinivasan, and P.L. Nunez, Pathintegral
evolution of chaos embedded in noise: Duffing neocortical
analog, Mathl. Computer Modelling
23 (3), 4353 (1996). 
22. 

L. Ingber, Pathintegral evolution of multivariate
systems with moderate noise, Phys. Rev. E
51 (2), 16161619 (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, 2186 (1900). 
26. 

M. C. Jensen, Some anomalous evidence regarding market
efficiency, an editorial introduction, J. Finan.
Econ. 6, 95101 (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, 225236 (1971). 
28. 

S. J. Taylor, Tests of the random walk hypothesis
against a pricetrend hypothesis, J. Finan. Quant.
Anal. 17, 3761 (1982). 
29. 

P. Brown, A. W. Kleidon, and T. A. Marsh, New evidence
on the nature of sizerelated anomalies in stock prices,
J. Fin. Econ. 12, 3356 (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), 40174060 (1991). 
32. 

K.S. Cheng, Quantization of a general dynamical system
by Feynman’s path integration formulation, J. Math.
Phys. 13, 17231726 (1972). 
33. 

H. Dekker, Functional integration and the
OnsagerMachlup Lagrangian for continuous Markov processes
in Riemannian geometries, Phys. Rev. A
19, 21022111 (1979). 
34. 

R. Graham, Pathintegral methods in nonequilibrium
thermodynamics and statistics, in Stochastic Processes in
Nonequilibrium Systems, (Edited by L. Garrido, P.
Seglar, and P.J. Shepherd), pp. 82138, 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, 419432 (1979). 
36. 

F. Langouche, D. Roekaerts, and E. Tirapegui, Short
derivation of Feynman Lagrangian for general diffusion
process, J. Phys. A
113, 449452 (1980). 
37. 

F. Langouche, D. Roekaerts, and E. Tirapegui,
Functional Integration and Semiclassical Expansions,
Reidel, Dordrecht, The Netherlands, (1982). 
38. 

M. RosaClot and S. Taddei, A path integral approach to
derivative security pricing: I. Formalism and analytic
results, INFN, Firenze, Italy, (1999). 
39. 

P. Hagedorn, NonLinear 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, 381408 (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,
SpringerVerlag, Berlin, Germany, (1983). 
44. 

R. C. Merton, An intertemporal capital asset pricing
model, Econometrica
41, 867887 (1973). 
45. 

L. Ingber, Statistical mechanics of neocortical
interactions: Stability and duration of the 7+−2 rule
of shorttermmemory capacity, Phys. Rev. A
31, 11831186 (1985). 
46. 

L. Ingber, Very fast simulated reannealing, Mathl.
Comput. Modelling
12 (8), 967973 (1989). 
47. 

S. Kirkpatrick, C.D. Gelatt, Jr., and M.P. Vecchi,
Optimization by simulated annealing, Science
220 (4598), 671680 (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), 721741 (1984). 
49. 

H. Szu and R. Hartley, Fast simulated annealing,
Phys. Lett. A
122 (34), 157162 (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), 2957 (1993). 
52. 

L. Ingber, Adaptive simulated annealing (ASA): Lessons
learned, Control and Cybernetics
25 (1), 3354 (1996). 
53. 

M.F. Wehner and W.G. Wolfer, Numerical evaluation of
pathintegral solutions to FokkerPlanck equations. I.,
Phys. Rev. A
27, 26632670 (1983). 
54. 

M.F. Wehner and W.G. Wolfer, Numerical evaluation of
pathintegral solutions to FokkerPlanck equations. II.
Restricted stochastic processes, Phys. Rev. A
28, 30033011 (1983). 
55. 

M.F. Wehner and W.G. Wolfer, Numerical evaluation of
path integral solutions to FokkerPlanck equations. III.
Time and functionally dependent coefficients, Phys. Rev.
A 35, 17951801 (1987). 
56. 

M. RosaClot 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 (56), 167182 (1999). 
59. 

E.G. Haug, The Complete Guide to Option Pricing
Formulas, McGrawHill, 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. 4366, Cambridge University Press, Cambridge,
UK, (1997). 
FIGURE CAPTIONS
FIG. 1.
(a)$F(S,{S}_{0},{S}_{\infty},x,y)$
for$x=1$
, the BlackScholes 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 shorttime probability distribution at
time$T=0.5$
years
for$x=1$
, the (truncated) BlackScholes distribution. The shorttime
probability distribution at
time$T=0.5$
years
for$x=0$
, the normal distribution. The shorttime probability
distribution at
time$T=0.5$
years
for$x=1$
The shorttime probability distribution at
time$T=0.5$
years
for$x=2$
FIG. 3. A twofactor 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$CR{R}_{European}$
,$CR{R}_{American}$
,$CR{R}_{variant}$
, 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$CR{R}_{American}$
differs somewhat from
the$CR{R}_{European}$
due to the discrete nature of the calculation. All CRR
calculations include averaging over 300 and 301 iterations
to minimize oscillatory errors.
Figure 1
Figure 2
Figure 3
Table 1
                                                             