Next Article in Journal
Deflecting an Asteroid on a Collision Course with Earth Using a Powered Swing-By Maneuver
Next Article in Special Issue
Long-Wave Anti-Plane Motion in a Pre-Stressed Compressible Elastic Laminate with One Fixed and One Free Face
Previous Article in Journal
Circuit Complexity from Supersymmetric Quantum Field Theory with Morse Function
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Fractional Brownian Motion: Estimating the Hurst Exponent via Variational Smoothing with Applications in Finance

1
Department of Computer Science, University of Verona, 37129 Verona, Italy
2
SMC Treviso, 31020 Villorba, Italy
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(8), 1657; https://doi.org/10.3390/sym14081657
Submission received: 2 July 2022 / Revised: 5 August 2022 / Accepted: 8 August 2022 / Published: 11 August 2022
(This article belongs to the Special Issue Asymmetric and Symmetric Study with PDE)

Abstract

:
Beginning with the basics of the Wiener process, we consider limitations characterizing the “Brownian approach” in analyzing real phenomena. This leads us to first consider the fractional Brownian motion (fBm)—also discussing the Wood–Chan fast algorithm to generate sample paths—to then focus on multi-fBm and methods to generate its trajectories. This is heavily linked to the Hurst exponent study, which we link to real data, firstly considering an absolute moment method, allowing us to obtain raw estimates, to then consider variational calculus approaches allowing to smooth it. The latter smoothing tool was tested in accuracy on synthetic data, comparing it with the exponential moving average method. Previous analyses and results were exploited to develop a forecasting procedure applied to the real data of foreign exchange rates from the Forex market.

1. Introduction

The properties of the Wiener process and the Gaussian distribution (appearing naturally in many contexts as a limiting distribution of the sum of random components, and the fact that it is relatively easy to numerically solve stochastic differential equations (SDEs) driven by the Brownian motion) are relevant elements in explaining why this process has been so heavily used to model several phenomena. Between the plethora of possible applications, great interest has been devoted to mathematical finance problems, where the renowned Black and Scholes model can be considered as the trigger of the modern analytical approach to finance. However, in recent years, the availability of large data sets has led to the identification properties that are common across a wide range of instruments, markets, and time periods, observed by independent studies and classified as “stylized facts” [1], offering empirical evidence that the geometric Brownian motion may not be a good candidate to describe the behaviors of financial assets [2]. Let us underline the following three: economic time series often exhibit cycles of all orders of magnitude, suggesting a strong interdependence between distant samples; the distribution of asset returns tends to be non-Gaussian, sharp-peaked, and heavy-tailed; the risky character of an asset is directly related to the smoothness of the observed price trajectories, so the local regularity of the sample paths generated by a model should attempt to reproduce their variations.
The fractional Brownian motion (fBm) can be considered as one of (if not the) first generalizations of the Wiener process. It incorporates long-range dependence and the possibility to tune the local regularity of its sample paths, both features being controlled with a constant parameter H ( 0 , 1 ) called the Hurst exponent. In the general case, the fBm is not a semimartingale. This fact entails two problems, which discouraged financial applications in their early years: the definition of proper stochastic calculus tools, based on it, and the fact that this model allows arbitrage opportunities. The former problem can be solved by exploiting Malliavin derivatives as well as the white noise theory [3]. Concerning the latter problem, since fBm is not a semimartingale, there can be no equivalent martingale measure for the price process; hence, the first fundamental theorem of asset pricing implies that there must be arbitrage. Rogers [4] and Cheridito [5] constructed examples where such situations can be seen clearly and proposed methods to exclude the arbitrage, such as restricting the class of admissible strategies or forbidding high-frequency trading (HFT), even if the HFT arbitrage is limited by the stock servers set up [6]; see also [7] for applications to option pricing.
The fBm was in turn generalized by the multi-fractional Brownian motion (multi-fBm or mBm for short), allowing to control the point-wise regularity of the paths of the process, i.e., the Hurst exponent becomes a function of time H = H ( t ) . This feature makes the process very flexible and capable of a more accurate description of reality. The mBm is able to reproduce many stylized empirical facts observed in the financial time series [1]. The mBm is able to harmonize the different approaches to the study of financial dynamics: efficient market theories, agent-based models, and behavioral finance. Many works in the literature show that the estimates for H ( t ) fluctuate around 0.5 [8,9], which is the only value compatible with the no-arbitrage condition. Moreover, the path traced by H ( t ) can capture the different aspects of the investor’s trading behavior: it can detect the trends ( H > 0.5 ), typical of “bull” and “bear” markets resulting from a gradual increase of confidence, and the mean-reversion ( H < 0.5 ), typical of the frenetic buy-and-sell activity during financial crises.

Contributions

The main contribution provided in the present paper consists of exploiting the aforementioned fBm and mBm stochastic processes within the financial sector, particularly using the smoothing method (see below) w.r.t. real Forex market data. We found a precise threshold to abide by for the smoothing method to be more accurate than the classical exponential smoothing, as well as to obtain accurate one-step forecasts simulating Forex data. In particular, the work is structured as follows: in the first two sections, we recall the definitions and properties of fBm and mBm, and the Wood–Chan algorithm to generate sample paths. Section 4 is devoted to the estimation of the Hurst exponent. First, we recall a method based on absolute moments to obtain a raw estimate, then we recall a technique based on variational calculus to smooth it. Section 5 is focused on testing the accuracy of the smoothing method on synthetic data and comparing it with the exponential moving average; finally, in Section 6, we apply all of the learned tools to real data from the Forex market to forecast the foreign exchange rates in single-step and multi-step frameworks. Let us underline that the ’goodness’ of the latter methods has also been highlighted by recent studies on multi-fractional processes in finance, see, e.g., [10,11,12].

2. Definitions and Main Properties

2.1. fBm

The standard fBm with the Hurst exponent H ( 0 , 1 ) on a probability space ( Ω , F , P ) is a stochastic process B H = { B t H } t R , such that
  • B 0 H = 0 .
  • It is a Gaussian process with B t H N ( 0 , | t | 2 H ) .
  • It has a covariance function
    E ( B t H B s H ) = 1 2 ( | t | 2 H + | s | 2 H | t s | 2 H ) .
