Apparatus and method of pricing financial derivatives

[0002] The pricing of an option of an asset is a fundamental problem of significant practical importance in today's financial markets. In 1973, a mathematician, Fischer Black, and an economist, Myron Scholes, devised one of the first mathematically accepted approaches for pricing what is known as a “European” option, which are options that can.only be exercised at its expiration date. What has become known as the Black-Scholes option formula was described first in “The pricing of options and corporate liabilities,” Journal of Political Economy81 (1973). The Black-Scholes option formula is presently of widespread use in financial markets all over the world. The price of such an option can be found by solving the Black-Scholes equation with the initial condition at expiration (i.e., the payoff of the option). The Black-Scholes equation is a reverse diffusion equation with parameters determined by the statistical characteristics of involved stocks and currencies such as risk free interest rate, holding cost or expected dividends, and volatility.

[0003] As an example, the Black-Scholes formula for the theoretical price of a vanilla European call option is:

[0004] C(S, t)=SN(d₁)−Ee^−r(T−t)N(d₂),  (1)

[0005] where the notation is fairly standard, as described by P. Wilmott, J. N. Dewynne and S. Howison, Option Pricing: Mathematical Models and Computation, Oxford Financial Press, Oxford (1993).

[0006] Unlike a European option, an “American” option gives the owner of the option the right of exercising the option before its expiration date. In recent years, the American option has become more prevalent than the European option. Due to the additional feature of early exercise, the pricing of an American option is generally considered to be more difficult than the pricing of a European option, especially when one considers “exotic” options which are variations and refinements of a basic American option. From a mathematical point of view, part of this difficulty is due to the price of American options obeying the Black-Scholes equation only in the region where it is statistically appropriate to hold the option rather than to exercise it immediately.

[0007] Therefore, in the case of American options, the above formula (1) and its analogs are no longer valid. In fact, as shown in P. Jaillet, D. Lamberton, and B. Lapeyre, “Variational inequalities and the pricing of American options,” Acta Applicandae Mathematicae21 (1990) 263-289, a rigorous mathematical model for pricing an American option is an infinite-dimensional free boundary problem. As such, there is in general no explicit formula or finite procedure for computing the exact price of an American option. As a result, various mathematical models have been devised in an attempt to approximate the price of such options.

[0008] However, the option prices computed from a mathematical model are of a theoretical nature. In computing these prices, various inputs are fed into the model and an algorithm produces an answer. In practice, the computed prices may not be consistent with the observed market prices, e.g., the prices on the trading floor. Ideally, these two sets of prices should coincide. However, such a result is difficult, if not impossible, using known models. Two principal reasons for this are: (i) the assumptions that lead to the construction of the mathematical model may not be realistic; and (ii) the inputs to the model are not correct.

[0009] Knowing the correct inputs is crucial to the success of any pricing model. Ideally the two sets of option prices, computed and observed, are within an acceptable range of one another. In general, the problem of computing a proper set of inputs to a forward pricing model so that the computed outputs obey certain prescribed criteria is called an “inverse” pricing problem. One such input which is crucial to the forward pricing model of American options is the implied volatility of an option. The volatility of an asset is an important input to an option pricing model and it is also an input that is most difficult to estimate. Black and Scholes assumed this parameter to be a constant when deriving their famous formula for the theoretical price of a vanilla European call option. In the above equation (1), for example, the two constants d₁and d₂contain the volatility parameter σ. As it is now well known, this parameter is, in general, not a constant; indeed it is a highly complicated function of several deterministic and random factors.

