Method and system for modeling financial markets and assets using fractal activity time

[0002] The invention relates generally to the field of mathematical modeling of random or stochastic systems, and more particularly to modeling of the behavior of financial markets and instruments.

[0003] Over the past thirty years many risky asset models have been proposed, but these are generally deficient in that they do not adequately explain key features from the statistical evidence related to the issues of heavy tails and long-range dependence (LRD) in financial time series. With regard to long-range dependence, there have been problems even with understanding the nature of the characteristic, let alone its assessment.

[0004] In particular, the Black-Scholes model for the price of a risky asset, which is widely used in the finance industry, makes two fundamental assumptions which are seriously flawed. They are:

[0005] (A) the differences of log(price) over time are “normally” distributed (i.e., the differences are a Gaussian process); and

[0006] (B) the differences of log(price) over different time periods are independent.

[0007] In order to elaborate on the flaws of the conventional model, it is helpful to introduce some notation. First, a value X_tis defined as log P_t−log P_t−1, where P_tis the asset price at time t. Then, according to assumption (A), the values {X_t, t=1,2, . . . } are normally distributed and, according to assumption (B), the values are independent.

[0008] The problem with assumption (A) can readily be seen by studying real data. For normally distributed data, values which are more than 5 standard deviations away from the mean are very rare, with a probability less than 10⁻⁶. However, empirical studies of actual risky asset data demonstrate that observations this far away from the mean are not uncommon. The data exhibit a pronounced leptokurtic distribution (much higher peaks and heavier tails than Gaussian). In fact, large values can often be observed out to 10 standard deviations. To put these observations into perspective, consider that a 10⁻⁶probability corresponds to a once-in-4000-years event for a normal (i.e., Gaussian) distribution of transaction data collected 5 days per week.

[0009] The problem with assumption (B) can readily be seen by examining lagged autocorrelations of the {X_t} process. While estimated values of corr(X_t, X_t+k) die out quickly as k increases, the values of corr(|X_t|,|X_t+k|) and corr(X_t², X_t+k²) remain significantly different from zero for long periods. In fact, non-zero values can be observed over periods of many years for some series. This has been convincingly illustrated by analysis of daily returns data from the S&P 500 index, the Nikkei index, foreign exchange rates (DM/US$), and Chevron Stock. See, e.g., Z. Ding & C. W. J. Granger, Modeling volatility persistence of speculative returns: A new approach,73 J. Enconometrics 185 (1996); C. W. J. Granger & Z. Ding, Varieties of long memory models,73 J. Econometrics 61 (1996). In the case of the absolute values from the S&P series, for example, the lagged autocorrelations are significantly different from zero for lags up to 2700, which is more than 10 years. The absolute values and squares of the returns exhibit LRD, classically defined (for stationary finite variance processes) as existing if Σ_k=0^∞ γ_kdiverges, γ_kbeing the autocorrelation at lag k. For a discussion of the definition and its extensions see C. C. Heyde and Y. Yang, On defining long-range dependence,34 J. Applied Prob. 939 (1997). For a comprehensive overview of the subject of LRD see J. Beran, Statistics for Long-Memory Processes,Chapman & Hall, New York (1994).

[0010] Yet, according to assumption (B), the functions corr(X_t, X_t+k), corr(|X_t|,|X_t+k|), and corr(X_t², X_t+k²) should all be zero for all k≧1. In other words, assumption (B) is belied by the fact that the corr(|X_t|,|X_t+k|) and corr(X_t², X_t+k²) functions exhibit strong dependence effects when applied to real, observed data. In addition, the actual data form a time series with high volatility and “intermittency” (clustering of extremes) quite unlike white noise. Failure to take these effects into account can adversely affect almost all probabilistic assessments involving price movements.

[0011] It is now widely appreciated that heavier tails than the normal are a necessary model feature for marginal distributions of returns. However, models with finite moments of all orders are still commonly advocated. Detailed comparisons using market index data have been given in S. R. Hurst, E. Platen, & S. R. Rachev, Subordinated Markov models: a comparison,Fin. Eng. Japanese Markets 97 (1997); S. R. Hurst & Platen, E., The marginal distribution of returns and volatility, inL₁—Statistical Procedures and Related Topics 301 (IMS Lecture Notes-Monograph Series, Vol. 31, IMS, Hayward, Calif., Y. Dodge ed., 1997). They favor a t-distribution with degrees of freedom v typically in the range 3-5. This, of course, implies an infinite k-th moment for k≧v.

[0012] The above-described feature of infinite k-th moment has been observed in actual financial data. For example, multiple tests on daily data from the well-known S&P 500 Index lead to the conclusion that the 4th moment is infinite. See, e.g., T. Mikosch, Heavy tails in finance,57 Bull. Internat. Statist. Inst. 109 (Book 2) (1997); Y. Yang, Long Range Dependence in the Study of Time Series with Finite or Infinite Variances,Ph.D. Dissertation, Columbia University (1998).

[0013] Most approaches to the marginal distribution of returns have involved the use of parametric families and the choice of a best fit within a family. Few authors have investigated general questions, such as whether particular moments are finite. Part of the difficulty has been a general supposition that tail behavior would be a power function, at least asymptotically, and that estimating the power index is problematic. Indeed, the usual method, based on the Hill estimator, is founded on order statistic theory for independent and identically distributed (hereinafter “iid”) samples and can suffer from major shortcomings outside that context. See for example, the “Hill Horror Plot” and the attendant discussion in S. Resnick, Heavy tail modeling and teletraffic data,25 Ann. Statist. 1805 (1997).

[0014] Accordingly, it is an object of the invention to provide a new mathematical model with better capability of predicting the behavior of random or stochastic processes such as the behavior of financial markets and instruments. The new model replaces the time parameter of the Black-Scholes model with a fractal activity time, which is a random activity time process having a heavy-tailed distribution, LRD, and stationary differences. The fractal activity time is standardized so that its mean value corresponds to the time parameter used in the Black-Scholes model.