Depending on the value of H, we have:
  • if H = 1 / 2 , then B H is the classical Brownian motion.
  • if H < 1 / 2 , the increments of B H are negatively correlated and the process exhibits short-range dependence, i.e., it is more likely that the trend is broken.
  • if H > 1 / 2 , the increments of B H are positively correlated and the process exhibits long-range dependence, i.e., it is more likely that the trend is preserved.
For the sake of completeness, let us briefly recall some properties of the fBm that will be useful to simulate its trajectories.
Denote with Y = { Y t } t R = { B t + 1 H B t H } t R the increment process, called fractional Gaussian noise (fGn). The covariance between two increments Y k and Y k + n is
ρ H ( n ) = 1 2 [ ( n + 1 ) 2 H + ( n 1 ) 2 H 2 n 2 H ] ,
a function depending on their relative distance n; hence, they cannot be independent in the case H 1 / 2 . In the limit of infinitely distant increments, i.e., as n
ρ H ( n ) H ( 2 H 1 ) n 2 H 2 ,
from which we can deduce the rate of decay of the dependence between two increments.
  • If H < 1 / 2 , the dependence decays exponentially and ρ H ( n ) < 0 , so the increments have opposite signs, and the particle (or price, etc.) tends to return, denoting a mean reversion. Such a process is called anti-persistent and is said to exhibit short-range dependence.
  • If H > 1 / 2 , the dependence decays slower than exponential and ρ H ( n ) > 0 , so the particle tends to insist on the same direction. Such a process is called persistent and is said to exhibit long-range dependence.
B H is a self-similar process since for every fixed H
B α t H = d α H B t H , α > 0 ,
i.e., for any positive stretching factor α , the rescaled process α H B α t H with time scale α t , is equal in distribution to the original process B t H .
B H has stationary increments. The distribution of any increment B t H B s H depends only on the length t s of the time interval. It is worth noting that the fBm is the only Gaussian self-similar process with stationary increments.

2.2. mBm

The mBm is a nonstationary-dependent Gaussian process generalizing the fBm by replacing the constant Hurst exponent H with a Hölder function of time. The mBm W H = { W H ( t ) } t 0 with Hurst exponent H ( t ) and scaling factor σ > 0 , is the only centered Gaussian process being zero at the origin and with the covariance function defined for s , t [ 0 , 1 ] by
E ( W H ( t ) W H ( s ) ) = σ 2 2 g ( H ( t ) , H ( s ) ) | t | H ( t ) + H ( s ) + | s | H ( t ) + H ( s ) | t s | H ( t ) + H ( s ) ,
with
g ( H ( t ) , H ( s ) ) = K ( 2 H ( t ) ) K ( 2 H ( s ) ) K ( H ( t ) + H ( s ) ) , K ( a ) = Γ ( 1 + a ) sin ( a π / 2 ) .
Considering an ordinary Brownian motion B, we can also write the mBm in integral form
W H ( t ) = σ π K ( 2 H ( t ) ) Γ ( H ( t ) + 1 / 2 ) ( 0 ( ( t s ) H ( t ) 1 / 2 ( s ) H ( t ) 1 / 2 ) d B ( s ) + 0 t ( t s ) H ( t ) 1 / 2 d B ( s ) ) .
Notice that as | t s | 0 , it is reasonable to think that | H t H s | 0 too, so in this case, g ( H ( t ) , H ( s ) ) = 1 , and from the covariance formula, we can deduce the variance for the mBm increments:
Var ( W H ( t ) W H ( s ) ) = E ( ( W H ( t ) W H ( s ) ) 2 ) = E ( W H 2 ( t ) ) + E ( W H 2 ( s ) ) 2 E ( W H ( t ) W H ( s ) ) σ 2 | t | 2 H ( t ) + σ 2 | s | 2 H ( s ) σ 2 | t | 2 H ( t ) σ 2 | s | 2 H ( t ) + σ 2 | t s | 2 H ( t ) = σ 2 | t s | 2 H ( t ) .
which we will use in Section 4 to estimate the Hurst exponent from the data.
Although globally the mBm is not a self-similar process, we can recover a local version of this property. Indeed, let ε > 0 , u R , t R + , and consider a fBm X with Hurst index H. Since it is a self-similar process with stationary increments, we have
X ( ε u ) = d ε H X ( u ) , X ( t + ε u ) X ( t ) = d X ( ε u ) ,
combining the two properties, we have
X ( t + ε u ) X ( t ) ε H = d X ( u ) .
Intuitively, in the general case where H only governs the self-similar property of a process in a neighborhood of t, we have to pass to the limit as ε goes to zero. It can be proved [13] that given a fBm B H ( t ) with Hurst index H ( t ) , then a mBm W H with Hurst exponent H is locally asymptotically self-similar at t R +
lim ε 0 + W H ( t + ε u ) W H ( t ) ε H ( t ) u R = d B H ( t ) ( u ) u R .
We can think of this property as the existence of a fBm “tangent” to the mBm at each time t where the mBm is defined, i.e., the function H can be considered constant in a sufficiently small neighborhood of any t.

3. Simulation of Sample Paths