[0010] Previous approaches for dealing with this difficult problem of unknown volatility are numerous and include: (i) statistical estimation methods based on historical data (such as J. Hull, Options, Futures, and Other Derivative Securities, Second Edition, Prentice Hall, New Jersey (1989), Section 10.4 and R. Gibson, Option Valuation: Analyzing and Pricing Standardized Option Contracts, McGraw-Hill, New York (1991), Section 1; (ii) mathematical models of stochastic volatilities such as those in J. Hull and A. White, “The pricing of options on assets with stochastic volatilities,” The Journal of Finance42 (1987) 281-300; H. Johnson and D. Shanno, “Option pricing when the variance is changing,” Journal of Financial and Quantitative Analysis22 (1987) 143-151; and (iii) implied volatilities based on observed option prices (suggested originally by H. A. Latant and R. J. Rendleman, “Standard deviations of stock price ratios implied in option prices,” The Journal of Finance31 (1976) 369-381 and empirically tested by S. Beckers, “Standard deviation implied in option prices as predictors of future stock price volatility” Journal of Banking and Finance5 (1981) 363-381).

[0011] The problem of computing implied volatilities of (European or American) options is an instance of an inverse problem that is the counterpart of the forward problem of pricing these options. Specifically, in the forward option pricing problem, a constant volatility parameter, along with other constants, such as the interest rate and asset dividend, is taken as an input to a mathematical model that produces the theoretical or modeled option prices. By using an incorrect volatility parameter in the forward pricing model, the computed option price is bound to deviate, often substantially, from the option price observed on the trading floor.

[0012] An apparatus for and method of determining the price of financial derivatives such as options. One preferred embodiment of the invention employs a discretized partial differential linear complementarity problem (PDLCP) based system to determine the forward pricing of financial instruments such as vanilla American options. In this embodiment, an optimization problem in the form of a mathematical program with equilibrium constraints (MPEC) is implemented to derive implied volatilities of the assets underlying the subject derivatives. The implied volatilities thus derived are used as inputs in the PDLCP-based system to accurately determine the forward pricing of the subject derivatives.

[0013] Many advantages, features, and applications of the invention will be apparent from the following detailed description of the invention which is provided in connection with the accompanying drawings in which:

[0014] FIG. 1 is a block diagram of a financial engine in accordance with a preferred embodiment of the invention;

[0015] FIG. 2 illustrates a flow chart which describes the operational flow of a preferred embodiment of the invention;

[0016] FIG. 3 illustrates a flow chart which describes a preferred embodiment of the invention that employs an implicit programming algorithm (IMPA);

[0017] FIG. 4 illustrates a flow chart which describes a preferred embodiment of the invention that employs a penalty interior point algorithm (PIPA);

[0018] FIG. 5 is a graph which depicts the volatility surface produced by the exemplary penalty interior point algorithm employed in accordance with a preferred embodiment of the invention;

[0019] FIG. 6 is a graph which depicts the volatility surface produced by the exemplary implicit programming algorithm employed in accordance with a preferred embodiment of the invention;

[0020] FIG. 7 is a graph which depicts the exemplary option values calculated using the volatility surface shown in FIGS. 5 and 6;

[0021] FIG. 8 is another graph which depicts the exemplary option values calculated using the volatility surface shown in FIGS. 5 and 6;

[0022] FIG. 9 is another graph which depicts the volatility surface produced by the exemplary, penalty interior point algorithm employed in accordance with a preferred embodiment of the invention;

[0023] FIG. 10 is another graph which depicts the volatility surface produced by the exemplary implicit programming algorithm employed in accordance with a preferred embodiment of the invention;

[0024] FIG. 11 is a graph which depicts the exemplary option values calculated using the volatitity surface shown in FIGS. 9 and 10;

[0025] FIG. 12 is another graph which depicts the exemplary option values calculated using the volatility surface shown in FIGS. 9 and 10;

[0026] FIG. 13 is another graph which depicts the volatility surface produced by the exemplary use of the penalty interior point method in accordance with a preferred embodiment of the invention;

[0027] FIG. 14 is another graph which depicts the volatility surface produced by the use of the implicit programming method in accordance with a preferred embodiment of the invention;

[0028] FIG. 15 is a graph which depicts the option values calculated using the exemplary volatility surface shown in FIGS. 13 and 14; and

[0029] FIG. 16 is another graph which depicts the option values calculated using the exemplary volatility surface shown in FIGS. 13 and 14.

[0030] Preferred embodiments and applications of the invention will now be described with reference to FIGS. 1-16. Other embodiments may be realized and structural or logical changes may be made to the disclosed embodiments without departing from the spirit or scope of the invention. Although the invention is particularly described as applied to the analysis, synthesis, and pricing of financial derivatives, particularly calculation and application to American options, it should be readily apparent that the invention may be applied to any other financial or economic instruments having the same or similar problems.

[0031] In a preferred embodiment, the invention is implemented in a financial instrument engine 100, as shown in FIG. 1, used to analyze financial and economic signals provided by the financial markets and, based upon user input data and commands, issue position statements or reports for use by the user (or other components or systems), as well as issuing control signals for use by automated systems to effect positions (e.g., increase, decrease, change, etc.) held by the user in the financial markets. The financial engine 100 may be a stand-alone computer hardware system, incorporated in (or distributed among) one or more locally or remotely located computer systems.

[0032] In a preferred embodiment, financial engine 100 is composed of a plurality of modules: data source 10, implied volatility calculator 12, financial instrument controller 14, forward pricing unit 16, and user interface device 18. In the illustrated embodiment, the modules are connected by a single transmission bus 13. (It should be understood that the illustration of bus 13 is merely representative of the various connectivity technologies available to those of ordinary skill in the art including single/multiple, wired, wireless, and other transmission mediums.) For illustrative purposes only, the modules depicted within financial engine 100 have been selected to illustrate a preferred embodiment of the invention that provides forward pricing of “vanilla” options (i.e., options limited to either buying or selling an underlying asset), particularly vanilla American options. The modules may be changed as needed for the specific financial instrument to be used (e.g., “exotic” options such as barrier options, digital options, compound options, Asian options, etc.).

[0033] In a preferred embodiment, data source 10 is used to provide external financial and economic data, signals, or other information to financial engine 100. Data source 100 may include one or more transmission links or connections (wired, wireless, etc.) to a variety of automatic quotation systems/services that provide current market data, and may also include its own information retrieval system(s). The information received by data source 100 is forwarded to one or more of the other modules depending on the particular information provided. Volatility calculator 12, which may be implemented using a computer program or algorithm (described below) executed by central processing unit (CPU) 120, and stored in memory 122, is used to calculate the implied volatility in price of the underlying asset of the American option.

[0034] Forward pricing unit 16, which may also be implemented using a computer program or algorithm (described below) stored in memory 162 and executed by CPU 160, receives information from data source 10, volatility calculator 12, and user interface device 18 to determine the forward price or value V_BS(S, t) (as a function of the price (S) of the underlying asset over time (t)). User interface device 18 is used to exchange information between the user and financial engine 100. Depending on the commands to be implemented (e.g., forward price calculation, implied volatility calculation, option exercise, etc.), the information may include the observed market price of the underlying asset, volatility (local, implied, etc.), interest rate (e.g., risk-free), maturity (e.g., expiration date), strike price, dividend yield, and other pertinent information. (This information may also be provided by (or output to) data source 10.)

[0035] Financial instrument controller 14 may exchange information from all of the modules in financial engine 100. Based on the information provided, financial instrument controller 14 can provide statements and reports to the user (e.g., containing pricing, volatility, position information, etc.) or to others through network 110, as well as provide a host of services based on the information such as controlling trading or allocation of funds in options (or other financial instruments), simulating market reactions based on input conditions, through signals output to network 110. (The depiction of network 110 is made to represent a variety of known networks and connected systems such as local or wide area networks, e.g., as a company intranet, the Internet, electronic communications network (ECNs), small order exchange systems (SOES), on-line brokers or other trading networks, etc.)

[0036] In order to accurately project or determine the forward price or value of financial or economic instruments such as American options, a preferred embodiment of the financial engine 100 utilizes an inverse pricing model to determine the inputs to be used by forward pricing unit 16 such that a more accurate forward theoretical price can result by solving the equation:

[0037] V_BS(S, t)=V_obs  (2)

[0038] where V_obsrepresents the value of the option observed in the market, for example, as reported by data source 10 or input to user interface device 18. In determining the forward option price, a preferred embodiment of the invention utilizes a calculation of the implied volatility of the underlying asset as determined by volatility calculator 12.

[0039] In a preferred embodiment, volatility calculator 12 includes CPU 120 that is used to execute a computer program or algorithm stored in memory 122. As will be discussed in detail below, in a preferred embodiment of the invention, the computer program calculates the implied volatilities of American options using a mathematical programming equilibrium constraints (MPEC) approach. A detailed explanation of the MPEC approach is given by Z. Q. Luo, J. S. Pang, and D. Ralph, in Mathematical Program with Equilibrium Constraints, Cambridge University Press, Cambridge (1996), which is incorporated herein by reference in its entirety.

[0040] Using the calculations made by volatility calculator 12, forward pricing unit 16 can project or determine the forward price or value of the option as a function of the price of the underlying asset and time. In a preferred embodiment, forward pricing unit 16 employs CPU 160 to execute a computer program or algorithm stored in memory 162. The computer program formulates the forward pricing task as a discretized linear complementarity problem (LCP), as described in “Option Pricing and Linear Complementarity,” Journal of Computational Finance2 (Fall 1998), by Jacqueline Huang and Jong-Shi Pang, which is incorporated herein by reference in its entirety. The basic framework of the forward pricing program is that of Black and Scholes, wherein the price of a risky asset is assumed to satisfy the following stochastic differential equation:

[0041] dS=μSdt+σ(S, t)SdW,

[0042] where S denotes the asset price that is a function of the time tε[0, T], with T>0 being the duration of a vanilla European option written on the asset, μ is the drift of the stochastic price process of the asset, dW is a standard Wiener process, and σ(S, t) is the volatility that is taken to be an unknown function of the pair (S, t), which defines the volatility surface.

[0043] In a preferred embodiment, the value V(S, t) of an American option should satisfy the following partial differential linear complementarity problem (PDLCP) for t in [0, T) and S in [0, ∞):

[0044] 0≦V(S, t)−Λ(S, t),

[0045] 0≧L_BS(V),

[0046] 0=[V(S, t)−Λ(S, t)]L_BS(V),  (3)

[0047] where ${ℒ}_{\mathrm{BS}}\equiv \frac{\partial }{\partial t}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^2}{\partial S^2}+\left(r-D_0\right)S\frac{\partial }{\partial S}-r,$