[0015] In accordance with the invention, it is possible to test for finiteness of moments based on asymptotics obtained via the use of the ergodic theorem if the assumption of stationarity of the returns is retained. In particular, if Σ_i=1ⁿ|X_i|^p/n converges as n→∞, then E|X|^p<∞, where E|X|^pdenotes the mean of |X|^p. If, on the other hand, max_1≦i≦n|X_i|^p/Σ_i=1ⁿ|X_i|^pstays away from zero as n→∞, then E|X|^p=∞. Plots as n increases for large data sets should give a clear picture.

[0016] In a method according to the invention, a random process is modeled by providing time-parameter data which corresponds to a fractal activity time process or an approximation thereof. The time-parameter data is used to determine whether the random process exhibits fractal scaling or self-similarity, in which case it is determined that the random process can be validly described using a model based on the time-parameter data.

[0017] Further objects, features, and advantages of the invention will become apparent from the following detailed description taken in conjunction with the accompanying figures showing illustrative embodiments of the invention.

[0018] The invention is explained in greater detail below by reference to the drawings, in which:

[0019] FIG. 1 is a flow diagram of an exemplary procedure for modeling the behavior of a random process in accordance with the invention;

[0020] FIG. 2 is a graph illustrating an autocorrelation function of typical financial data;

[0021] FIG. 3ais a graph illustrating the behavior of an exemplary fractal activity time in accordance with the invention, plotted with respect to a conventional time parameter;

[0022] FIG. 3bis a graph illustrating the behavior of a {T_t−t} process, plotted with respect to a conventional time parameter t, where T_tis an exemplary fractal activity time in accordance with the invention;

[0023] FIG. 3cis a graph illustrating the probability densities of a {T₁−t} process and a standard Gaussian process, where T₁is a scaled generalized inverse Gaussian distribution;

[0024] FIG. 4 is a graph illustrating an exemplary moment function of typical financial data;

[0025] FIGS. 5a-fare graphs illustrating the linearity of a log moment function generated by a method according to the invention;

[0026] FIGS. 6a-iare graphs illustrating the linearity of another function generated by a method according to the invention;

[0027] FIG. 7 is a pictorial diagram of an exemplary computer system suitable for use in practicing the invention; and

[0028] FIG. 8 is a block diagram illustrating a computer section suitable for use as a component of the computer system of FIG. 7.

[0029] The New Model

[0030] According to the commonly used Black-Scholes model, the price P_tof a risky asset at time t is:

[0031] P_t=P₀exp[μt+σW(t)]  (1)

[0032] where P₀is the price at t=0, σ²>0, μ and σ are fixed constants, and W(t) is a standard Brownian motion function with respect to time t. Under this model, the corresponding log returns

[0033] X_t=log P_t−log P_t−1=μ+σ(W(t)−W(t−1))  (2)

[0034] are iid Gaussian with mean μ and variance σ².

[0035] The Black-Scholes model usually does not accurately represent the behavior of actual log returns data, since actual data tend to exhibit a pronounced leptokurtic distribution (much higher peaks and heavier tails than Gaussian). In addition, The Black-Scholes model fails to account for the strong dependence and the high volatility and intermittency of actual data, as discussed above.

[0036] In order to remedy the aforementioned shortcomings, the present invention employs a new (fractal activity time) model which replaces the Black-Scholes model by

[0037] X_t=log P_t−log P_t−1=μ+σ(W(T_t)−W(T_t−1))  (3)

[0038] where {T_t} is a (random) activity time process having a heavy tailed distribution, in order to avoid the above-described problem with assumption (A), and LRD in order to avoid the above-described problem with assumption (B). The above use of the term “heavy tailed distribution” means that {T_t−t} has, for each fixed t, a distribution whose tail probabilities decrease faster than exponentially, for example such that log P(|T_t−t|>x) behaves asymptotically like −α log x for large x, where α is a positive constant. The {T_t} process is assumed to be independent of W(t) and to have stationary differences, i.e., the mean value of (T_t−T_t−1) does not vary with respect to t, and the distribution of {T_t+s−T_t} depends only on s, not on t. In addition, the {T_t} process is standardized so that its mean value (“ET_t”) equals t. The “log” function in Equations (1), (2), and (3) can represent a logarithm having a base of e (i.e., a “natural logarithm”), or a logarithm having a different base (e.g., 10).

[0039] It is also to be noted that:

[0040] X_t=μ+σ(W(T_t)−W(T_t−1))μ+σ(T_t−T_t−1)^1/2W(1).  (4)

[0041] where “d” denotes equality of finite dimensional distributions, i.e., for any finite set (t₁,t₂, . . . ,t_k), the set (X_t,X_t, . . . ,X_t)has the same distribution as

[0042] (μ+σ(T_t−T_t−1)^1/2W(1), μ+σ(T_t−T_t−1)^1/2W(1), . . . , μ+σ(T_t−T_t−1)^1/2W(1)).

[0043] The heavy tails of the X_tcome from those of τ_t=T_t−T_t−1.

[0044] Also, for k=1, 2, . . . and centered variables,

[0045] cov(X_t,X_t+k)=0  (5a)

[0046] cov(|X_t|,|X_t+k)=σ²cov(τ_t^1/2,τ_t+k^1/2)  (5b)

[0047] cov(X_t²,X_t+k²)=3σ⁴cov(τ_t,τ_t+k)  (5c)

[0048] where Equation (5c) is defined only if Eτ_t²<∞. Thus, the LRD of {|X_t|} and {X_t²} follows from that of {τ_t^1/2} and {τ_t}, respectively. Other features of the model are conditional heteroscedasity,

[0049] var(X_t|F_t−1)=σ²E(τ_t|F_t−1),  (6)

[0050] and leptokurtosis

[0051] Kurtosis(X_t)=3(1+var τ_t)>3  (7)

[0052] if the kurtosis is finite. The kurtosis in Equation (7) is infinite if var τ_tis infinite (i.e. if Eτ_t²<∞).

[0053] The above-described model of the present invention typically works well for return periods of a day or more. Short term effects, such as leverage, which would otherwise preferably be factored in, can then be safely ignored. For shorter return periods, a leverage parameter λ can be used to adjust the model as follows:

[0054] X_t=μ+σ(W(T_t)−W(T_t−1))−λ(T_t−2T_t−1+T_t−2),  (8)

[0055] where the time steps of 1 no longer represent days, but some shorter time interval.

[0056] LRD is often associated with fractal or multifractal scaling, which suggests a relationship of the form

[0057] T_ct−ct(c)(T_t−t)  (9)

[0058] for positive c, where M and T are independent random functions. In Equation (9), “d” denotes equality of finite dimensional distributions, as discussed above. That is, for any positive c, and any finite set (t₁, t₂, . . . , t_k), (T_ct−ct₁, T_ct−ct₂, . . . , T_ct−ct_k) has the same distribution as M(c)(T_t−t₁, T_t−t₂, . . . , T_tt_k). Multifractal scaling means that the time axis of T_ct−ct can be stretched or contracted, by changing c, without changing the distribution of T_t−t. Self-similar (i.e., single fractal) behavior is a special case which occurs if M(c) is of the form c^κ, where 0<κ<1. These behavioral properties have been found to be quite common for risky asset time series. Accordingly, the fractal activity time model of the present invention provides greater accuracy where financial market data on the price of a risky asset exhibits fractal scaling and/or self similarity behavior,

[0059] The presence of multifractal or self-similar behavior can be tested by an algorithm in accordance with the invention, as described below.

[0060] The Algorithm

[0061] Suppose, for example, that the price P_tof a risky asset is observed at intervals of time each having a size of δ, where the observations take place over a fixed time period [0, τ]. In accordance with the invention, an algorithm can then be used to analyze empirical data such as {X_kδ, k=0, 1, . . . , N}, where X_kδ=log P_kδ−log P_(k−1)δ and τ=Nδ. An example of such an algorithm is described below, with reference to the flow diagram of FIG. 1. Each of the steps depicted in FIG. 1 can be carried out in a manner known to those skilled in the art, by interacting with a computer apparatus such as illustrated in FIGS. 7 and 8.

[0062] In Step 102 of the flow diagram, the algorithm initially checks for failure of the Black-Scholes model by analyzing the autocorrelation functions of {X²_kδ} and {|X_kδ|}. The autocorrelation functions are calculated from the empirical data X_kδ. For example, if the data {X²_kδ, k=0, 1, . . . , N} is put into a vector x, the acf(x,n) command of the “Splus” computer programming language can be used to plot the autocorrelation function of X²_kδ to lag n. The “acf” function can similarly be used to calculate the autocorrelation of |X_kδ|.

[0063] In Step 104 the autocorrelations calculated in Step 102 are tested as to determine whether they die away slowly. If the autocorrelations do not die away to zero for hundreds of observations, strong LRD is present in the data, and the algorithm therefore proceeds to Step 108. On the other hand, if the autocorrelations vanish relatively quickly, the algorithm is terminated by proceeding to Step 150. FIG. 2 illustrates an example of an autocorrelation function (ACF) of {|X_kδ} versus lag, derived from actual S&P 500 Index data, in which LRD is evident.

[0064] It is to be noted that the X_kδ are obtained from Equation (3) by replacing the time intervals of 1 by those of δ. That is,

[0065] X_kδ=log P_kδ−log P_(k−1)δ=μδ+σ(W(T_kδ)−W(T_(k−1)δ)).  (10)

[0066] Then, since E(W(T_t))=0 for all t, it follows that EX_kδ=μδ. Also since

[0067] W(T_kδ)−W(T_(k−1)δ)(T_kδ−T_(k−1)δ)^1/2W(1)  (11)

[0068] and T_kδ−T_(k−1)δd T_δ, it follows that

[0069] varX_kδ=σ²ET_δW(1)²=σ²δ,  (12)

[0070] since T_δ is independent of W(1) with ET_δ=δ and EW(1)²=1.

[0071] Since, as discussed above, for each k, EX_kδ=μδ and var X_kδ=σ²δ, it can be advantageous to standardize the array X_kδ so that it is unnecessary to further account for the unknown parameters μ and σ. Accordingly, in Step 108 of the flow diagram of FIG. 1, the algorithm generates a standardized array {Z_kδ, k=0, 1, . . . , N}, where:

[0072] Z_kδ=(X_kδ−{circumflex over (μ)})/{circumflex over (σ)},  (13)

[0073] and

[0074] $\begin{array}{cc}\hat{\mu }=\sum _{k=0}^N{X}_{k\delta }/\left(N+1\right),{\hat{\sigma }}^2=\sum _{k=0}^N{\left(X_{k\delta }-\hat{\mu }\right)}^2/\left(N+1\right).& \left(14\right)\end{array}$

[0075] {circumflex over (μ)} is an estimate of the mean of X_kδ, and {circumflex over (σ)}²is an estimate of the variance of X_kδ, where both are calculated from the empirical data.

[0076] It is then possible to produce an approximation of the activity time process {T_t}, based on the empirical data. The approximation is based upon Ito's formula, which is well-known in the field of stochastic calculus. See, e.g., R. J. Elliott, Stochastic Calculus and Applications,Springer, New York (1982). Under Ito's formula: $\begin{array}{cc}{\sigma }^2{T}_t={\int }_0^t{(d{\log P}_S)}^2,& \left(15\right)\end{array}$

[0077] the right hand side being called the quadratic variation. This is a standard result in stochastic calculus. The above equation can be put in discrete (as opposed to continuous) form by replacing the ${\int }_0^t{(d{\log P}_S)}^2\mathrm{with}\sum _{k-0}^N{\left({\log P}_{k\delta }-\log P_{\left(k-1\right)\delta }-\mathrm{\mu \delta }\right)}^2,$

[0078] which equals $\sum _{k-0}^N{\left(X_{k\delta }-\mathrm{\mu \delta }\right)}^2.$

[0079] T_tis replaced by T_nδ to yield the discrete approximation: $\begin{array}{cc}{\sigma }^2{T}_{n\delta }\simeq \sum _{k=0}^N{\left(X_{k\delta }-\mathrm{\mu \delta }\right)}^2.& \left(16\right)\end{array}$