Several methods have been proposed to simulate the trajectories of the fBm. The first group of techniques exploits the properties of the fBm, while a second one aims to obtain a sample of a fBm indirectly via accumulated sums of a fGn. We consider the Wood–Chan method [14,15] as an efficient solution belonging to the latter group.
In particular, let us consider a Gaussian vector X N N ( 0 , G ) , whose N components X j = B j H B j 1 H , j = 1 , , N , represent the fBm increments. The idea of the method is to find the square root of G, in the sense that G = S S T for some square matrix S, by embedding G in a circulant symmetric matrix of size M = 2 ( N 1 )
C = ( c 0 c 1 c M 1 c M 1 c 0 c M 2   c 1 c 2 c 0 ) , c i = ρ H ( i ) , if 0 i N 1 , ρ H ( M i ) , if N i M 1 ,
where ρ H ( n ) is the covariance of the fGn (1). The peculiar feature of circulant matrices is that each row is shifted one element to the right relative to the preceding row; hence, only the first row is needed to construct a circulant matrix, which is a particular kind of Toeplitz matrix. Moreover, they are diagonalized by a discrete Fourier transform.
Recall that our aim is to generate the vector X containing the increments of the fBm. Suppose C is positive definite, since C is also symmetric, it is similar to the diagonal matrix containing its eigenvalues; it can be decomposed as
C = Q Λ Q ,
where Q is defined by
q j k = 1 M exp ( 2 π i j k M ) , j , k = 0 , , M 1 ,
and it is a unitary matrix, i.e., Q Q = Q Q = I , and Λ = diag ( λ 0 , , λ M 1 ) where
λ k = j = 0 M 1 c j exp ( 2 π i j k M )
are the eigenvalues of C. Now, define
Y = Q Λ 1 / 2 Q Z ,
where Z is a vector of M independent N ( 0 , 1 ) . Then Y N M ( 0 , C ) and its subvector ( Y 0 , , Y N 1 ) T N N ( 0 , G ) has the same distribution of the fGn vector X.
The Wood–Chan algorithm is executed through the following steps:
  • Compute the eigenvalues of C with the fast Fourier transform (FFT) of its first row.
  • Compute W = Q Z as a complex standard normal random vector: generate W 0 , W M / 2 N ( 0 , 1 ) , while for 0 < j < M / 2 , generate two independent U j , V j N ( 0 , 1 ) ; then compute
    W j = 1 2 ( U j + i V j ) , W M j = 1 2 ( U j i V j ) .
  • A sample of fGn is obtained by computing the FFT of Y = Q Λ 1 / 2 W
    X ( k ) = 1 M j = 0 M 1 λ j W j exp ( 2 π i j k M ) , k = 0 , , N 1 .
  • A sample of fBm at times k is given by the cumulative sum of X. Thanks to the self-similarity of fBm, the rescaled fBm on an interval [ 0 , T ] partitioned in N + 1 points with distance T / N , is obtained by multiplying the cumulative sum of X by ( T / N ) H . Finally, we have to set the first value to 0, since fBm is a process starting at 0.
In order to use all of the speeds of the FFT, N has to be chosen in such a way that M is a power of 2, i.e., N = 2 a + 1 , with a positive integer, so that M = 2 a + 1 . Examples of fBm paths generated with the Wood–Chan method are shown in Figure 1.
Exploiting the local self-similarity, we can generate a mBm sample path through the simulation of many fBm sample paths, one for each discretization point. Let k = i / N where 1 i N and N is the size of the sample. Given a function t H ( t ) , for each value of H ( k ) , an fBm B H ( k ) is generated. The mBm W H is then obtained by setting
W H ( k ) = B H ( k ) ( k ) ,
i.e., the k-th value of an mBm path is given by the k-th value of the fBm generated using the k-th value of the discretized Hurst function. Examples of mBm paths generated with this approach are shown in Figure 2.

4. Estimation of the Hurst Exponent

As we have seen within previous sections, processes underlying the generation of fBm and mBm sample paths rely only on the associated Hurst exponent. Therefore, the problem of accurately estimating such a quantity is of great interest. In the literature, most algorithms for estimating the Hurst exponents are based on a rolling window technique, where in each window a local Hurst exponent is estimated, e.g., the rescaled range analysis (R/S), originally introduced by Hurst himself while studying Nile river floods, and the detrended fluctuation analysis (DFA). Alternatively, there is a group of estimators that stems from the absolute moments of a Gaussian random variable. In what follows, we recall an estimator based on absolute moments of order 2.
Assume X is a discretized mBm of length N in the time interval [ 0 , 1 ] , with spacing equal to δ = 1 / ( N 1 ) , and that X locally behaves like a fBm within a window of length ε < N . Recalling the variance of an increment (3), we have
X ( j + 1 ) X ( j ) N ( 0 , σ 2 ( 1 N 1 ) 2 H ( t ) ) , t = ε + 1 , , N , j = t ε , , t 1 .
To estimate the function H, we follow the guidelines of the seminal work by Peltier and Véhel [16]. The absolute moment of order 2 of a normally distributed variable Y with zero mean and variance σ 2 is E ( Y 2 ) = 2 Γ ( 3 / 2 ) Γ ( 1 / 2 ) σ 2 . Define the quantity
M 2 ( t ) = 1 ε i = 0 ε 1 | X ( t i δ ) X ( t ( i + 1 ) δ ) | 2 ,
then using (4)
E ( M 2 ( t ) ) = 2 Γ ( 3 / 2 ) Γ ( 1 / 2 ) σ 2 ( 1 N 1 ) 2 H ( t ) .
It can be proved that M 2 ( t ) / E ( M 2 ( t ) ) p 1 as ε [2], so by equating their formulas, we can extract an estimate of H ( t ) . First, we need to get rid of the unknown σ . Define the quantity M in the same window as before, but with a halved resolution
M 2 ( t ) = 2 ε i = 0 ε / 2 1 | X ( t 2 i δ ) X ( t 2 ( i + 1 ) δ ) | 2 .
As N , their ratio can be written as
M 2 ( t ) M 2 ( t ) = 2 σ 2 ( ( N 1 ) / 2 ) 2 H ( t ) Γ ( 3 / 2 ) / Γ ( 1 / 2 ) 2 σ 2 ( N 1 ) 2 H ( t ) Γ ( 3 / 2 ) / Γ ( 1 / 2 ) = 2 2 H ( t ) ,
from which we can isolate H
H ^ ( t ) = 1 2 log 2 ( M 2 ( t ) M 2 ( t ) ) .

Smoothing