[0048] is the Black-Scholes partial differential operator, with r being the constant interest rate of a risk-free asset (e.g., U.S. government security), D₀being the constant dividend rate of the asset, and Λ(S, t) being the payoff function when early exercise occurs. Preferably, boundary values of V(S, t) at S=0 and S=∝, and terminal values of V(S, t) at t=T, the expiry time, are assumed.

[0049] In a preferred embodiment, the financial instrument engine 100 employs a finite difference scheme for approximating the partial differential equation (3). Specifically, the state variable S is truncated to a finite range [0, NδS], where N is a positive integer and δS>0 is the step size that will be used for discretizing the partial differentiation with respect to S. Preferably, a time step δt>0 is used for which M≡T/δt is an integer for, discretizing the partial differentiation with respect to t. In the program, the following approximations may be employed: $\begin{array}{c}\frac{\partial V}{\partial t}\left(S,t\right)=\frac{V\left(S,t+\delta t\right)-V\left(S,t\right)}{\delta t}+O\left(\delta t\right),\\ \frac{{\partial }^{2}V}{\partial S^2}\left(S,t\right)=\theta \frac{V\left(S+\delta S,t\right)-2V\left(S,t\right)+V\left(S-\delta S,t\right)}{{\left(\delta S\right)}^2}+\\ \left(1-\theta \right)\frac{\begin{array}{c}V\left(S+\delta S,t+\delta t\right)-2V\left(S,t+\delta t\right)+\\ V\left(S-\delta S,t+\delta t\right)\end{array}}{{\left(\delta S\right)}^2}+O\left({\left(\delta S\right)}^2\right),\\ \frac{\partial V}{\partial S}\left(S,t\right)=\theta \frac{V\left(S+\delta S,t\right)-V\left(S-\delta S,t\right)}{2\delta S}+\\ \left(1-\theta \right)\frac{V\left(S+\delta S,t+\delta t\right)-V\left(S-\delta S,t+\delta t\right)}{2\delta S}+O\left({\left(\delta S\right)}^2\right),\end{array}$

[0050] where a function 𝑓(x) is said to be O(g(x)) if $\underset{g\left(x\right)->0}{\lim \sup }\frac{\left|f\left(x\right)\right|}{\left|g\left(x\right)\right|}<\infty$

[0051] and θε[0, 1] is a given parameter whose specialization yields the explicit approximation (θ=0), the implicit approximation (θ=1), and the Crank-Nicolson approximation (θ=½). Based on the above finite difference scheme and the initial and boundary values of V(S, t), the financial engine will derive the unknown option and volatility values at the grid points (ηδS, mδt), for n=1, 2, . . . , N−1 and m=0, 1, 2, . . . M−1, in the (state, time)-product space. Let

[0052] V_mn≡V(nδS, mδt) and σ_mn≡σ(nδS, mδt)

[0053] denote these unknown values. The boundary values

[0054] V_m0≡V(0, mδt), V_mN≡V(NδS, mδt), V_Mn≡V(nδS, T),

[0055] for m=0, . . . , M−1 and n=1, . . . , N−1 are all given. We also let

[0056] Λ_mn≡Λ(nδS, mδt)

[0057] denote the payoff function evaluated at the grid point (ηδS, mδt). For each m, let V_m, σ_mand Λ_m, denote, respectively, the (N−1)-vectors (V_mn)_n−1^N−1, (σ_mn)_n−1^N−1and (Λ_mn)_n−1^N−1. Along with suitable boundary conditions, the PDLCP (3) is approximated by the following (M−1) LCPs, each of order (N−1): at each time t=mδt for m=M−1, M−2, . . . , 1, 0,

[0058] 0≦V_m−Λ_m⊥b_m(σ_m)+Q(σ_m)V_m+N(σ_m)V_m+1≧0  (4)

[0059] where ⊥ is the notation for “perpendicular to” and for an arbitrary (N−1)-vector ω,

[0060] Q(ω)≡(δt⁻¹+r)I_N−1+L_θ(ω), N(σ)≡−δt⁻¹I_N−1+L_1−θ(ω), $b_m\left(\omega \right)\equiv \frac{1}{2}\left(\begin{array}{c}\frac{1}{2}\left[\theta V_{\mathrm{m0}}+\left(1-\theta \right)V_{\left(m+1\right)0}\right]\left(r-D_0-{\omega }_1\right){\omega }_1\\ 0\\ ⋮\\ 0\\ -\frac{1}{2}\left[\theta V_{\mathrm{mN}}+\left(1+\theta \right)V_{\left(m+1\right)N}\right]\left(r-D_0+{\omega }_{N-1}\right){\omega }_{N-1},\end{array}\right)$