[0080] This equation can be put in standardized form such that: $\begin{array}{cc}{T}_{n\delta }\simeq \sum _{k=0}^N{Z}_{k\delta }^2.& \left(17\right)\end{array}$

[0081] As the interval δ approaches zero, the discrete approximation becomes exact.

[0082] To produce the approximation of {T_t}, the algorithm, in Step 110 of FIG. 1, generates an array {{circumflex over (T)}_nδ,n=1,2, . . . ,N} where ${\hat{T}}_{n\delta }=\sum _{k-1}^n{Z}_{k\delta }^2.$

[0083] As discussed above, the empirical quadratic variation process {{circumflex over (T)}_nδ} serves as an approximation of the {T_nδ} process. Convergence should take place under broad conditions as δ decreases. For illustration, FIG. 3ashows an example of a relationship between {circumflex over (T)}_tand t for a set of actual data, and FIG. 3bshows a corresponding plot of {{circumflex over (T)}_t−t} versus t.

[0084] It is noted that the empirical quadratic variation process {{circumflex over (T)}_nδ} approximates the (unobserved) activity time process {T_nδ} in the following sense. As the time interval δ between observations decreases, the {circumflex over (T)} process converges to the T process with probability 1. That is, max_0<t<τ|{circumflex over (T)}_t−T_t|→0 with probability one for any fixed τ>0. The approximation is satisfactory even for daily data as only broad features, like averages of the {circumflex over (T)} are used, and not fine detail.

[0085] Next, using E|{circumflex over (T)}_δ|^pas an approximation of E|T_δ|^p, the algorithm, in Step 112 of FIG. 1, determines the range of p values for which E|T_δ|^p<∞, or equivalently, E|X_δ|^2p<∞, where {circumflex over (T)}_δ denotes T_nδ|n=1, T_δdenotes T_nδ|n=1, and X_δ denotes X_nδ|n=1. In other words, the algorithm determines the range of p values for which the p-th moment of {T_nδ} and the 2p-th moment of {X_nδ} are finite. Since the process has stationary differences, T_δ has the same distrbution as T_kδ−T_(k−1)δ for all k. Similarly for Xδl . Empirical studies suggest that the maximum value of p is ordinarily less than 3 and sometimes less than 2, but generally greater than 1.5. For example, if the data Z_kδ is put into a vector Z, a routine written in Splus can calculate mean((abs(Z))^r) for a range of r values (e.g., at intervals of 0.25), where r=2p, until the resulting function begins to increase steeply, indicating that the limit of the range of convergence of both |Z_δ|^pand |X_δ|^2pis being approached. The results of such a procedure are illustrated in FIG. 4, which shows a plot of “moment” E|Z_δ|^rversus “index” r, based on S&P 500 Index data. The plot shows the range of values of p for which convergence takes place for an empirical quadratic variation process calculated from the aforementioned exemplary set of actual index data. The plot suggests that for r<4, the value of the moment E|Z_δ|^ris finite, meaning that the moment converges. However, for r≧4, the moment becomes extremely large, suggesting that it approaches ∞ for such larger values of r.

[0086] Next, in Step 114 of the flow diagram of FIG. 1, the algorithm computes a set of values of the function $\begin{array}{cc}{\hat{S}}_q\left(m\delta \right)=\sum _{k=1}^{\left[N/m\right]}{\left|{\hat{T}}_{\mathrm{km}\delta }-{\hat{T}}_{\left(k-1\right)m\delta }-m\delta \right|}^q/N,& \left(18\right)\end{array}$

[0087] where m is an integer and “[N/m]” denotes the integer part of N/m. In Step 114, log Ŝ_q(mδ) is computed for a wide range of q and m values. Values of q are chosen such that 0<q≦p, where p has been determined from Step 112 to ensure that E|T_δ|^p<∞. For example, values of q can be selected at intervals of 0.25. The selected values of m can be, for example, 2, 2², . . . , 2¹⁰.

[0088] If the random process being analyzed obeys Equation (9), and if $\begin{array}{cc}{S}_q\left(m\delta \right)=\sum _{k=1}^N|T_{\mathrm{km}\delta }-T_{\left(k-1\right)m\delta }-m\delta {|}^q/N& \left(19\right)\end{array}$

[0089] then

[0090] ES_q(mδ)=E|T_mδ−mδ|^q=E|M(m)|^qE|T_δ−δ|^q.  (20)

[0091] Furthermore, Equation (9) implies that

[0092] log E|M(m)|^q=log E|M(e)|^qlog m.  (21)

[0093] In addition, the ergodic theorem widely used in probability states that if {W_k, k=1, 2, . . . } is a stationary process, and E|W₁|<∞, then $\begin{array}{cc}{n}^{-1}\sum _{k=1}^n{W}_k\to {\mathrm{EW}}_1& \left(22\right)\end{array}$

[0094] almost surely (hereinafter “a.s.”, which indicates an outcome with probability=1) as n→∞.

[0095] Moreover, for large enough sample sizes, log Ŝ_q(mδ) approaches, and can therefore be used to approximate, log ES_q(mδ).

[0096] Therefore, using the ergodic theorem, $\begin{array}{cc}\log {\hat{S}}_q\left(m\delta \right)\stackrel{a.s.}{\to }\log {\mathrm{ES}}_q\left(m\delta \right)=\log E{\left|M\left(e\right)\right|}^q\log m+\log E{\left|T_{\delta }-\delta \right|}^q& \left(23\right)\end{array}$

[0097] which is linear with respect to logm as N→∞. Accordingly, the algorithm, in Step 116, checks for linearity of log S_q(mδ) with respect to log m using the previously calculated values of Ŝ_q(mδ) as an approximation for S_q(mδ). Linearity can often be checked by visual observation on a computer display or a hard copy of a graphical plot of log S_q(mδ) versus log m. If linearity is found in Step 118, this indicates that a multifractal model can properly be used to describe {T_t−t} and the algorithm proceeds to Step 120. Multifractal behavior means that scaling or magnifying the time axis of the {T_t−t} function will produce a function with an identically-shaped probability density distribution. If linearity is not found in Step 118, the algorithm proceeds to Step 150 and terminates.