By definition, the Hurst exponent of a mBm is a continuous Hölder function; hence, it is reasonable to assume a smooth evolution of H over time. In Figure 3, we see that our estimate of H is raw and erratic, so we look for a way to eliminate short-term fluctuations and highlight long-term trends and patterns. In the following, we recall a method based on variational calculus, introduced by Garcin [17], in which we can think of smoothing H ^ as a trade-off problem between a fidelity term and a smoothness penalty.
Let L ( t , h ( t ) , h ˙ ( t ) ) = ( h ( t ) H ^ ( t ) ) 2 + λ h ˙ ( t ) 2 , with λ 0 quantifying the smoothness. We aim to minimize
H = argmin h C 2 ( [ 0 , 1 ] ; R ) 0 1 L ( t , h ( t ) , h ˙ ( t ) ) d t .
Let H ^ : [ 0 , 1 ] R be C 2 since h C 2 , also L C 2 , so we can apply the Euler–Lagrange equation
h L ( t , h ( t ) , h ˙ ( t ) ) d d t h L ( t , h ( t ) , h ˙ ( t ) ) = 0 .
The associated linear ODE can be solved with the method of variation of constants:
H ( t ) = e t / λ ( A 1 2 λ 0 t e s / λ H ^ ( s ) d s ) + e t / λ ( B + 1 2 λ 0 t e s / λ H ^ ( s ) d s ) .
The parameters A and B have to be adjusted in order to fit the minimal quadratic distance between H and H ^ . In a discrete framework, let 0 = t 1 < < t n = 1 and H ^ ( t 1 ) , , H ^ ( t n ) , and define the operator Φ d , such that
Φ d H ^ ( t i ) = e t i / λ 2 λ j = 0 i e t j / λ H ^ ( t j ) t j + 1 t j 1 2 + e t i / λ 2 λ j = 0 i e t j / λ H ^ ( t j ) t j + 1 t j 1 2 .
We look for the unique function H : { t 1 , , t n } R of the form
H : t Φ d H ^ ( t ) + A e t / λ + B e t / λ ,
such that i = 1 n ( H ( t i ) H ^ ( t i ) ) 2 is minimal. This is an ordinary least-square problem of the form Y = X a + u , where u is a residual, a = ( A , B ) T , Y = ( ( 1 Φ d ) H ^ ( t 1 ) , , ( 1 Φ d ) H ^ ( t n ) ) T and
X = ( e t 1 / λ e t 1 / λ e t n / λ e t n / λ ) .
The solution is provided by the estimate a ^ = ( X T X ) 1 X T Y , for any finite λ > 0 and n > 1 . After some computations, we find that the function H we are looking for is characterized by
( A B ) = 1 S 1 S 2 n 2 ( S 2 Z 1 n Z 2 S 1 Z 2 n Z 1 ) ,
with
S 1 = i = 1 n e 2 t i / λ , S 2 = i = 1 n e 2 t i / λ , Z 1 = i = 1 n e t i / λ ( H ^ ( t i ) Φ d H ^ ( t i ) ) , Z 2 = i = 1 n e t i / λ ( H ^ ( t i ) Φ d H ^ ( t i ) ) .
In the following simulations, we will use exponential smoothing as the baseline method to compare the performances of the variational smoothing. Given a sequence { x t } t 0 , and a smoothing factor 0 α 1 , the exponential moving average is given by the recursive formula
s t = x 0 , t = 0 , α x t + ( 1 α ) x t 1 , t > 0 .
Notice that, contrary to the variational method, here, the smoothing decreases as the smoothing factor increases.

5. Accuracy Test

Consider the time interval [ 0 , 1 ] discretized in N points; we use three accuracy metrics to compare the two smoothing methods studied in the previous Section.
  • The mean integrated squared error, representing the average accuracy of the method:
    R 0 = 1 N i = 1 N [ H ( i N ) H ( i N ) ] 2 .
  • The average bias at t = 1 , representing the ability of the technique not to create artificial gaps at the border:
    R 1 = H ( 1 ) H ( 1 ) .
  • The average quadratic slope, determining to which extent the erratic estimate has been smoothed:
    R 2 = N 2 N 1 i = 2 N [ H ( i N ) H ( i 1 N ) ] 2 .
Generate a mBm trajectory on the time interval [ 0 , 1 ] discretized in N = 2 8 points, using H ( t ) = 0.2 + 0.6 ( e t 2 1 ) / ( e 1 ) . Compute H ^ on a window of ε = 26 observations, and the smoothed estimates using λ = 0.02 and α = 0.04 , chosen to obtain overall good accuracy. The four curves are plotted in Figure 4.
The comparison of the two smoothing methods follows these steps:
  • Generate 100 mBm trajectories and compute the associated raw estimates.
  • Define a set of 200 logarithmically spaced (the variational smoothing method is sensitive to small values of λ ) parameters λ between decades 10 3 and 10 3 , and a set of 200 evenly spaced parameters α between 0.02 and 0.07. Compute the smooth estimates of the raw ones from the previous step.
  • For each smooth estimate, compute the three accuracy metrics.
  • For each parameter, compute the average metrics and visualize the results with plots: R2 vs. R0 and R2 vs. R1.
From Figure 4, we observe that the variational method smooths the raw Hurst exponent more accurately, since (almost all of the time) it is closer to H than the exponential smoothing estimate. In particular, the exponential moving average slightly underestimates the raw estimate after the initial part, both for low and high t. The method mitigates a value with the preceding ones, which in this case are mostly lower since the exponential Hurst exponent monotonically increases over time. This behavior seems more evident at the boundaries of the time interval, while the smoothing in the middle part is more consistent with the raw estimate.
The superior accuracy of the variational smoothing over the exponential moving average can also be observed when comparing the three accuracy metrics in Figure 5 and Figure 6. We first note that, in all cases, the average quadratic slopes R 2 vary monotonically within the range of the parameters. Moreover, for a fixed R 2 , the value of R 0 is always closer to zero for the variational smoothing; however, as the smoothing increases, the value of R 0 decreases only up to a certain value of λ , which we found to be approximately 0.08. Furthermore, R 1 is lower than the one for exponential smoothing only if λ < 0.06 . From these observations, it follows that a good trade-off is to choose a small value for λ , and we found a precise threshold to abide.
One more thing to notice is that the average bias in t = 1 is negative for both smoothing methods and for all of the values of the parameters. Since its formula is simply R 1 = H ( 1 ) H ( 1 ) , this means that the techniques underestimate H ( 1 ) , whatever the parameter, and such behavior is probably inherent in the growing feature of the exponential Hurst function used in the simulation. Finally, observe that, in general, the bias does not asymptotically approach zero, because even if the smooth estimates converge toward the raw series for decreasing λ and increasing α , this does not imply that the raw series converges toward the real Hurst function, since the raw estimate is independent of the choice of parameters λ and α .

6. Simulation of Financial Data