[0061] with I_N−1being the identity matrix of order N−1 and L_α(ω) being the (N−1)×(N−1) tridiagonal matrix whose entries are given by: for i, j=1, 2, . . . , N−1, ${\left(L_{\alpha }\left(\omega \right)\right)}_{\mathrm{ij}}\equiv \{\begin{array}{cc}-\frac{\alpha }{2}{i}^2{\omega }_i^2+\frac{\alpha }{2}i\left(r-D_0\right)& \mathrm{if}j=i-1\\ \alpha i^2{\omega }_i^2& \mathrm{if}j=i\\ -\frac{\alpha }{2}{i}^2{\omega }_i^2-\frac{\alpha }{2}i\left(r-D_0\right)& \mathrm{if}j=i+1\\ 0& \mathrm{otherwise}.\end{array}$

[0062] Provided that $\begin{array}{cc}\frac{1}{\delta t}+r\ge N\left(r-D_0\right),& \left(5\right)\end{array}$

[0063] the matrix L_α(ω), and hence Q(σ_m) is strictly diagonally dominant, thus positive definite, for all vectors ω and scalars αε[0, 1]. It is assumed that the condition of equation (5) on δt is satisfied.

[0064] To solve the forward pricing problem where the volatilities σ_mnare all given, the LCPs (4) are time stepped starting with m=M−1; since V_Mis known, by solving the LCP at time t=(M−1)δt, a unique solution V_M−1can be obtained. Proceeding backward in time, a set of discrete option prices V_MNcan be computed, that depend on the input volatilities σ_mn. In order to better illustrate the details of this dependence the M LCPs (4) are written as an aggregate LCP of size M(N−1). Indeed, define the M(N−1)×M(N−1) matrix: $A\left(\sigma \right)\equiv \left[\begin{array}{ccccc}Q\left({\sigma }_0\right)& N\left({\sigma }_0\right)& & & \\ & Q\left({\sigma }_1\right)& N\left({\sigma }_1\right)& & \\ & & ⋰& ⋰& \\ & & & Q\left({\sigma }_{M-2}\right)& N\left({\sigma }_{M-2}\right)\\ & & & & Q\left({\sigma }_{M-1}\right)\end{array}\right]$

[0065] and the M(N−1)-vector $b\left(\sigma \right)\equiv \left[\begin{array}{c}{b}_1\left({\sigma }_1\right)\\ ⋮\\ b_M\left({\sigma }_M\right)\end{array}\right]$

[0066] where σ is the M(N−1)-vector whose entries are the discretized volatilities σ_mnfor m=0, . . . , M−1 and n=1, . . . , N−1: $\sigma \equiv \left(\begin{array}{c}{\sigma }_0\\ ⋮\\ {\sigma }_{M-1}\end{array}\right).$

[0067] Similarly, we define the aggregate M(N−1)-vectors of option prices and payoffs: $V\equiv \left(\begin{array}{c}{V}_1\\ ⋮\\ V_M\end{array}\right)\mathrm{and}\Lambda \equiv \left(\begin{array}{c}{\Lambda }_1\\ ⋮\\ {\Lambda }_M\end{array}\right).$

[0068] The above time-stepping scheme for computing the forward prices of a vanilla American option can be summarized as the following LCP of size M(N−1):

[0069] 0≦V−Λ⊥b(σ)+A(σ)V≧0.  (6)

[0070] Under the condition (5), A(σ), preferably being a block upper triangular matrix with positive definitive diagonal blocks, is a P-matrix, albeit A(σ) is not symmetric or positive definite. A comprehensive treatment of P-matrices and the LCP is described in R. W. Cottle, J. S. Pang, and R. S. Stone, The Linear Complementarity Problem, Academic Press, Boston (1992), which is incorporated herein by reference in its entirety. A(σ) depends on the finite-difference scheme preferably being used to discretize the Black-Scholes partial differential operator L_BSand is independent of the option, i.e., the payoff function and the expiration. Presumably, if a different discretization of this operator is used, a different matrix is obtained. The embodiment presented below applies to other discretization schemes as well, provided that the resulting A(σ) is a P-matrix. In contrast to A(σ), the vectors Λ and b(σ) are both dependent on the option.

[0071] In another embodiment of the invention, multiple options written on the same underlying asset may be present. In such a case the LCP (6) can be embedded in a larger system in accordance with the preferred embodiment of the invention. Specifically, suppose that there are K American options each being characterized by its payoff function Λ^k(S, t) at exercise. The prices of these options are all calculated under the same volatility function σ(S, t) of the asset. Thus for each k=1, . . . , K, the LCP:

[0072] 0≦V^k−Λ^k⊥b^k(σ)+A(σ)V^k≧0,

[0073] is solved to obtain the discretized option prices V^k_mnof type k. Preferably, all these LCPs may be defined by the same matrix A(σ). Concatenating these K LCPs, the final LCP formulation is determined for the forward pricing problem of multiple American options on a single asset:

[0074] 0≦x−p⊥q(σ)+M(σ)x≧0,  (7)

[0075] where $x\equiv {\left(V^k\right)}_{k=1}^K,p\equiv {\left({\Lambda }^k\right)}_{k=1}^K,$

[0076] are, respectively, the KM(N−1)-dimensional vectors of unknown option prices and known payoffs at the discretized grid points, $q\left(\sigma \right)\equiv {\left(b^k\left(\sigma \right)\right)}_{k=1}^K,$

[0077] is the KM(N−1)-dimensional vector that contains the given initial and boundary values of the options, and M(σ) is the KM(N−1)×KM(N−1) block diagonal matrix all of whose K diagonal blocks are equal to the M(N−1)×M(N−1) matrix A(σ).

[0078] Thus, the forward option pricing problem is implemented by the program as the LCP (7), referred to herein as the “forward option LCP.” The forward option LCP preferably requires the discretized volatility matrix σ as an input, while the output from the forward option LCP yields the discrete option prices V^k(nδS, mδt) of K options. By the P-property of the matrix M(σ), the forward option LCP has a unique solution x(σ) that depends on the input volatilities. This option function x(σ) is only implicitly known; for any given σ, x(σ) can be evaluated by solving the forward option LCP.

[0079] The optimization problem overcome by financial engine 100 of the joint computation of the volatility matrix σ and the option price vector x(σ) is subject to a set of criteria that is the result of some input or detected market considerations. Preferably, the constraints on the unknown volatilities are modeled by the set Γ⊂R^M(N−1)in the computer program stored in memory 122. The constraints preferably include, for example, upper and lower bounds on σ and may include any other constraints known or desired in the art. The objective is described by the function θ(σ, x). A variety of objective functions may serve as the fundamental basis for the computations. θ(σ, x) may be, for example, a standard least-squares deviation from observed option prices and historical volatilities. More precisely, a subset of option prices V^obs,k_mnwith (k, m, n) belonging to a subset O of {1, . . . , K}×{0, . . . , M−1}×{1, . . . , N−1} and a subset of volatilities σ^his_mnwith (m, n) belonging to a subset S of {0, . . . , M−1}×{1, . . . , N−1} are given (e.g., the observed market prices and historical volatilities). Knowing (σ, x) allows the computed values to be minimally deviated from the given values. Mathematically, the objective function is then represented by: $\theta \left(\sigma ,x\right)\equiv \sum _{\left(k,m,n\right)\in O}{\left(V_{\mathrm{mn}}^k=V_{\mathrm{mn}}^{\mathrm{obs},k}\right)}^2+\sum _{\left(m,n\right)\in S}{\left({\sigma }_{\mathrm{mn}}-{\sigma }_{\mathrm{mn}}^{\mathrm{his}}\right)}^2,$

[0080] where V^k_mnis the theoretical option price (collectively, these prices are the components of x). Another example of θ would be a function similar to the one used in T. Coleman, Y. Li and A. Verma, “Reconstructing the unknown local volatility function,” Risk(1998), which is incorporated herein by reference in its entirety.

[0081] In a preferred embodiment, financial engine 100 computes a discretization of the volatility surface so that a prescribed objective function is minimized subject to given constraints on the unknown volatilities and option prices that includes a measure of the smoothness of the discretized volatility surface. Variations of these functions and/or other realistic objective functions are amenable to the same implementation.

[0082] In accordance with a preferred embodiment of the invention, with the set Γ and the objective function θ given, the implied volatility problem of American options is implemented by the computer program as the following constrained optimization problem: compute (σ, x) to

[0083] minimize θ(σ, x)

[0084] subject to σεΓ

[0085] and 0≦x−p⊥q(σ)+M(σ)x≧0.  (8)

[0086] In terms of the option function x(σ), this optimization problem may be rewritten as an implicit program in the variable σ alone:

[0087] minimize φ(σ)≡θ(x(σ), σ)

[0088] subject to σεΓ.

[0089] The terms σ and x are vectors in the Euclidean spaces R^M(N−1)and R^KM(N−1), respectively. Let

[0090] F(σ, x)≡q(σ)+M(σ)x

[0091] be the function defining the option LCP (7). The partial Jacobian matrix of F(σ, x) with respect to a is denoted J_σF(σ, x); since F(σ, x) is linear in x, the partial Jacobian matrix of F(σ, x) with respect to x is equal to the matrix M(σ). For an arbitrary σ, three index sets can be defined associated with the solution x(σ):

[0092] α(σ)≡{i:(x(σ)−p)_i>0=(q(σ)+M(σ)x)_i}

[0093] β(σ)≡{i:(x(σ)−p)_i=0=(q(σ)+M(σ)x)_i}

[0094] γ(σ)≡{i:(x(σ)−p)_i<0=(q(σ)+M(σ)x)_i}.

[0095] These index sets play an important role in the local properties of x(σ) when σ undergoes small perturbations. The particular case where β(σ), called the degenerate set, is the empty set is particularly noteworthy. This case corresponds to the solution x(σ) being nondegenerate. As seen from the theorem below, x(.) is then locally smooth around this value σ. Similar to the “big O” notation, 𝑓(x)=o(g(x)) can be written if $\underset{g\left(x\right)->0}{\lim }\frac{f\left(x\right)}{g\left(x\right)}=0.$

[0096] The following result summarizes the key properties of the option function x(σ).

[0097] Assuming that the condition (5) holds, the option function x(σ) is Lipschitz continuous and directionally differentiable in its argument σεR^M(N−1). The directional derivative of x(σ) at any σεR^M(N−1)along any direction dσεR^M(N−1), denoted x′(σ; dσ), is the unique solution dxεR^KM(N−1)to the following mixed LCP: $\begin{array}{cc}\left\{\begin{array}{c}{\left(J_{\sigma }F\left(\sigma ,x\left(\sigma \right)\right)d\sigma +M\left(\sigma \right)dx\right)}_i=0,\forall i\in \alpha \left(\sigma \right),\\ 0\le {\left(dx\right)}_i\perp {\left(J_{\sigma }F\left(\sigma ,x\left(\sigma \right)\right)d\sigma +M\left(\sigma \right)dx\right)}_i\ge 0,\\ \forall i\in \beta \left(\sigma \right),\\ {\left(dx\right)}_i=0,\forall i\in \gamma \left(\sigma \right).\end{array}\right\}& \left(10\right)\end{array}$

[0098] Furthermore,

[0099] x=(σ+dσ)=x(σ)+x′(σ; dσ)+o(∥dσ∥).  (11)

[0100] Finally, if β(σ) is empty, then x is Fréchet differentiable at σ.

[0101] In a preferred embodiment, the implied volatility problem as formulated as either (8) or (9) preferably has a solution if the objective function θ(σ, x) is continuous and the feasible volatility region Γ is compact.

[0102] Thus, in a preferred embodiment where the objective function θ(σ, x) is continuous and Γ is compact, then an optimal solution to (8) exists.

[0103] A consequence of programming financial engine 100 in this manner is that it provides a computationally effective way of approximately updating a set of option values when the volatilities undergo small perturbations, without recomputing the exact values from scratch. Indeed, the formula (11) yields

[0104] x(σ′) approx. equal to ∝x(σ)+x′(σ; dσ), where dσ≡σ′−σ

[0105] preferably with σ′ being a sufficiently small perturbation of σ. Thus, financial engine 100 can compute an approximation to x(σ′) from x(σ) by simply computing the directional derivative x′(σ, dσ). In general, the computation of this derivative involves the solution of a mixed LCP of reduced size. The latter problem may be of reduced size because the last equation in (10) allows the γ-components of dx to be fixed at zero. For example, when x(σ) is a nondegenerate solution (i.e., when β(σ) is empty), the mixed LCP (10) reduces to a single system of linear equations involving only the variables dx_ifor iεα(σ). In the general case, the mixed LCP (10) is preferably solved as stated or converted into a standard LCP involving only the variables dx_ifor iεβ(σ), by using the first equation in (10) to eliminate the variables dx_ifor all iεα(σ).

[0106] Since the option function x(σ) is preferably not a smooth function, (9) is a nonsmooth optimization problem in general. Nevertheless, if the original objective function θ(σ, x) is smooth, then the composite objective function φ(σ) is “B(ouligand)-smooth,” meaning that it is locally Lipschitz continuous and directionally differentiable. As such, a descent method can be applied in solving the problem (9) in accordance with a preferred embodiment of the invention. Furthermore, the “B-derivative” of φ can be used to describe the stationarity conditions of the program (9).

[0107] The operation of the financial engine 100 in accordance with a preferred embodiment of the invention can be used to implement the method of pricing financial derivatives such as vanilla American options shown in FIG. 2. In the retrieval step S20, for example, both data source 10 and user interface device 18 can be used to retrieve data to be input into one or more of the modules contained in financial engine 100. In one embodiment, for example, initial parameters such as the price of the underlying asset (e.g., stock price), risk-free interest rate, maturity date of derivative, dividend yield, and initial volatility parameter can be input to forward pricing unit 16 through data source 10 and/or user interface device 18. Based on this information, initial option pricing data can be determined and output by forward pricing unit 16 for use by other modules in financial engine 100, including volatility calculator 12 (as represented by the broken line path in FIG. 2). Volatility calculator 12, for example, may utilize the theoretically calculated option price for comparison with observed option prices input from data source 10 (or user interface device 18) in optimizing its calculations (in step S22) of the volatility implied from the market, as will be discussed in detail below.

[0108] Once optimized, the implied volatilities calculated in S22 are used as inputs by forward pricing unit 16, together with additional inputs to derive the forward price of the subject derivative (e.g., vanilla American option). The forward pricing information can then be included in a report output to the user via, for example, user interface device 18, or via network 110 (in step S26). The may include the volatility calculation itself, the forward pricing data, additional calculations based on the calculated data (e.g., positions taken, projected, predicted, simulated, etc.), and other pertinent data. The calculated volatility and pricing information may also be used in effecting trading by issuing control signals through financial instrument controller 14 or the like (step S26) to trading systems (represented by network 110) based on the calculations made in steps S22 and S24.

[0109] The calculations made within step S22 may be performed, in accordance with a preferred embodiment of the invention, through the use of computer program or algorithm referred to herein as “Implicit Programming Algorithm (IMPA)” of the invention. In a preferred embodiment, IMPA is used by financial engine 100 to solve the implied volatility problem based on the optimization formulation (9). For illustration purposes, the feasible volatility region Γ is assumed to be a relatively simple set such as a polyhedron; as a result, the directional subprograms in IMPA can be easily solved.

[0110] With φ(σ) being directionally differentiable, an iterative descent algorithm of the sequential quadratic programming kind is employed to minimize the formulation (9). In describing IMPA we let σ^vεΓ be a given volatility vector that is not stationary for (9). Write x^v≡x(σ^v) and

[0111] α_v≡α(σ^v), β_v≡β(σ^v), and γ_v≡γ(σ^v).

[0112] A descent direction of the objective function φ(σ) at σ^vcan be generated by solving the following minimization problem: for any symmetric positive definite matrix

[0113] Q_vεR^{M(N−1)×M(N−1)},

[0114] minimize φ′(σ^v; dσ)+½dσ^TQ_vdσ

[0115] subject to σ^v+dσεΓ.

[0116] Since

[0117] φ′(σ^v; dσ)=∇_σφ(σ^v, x^v)^Tdσ+∇_xφ(σ^v, x^v)^Tx′(σ^v; dσ),

[0118] the above directional minimization problem is equivalent to:

[0119] minimize ∇_σφ(σ^v, x^v)^Tdσ+∇_xφ(σ^v, x^v)^Tdx+½dσ^TQ_vdσ

[0120] subject to σ^v+dσεΓ

[0121] (J_σF(σ^v, x^v)dσ+M(σ^v)dx)_i=0, ∀iεα_v,

[0122] 0≦(dx)_i⊥(J_σF(σ^v, x^v)dσ+M(σ^v)dx)_i≧0, ∀iεβ_v,

[0123] (dx)_i=0, ∀iεγ_v,  (12)

[0124] where the directional derivative x′(σ^v; dσ) has been substituted. The latter optimization problem is another MPEC whose solution is not trivial. A preferred embodiment provides a solution which simplifies directional computation while providing accurate results for implied volatility inputs. The method of solving the optimization problem is an implicit programming algorithm (IMPA) formulated in accordance with the invention. The following additional notation is used in the description of the use of IMPA in the preferred embodiment illustrated in FIG. 3:

[0125] φ_v≡φ(σ^v, x^v), dφ_σ^v≡∇_σφ(σ^v, x^v), dφ_x^v≡∇_xφ(σ^v, x^v), and J_σF^v≡J_σF(σ^v, x^v).

[0126] An initialization step S30 initializes parameters ρ, γε(0, 1) as given scalars and Q₀as a given symmetric positive definite matrix of order M(N−1). σ⁰εΓ is also given, and v=0 is set.

[0127] A direction generation step S32 allows β′ to be an arbitrary subset of β_v, and solves the convex quadratic program

[0128] minimize (dφ_σ^v)^Tdσ+(dφ_x^v)^Tdx+½dσ^TQ_vdσ

[0129] subject to σ^v+dσεΓ

[0130] (J_σF^vdσ+M^vdx)_i=0, ∀iεα_v∪β′,

[0131] (dx)_i=0, ∀iεγ_v∪(β_v\β′),  (13)

[0132] to obtain the search direction dσ and the auxiliary vector dx, both of which is preferably necessarily unique.

[0133] In the termination test step S34 a check is made to determine if the computed implied volatility is within acceptable values, e.g., if

[0134] ∥dσ∥≦prescribed tolerance,

[0135] IMPA is stopped, as shown in step S36, and the pair σ^v, x^vis taken to be the desired approximate solution of(9). If

[0136] ∥dσ∥>prescribed tolerance

[0137] IMPA is continued with step S38. The step size determination step S38 requires that with

[0138] σ^v(τ)≡σ^v+τdσ and x^v(τ)≡x(σ^v(τ)), ∀τε[0, 1],  (14)

[0139] set

[0140] τ_v≡max(0.01, ρ^l)

[0141] where l_vis the smallest nonnegative integer l such that with τ=ρ^l,

[0142] φ(σ^v(τ), x^v(τ)−φ_v≦γτ[(dφ_σ^v)^Tdσ+(dφ_x^v)^Tdx].

[0143] Set σ^p+1≡σ^v(T_v). A symmetric positive definite matrix Q_v+1is chosen and v is replaced by v+1. IMPA then returns to step S30.

[0144] The choice of the matrices {Q} could presumably affect the practical performance of the algorithm. Preferably, the matrices {Q} are chosen to reflect some second-order information of the objective function φ(σ, x). In a preferred embodiment, a lower bound of 0.01 on the step size is used. The step size determination is preferably the well-known Armijo inexact line search rule in standard unconstrained optimization algorithms, although other known rules may be used. A convergence analysis of the algorithm without the modification of the directional step may be used as described in Z. Q. Luo, J. S. Pang, and D. Ralph, Mathematical Program with Equilibrium Constraints, Cambridge University Press, Cambridge (1996), which is incorporated herein by reference in its entirety.

[0145] In another embodiment of the invention, a computer program or algorithm for solving the implied volatility problem based upon formulation (8) is provided, which is referred to herein as the “Penalty Interior Point Algorithm (PIPA),” and executed by financial engine 100. For illustration purposes, the feasible volatility region Γ is assumed to be a relatively simple set such as a polyhedron; as a result, the directional subprograms in PIPA can be easily solved. Unlike IMPA (described above) that originates from an iterative descent algorithm for solving smooth nonlinear programs, PIPA originates from a penalty method coupled with an interior point routine to deal with complementarity constraint in the other formulation (8) of the implied volatility problem.

[0146] In a preferred embodiment, PIPA is executed by volatility calculator 12 to operate on the following reformulation of (8):

[0147] minimize θ(σ, x)

[0148] subject to σεΓ

[0149] w−q(σ)−M(σ)x=0

[0150] w◯(x−p)=0

[0151] (w, x−p)≧0

[0152] where ◯ denotes the Hadamard operator on two vectors; that is, w◯x is the vector whose i-th component is equal to the product of the i-th components of w and x. The following additional notation is useful in describe the use of PIPA in a preferred embodiment of the invention, as shown in FIG. 4.

[0153] G(σ, x, w)≡w−q(σ)−M(σ)x=w−F(σ, x)

[0154] In addition, the residual function of the complementarity conditions is preferably defined as:

[0155] r(σ, x, w)≡(G(σ, x, w))^TG(σ, x, w)+(x−p)^Tw

[0156] The penalized objective function is preferably defined as:

[0157] P_c(σ, x, w)≡θ(σ, x)+cr(σ, x, w),

[0158] where c>0 is a penalty scalar to be adjusted. For a given vector x, diag(x) is assumed to be the diagonal matrix whose diagonal entries are the components of this vector.

[0159] PIPA is an iterative algorithm that generates a sequence of iterates {(σ^v, x^v, w^v)} preferably satisfying the following conditions for all v:

[0160] a. (feasible volatilities) σ^vεΓ;

[0161] b. (positivity of state variables) x^v−p>0 and w^v>0; and

[0162] c. (centrality condition) (x^v−p)◯w^v≧ηg_v1_M(N−1), where ηε(0, 1) is a given scalar, 1_M(N−1)is the M(N−1)-vector of all ones, and $g_v\equiv \frac{{\left(x^v-p\right)}^T\left(w^v\right)}{M\left(N-1\right)}$

[0163] is the average “complementary gap” between x ^vand w^v.

[0164] While maintaining these conditions, PIPA as executed by volatility calculator 12 attempts to decrease the objective function θ(σ, x) by reducing the penalty function Pc(σ, x, w). In a preferred embodiment, this can be accomplished via an Armijo inexact line search using a properly defined search direction and a proper choice of the penalty scalar c. PIPA drives the residual r(σ, x, w) toward zero, thereby achieving feasibility to the MPEC, and obtains a satisfactory “minimum” value of the objective θ(σ, x).

[0165] At the beginning of each iteration v, a triple σ^v, x^v, w^vsatisfying the above conditions is preferably given. Consistent with the notation used so far, the triple can be written:

[0166] G^v≡G(σ^v, x^v, w^v) and r_v≡r(σ^v, x^v, w^v).

[0167] At this triple, a Newton linearization step is applied to the equation:

[0168] w−q(σ)−M(σ)x=0

[0169] and also to the following perturbed complementarity equation:

[0170] w◯(w−q)=κg_v1_M(N−1),

[0171] where κ is a given constant in the interval (0, 1). The resulting linear equations are preferably used to define a directional quadratic program whose solution yields a search direction along which the aforementioned line search is carried out.

[0172] As shown in FIG. 4, running PIPA calculator 12 initializes (8) in initialization step S40 by letting ρ,γ,κ₀, and η be given scalars in the interval (0, 1). In addition, c_−≧1 and ζ>0 are given scalars, and H_ois a given symmetric positive (semi) definite matrix of order 2M(N−1). (σ⁰, x⁰, w⁰) is a given triple satisfying the three conditions (a, b, c), while v is set to 0.

[0173] In the direction generation step S42 is applied to solve the convex quadratic program $\mathrm{minimize}{\left(d{\varphi }_{\sigma }^v\right)}^{T}d\sigma +{\left(d{\varphi }_x^v\right)}^{T}dx+\frac{1}{2}{\left(\begin{array}{c}d\sigma \\ dx\end{array}\right)}^T{H}_v\left(\begin{array}{c}d\sigma \\ dx\end{array}\right)$  subject to σ^v+dσεΓ

[0174] G^v+dw−J_σF^vdσ−M^vdx=0

[0175] (x^v−p)◯w^v+diag(x^v)dw+diag(w^v)dx=κ_vg_v1_M(N−1)

[0176] ∥dσ∥_∞≦ζr_v  (15)

[0177] to obtain the search triple (dσ, dx, dw).

[0178] PIPA then applies a termination test in step S44 to determine if the computed implied volatility is within acceptable values. If

[0179] ∥(dσ, dx, dw)∥+r_v≦prescribed tolerance,

[0180] PIPA as shown in step S46, and the triple (σ^v, x^v) is taken to be a desired approximate solution of (8).

[0181] If the termination condition is not met, the PIPA of the invention applies a penalty update in step S47. In this step s_v≧1 is taken as the smallest integer s≧1 such that

[0182] (dφ_σ^v)^Tdσ+(dφ_σ^v)^Tdx−c_v−1^s[2(G^v)^TG^v+(1−κ_v)(x^v)^Tw^v]<−r_v.

[0183] Set c_v≡c_v−1^sv.

[0184] In step S48 PIPA determines the next step size, preferably using centrality, with $\left(\begin{array}{c}{\sigma }^v\left(\tau \right)\\ x^v\left(\tau \right)\\ w^v\left(\tau \right)\end{array}\right)\equiv \left(\begin{array}{c}{\sigma }^v\\ x^v\\ w^v\end{array}\right)+\tau \left(\begin{array}{c}d\sigma \\ dx\\ dw\end{array}\right)$

[0185] determine the largest τε(0, 1] such that

[0186] x^v(τ)◯w^v(τ)≧ηg_v(τ)1_M(N−1),

[0187] where $g_v\left(\tau \right)=\frac{{\left(x^v\left(\tau \right)-p\right)}^T\left(w^v\left(\tau \right)\right)}{M\left(N-1\right)}.$

[0188] Set ${\tau }_v^1\equiv 0.9999\tau .$

[0189] In step s49, calculator 12 preferably determines step size using an Armijo line search by setting

[0190] τ_v≡τ_v′ρ^l

[0191] where l_vis the smallest nonnegative integer l such that with $\tau ={\tau }_v^1{\rho }^l,$  P_c(σ^v(τ), x^v(τ), w^v(τ))−P_c(σ^v, x^v, w^v)< $\gamma \tau \left[{\left(d{\varphi }_{\sigma }^v\right)}^{T}d\sigma +{\left(d{\varphi }_x^v\right)}^T\mathrm{dx}-c_{v-1}^s\left(2{\left(G^v\right)}^T{G}^v+\left(1-{\kappa }_v\right){\left(x^v\right)}^T{w}^v\right)\right].$

[0192] Set (σ^v+1, x^v+1, w^v+1)≡(σ^v(τ_v), x^v(τ_v), w^v(τ_v)). Choose a scalar κ_v+1ε(0, κ_v] and a symmetric positive (semi) definite matrix H_v+1. Return to Step 1 with v replaced by v+1.

[0193] It should be noted that the program or algorithm such as IMPA or PIPA, executed by calculator 12 can be written in different computer languages for different computer modules or systems. The programs can be stored on a hard drive, floppy disk, CD-ROM or other permanent or semi-permanent storage medium, as well as downloaded, for example, from a server computer, or the Internet, or transmitted as a data signal over a communications path, e.g., over the Internet or on a transmitted carrier wave, as is known in the art.

[0194] The following examples were performed through a computer program written using MATLAB® to illustrate the numerical performance of algorithms such as IMPA and PIPA, in accordance with preferred embodiments of the invention. The inverse pricing of Kε{2, 11, 14} American put options as examples in order to determine the discretized volatilities implied by certain given option values. Thus for k=1, . . . , K, the payoff function of option k is given by:

[0195] Λ^K(S, t)≡max(S−E_k, 0)

[0196] where E_kis the strike price of the option at its expiry date T_k. For each option k, a set of discretized American option prices V_n^obs,k≡V_n^obs,k(nδS, 0), representing the set of currently observable option prices, is generated by solving the single-option LCP (4) using a constant volatility of σ₀=0.04 and the following parameters: N≡16, δS=1, T=T_k, δt=0.125, r=0.05, and D₀=0.02.

[0197] The constant σ₀is not known to IMPA or PIPA; instead, we use the values V_n^obs,kwere used to define the following least-squares objective function: $\theta \left(x\right)\equiv \frac{1}{2}\sum _{k=1}^K\sum _{n\in N_k}{\left(V^k\left(n\delta S,0\right)-V_n^{\mathrm{obs},k}\right)}^2.$

[0198] The index set N_kis a singleton in all the runs; it consists of {n′} such that n′δS≡S₀the current asset price. With the objective of minimizing θ(x), the following examples show the use of IMPA or PIPA in obtaining an objective value of θ close to zero. In addition, Γ is taken to be a rectangular box, meaning that the unknown volatility vector σ is subject to simple upper and lower bounds.

[0199] In implementing IMPA in accordance with a preferred embodiment, the major computations in each iteration were (i) the direction generation step which involves solving the convex quadratic program (13), and (ii) the evaluation of x^v(τ) for various step sizes τ=ρl^v. By using the equality constraints in (13), the variables dx_ifor iεx_vUβ′ can be solved in terms of dσ, thereby converting (13) into a strictly convex quadratic program in the variable dσ alone. In the experiments below, each Qv is chosen to be the identity matrix to simplify the resulting quadratic program.

[0200] The evaluation of x^v(τ) involves solving the forward option LCP (7) corresponding to σ^v(τ) This LCP decomposes into K independent linear complementarity subproblems each corresponding to a given option. In turn, as mentioned above, each single-option LCP can be solved backward in time. Thus, the evaluation of x^v(τ) amounts to solving KM time-stepped LCPs each of the form (4). In principle, there are many methods applicable for solving the latter LCP (4). For example, when N (the number of discretized asset prices) is large, an interior point algorithm is the preferred method for solving the LCP (4). Thus, this interior point algorithm is chosen as the principal tool for evaluating x^v(τ).

[0201] The main computation in each iteration of PIPA is the solution of the convex quadratic program (15). Unlike IMPA, there is no need to solve herein any forward option LCP. As in IMPA, the program (15) can be converted into one in the variable do, alone. In both cases, the highly special structure of the matrices involved is exploited to facilitate the computations.

[0202] The termination test for IMPA is

[0203] min(∥dσ∥, θ(x^v, σ^v))≦10⁻⁶  (16)

[0204] A maximum number of 80 iterations has also been imposed. Since IMPA maintains the feasibility of the pair (σ^v, x^v) to the MPEC (8), if IMPA terminates after 80 iterations without satisfying (16), an implied volatility vector σ^vεΓ is obtained with associated American option prices x^vthat are deemed satisfactory. The resulting objective values θ(x^v) at termination of IMPA in all runs were invariably very small.

[0205] The termination test for PIPA is

[0206] ∥(dσ, dx, dw)∥+r_v≦10⁻⁸.  (17)

[0207] At termination, the iterate (σ^v, x^v) produced by PIPA is an approximate feasible solution to the MPEC (8) with the feasibility accuracy being less than 10⁻⁸. The objective value θ(x^v) is reported in each run.

[0208] Two observed put option values are given (K=2):

[0209] The current asset price is S₀=8. An initial value σ_init=0.255 is used to start both algorithms. The discretized volatility surfaces computed by PIPA and IMPA are shown in FIGS. 5 and 6, respectively. The outputs from the two algorithms are summarized in Table 1 below. Note that the objective value obtained by PIPA is larger than that obtained by IMPA even though the residual is quite small. The two surfaces in FIG. 5 and FIG. 6 have similar shape, but are not identical.

[0210] Since PIPA and IMPA did not produce identical surfaces and neither recovered the original constant volatility surface used to generate the observed option data, as a test to determine how close the option prices calculated from given volatility are being reproduced, the volatility surfaces obtained by PIPA and IMPA are used now as parameters in the forward problem to re-calculate the option values. For given strike prices E=8 and E=9, the three set of option values were computed and plotted using: 1) σ(S, t)=0.4, 2) σ(S, t)=σ_PIPA(S, t), and 3) σ(S, t)=σ_IMPA(S, t), shown in FIGS. 7 and 8. It is apparent that the IMPA of the invention reproduces the option prices V_obsrather closely.