[0098] A self-similar (i.e., single fractal) model can be used in the special case where, in Equation (9), M(c)≡c^κ for all c and some κ (i.e., where 0<κ<1). In such a case,

[0099] log E|M(m)|^q=qκ log m  (24)

[0100] and Equation (23) becomes $\begin{array}{cc}\log {\hat{S}}_q\left(m\delta \right)\stackrel{a.s.}{\to }\log {\mathrm{ES}}_q\left(m\delta \right)=\mathrm{qK}\log m+\log E{\left|T_{\delta }-\delta \right|}^q& \left(25\right)\end{array}$

[0101] Accordingly, the algorithm, in Step 122, checks log S_q(mδ)−log E|T_δ−δ|^qfor linearity with respect to q, again using the previously calculated values of Ŝ_q(mδ) as approximations for S_q(mδ). In Step 122a log E|T_δ−δ|^qis estimated by $\log \left(\sum _{k-1}^N{\left|Z_{k\delta }^2-\delta \right|}^q/N\right).$

[0102] Linearity can again often be checked by eye using a computer display or a hard copy of the graphical representation of the relation set forth in Step 112. If linearity is found in the test of Step 124, then a self-similar model satisfactorily represents the process {T_t−t}, and the algorithm proceeds to Step 126. If linearity is not found in this test of Step 124, the algorithm proceeds to Step 150 and terminates.

[0103] An exemplary interactive computing procedure for carrying out the exemplary algorithm represented by the flow diagram 100 in FIG. 1 is provided in the Appendix hereto. The computing steps of the procedure are expressed in the Splus programming language, but other programming languages can also be used.

[0104] FIGS. 5a-fand 6a-ishow results generated by an algorithm in accordance with the invention, using as input 20 years of daily data from the S & P 500 Index. To generate these results, the S&P 500 Index data, represented by the vector {P_kδ, k=0, 1, . . . , N}, were used to calculate the vector X_kδ=log P_kδ−log P_(k−1)δ. The vector {X_kδ} was then used to generate the vector Z_kδ=(X_kδ−{circumflex over (μ)})/{circumflex over (σ)} as discussed above with respect to the exemplary algorithm of the invention. The vector {Z_kδ} was then used to generate the vector {{circumflex over (T)}_hδ,n=1,2, . . . ,N}, where ${\hat{T}}_{n\delta }=\sum _{k-1}^n{Z}_{\mathrm{k\delta }}^2.$

[0105] The vector {{circumflex over (T)}_nδ} was then used to calculate the value of the function ${\hat{S}}_q\left(m\delta \right)=\sum _{k=1}^{\left[N/M\right]}{\left|{\hat{T}}_{\mathrm{km}\delta }-{\hat{T}}_{\left(k-1\right)m\delta }-m\delta \right|}^q/N$

[0106] for various values of q and m. The resulting values of the function Ŝ_q(mδ) were used to approximate S_q(mδ) in order to calculate log S_q(mδ). In addition, the vector {{circumflex over (T)}_nδ} was used to approximate {T_nδ} in order to calculate log S_q(mδ)−log E|T_δ−δ|^q. The resulting values of log S_q(mδ) and log S_q(mδ)−log E|T_δ−δ|^qwere used to produce the plots in FIGS. 5a-fand 6a-i.

[0107] The results are provided in order to illustrate the effectiveness of the algorithm. In FIGS. 5a-f,log S_q(mδ) is plotted on the horizontal axis, and log m is plotted on the vertical axis. In FIGS. 6a-i,q is plotted on the horizontal axis and log S_q(mδ)−log E|T_δ−δ|^qis plotted on the vertical axis. The data in the figures, which clearly illustrate the linearity of log S_q(mδ) with respect to log m (FIGS. 5a-f), and the linearity of log S_q(mδ)−log E|T_δ−δ|^qwith respect to q (FIGS. 6a-i), demonstrate that a self-similar model with κ=0.8 for {T_t−t} provides a good representation of an exemplary random process in real life. The linearity of the above-described plots indicates that the value of the S&P 500 Index can be modeled, in accordance with the invention, by the function

[0108] P_t=P_oexp[μt+σW(T_t)]  (26)

[0109] where P_ois the value of the Index at t=0, μ and σ>0 are fixed constants, W(t) is a standard Brownian motion function, and T_tis a random activity time process which is approximated by the vector {{circumflex over (T)}_nδ} that is generated as described above. In other words, the S&P 500 index can be modeled by the discrete function

[0110] P_nδ=P_oexp[μnδ+σW({circumflex over (T)}_nδ)]  (27)

[0111] It has been found from other similar investigations that good self-similar fits can be obtained for a wide range of time series data. The computer implemented steps of an algorithm in accordance with the invention can, optionally, be implemented in Splus, C, C+, C++, BASIC, FORTRAN, or any other suitable computer programming language. The models in accordance with the invention have characteristics which differ substantially from those of the Black-Scholes model. They can provide significant improvement in the accuracy of assessed risk calculations involving financial data.

[0112] It is to be noted that the increments of the process {T_t−t} can easily be shown to be LRD (SAV) for κ>1/p if E|T₁−1|^p<∞ for some p, where 0<p<2. C. C. Heyde & Y. Yang, On defining long-range dependance,34 J. Applied Prob. 939 (1997). The notation “LRD (SAV)” refers to a sample Allen variance based definition of LRD which copes with processes having possible infinite variance, such as is the case here, but which essentially reduces to the classical definition in the finite variance case. Such characteristics readily explain the empirically observed strong dependence.

[0113] The present invention can, for example, be used to valuate a European Call Option for a traded security where the time parameter data generated from the price data of the security over time exhibits self similarity, such that the security price can be modeled as a random process in accordance with the present invention so that a model of the price of the security in accordance with the present invention affords a higher degree of accuracy than the prior art. The market price of the security at time t is denoted by P_t. A European Call Option for such a security gives a holder of the option the right, but not the obligation, to buy the security when it matures at time t₁for the exercise price K. The value of such an option C(t,t₁,K) at time t<t₁may be calculated under the fractal activity time model assuming a fixed interest rate r for bonds.