In this section, we recall the procedure by Garcin [17] to forecast the trend of particular financial assets, the foreign exchange rates. Let { X 1 , , X n } be the time series of log prices of a financial asset driven by a mBm. Let 0 , δ , , T δ be the equispaced times, and τ < T define the length of a sample in which we smooth the raw estimate H ^ . For t { 0 , , T τ } , we make successive out-of-sample predictions for X ( ( t + τ ) δ ) , based only on the knowledge of { X ( s δ ) } t s t + τ 1 .
  • Compute the smooth version { H ( s δ ) } t s t + τ 1 of the raw estimate { H ^ ( s δ ) } t s t + τ 1 .
  • Forecast H ( ( t + τ ) δ ) by means of a second-order Taylor expansion. Consider a smooth function f on an equispaced grid with step h = 1 , and approximate the first derivative with the first order backward difference and the second derivative with the second order backward difference:
    f ( x + 1 ) f ( x ) + f ( x ) + 1 2 f ( x ) f ( x ) + ( f ( x ) f ( x 1 ) ) + 1 2 ( f ( x ) 2 f ( x 1 ) + f ( t 2 ) ) = 5 2 f ( x ) 2 f ( x 1 ) + 1 2 f ( x 2 ) .
    Hence, applying the same rule for the smooth Hurst exponents, we obtain
    H ( ( t + τ ) δ ) = 5 2 H ( ( t + τ 1 ) δ ) 2 H ( ( t + τ 2 ) δ ) + 1 2 H ( ( t + τ 3 ) δ ) .
  • Estimate the scale parameter σ by Equations (5) and (6) in the case of the absolute moment of order 1 and solving for σ
    σ = N H ( t ) 1 ε π 2 i = 0 ε 1 | X ( t i δ ) X ( t ( i + 1 ) δ ) | .
    The raw estimation H ^ and its smoothing H are distinct steps of the estimation procedure, but it is reasonable to assume that the time steps in both computations are the same. Therefore, we can write the estimate of the scale parameter as
    σ ^ = 1 τ π 2 s = t t + τ 1 | X ( s δ ) X ( ( s 1 ) δ ) | δ H ( s δ ) .
  • Estimate the underlying standard Bm W that appears in the definition of the mBm (2). Firstly, notice that in this case t 0 , so the first integral in the definition vanishes. Moreover, since at the origin t of the time interval the mBm has to be 0, we translate X by X ( t δ ) , obtaining
    Γ ( H ( s δ ) + 1 / 2 ) σ Γ ( 2 H ( s δ ) + 1 ) sin ( H ( s δ ) π ) ( X ( s δ ) X ( t δ ) ) = 0 s ( ( s u ) δ ) H ( s δ ) 1 / 2 d W ( u δ ) .
    The integral on the rhs, call it I, can be approximated using the Riemann sum as
    I i = 1 n ( ( s u i 1 ) δ ) H ( s δ ) 1 / 2 ( W ( u i δ ) W ( u i 1 δ ) ) .
    Moving the first δ on the other side and defining
    Π ( s δ ) = Γ ( H ( s δ ) + 1 / 2 ) σ ^ δ H ( s δ ) 1 / 2 Γ ( 2 H ( s δ ) + 1 ) sin ( H ( s δ ) π ) ,
    we can write
    Π ( s δ ) ( X ( s δ ) X ( t δ ) ) i = 1 n ( s u i 1 ) H ( s δ ) 1 / 2 ( W ( u i δ ) W ( u i 1 δ ) ) .
    Now, since ( s u i 1 ) equals 1 when i = n , we can move the term W ( s δ ) W ( ( s 1 ) δ ) out from the series, obtaining
    W ^ ( s δ ) W ^ ( ( s 1 ) δ ) + ( X ( s δ ) X ( t δ ) ) Π ( s δ ) 1 s > t u = t s 1 [ W ( u δ ) W ( ( u 1 ) δ ) ] ( s u + 1 ) H ( s δ ) 1 / 2 .
  • Recall that X is a log price, i.e., X ( t δ ) = log P ( t δ ) , where P is the price of the financial asset; thus, its increment can be written as
    X ( ( t + τ ) δ ) X ( ( t + τ 1 ) δ ) = log P ( ( t + τ ) δ ) P ( ( t + τ 1 ) δ ) = r ( ( t + τ ) δ ) ,
    where r is known as log return. We will exploit this relation to compute the forecasted log price X ( ( t + τ ) δ ) . However, first, we have to derive the formula for the log return. Consider the increment of the mBm at time ( t + τ ) δ , then add and subtract X ( t δ ) and rearrange the terms in order to reach a form resembling the left-hand side of Equation (7) (with the term Π ( s δ ) moved to the right-hand side)
    X ( ( t + τ ) δ ) X ( t δ ) X ( ( t + τ 1 ) δ ) X ( t δ ) .
    Applying (7) to the two terms and combining the two series, we arrive at the formula for the log price
    r ( ( t + τ ) δ ) = s = t t + τ 1 ( W ^ ( s δ ) W ^ ( ( s 1 ) δ ) ) ( ( t + τ s + 1 ) H ( ( t + τ ) δ ) 1 / 2 Π ( ( t + τ ) δ ) ( t + τ s ) H ( ( t + τ 1 ) δ ) 1 / 2 Π ( ( t + τ 1 ) δ ) ) .
Notice that each time this procedure lets us compute a single future price; therefore, we must apply it several times until we reach the horizon T.

6.1. One-Step Forecasting

