A Monte Carlo Approach to Bitcoin Price Prediction with Fractional Ornstein–Uhlenbeck Lévy Process

: Since its inception in 2009, Bitcoin has increasingly gained main stream attention from the general population to institutional investors. Several models, from GARCH type to jump-diffusion type, have been developed to dynamically capture the price movement of this highly volatile asset. While ﬁtting the Gaussian and the Generalized Hyperbolic and the Normal Inverse Gaussian (NIG) distributions to log-returns of Bitcoin, NIG distribution appears to provide the best ﬁt. The time-varying Hurst parameter for Bitcoin price reveals periods of randomness and mean-reverting type of behaviour, motivating the study in this paper through fractional Ornstein–Uhlenbeck driven by a Normal Inverse Gaussian Lévy process. Features such as long-range memory are jump diffusion processes that are well captured with this model. The results present a 95% prediction for the price of Bitcoin for some speciﬁc dates. This study contributes to the literature of Bitcoin price forecasts that are useful for Bitcoin options traders.


Introduction
Cryptocurrency is a digital currency that uses advanced encryption techniques known as cryptography to secure financial transactions, control the creation of additional units, and verify the transfer of assets. It uses decentralized control by using distributed ledger technology called Blockchain, which serves as a public financial transaction database. Bitcoin is the first decentralized cryptocurrency introduced in 2009 by a person or a group of persons under the pseudonym Sathoshi Nakamoto [1]. Many other alternatives of Bitcoin have been created. On Coinmarketcap (https://coinmarketcap.com/, accessed on 25 February 2022), there are 2322 tradable cryptocurrencies with a total market capitalization of USD 345 billions at the time of this paper (see [2,3] for further background on the Bitcoin and its technology). The cryptocurrency market is known to be highly volatile [4][5][6][7] due, on the one hand, to its sensibility to new information, whether fundamental or speculative [8], since it does not rely on the stabilizing policy of a central bank. On the other hand, the relative liquidity of the market with no official market makers makes it fundamentally fragile to large trading volumes and to market imperfections; thus, it is more prone to large swings than other traded assets (see [9]). In such environments, suitable forecasting models for cryptocurrencies have to reflect all these underlying features. Bitcoin prices have been already modeled and predicted by means of a noncausal autoregressive process [10], fractional geometric Brownian motion [11], and machine learning [12,13] to name the few. Andrew and Christian [10] applied the noncausal autoregressive process with Cauchy errors to model and predict Bitcoin prices. Their results display episodes of local trends, which can be modeled and interpreted as speculative bubbles. Tarnopolski [11] used the fractional geometric Brownian motion to account for the long-memory or long-term dependence of Bitcoin (BTC), manifesting itself through a Hurst exponent H > 0.5, in order to predict future BTC/USD price. To achieve this, a Monte Carlo simulation with 10 4 geometric fractional Brownian motion realizations was performed as extensions of historical data, and the accuracy of statistical inferences is taken to be 10%. Tarnopolski [11] concluded the prediction that the most probable Bitcoin prices at the beginning of 2018 will be 6358 USD. Jalali and Heidari [14] predicted the price of Bitcoin and changes therein using a first-order Grey model (GM (1,1)). Other studies on cryptocurrency prediction using variants of Grey models can be found in Gatabazi et al. [15][16][17]. However, the grey system theory prediction works better only with small size datasets, as the error of prediction will increase as the dataset becomes larger (Wu et al. [18]). Moreover, small size dataset does not exhibit some of the important characteristics of the data, such as long range dependence as well as the size and frequency of jumps. McNally et al. ascertained information on what accuracy the direction of Bitcoin price can be predicted. This task is achieved with varying degrees of success through the implementation of a Bayesian optimized recurrent neural network (RNN) and a Long Short Term Memory (LSTM) network. They showed that these non-linear deep learning methods outperform the popular ARIMA model for time series forecasting. As shown by Jalali and Heidari [14], these neural network models are sensitive to input variables, given that the Bitcoin price depends on different inputs with complex behaviors, as illustrated by Chen et al. [19] and Kristoufek and Ladislav [20]. Based on new technologies, economic policies, and cultural behaviors, these inputs may change. Motivated by these shortcomings, we introduce a mean-reverting fractional jump-diffusion stochastic model in this study to predict Bitcoin exchange rates.
Since Bitcoin exchange rates not only exhibit long-range dependence [21], they also present some periods of randomness and mean-reverting feature. Moreover, they exhibit abrupt and discontinuous changes in prices known as jumps, which are assumed to follow some probability law [9]. We present in this paper a fractional Ornstein-Uhlenbeck driven by a Normal Inverse Gaussian Lévy process for short-term predictions of Bitcoin exchange rates.
The rest of the paper is organized as follows: Section 2 introduces the methodology used, Section 3 discusses the empirical findings, and the last Section concludes the study.

Materials and Methods
In this section, the basic structure of the fractional Ornstein-Uhlenbeck Lévy (fOUL) model is described where the Bitcoin price process at any time point t, {X t } t≥0 is modeled as an exponential of a Lévy process. The model augments two continuous-time random processes: a fractional Brownian process and a compound Poisson jump process, which also generates two random sources. Fractional Brownian motion could be divided into three different categories: The process is likely to reverse the trend over the time frame considered.

2.
The process is random in which knowing one data point does not provide insight into predicting future data points.

3.
The process is persistent in the sense that a future data point is likely to be similar to a data point preceding it.
This is of particular importance in examining long-range dependence in processes. The below subsection presents some details of the a fractional Brownian motion.
The increments (b H (k)) follow a standard normal distribution with zero mean and variance The corresponding autocovariance function γ H (·) is of the following form.
, all the covariances are zero (except for k = 0), implying independence. This is in agreement with the standard Brownian motion, which has independent increments. Using the Taylor expansion at the origin of the function h( In particular, we have the following.
Therefore, the following is the case.
Hence, only in the case of H > 1 2 , fractional Brownian motions display long-memory dependence. For a fixed real number H ∈ ( 1 2 , 1), let {B H t , t ≥ 0} be a fractional Brownian motion (fBm) with Hurst parameter H on a complete probability space (Ω, F , P). The integral with respect to B H in the generalized sense of Lebesgue-Stieltjes is defined using the fractional derivatives for a < b and α ∈ (0, 1) as follows.
From the Hölder continuity of B H , it follows that, for α ∈ (1 − H, 1), The tracking of the variation of the Hurst parameter H Figure 1 suggests a meanreverting process for the data (since H < 1/2 for a large part of the time period). This justifies our choice for a fractional Ornstein-Uhlenbeck Lévy process. Let us start by defining a Lévy process. A Lévy process L = (L t ) t≥0 can be decomposed into a Brownian motion W = (W t ) t≥0 and a pure jump process Z = (Z t ) t≥0 independent from the Brownian motion as follows.
The building block of the pure jump process Z is the Poisson process. While fitting the Normal, Hyperbolic, Generalized Hyperbolic, and Normal Inverse Gaussian (NIG) distributions to the log returns of Bitcoin data (see Figure 2), NIG displays the best fit. We consider, in this work, the fractional Ornstein-Uhlenbeck Lévy process where the jump part is modeled by an NIG process. Before presenting this model, it is important to recall what a Levy process as well as the Normal Inverse Gaussian distribution are.

Lévy Process 2.2.1. Infinitely Divisible Distributions
A random variable X is called infinitely divisible, if for each positive integer n ≥ 2, there exist independent and identically distributed (i.i.d) random variables X 1 , X 2 , · · · , X n such that the distribution of Y is the same as the distribution of ∑ n k=1 X k and has a characteristic function. In the literature, the characteristic function of the one-dimensional infinitely divisible distribution is generalized by the Lévy-Khinchin-Kolmogorov formula given below.
Let us start by defining the cumulative generating function of the random vector x = (x 1 , · · · , x n ) ∈ R n , which determines the Lévy measure.
The p th order joint cumulant of x is given by the following: where E is the expectation operator, and the summation extends over all partitions (p 1 , · · · , p k ), k = 1, · · · , n, of (1, · · · , n). Note that cum(x) is the coefficient of i n ξ 1 · · · ξ n in the Taylor series expansion about the origin of the following function, known as the cumulant generating function of the random vector x.
If function C(ξ, x) has the following Lévy-Khinchin-Kolmogorov form: then the random variable x follows an infinitely divisible distribution, where the following is the case: a, b ∈ R, b ≥ 0, and Q is a Radon measure on R with the following: and where function φ in Equation (4) is of the following form.
The measure Q in Equation (4) is called the Lévy measure of x. The variable a is referred to as the center or drift and determines the location, and b is the variance. The Lévy triplet (a, b, Q) determines x uniquely. If b = 0, then the distribution is referred to as a purely non-Gaussian distribution.

Continuous-Time Stochastic Processes
Using infinitely divisible distributions, one can define a continuous sequence of independently infinitely divisible random variables, which will be referred to as the continuoustime stochastic processes. The two basic classes of continuous-time stochastic processes are Brownian motion and the Poisson process. The Poisson process generated by the Poisson distribution is the building block of pure jump processes. Both processes are fundamentally different concerning their path properties and the fact that they belong to the larger class of Lévy processes (for more details about Lévy processes, see [5]. Cont and Tankov [6] provide details of Lévy processes with applications to finance). A stochastic process X = (X t ) t≥0 is a family of R-valued random variables X t with parameter t ≥ 0, defined on the sample space Ω.
A real-valued stochastic process L = {L t , t ≥ 0} defined in a complete probability space (Ω, F , P) is called a Lévy process if the following is the case: 1.
L t+s − L t has the same distribution as L s for all s, t ≥ 0 (stationary increments); 3.
L is stochastically continuous; that is, for all t ≥ 0 and a > 0.

4.
The paths t → L t are right-continuous with left limits (cadlag-continue á droite et limite á gauche).

Normal Inverse Gaussian
An N-vlued stochastic process N = (N t ) t≥0 is called a Poisson process with intensity λ ∈ (0, ∞) if N satisfies (i)-(iii) and if the following is the case.
We can define more general Lévy processes as follows: for a Poisson process (N t ) t≥0 and independently identically distributed Y k , and with k ≥ 1 being the jump size distribution f (.). Such processes are called compound Poisson processes. The characteristic function of C t is given by the following.
If f is given by the probability density function of the normal distribution, then C = (C t ) t≥0 is referred to as a jump diffusion process. The total intensity of Poisson processes with jump sizes in the interval [x, x + dx] is determined by density λ f (dx).
Consider a process X s = (X s t ) t≥0 with X s t = sN , where s is a given a real number is the Poisson process with intensity λ(s). The number s represents the jump size, and intensity λ(s) is the expected number of jumps with size s in the unit time.
Let S = {s i ∈ R : s i = 0, i = 1, 2, · · · } be a subset of jump sizes, λ(s i ) > 0 for all s i ∈ S. Consider a process Y = (Y t ) t≥0 defined by the following.
For process Y, function Q defined for any subset V of S to be Q(V) = ∑ s i ∈V λ(s i ) represents the expected number of jumps with size s i ∈ V in the unit time interval. Thus, the map Q from a subset of R to R + is a measure. Using Q, we can obtain an extended process Y such that the characteristic function of Y t is given by the following: where γ ∈ R, and Q is referred to as a Lévy measure, which captures both the discontinuities rate of occurrence and their overall size. The process Y will be called a pure jump process if If t = 1, then we have a purely non-Gaussian infinitely divisible random variable. Thus, there is a one-to-one correspondence between a purely non-Gaussian infinitely divisible random variable and a pure jump process. If a pure jump process T = (T t ) t≥0 is nondecreasing (that is, T t ≥ 0 almost surely for t > 0, and T t ≥ T l almost surely for l ≤ t), then process T is called a subordinator or intrinsic time process. The Poisson, gamma, and inverse Gaussian processes are nondecreasing; hence, they are subordinators. By considering the inverse Gaussian process T = (T t ) t≥0 as the subordinator of Brownian motion W = (W t ) t≥0 , we obtain the Normal Inverse Gaussian (NIG) process L = (L t ) t≥0 with the following.
We can now present the model to be used in this study.

Fractional Ornstein-Uhlenbeck Lévy Process
The fractional Ornstein-Uhlenbeck Lévy (fOUL) process X = {X t } t≥0 is a solution of the following stochastic differential equation: where α > 0 is a parameter, and L H = (L H ) t≥0 ) is a fractional NIG Lévy process with Hurst parameter H.
By applying Ito's formula to function f (t, X t ) = X t e αt , we obtain the following.
Thus, the solution of the Equation (7) is given by the following.
In practice, as the price is observed at fixed times 0 = t 0 < t 1 < t 2 < · · · < t n = T, with constant ∆t = t k+1 − t k , the solution (8) can be discretized by

Results
The data comes from the Coinmarketcap website (https://coinmarketcap.com/, accessed on 25 February 2022) and covers a period of 5 years starting from 28 April 2013 to 18 July 2019 with a total of 2283 data points (trading days). The evolution of closing prices S t and the corresponding time paths of logarithmic returns r t series over the study period are depicted in Figure 3. As illustrated, the daily log-returns lie between [−0.27, 0.36] over the observed period. The plots confirm the presence of volatility clustering and jumps. Figure 4 depicts a Monte Carlo simulation with 1000 trajectories at equidistant time step t ∈ (0, ∆t, 2∆t, · · · , T) where the time step is set at ∆t = 1 and T = 30. Figure 4 depicts a Monte Carlo simulation with 1000 trajectories at equidistant time step t ∈ (0, ∆t, 2∆t, · · · , T) where the time step is set at ∆t = 1 and T = 30. This technique is useful in the estimation of the conditional value at risk(CVaR) and most notably in pricing derivatives such as over the counter (OTC) options where market prices are not observable. Other applications have been found in quantifying the extent of mispricings in vanilla options. Figure 5 is a graphical display of histograms for an empirical distribution superimposed with the best fitting density distribution, which has been taken at specified time points:

Discussion
On one hand, Figure 2 reveals that the Normal Inverse Gaussian (NIG) distributions fits better to the log-returns of data as compared to the Gaussian distribution, Hyperbolic distribution, and Generalized Hyperbolic distribution. This indicates that Bitcoin data are non-Gaussian, and any prediction model with only a Gaussian process may be misleading. Moreover, a stochastic model driven by a Lévy process is more likely to bring reliable forecast results.
On the other hand, given the self-similarity property of Bitcoin, a fractional Brownian motion model will be convenient to simulate its future prices. Tracing the variation of Hurst parameter H estimated from Bitcoin data over a one year (entire 2018) period (see Figure 1) reveals that H remains below 1/2 during almost the entire second half of 2018. Thus, we choose the mean-reverting fractional stochastic model, such as the fractional Ornstein-Uhlenbeck stochastic driven by NIG processes for the prediction of Bitcoin's future prices. Table 1 is a percentile table that quantifies the risk of the Bitcoin price being above or below a certain level. The economic interpretation of the p-th percentile is the value such that P(X t ≤ x t ) = p α α ∈ (0.25, 0.50, 0.75, 0.90, 0.95, 0.99). The 50th percentile is also called the median, whereas the 25th and 75th percentiles are called quantiles. Using the closing price on 19 July, we can easily forecast future prices at specific time intervals: t ∈ (∆t = 10 day (29July), ∆t = 20 day (8Aug) and ∆t = 30 day (18Aug)).

Conclusions
This paper assess bitcoin exchange rate dynamics by using a time-variant Hurst parameter and forecasts its future value by using a fractional Ornstein-Uhlenbeck driven by a Levy process. The time-dependent Hurst parameter H reveals that Bitcoin displays various memory dependency across the study period: mean-revertingness (H < 1 2 ), randomness (H = 1 2 ), and persistency (H > 1 2 ). While fitting the Gaussian, the Hyperbolic, the Generalized Hyperbolic, and the Normal Inverse Gaussian (NIG) distributions to logreturns of bitcoin, the NIG distribution appears to provide the best fit. This justifies our selection of a stochastic model driven by NIG process. As a result, a Monte Carlo simulation through a fractional Ornstein-Uhlenbeck driven by the NIG process predicts that there is a 95% chance for Bitcoin prices to be less that USD 11,881, USD 12,245, and USD 13,157 on Funding: This research received no external funding .

Data Availability Statement:
The data used for this study can be obtained from the authors upon request or visit Coinmarketcap: https://coinmarketcap.com/, accessed on 25 February 2022.