[0114] As an initial step in valuating the European Call Option, it is necessary to generate an array {{circumflex over (T)}_nδ, n=1,2, . . . ,N} by using an algorithm in accordance with the invention to process an array of price data {P_kδ} representing the history of the price P_tof the underlying security. For example, an algorithm such as illustrated by the flow diagram 100 of FIG. 1 can be used to generate {{circumflex over (T)}_nδ}, which can be used to approximate the fractal activity time process {T_t} in the subsequent calculations. In addition, the algorithm calculates the fractal scaling parameter κ of Equation (24).

[0115] It is to be noted that the activity time process {T_t} is independent of the {W} process in the model. Then, conditionally on a fixed outcome for the {T_t} process, the fractal activity time model can be written in geometric Brownian motion form as

[0116] log P_t=μt+σW(T_t)  (28)

[0117] which differs from the Black-Scholes formulation in that T_tis a function of t. We initially assume that T_tis non-random with respect to t.

[0118] The standard Black-Scholes formula calculation gives the result in the risk-neutral (i.e., arbitrage is not allowed) valuation form that

[0119] E(C(t,t₁,K))=e^−r(t−t){tilde over (E)}((P_t−K)⁺|  (29)

[0120] where {tilde over (E)} is the expectation with respect to a new measure for the discounted price, x⁺=max(0,x), where x=P_t−K, and the conditioning is on the past history _tof the security price up to time t. This calculation applies unchanged to the fractal activity time model conditioned on fixed (i.e., non-random, with known history) {T_t}, which gives

[0121] E(C(t,t₁,K)|{T})=e^−r(t−t){tilde over (E)}((P_t−K)⁺|_t)  (30)

[0122] However, under this new measure for the fractal time model, still conditioned on fixed {T_t}, there is W* (t) given by $\begin{array}{cc}\sigma W^{*}\left(T_t\right)=\sigma W\left(T_t\right)+\left(\mu -r\right)t-\frac{1}{2}{\sigma }^2{T}_t,& \left(31\right)\end{array}$

[0123] which is standard Brownian motion, so that $\begin{array}{cc}\log P_{t_1}=\log P_t+r\left(t_1-t\right)+\frac{1}{2}{\sigma }^2\left(T_{t_1}-T_t\right)+\sigma \left(W^{*}\left(T_{t_1}\right)-W^{*}\left(T_t\right)\right)& \left(32\right)\end{array}$

[0124] and W* (T_t)−W* (T_t) is independent of _t, with the distribution of W* (T_t−T_t). Then,

[0125] e^−r(t−t){tilde over (E)}(P_t−K)⁺|_t)=e^−r(t−t)E((P_te^Z−K)⁺|P_t)  (33)

[0126] where Z has the normal distribution with mean r(t₁−t)+½σ²(T_t−T_t) and variance σ²(T_t−T_t). Then rewriting equation (33) explicitly in terms of the normal distribution for Z, there is now obtained the expression $\begin{array}{cc}E\left({\left(P_t{e}^Z-K\right)}^{+}|P_t\right)=\frac{1}{\sigma \sqrt{2\left(T_{t_1}-T_t\right)\pi }}{\int }_{z>-\log \left(P_t/K\right)}\left(P_t{e}^z-K\right)\exp \left\{-\frac{{\left(z-r\left(t_1-t\right)-\frac{1}{2}{\sigma }^2\left(T_{t_1}-T_t\right)\right)}^2}{2{\sigma }^2\left(T_{t_1}-T_t\right)}\right\}dz,& \left(34\right)\end{array}$

[0127] which transforms to give, still conditioned on fixed (non-random) {T},