We apply the forecasting procedure to several foreign-exchange rates versus the Euro (USD/EUR, GBP/EUR, etc.). Prices were observed every 15 min from 7 March to 7 September 2016 (ca. 12900 observations). For each series, the raw estimate H ^ of the Hurst exponent is computed with a window of ε = 26 observations. Recall that this formula does not introduce unfairness in the forecast since it is only based on preceding observations. We compute the smooth estimate H on a sample of τ = 48 raw estimates (12 h) using λ = 0.02 .
In Figure 7, notice that H may vary significantly even in a small time span of 12 h: a strong justification to use the mBm to analyze financial time series.
As a measure to check the goodness of the results, we compute the hit ratio, which is the proportion of times in which the trend is correctly forecasted. Results are shown in Table 1. All of the hit ratios are above 50%, but while they are near to an average of 59.3% in the case of the raw series, they are considerably better in the case of the variational smoothing, with an average of 71.6%.
Let us also consider the behavior of the hit ratio conditionally to the Hurst exponent. In particular, for each currency, we consider the time series of the forecasted Hurst exponents computed during the forecasting procedure, then sort it in ascending order. Therefore, the first values are the ones near to 0, the last ones being closer to 1. Then we split the interval [ 0 , 1 ] into 20 subintervals of length 0.05 and compute the hit ratio for each one of them, which counts how many forecasted Hurst exponents lie in each interval and for how many of them the forecast is good, i.e., for how many of them the trend of the corresponding forecasted exchange rate has the same sign as the trend of the real exchange rate.
The results of this procedure are shown in Figure 8. We underline that value 0.5 of the Hurst exponent basically divides the plots into two parts. For values lower than 0.5, the hit ratios are usually low, while for values higher than 0.5, the hit ratios are high almost every time. Since such scenarios happen for all of the currencies, we can reasonably state that if the estimates for the Hurst exponents are higher than the threshold, the predictions are more reliable, as indeed confirmed in the literature, see, e.g., [18].
A possible explanation for this phenomenon relies on the fBm properties. Indeed, this process is characterized by long memory when the Hurst index is greater than 0.5, contrary to lower values. Moreover, for some of the currencies, we can notice that the hit ratio decreases as the Hurst exponent approaches value 1. Such high values of the Hurst exponent correspond to a process with a very long memory; hence, for a forecast to be reliable, it should be computed using many observations in the past; that is, we should use a very large window in the forecasting procedure.
In Table 2, we report the most common measures of errors e i = y i y ^ i between the observed values y i and the predicted values y ^ i . Namely, the mean absolute error MAE = 1 n | e i | , the mean-square error MSE = 1 n e i 2 , the root-mean-square error RMSE = MSE , and the relative error η = e 2 / y 2 . The low orders of magnitude of the errors are determined by two concurrences. Firstly, recall that the prices in the time series are observed every 15 min; thus, they vary by a small amount at each time step. For example, in the case of USD/EUR, the average variation between adjacent prices is 3.3 × 10−4, and the standard deviation of the prices is 4.1 × 10−4. Secondly, recall that the method uses a window of past prices to predict only one future price; hence, it is reasonable to think that the forecasted value will not be so far away from the current price, as seen in Figure 9.
For the sake of completeness, let us underline that a limited time window may not allow for robust results. This could occur for many reasons, such as, e.g., market manipulations or a financial crisis. Therefore, testing on different time windows allows for improving the reliability of the results. Analogously to what we have previously done, we apply our forecasting procedure to both GBP/EUR and USD/EUR time series, but consider a significantly larger window, namely from 2016 to 2021. Prices are observed every 15 min from 3 January 2016 to 31 December 2021, for a total of about 153,000 samples. The raw estimate H ^ of the Hurst exponent is computed exploiting a window of ε = 26 observations, while the smooth estimate H is derived considering a sample of τ = 48 raw estimates (12 h) using λ = 0.02 .
As before, in Table 3 we show the hit ratio to measure the goodness of the results. Note that, differently from what was already reported, this time we have one value for each year of analysis. Again, all of the hit ratios are above 50%, and are considerably better in the case of variational smoothing, with percentages similar to what we obtained in the previous test.
In Figure 10, we consider the behavior of the hit ratio conditionally to the Hurst exponent. Moreover, in this simulation, the value H = 0.5 divides the plots into two parts, with low hit ratios on the left and high values on the right. Results confirm that predictions are more reliable if the estimates for the Hurst exponents are higher than the threshold, but with a tendency to decrease in the case of very long memory, i.e., as H 1 .
Finally, in Table 4, we present the measures of errors between the observed and predicted values for the time series USD/EUR. Errors are similar to the ones obtained in the previous test, with smaller values, especially for the year 2019.

6.2. Multi-Step Forecasting

In the previous test, each future price was computed using a window of past observations. In the following, we attempt a more challenging task, using past predictions to compute new ones; in particular, using a sliding window containing both observations and predictions. There are two main possibilities to run such a procedure: (1) train a model on past real data to then predict multiple future data at once, which is the typical machine learning approach; (2) adopt a multi-stage prediction solution that consists of computing one future value at a time, using the new data to compute the following ones.
Adopting the same parameters as in the previous test, we implemented the latter approach to compute 24 future prices (6 h). Figure 11 depicts both the real and the forecasted prices for the USD/EUR in a time frame of about 700 time steps, with the smaller range shown in Figure 9. This time, the model was not able to capture the price variation, but rather it acted as a sort of piecewise approximation of the path, where changes in height correspond to the start of a new multi-step window (marked in the figure with black circles).
With respect to the one-step forecasting, the errors MAE = 1.3e-03 and MSE = 5.6e-06 are bigger by two orders of magnitude, while RMSE = 2.4e-03 and the relative error = 2.1e-03 are about twenty times bigger. Moreover, the hit ratio has drastically worsened, dropping at 51.1%.
The poor performance of the multi-stage approach is likely due to two factors:
  • The fast variation of H in time (as seen in Figure 7), which seems to be hard to catch with this approach.
  • A flaw inherent in the very nature of the multi-stage prediction: the propagation of errors. This approach uses past predicted values to compute new ones, so previous errors are propagated into future predictions.

7. Conclusions

In this work, we studied a method based on the absolute moment of order 2 to obtain a raw estimate of the mBm Hurst exponent. To reduce the noise of the estimate and move closer to the real Hurst function, we studied a smoothing method based on variational calculus and compared it with exponential smoothing, which turned out to be (globally) less accurate. In particular, we found a threshold of 0.06 for the penalizing parameter λ . Below this value, the variational smoothing is always more accurate than the exponential one. We believe this finding represents a step forward in the understanding of financial dynamics, offering a robust tool that can be possibly used as a teammate of other approaches based on more classical stochastic process frameworks.
After deploying all of the necessary tools, we also considered a procedure to forecast foreign exchange rates. We applied it in both a one-step and a multi-step approach, assessing the goodness of the results with the hit ratio and other widespread measures of errors. Analyzing the relationship between the hit ratio and the local Hurst exponent, it turned out that in the presence of long-range dependence, in particular for 0.5 < H < 0.9 , making some predictions on the Forex market is quite promising since the hit ratio is usually higher than 50%. However, the same does not apply to the multi-stage prediction test, in which the model was not able to capture the variation of the prices.
Future research will be based on refinements of the formulas we derived for the raw estimate of H, as well as for the forecast value H ( ( t + τ ) δ ) . Analogously, we aim at improving the results we obtained for the multi-step forecasting by implementing optimization algorithms, e.g., the stochastic gradient descent to find the model’s parameters, minimizing the approximation error. Finally, it would be interesting to replace the time-varying Hurst exponent with a stochastic process, which is to study the multi-fractional process with a random exponent; hence, generalizing the mBm approach.