[0211] Eleven observed put option values are given (K=11):

[0212] The corresponding parameters are the same as those given in Example 1.

[0213] The discretized volatility surfaces computed by PIPA and IMPA are shown in FIGS. 9 and 10, respectively. Again, the re-calculated option values were plotted for strike prices E=8 and E=9, as shown in FIGS. 11 and 12. In this example, both. IMPA and PIPA reproduce the observed option curve very closely.

[0214] Fourteen observed put option values are given (K=14):

[0215] Following the procedure of Example 1, the following results were obtained:

[0216] The discretized volatility surfaces computed by PIPA and IMPA are shown in FIGS. 13 and 14, respectively. The re-calculated option values for strike prices E=8 and E=9 are plotted in FIGS. 15 and 16.

[0217] In all three examples, the volatility surface calculated by PIPA is lower than that calculated by IMPA. This is also reflected in the plots of the option values: the option values calculated from σ_PIPA(S, t) are below those calculated from σ₀=0.4 and σ_IMPA.

[0218] Again taking the option values calculated from the constant volatility (σ₀=0.4) ${\tilde{V}}_n^{\mathrm{obs},k}\equiv V_n^{\mathrm{obs},k}+0.0707$

[0219] to be the observed values, each observed values is perturbed by a small amount; i.e.