[0128] E(C(t,t₁,K)|{T})=P_t[Φ(d(t₁,t))]−Kexp(r(t−t₁))[Φ(d(t−t₁)−σ(T_t−T_t)^1/2],  (35)

[0129] where $\begin{array}{cc}d\left(t_1,t\right)=\frac{\log \left(P_t/K\right)+r\left(t_1-t\right)+\frac{1}{2}{\sigma }^2\left(T_{t_1}-T_t\right)}{{\sigma \left(T_{t_1}-T_t\right)}^{\frac{1}{2}}}& \left(36\right)\end{array}$

[0130] and $\Phi \left(x\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^x{e}^{-\frac{1}{2}{u}^2}du$

[0131] is the distribution function of the standard normal. Thus, unconditionally (i.e., regardless of whether T_t, is a random or non-random function of t),

[0132] C(t,t₁,K)=P_tE[Φ(d(t₁,t))]−Kexp(r(t−t₁))E[Φ(d(t₁,t)−σ(T_t−T_t)^1/2]  (37)

[0133] with the expectation “E” being taken over the distribution T_t−T_t.

[0134] In the special case of the geometric Brownian motion model and the Black-Scholes formula, T_t−T_t=t₁−t with probability one, where t₁−t is constant, and the expectation disappears. However, in general T_t−T_tis random.

[0135] In the case of the fractal activity time model of the present invention, since {T_t} has stationary differences, it is known that

[0136] T_t−T_tT_t−tt₁−t+(t₁−t)^κ(T₁−1)  (38)

[0137] where d denotes equality in distribution and κ is the fractal scaling parameter.

[0138] For calculation of option valuations an advantageous choice for a distribution for T_tis a scaled generalized inverse Gaussian distribution whose probability density is of the form $\begin{array}{cc}{g}_v\left(x\right)=\frac{{\left(\left(v-2\right)/2\right)}^{v/2}}{\Gamma \left(v/2\right)}{x}^{-v/2-1}{e}^{-\left(v-2\right)/2x},x>0,v>2.& \left(39\right)\end{array}$

[0139] where Γ is the gamma function $\Gamma \left(y\right)={\int }_0^{\infty }{u}^{y-1}{e}^{-u}du,$

[0140] v is a fixed parameter to be estimated from the data, and x is a dummy variable defining the function over the interval (0, ∞).

[0141] This choice of distribution gives rescaled t distributions for the one dimensional distribution of the log price which has been shown to fit the data in various empirical studies. (See, e.g., S. R. Hurst, E. Platen, and S. T. Rachev, “Subordinated Market Index Models: A Comparison,” Financial Engineering and the Japanese Markets 4 (1997), pp. 97-104. The parameter v is fitted from data in any particular case. If {X_t=log P_t−log P_t−1}, then the X_thas a rescaled t distribution with v degrees of freedom. If N_x=number of |X_t|>x is regressed on logx then the slope estimates −v. For the model of the S&P 500 Index data of Example 1, the estimated value of v is 3.4. An exemplary distribution for T₁is the inverse Gaussian given in Equation (39). The probability densities of (T₁−1) with v=3.4 are illustrated in the graph of FIG. 3c.For comparison, FIG. 3calso includes a plot of the probability densities of a standard Gaussian distribution.

[0142] Calculations of the value C(t, t₁, K) are most conveniently done by expressing it in terms of the “moneyness” x defined by

[0143] P_te^r(t−t)=Ke^x.  (40)

[0144] Then, using (37) and (38), $\begin{array}{cc}C\left(t,t_1,K\right)=P_t\left[E\right(\Phi \left(\frac{x+\frac{1}{2}{\sigma }^2{T}_{t_1-t}}{\sigma T_{t_1-t}^{1/2}}\right)-e^{-x}E\left(\Phi \left(\frac{x-\frac{1}{2}{\sigma }^2{T}_{t_1-t}}{\sigma T_{t_1-t}^{1/2}}\right)\right].& \left(41\right)\end{array}$

[0145] Also, again using (38), T_t−thas the same distribution as

[0146] t₁−t+(t₁−t)^κ(T₁−1)=t₁−t−(t₁−t)^κ+(T₁−t)^κT₁.  (42)

[0147] Now, if g^v(x) is the probability density function (39) for T₁, then $\begin{array}{cc}E\left(\Phi \left(\frac{x\pm \frac{1}{2}{\sigma }^2{T}_{t_1-t}}{\sigma T_{t_1-t}^{1/2}}\right)\right)={\int }_0^{\infty }\Phi \left(\frac{x\pm \frac{1}{2}{\sigma }^2\left(t_1-t-{\left(t_1-t\right)}^{\kappa }+{\left(t_1-t\right)}^{\kappa }x\right)}{{\sigma \left(t_1-t-{\left(t_1-t\right)}^{\kappa }+{\left(t_1-t\right)}^{\kappa }x\right)}^{1/2}}\right)g_v\left(x\right)dx,& \left(43\right)\end{array}$

[0148] where κ is the fractal scaling parameter of Equations (24) and (25). The factal scaling parameter κ is calculated from the {T_nδ} data using, for example, the algorithm of FIG. 1. This integration can be done in any particular case using the Splus programming language.

[0149] To give a numerical illustration, the case of the S&P 500 Index data is used, for which the estimated values of the parameters are J=0.01, κ=0.8, v=3.4. It is also assumed that t₁−t=90 (i.e. a 3 month time interval). Then, in the fundamental case of x=0 (called “at the moneyness”) $\begin{array}{cc}C\left(t,t+90,K\right)\left(0\right)=P_t\left(2E\left(\Phi \left(\frac{1}{2}\sigma T_{90}^{1/2}\right)\right)-1\right).& \left(44\right)\end{array}$

[0150] This gives (0.0362) P_t, in contrast to the Black-Scholes valuation, $\begin{array}{cc}{P}_t\left(2\Phi \left(\frac{1}{2}\sigma \sqrt{90}\right)-1\right)=\left(0.0378\right))P_t,& \left(45\right)\end{array}$

[0151] resulting in roughly a 5% difference. Accordingly, if the analysts driving the market price of the option employ the conventional Black-Scholes model, then the new method of the invention will provide the user of the new method with the important information that the option is overvalued by 5% and should therefore be sold, not bought, at the market price.

[0152] The present invention can also be applied, for example, to the valuation of a European put option for a traded financial security where the time parameter data generated from the data of the security over time exhibits self similarity. A European put option gives the holder the right, but not the obligation, to sell the security for an exercise price K when the option matures at time t₁. If arbitrage is not permitted and there is a fixed interest rate r for bonds, the price P(t,t₁,K) of such a European put option may be obtained from the relation

[0153] P_t+P(t,t₁,K)−C(t,t₁,K)=Ke^−r(t−t)  (46)

[0154] where C(t,t₁,K) is the corresponding European call option value for the same security calculated in accordance with Example 2 hereinabove, and P_tis the price of the underlying traded asset at time t.

[0155] The present invention can also be applied, for example, to valuate an American call option for a traded security where the time parameter data generated from the price data of the security over time exhibits self similarity, such that the security price can be modeled as a random process in accordance with the present invention. An American option gives the holder the right, but not the obligation, to purchase the security at an exercise price K at any time t before an expiration time t₁. It is well known that it is suboptimal to exercise an American call option before its expiration time t₁. Therefore, the value of an American call option is the same as that of the corresponding European call option described in Example 2 hereinabove.

[0156] The invention can also be applied, for example, to valuate an American put option for a traded security where the time parameter data generated from the price data of the security over time exhibits self similarity, such that the security price can be modeled as a random process in accordance with the present invention. An American put option gives the holder the right to sell the security for an exercise price K at any time t before an expiration time t₁. No closed form expression is available for the value of an American put option in the case of the geometric Brownian motion model. Results must be obtained by numerical methods. A simple adaptation of the binomial tree method may be used to obtain approximate valuations using the fractal activity time model of the present invention.

[0157] The present invention is likewise applicable, for example, to the valuation of so-called exotic options for traded securities where the time parameter data generated from the price data of the security over time exhibits self similarity, such that the security price can be modeled as a random process in accordance with the present invention. Various exotic options are easily valued using risk-neutral valuation principles and the European call option valuation within the classical Black-Scholes framework. Examples of exotic options are collar, break forward, range forward, forward start, chooser, digital, barrier and lookback options, as described, for example, in Chapter 9 of the book by M. Musiela and M. Rutwowski, entitled Martingale Methods in Financial Modelling,Springer, Berlin, 1998. The valuations described are easily replaced by new versions calculated under the fractal activity time model of the present invention.

[0158] Other exotic options such as Asian and basket options require numerical methods for their calculation in the case of geometric Brownian motion model. Valuations can be done under the fractal activity time model of the present invention using Monte Carlo methods.

[0159] Computer Implementation of Methods According to the Invention

[0160] It will be appreciated by those skilled in the art that the exemplary algorithm represented by the flow diagram of FIG. 1 can be implemented through interactive computing on various standard computer hardware platforms operating under the control of suitable software implementing the interactive computing of the algorithm defined by the flow diagram of FIG. 1. In some cases, dedicated computer hardware, such as a peripheral card which resides on the bus of a standard personal computer, may enhance the operational efficiency of the above methods.

[0161] FIGS. 7 and 8 illustrate typical computer hardware suitable for practicing the present invention. The exemplary computer system of FIG. 7 includes a computer section 710, a display 720, a keyboard 730, and a communications peripheral device 740, such as a modem. The system may also include other input devices such as a CD ROM drive 750 for reading a CD ROM 700 and a printer 760. The computer system generally includes one or more disk drives 770 which can read and write to computer-readable media, such as magnetic media (i.e., diskettes or hard disks) for storing data and software. The data and software can include, e.g., instructions to direct the computer system to perform the steps of the invention. While not shown, other input devices, such as a digital pointer (e.g., a “mouse”) and the like may also be included.

[0162] FIG. 8 is a functional block diagram which further illustrates the computer section 710. The computer section 710 typically includes a processing unit 810, control logic 820 and a memory unit 830. Preferably, computer section 710 can also include a timer 850 and input/output ports 840. The computer section 710 can also include a co-processor 860, depending on the microprocessor used in the processing unit. Control logic 820 provides, in conjunction with processing unit 810, the control necessary to handle communications between memory unit 830 and input/output ports 840. Timer 850 provides a timing reference signal for processing unit 810 and control logic 820. Co-processor 860 provides an enhanced ability to quickly perform numerous and/or complex computations, such as those involving large data sets.

[0163] Memory unit 830 may include different types of memory, such as volatile and non-volatile memory and read-only and programmable read-only memory. For example, as shown in FIG. 8, memory unit 830 may include read-only memory (ROM) 831, electrically erasable programmable read-only memory (EEPROM) 832, and random-access memory (RAM) 833. Different computer processors, memory configurations, data structures and the like can be used to practice the present invention, and the invention is not limited to a specific hardware or software platform.

[0164] An example of a programming language which is well-suited for practicing the invention is Splus, version 3.4 or 5.0, which can run on a Unix operating system. Splus was originally developed by A. T. & T. Bell Laboratories, and is currently marketed by MathSoft. An Splus program can, for example, be run on a Sun Sparkstation or an equivalent computer.

[0165] Although the present invention has been described by way of specific exemplary embodiments, it should be understood that various changes, substitutions and alterations to the disclosed embodiments will be apparent to those skilled in the art, and can be made without departing from the spirit and scope of the invention as set forth in the appended claims.

[0166] Note that this exemplary interactive computing procedure utilizes human judgment at various steps. It is preferable to use at least several thousand observations of the random process being modeled. For financial processes the return data is typically daily, so that for this example δ=1. An exemplary procedure is provided as follows:

[0167] Return data {P_k} is placed into a vector u.

[0168] Calculate v which is log u (gives vector of the log P_k).

[0169] Calculate x which is diff(v, lag=1) (gives vector of the X_t=log P_k−log P_(k−1)).

[0170] acf(x,n) generates autocovariances of elements in vector x out to lag n and plots these. This is done for the vectors x, abs(x) and various n values, say n=50 or n=100. The autocovariances for x should die out very quickly, while those for abs(x) can be expected to be significantly different from zero even for n=50 or more.

[0171] Generate a standardized array z which is ((x)−mean(x))/(var(x))^1/2.

[0172] Calculate mean((abs(z))^r) for a range of r values, say at intervals of 0.25 until the function begins to increase steeply. (This indicates where E|X₁|^rceases to be finite.)

[0173] Calculate w which is cumsum (z²) (Generates the empirical activity time process {circumflex over (T)} which can be plotted using ts.plot(w).)

[0174] Calculate the Ŝ_q(m). The vector of {circumflex over (T)}_k−{circumflex over (T)}_(k−m)is given by um which is diff(w, lag=m). Then take vm which is abs((um)−m) and Ŝ_q(m) is given by mean((vm)^q). For m values, use, for example, m=5.2^j, j=0, 1, . . . , 7.

[0175] Plot log Ŝ_q(m) against log m using the plot (x, y) function in Splus for the range of q values for which T, has finite absolute moments as determined earlier. Check the plots for linearity of the plots.

[0176] Estimate the E|T₁−1|^qusing mean(abs(u1−1)^q). Plot log Ŝ_q(m)−log E|{circumflex over (T)}₁−1|^qfor the q values used in the previous step. Check for linearity of the plots and constancy of estimated {circumflex over (K)}.

[0177] A summary of an exemplary code sequence for implementing the computing portions of the above procedure is provided as follows:

[0178] v<−log(u)

[0179] x<−diff(v, lag=1)

[0180] acf(x, 200)

[0181] z<−((x)−mean(x))/(var(x)){circumflex over ( )}(1/2)

[0182] mean((abs(z){circumflex over ( )}r)

[0183] w<−cumsum(z{circumflex over ( )}2)

[0184] um<−diff(w, lag=m)

[0185] vm<−abs((um)−m)

[0186] log(mean(vm){circumflex over ( )}q)

[0187] mean((abs((u1)−1)){circumflex over ( )}q)

[0188] plot(x, y) (various values of x, y)