Author Contributions

All authors contributed equally to the presenting research and preparing the publication. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Forex market data provided by “Centro Interdipartimentale di Documentazione Economica” (C.I.D.E.) resources, via University of Verona access.

Acknowledgments

The authors would like to thank the “Centro Interdipartimentale di Documentazione Economica” (C.I.D.E.) for providing the Forex data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cont, R. Empirical properties of asset returns: Stylized facts and statistical issues. Quant. Financ. 2001, 1, 223–236. [Google Scholar] [CrossRef]
  2. Bianchi, S. Pathwise identification of the memory function of multifractional Brownian motion with application to finance. Int. J. Theor. Appl. Financ. 2005, 8, 255–281. [Google Scholar] [CrossRef]
  3. Biagini, F.; Sulem, A.; Øksendal, B.; Wallner, N.N. An introduction to white-noise theory and Malliavin calculus for fractional Brownian motion. Proc. R. Soc. 2004, 460, 347–372. [Google Scholar] [CrossRef]
  4. Rogers, L.C.G. Arbitrage with fractional Brownian motion. Math. Financ. 1997, 7, 95–105. [Google Scholar] [CrossRef] [Green Version]
  5. Cheridito, P. Arbitrage in fractional Brownian motion models. Financ. Stoch. 2003, 7, 533–553. [Google Scholar] [CrossRef]
  6. Staszkiewicz, P.; Łosiewicz-Dniestrzańska, E.; Grygiel-Tomaszewska, A. European financial institution physical geolocation and the high-frequency trading potential. In The Digitalization of Financial Markets; Marszk, A., Lechman, E., Eds.; Routledge: London, UK, 2021; pp. 19–37. [Google Scholar]
  7. Bezborodov, V.; Di Persio, L.; Mishura, Y. Option Pricing with Fractional Stochastic Volatility and Discontinuous Payoff Function of Polynomial Growth. Methodol. Comput. Appl. Probab. 2019, 21, 331–366. [Google Scholar] [CrossRef] [Green Version]
  8. Alvarez-Ramirez, J.; Alvarez, J.; Rodriguez, E.; Fernandez-Anaya, G. Time-Varying Hurst Exponent for US Stock Markets. Phys. A 2008, 387, 6159–6169. [Google Scholar] [CrossRef]
  9. Jiang, J.; Ma, K.; Cai, X. Non-linear characteristics and long-range correlations in Asian stock markets. Phys. A 2007, 378, 399–407. [Google Scholar] [CrossRef]
  10. Bianchi, S.; Pianese, A. Multifractional Processes in Finance. SSRN Electron. J. 2014, 5, 1–22. [Google Scholar] [CrossRef]
  11. Corlay, S.; Lebovits, J.; Vehel, J.L. Multifractional Stochastic volatility models. Math. Financ. 2014, 24, 364–402. [Google Scholar] [CrossRef]
  12. Park, J. Estimating the DJI Series by Multifractional Brownian Motion; Econometric Modeling: International Financial Markets; SSRN: Rochester, NY, USA, 2020. [Google Scholar]
  13. Benassi, A.; Jaffard, S.; Roux, D. Elliptic Gaussian Random Processes. Rev. Mat. Iberoam. 1997, 13, 19–90. [Google Scholar] [CrossRef] [Green Version]
  14. Coeurjolly, J.-F. Simulation and identification of the fractional Brownian motion: A bibliographical and comparative study. J. Stat. Softw. 2000, 5, 1–53. [Google Scholar] [CrossRef] [Green Version]
  15. Wood, A.; Chan, G. Simulation of stationary Gaussian processes in [0, 1]d. J. Comput. Graph. Stat. 1994, 3, 409–432. [Google Scholar] [CrossRef]
  16. Peltier, R.F.; Lévy Véhel, J. Multifractional Brownian Motion: Definition and Preliminary Results; Rapport de recherche de l’INRIA: Paris, France, 1995; p. 2645. [Google Scholar]
  17. Garcin, M. Estimation of time-dependent Hurst exponents with variational smoothing and application to forecasting foreign exchange rates. Phys. Stat. Mech. Its Appl. 2017, 483, 462–479. [Google Scholar] [CrossRef]
  18. Mitra, S.K. Is Hurst exponent value useful in forecasting financial time series? Asian Soc. Sci. 2012, 8, 111–120. [Google Scholar] [CrossRef]