[0220] Thus, the new objective function is given by $\begin{array}{c}\theta \left(x\right)=\frac{1}{2}\sum _{k=1}^K\sum _{n\in N_k}{\left(V^k\left(n\delta S,0\right)-{\tilde{V}}_n^{\mathrm{obs},k}\right)}^2\\ =\frac{1}{2}\sum _{k=1}^K\sum _{n\in N_k}{\left(V^k\left(n\delta S,0\right)-{\tilde{V}}_n^{\mathrm{obs},k}-0.0707\right)}^2\end{array}$

[0221] Clearly, σ(S, t)=σ₀for all (S, t) is a feasible solution to the optimization problem (9). This gives an upper bound on the optimal objective value:

[0222] θ_optimal≦K(0.0707)²≡UB.

[0223] Examples 1-3 were repeated using this new objective function to investigate how well the IMPA and PIPA of the invention perform in relation to this upper bound value UB, starting at σ_init=0.255. The results are summarized below (with the “IMPA” results achieved after 80 iterations):

[0224] While preferred embodiments of the invention have been described and illustrated, it should be apparent that many modifications to the embodiments and implementations of the invention can be made without departing from the spirit or scope of the invention. For example, while only vanilla American options are explained in the interest of simplicity, the same general approach can be applied to computing volatilities implied by exotic American options and/or American options with transaction costs ,as well as the inverse pricing of other financial derivatives not described herein. The modules illustrated in FIG. 1 as making up financial engine 100 may be one or more hardware, software, or hybrid components residing in (or distributed among) one or more local or remote computer systems. Although the modules are shown as physically separated components, it should be readily apparent that the modules may be combined or further separated into a variety of different components, sharing different resources (including processing units, memory, clock devices, software routines, etc.) as required for the particular implementation of the embodiment. Indeed, even a single general purpose computer executing a computer program to produce the functionality described herein may be utilized to implement the illustrated embodiments. User interface device 18 may be any device used to input and/or output information during an exchange of information between user and financial engine 100. The interface device 18 may be implemented as a graphical user interface (GUI) containing a display or the like, or may be a link to other user input/output devices known in the art.

[0225] In addition, memory units 122, 162 described herein as used respectively in volatility calculator 12 and forward pricing unit 16 may be any one or more of the known storage devices (e.g., Random Access Memory (RAM), Read Only Memory (ROM), hard disk drive (HDD), floppy drive, zip drive, compact disk-ROM, DVD, bubble memory, etc.), and may also be memory devices embedded within CPUs 120, 162, or shared with one or more of the other components within financial engine 100. The computer programs or algorithms described herein may easily be configured as one or more hardware modules, and the hardware module shown may easily be configured as one or more software modules without departing from the invention. Although the illustrated embodiments may have described the use the class of iterative algorithms in implementing the financial engine to approach the linear complementarity issues forming part of the option pricing problem, it should be readily apparent that other known algorithms such as the class of pivotal algorithms may also be used in any implementation of the invention. Accordingly, the invention is not limited by the foregoing description or drawings, but is only limited by the scope of the appended claims.