Figure 1. fBm sample paths of length N = 2 9 for different values of the Hurst exponent.
Figure 1. fBm sample paths of length N = 2 9 for different values of the Hurst exponent.
Symmetry 14 01657 g001
Figure 2. mBm trajectories (blue) of length N = 2 9 generated, respectively, with the Hurst exponent (red) H ( t ) = 0.6 t + 0.2 , t [ 0 , 1 ] and H ( t ) = sin ( t ) / 4 t + 1 / 2 , t [ 0 , 4 π ] .
Figure 2. mBm trajectories (blue) of length N = 2 9 generated, respectively, with the Hurst exponent (red) H ( t ) = 0.6 t + 0.2 , t [ 0 , 1 ] and H ( t ) = sin ( t ) / 4 t + 1 / 2 , t [ 0 , 4 π ] .
Symmetry 14 01657 g002
Figure 3. Raw estimate (blue) of a mBm generated with an exponential Hurst function (red) H ( t ) = 0.2 + 0.6 ( e t 2 1 ) / ( e 1 ) , with N = 2 8 observations and a window of length ε = 26 .
Figure 3. Raw estimate (blue) of a mBm generated with an exponential Hurst function (red) H ( t ) = 0.2 + 0.6 ( e t 2 1 ) / ( e 1 ) , with N = 2 8 observations and a window of length ε = 26 .
Symmetry 14 01657 g003
Figure 4. Real Hurst exponent and its estimates.
Figure 4. Real Hurst exponent and its estimates.
Symmetry 14 01657 g004
Figure 5. Averaged accuracy metrics R 2 and R 0 for the variational and exponential smoothing.
Figure 5. Averaged accuracy metrics R 2 and R 0 for the variational and exponential smoothing.
Symmetry 14 01657 g005
Figure 6. Averaged accuracy metrics R 2 and R 1 for the variational and exponential smoothing.
Figure 6. Averaged accuracy metrics R 2 and R 1 for the variational and exponential smoothing.
Symmetry 14 01657 g006
Figure 7. Raw Hurst exponent estimate (blue) and smoothed estimate (red) for the USD/EUR exchange rate, from 20:30 on 6 April to 8:30 on 7 April 2016.
Figure 7. Raw Hurst exponent estimate (blue) and smoothed estimate (red) for the USD/EUR exchange rate, from 20:30 on 6 April to 8:30 on 7 April 2016.
Symmetry 14 01657 g007
Figure 8. USD/EUR. The frequency with which the model correctly forecasts the sign of the return (i.e., the trend), as a function of the local Hurst exponent. Data from 7 March to 7 September 2016.
Figure 8. USD/EUR. The frequency with which the model correctly forecasts the sign of the return (i.e., the trend), as a function of the local Hurst exponent. Data from 7 March to 7 September 2016.
Symmetry 14 01657 g008
Figure 9. USD/EUR. Comparison between the forecasted and real prices.
Figure 9. USD/EUR. Comparison between the forecasted and real prices.
Symmetry 14 01657 g009
Figure 10. USD/EUR. The frequency with which the model correctly forecasts the sign of the return (i.e., the trend), as a function of the local Hurst exponent.
Figure 10. USD/EUR. The frequency with which the model correctly forecasts the sign of the return (i.e., the trend), as a function of the local Hurst exponent.
Symmetry 14 01657 g010
Figure 11. USD/EUR. Comparison between the forecasted and real prices.
Figure 11. USD/EUR. Comparison between the forecasted and real prices.
Symmetry 14 01657 g011
Table 1. Values of the hit ratio for each currency used in the tests, computed for the raw series (column H ^ ) and for the smoothed estimate (column H ). Data from 7 March to 7 September 2016.
Table 1. Values of the hit ratio for each currency used in the tests, computed for the raw series (column H ^ ) and for the smoothed estimate (column H ). Data from 7 March to 7 September 2016.
Currency vs. Euro H ^ H
AUD/EUR59.5%72.3%
CHF/EUR57.5%68.5%
GBP/EUR59.6%71.8%
JPY/EUR60.4%73.4%
SEK/EUR58.1%70.2%
SGD/EUR58.5%70.5%
USD/EUR59.8%73.3%
Table 2. Various measures of the differences between the observed values and values predicted by the model incorporating the variational smoothing method. Data from 7 March to 7 September 2016.
Table 2. Various measures of the differences between the observed values and values predicted by the model incorporating the variational smoothing method. Data from 7 March to 7 September 2016.
Currency vs. EuroMAEMSERMSERelative Error
AUD/EUR1.0 × 10−41.7 × 10−74.2 × 10−42.8 × 10−4
CHF/EUR3.2 × 10−55.1 × 10−97.1 × 10−56.5 × 10−5
GBP/EUR5.7 × 10−51.8 × 10−81.3 × 10−41.7 × 10−4
JPY/EUR9.3 × 10−33.8 × 10−42.0 × 10−21.6 × 10−4
SEK/EUR3.3 × 10−44.7 × 10−76.9 × 10−47.4 × 10−5
SGD/EUR5.7 × 10−51.1 × 10−81.1 × 10−47.0 × 10−5
USD/EUR5.3 × 10−51.3 × 10−81.1 × 10−41.0 × 10−4
Table 3. Yearly values of the hit ratios for GBP and USD, computed for the raw series (column H ^ ) and the smoothed estimate (column H ).
Table 3. Yearly values of the hit ratios for GBP and USD, computed for the raw series (column H ^ ) and the smoothed estimate (column H ).
YearGBP/EURUSD/EUR
H ^ H H ^ H
201659.4%71.8%60.2%73.5%
201758.3%70.5%61.0%75.0%
201857.5%69.2%60.4%73.9%
201958.4%70.7%59.6%72.6%
202059.8%71.8%60.7%74.9%
202159.6%71.7%61.2%74.9%
Table 4. USD/EUR. Yearly measures of the differences between the observed values and the values predicted by the model incorporating the variational smoothing method.
Table 4. USD/EUR. Yearly measures of the differences between the observed values and the values predicted by the model incorporating the variational smoothing method.
YearMAEMSERMSERelative Error
20166.0 × 10−51.4 × 10−81.2 × 10−41.1 × 10−4
20175.3 × 10−51.2 × 10−81.1 × 10−49.7 × 10−5
20185.7 × 10−59.8 × 10−99.9 × 10−58.4 × 10−5
20193.6 × 10−54.5 × 10−96.7 × 10−56.0 × 10−5
20205.8 × 10−51.0 × 10−81.0 × 10−48.9 × 10−5
20214.6 × 10−58.0 × 10−99.0 × 10−57.6 × 10−5
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Di Persio, L.; Turatta, G. Multi-Fractional Brownian Motion: Estimating the Hurst Exponent via Variational Smoothing with Applications in Finance. Symmetry 2022, 14, 1657. https://doi.org/10.3390/sym14081657

AMA Style

Di Persio L, Turatta G. Multi-Fractional Brownian Motion: Estimating the Hurst Exponent via Variational Smoothing with Applications in Finance. Symmetry. 2022; 14(8):1657. https://doi.org/10.3390/sym14081657

Chicago/Turabian Style

Di Persio, Luca, and Gianni Turatta. 2022. "Multi-Fractional Brownian Motion: Estimating the Hurst Exponent via Variational Smoothing with Applications in Finance" Symmetry 14, no. 8: 1657. https://doi.org/10.3390/sym14081657

APA Style

Di Persio, L., & Turatta, G. (2022). Multi-Fractional Brownian Motion: Estimating the Hurst Exponent via Variational Smoothing with Applications in Finance. Symmetry, 14(8), 1657. https://doi.org/10.3390/sym14081657

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop