Bayesian Joint Estimation of the Hurst Parameter and Volatility with Applications to Fractional Option Pricing

1. Introduction

The pricing of financial derivatives remains a central topic in mathematical finance, with the Black–Scholes model representing a fundamental breakthrough in the field [7,25]. By modeling asset prices through geometric Brownian motion under a no-arbitrage framework, the model provides a closed-form solution for European option pricing and has become a cornerstone of modern financial theory. Despite its widespread use, empirical studies have shown that financial time series often exhibit features that are inconsistent with the model assumptions, including long-range dependence, volatility clustering, and heavy-tailed return distributions.

Fractional Brownian motion directly addresses these shortcomings by introducing memory into the stochastic dynamics of asset prices. Fractional Brownian motion (fBm) was formally introduced by Mandelbrot Van Ness [22], who defined it as a generalization of standard Brownian motion with stationary, correlated increments. Its relevance to financial time series—particularly the presence of long-range dependence—was subsequently argued by Mandelbrot [21]. Fractional Brownian motion incorporates a Hurst parameter that characterizes long-range dependence in stochastic processes. In this framework, the Hurst parameter $H\in (0,1)$ governs the degree of dependence in the process, while $\sigma >0$ denotes the volatility. In particular, when the Hurst parameter exceeds one-half, the process exhibits persistence, reflecting the tendency of financial returns to display long memory. This has led to the development of fractional Black–Scholes models, which extend the classical framework by incorporating memory effects and providing a more flexible representation of asset dynamics [6]. More recent developments in financial modeling have also emphasized rough volatility and fractional dynamics as realistic representations of market behavior, particularly at high frequencies [3,14]. Several studies have proposed explicit formulations of such models based on fractional Brownian motion, illustrating how the Hurst parameter influences option pricing and volatility dynamics [10,16,27].

Early studies also investigated option pricing in markets driven by fractional Brownian motion and demonstrated the potential impact of long-range dependence on derivative valuation [28]. Notwithstanding these advances, the problem of jointly estimating the Hurst parameter H and the volatility $\sigma$ within a unified statistical framework has received comparatively little attention. In particular, Bayesian approaches that simultaneously quantify uncertainty in both parameters remain underdeveloped, despite their natural suitability for this inference problem.

While fractional models offer improved realism, their practical implementation introduces significant statistical challenges. In particular, the Hurst index and volatility must be inferred from observed data, and their joint behavior directly influences pricing outcomes. Many existing approaches rely on point estimates or marginal intervals, which may fail to adequately capture the dependence between these quantities. As a result, estimation error can be misrepresented, potentially leading to inaccurate assessments in financial applications. From a practical perspective, these challenges are especially important in decision-making, where even small inaccuracies can lead to meaningful differences in option values and affect hedging strategies, risk management, and portfolio allocation.

Accurate estimation of the Hurst parameter has been a central challenge in the application of fractional models to financial data. The statistical theory of long-range dependence has been extensively developed in both Gaussian and non-Gaussian settings, with applications across telecommunications, hydrology, economics, and finance [4,12,13,33]. These developments have established the Hurst parameter as a key measure of persistence and dependence in stochastic processes.

A wide range of methods has been proposed for estimating long-range dependence, each with different trade-offs between bias, efficiency, and computational complexity. Early work by Beran [4] provided a systematic treatment of statistical methods for long-memory processes, while subsequent studies highlighted the sensitivity of classical estimators to short-range dependence and finite-sample effects [34]. Semi-parametric approaches, such as log-periodogram regression and local Whittle estimation, have been widely used because of their attractive asymptotic properties, although they often require relatively large sample sizes to achieve reliable performance [17,31]. Alternative approaches include wavelet-based estimators, which exploit the scale invariance of fractional processes and have been shown to improve estimation accuracy in certain settings [1,2]. A broader overview of commonly used Hurst exponent estimation procedures and their practical characteristics is provided by Zhang et al. [37].

Despite these developments, much of the existing literature focuses primarily on point estimation and does not fully account for uncertainty in the estimated parameters. This limitation is particularly important in financial applications, where small errors in estimating the Hurst parameter and volatility can propagate through pricing formulas and lead to substantial differences in option valuation. These considerations motivate the use of Bayesian approaches, which provide full posterior distributions and allow uncertainty in model parameters to be carried directly into subsequent inference and pricing decisions.

More recently, Bayesian approaches have been proposed as a flexible framework for estimation in long-memory models, allowing for the incorporation of prior information and providing full posterior distributions [5,8,11,20,35,36]. These methods offer a natural way to quantify uncertainty through posterior distributions. In many applications, however, inference is primarily reported through marginal posterior summaries such as posterior means and credible intervals [9,19,20].

At the same time, maximum likelihood and the method of moments estimators are known to exhibit systematic bias when applied to fractional stochastic differential equation models, particularly for values of the Hurst parameter away from 0.5. This issue has been documented in both parametric and moment-based settings [29,32], where it is shown that both maximum likelihood and method of moments estimators require bias correction for reliable practical use.

Recent work has also considered Bayesian modeling of static and dynamic Hurst parameters under stochastic volatility frameworks, further highlighting the importance of jointly modeling persistence and volatility in long-memory processes [35]. Despite these developments, the joint dependence between key parameters, such as the Hurst parameter and volatility, is often not fully captured, even though this dependence has a direct impact on derived quantities such as option prices. In particular, while bias-corrected estimators improve point estimation, they do not address the propagation of joint parameter uncertainty into financial quantities. This limitation is especially important in financial applications, where interactions between parameters can significantly influence pricing outcomes. These observations highlight the need for methods that explicitly account for joint parameter uncertainty in fractional models.

Bayesian inference provides a coherent framework for addressing these challenges by combining prior information with observed data to produce a full posterior distribution over model parameters [26]. A key advantage of the Bayesian approach is its ability to capture dependence between parameters through joint inference. In this context, joint highest-density regions (HDRs) offer a principled method for summarizing uncertainty in a multidimensional setting, providing a more accurate representation than separate marginal intervals.

Motivated by these considerations, this paper focuses on Bayesian estimation of the Hurst parameter and volatility within a fractional Gaussian noise framework. Using a simulation-based approach, we construct empirical joint highest-density regions for $(H,\sigma )$ and evaluate their repeated-sampling performance. The results reveal strong posterior dependence between parameters and highlight the limitations of marginal interval-based inference.

To assess the practical implications of parameter uncertainty, posterior samples are propagated through the fractional Black–Scholes pricing formula to obtain distributions of option prices. This approach allows for the construction of credible intervals that reflect joint parameter uncertainty, providing a more informative basis for decision-making compared to traditional point estimates.

In this paper, we address this gap by developing a Bayesian framework for the joint estimation of H and $\sigma$ from discretely observed asset price data following a fractional Black–Scholes model. We construct prior distributions informed by the theoretical constraints on each parameter, derive the corresponding posterior, and employ Markov chain Monte Carlo (MCMC) methods for inference. The resulting estimates are then applied to fractional option pricing, where we demonstrate that incorporating posterior uncertainty in H and $\sigma$ leads to materially different option prices compared to approaches that treat these parameters as fixed or known. The remainder of the paper is organized as follows. Section 2 presents the fractional Black–Scholes framework and the Bayesian methodology used to jointly estimate H and $\sigma$, including prior specification, posterior computation, and the construction of joint highest-density regions. Section 3 reports a simulation study evaluating the performance of the proposed estimators, with emphasis on coverage, posterior dependence, and estimation accuracy. Section 4 provides an empirical application in which posterior parameter draws are propagated through the fractional Black–Scholes pricing formula to quantify option pricing uncertainty. Section 5 concludes the paper and discusses limitations and directions for future research.

2. Methodology

This section describes the simulation framework used to study Bayesian estimation of the Hurst parameter H and the scale parameter $\sigma$ under fractional Gaussian noise (fGn). We construct joint highest-density regions (HDRs) for $(H,\sigma )$ and assess repeated-sampling coverage with respect to the true pair $(H_c,{\sigma }_0)$.

2.1. Data-Generating Mechanism

Let $d={(d_1,\dots ,d_n)}^{\top }$ denote an observed fGn sample of length n. For a fixed Hurst parameter $H\in (1/2,1)$, the lag-k autocovariance of unit-variance fGn is

$${\rho }_H\left(k\right)={\displaystyle \frac{1}{2}}\left({|k+1|}^{2H}-{2\left|k\right|}^{2H}+{|k-1|}^{2H}\right),k=0,1,2,\dots .$$

(1)

Since ${\rho }_H\left(0\right)=1$, ${\rho }_H\left(k\right)$ also equals the lag-k autocorrelation, as derived from fractional Brownian motion [4,22].

Using (1), the $n\times n$ Toeplitz correlation matrix $R\left(H\right)$ is defined by

$$R\left(H\right)=\left(\begin{array}{cccc}{\rho }_H\left(0\right)& {\rho }_H\left(1\right)& \dots & {\rho }_H(n-1)\\ {\rho }_H\left(1\right)& {\rho }_H\left(0\right)& \dots & {\rho }_H(n-2)\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }_H(n-1)& {\rho }_H(n-2)& \dots & {\rho }_H\left(0\right)\end{array}\right).$$

(2)

This Toeplitz structure is standard in long-memory Gaussian processes [4] and provides the basis for data generation in our simulation study. For fixed true values $({\mu }_0,{\sigma }_0,H_c)$, the simulated sample is generated from the corresponding Gaussian model using this covariance structure.

$$d\sim N\left({\mu }_0{\mathbf{1}}_n,{\sigma }_0^{2}R\left(H_c\right)\right),$$

(3)

where ${\mathbf{1}}_n$ denotes the n-vector of ones.

To simulate exactly from (3), we compute the Cholesky factorization

$$R\left(H_c\right)=U^{\top }U,$$

(4)

generate $z\sim N(0,I_n)$, and set

$$d={\mu }_0{\mathbf{1}}_n+{\sigma }_0{U}^{\top }z.$$

(5)

This yields an exact Gaussian sample with mean ${\mu }_0{\mathbf{1}}_n$ and covariance ${\sigma }_0^{2}R\left(H_c\right)$. Cholesky factorization provides a numerically stable and computationally efficient approach for simulating Gaussian vectors with positive-definite covariance matrices [15].

2.2. Bayesian Model Specification

Conditional on $(\mu ,\sigma ,H)$, the data are modeled as

$$d\mid \mu ,\sigma ,H\sim N\left(\mu {\mathbf{1}}_n,{\sigma }^{2}R\left(H\right)\right),$$

(6)

with $H\in (H_L,H_U)$, where in our implementation $H_L=0.5$ and $H_U=1$.

The prior distributions are chosen as

$$\begin{array}{cc}\hfill \mu & \sim N({\mu }_0^{*},{\tau }_{\mu }^2),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \sigma & \sim \mathrm{Lognormal}(m_{\sigma },s_{\sigma }^2),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill H& \sim \mathrm{Uniform}(H_L,H_U),\hfill \end{array}$$

(9)

where ${\mu }_0^{*}$, ${\tau }_{\mu }^2$, $m_{\sigma }$, and $s_{\sigma }^2$ are fixed hyperparameters, following standard Bayesian modeling practice [26].

Hence, the posterior density is

$$\pi (\mu ,\sigma ,H\mid d)\propto L(d\mid \mu ,\sigma ,H)\pi \left(\mu \right)\pi \left(\sigma \right)\pi \left(H\right),$$

(10)

where $L(d\mid \mu ,\sigma ,H)$ is the multivariate normal likelihood induced by (6).

2.3. Posterior Computation

Under (6), the log-likelihood is

$$\begin{array}{cc}\hfill logL(d\mid \mu ,\sigma ,H)=-{\displaystyle \frac{1}{2}}[& nlog\left(2\pi \right)+nlog\left({\sigma }^2\right)+log\left|R\left(H\right)\right|\hfill \\ +{\displaystyle \frac{1}{{\sigma }^2}}{(d-\mu {\mathbf{1}}_n)}^{\top }R{\left(H\right)}^{-1}(d-\mu {\mathbf{1}}_n)].\hfill \end{array}$$

(11)

which follows from the multivariate Gaussian likelihood [4]. To compute this efficiently, we use the Cholesky factorization

$$R\left(H\right)=C{\left(H\right)}^{\top }C\left(H\right),$$

(12)

so that

$$log\left|R\left(H\right)\right|=2{\displaystyle \sum _{i=1}^n}log{C}_{ii}\left(H\right),$$

(13)

and the quadratic form is evaluated through triangular solves rather than direct matrix inversion, which improves numerical stability and computational efficiency.

Combining (11) with the priors (7)–(9), the log-posterior is

$$log\pi (\mu ,\sigma ,H\mid d)=logL(d\mid \mu ,\sigma ,H)+log\pi \left(\mu \right)+log\pi \left(\sigma \right)+log\pi \left(H\right)+\mathrm{constant}.$$

(14)

2.4. Bayesian Implementation

To facilitate efficient MCMC sampling and avoid boundary issues, we work with transformed parameters [9]. For the volatility parameter $\sigma >0$, we define $\eta =log\sigma$. For the Hurst parameter $H\in (H_L,H_U)$, we apply a logit transformation

$$z_H=log\left({\displaystyle \frac{H-H_L}{H_U-H_L-(H-H_L)}}\right),$$

(15)

which maps H to the real line. The corresponding Jacobian adjustments are included in the posterior target density.

Posterior inference is carried out using a Metropolis-within-Gibbs scheme with random-walk proposals [9]. To improve sampling efficiency, we incorporate a joint block update for $(H,\sigma )$, motivated by their strong posterior dependence. The proposal covariance is estimated from a preliminary pilot run and used to construct a bivariate Gaussian random-walk proposal aligned with the posterior geometry. This substantially improves mixing in regions of high posterior curvature.

2.5. Joint Highest-Density Region Coverage

Following the posterior sampling step, we evaluate the transformed posterior log-density at each draw as

$$h^{\left(m\right)}=log\tilde{\pi }({\mu }^{\left(m\right)},{\eta }^{\left(m\right)},z_H^{\left(m\right)}\mid d),$$

(16)

which is used to construct the joint highest-density regions.

The transformed posterior log-density is evaluated at each retained posterior draw $m=1,\dots ,M$, producing a set of posterior heights ${\left\{h^{\left(m\right)}\right\}}_{m=1}^M$. These values are ranked from largest to smallest. For a target probability level $\alpha \in \{0.90,0.95\}$, we retain the top $\alpha M$ posterior draws following the highest-density region (HDR) construction of [18]. The retained subset defines an empirical joint HDR for $(H,\sigma )$:

$${\mathcal{R}}_{\alpha }=\left\{(H^{\left(m\right)},{\sigma }^{\left(m\right)}):log\tilde{\pi }({\mu }^{\left(m\right)},{\eta }^{\left(m\right)},z_H^{\left(m\right)}\mid d)\ge c_{\alpha }\right\},$$

(17)

where $c_{\alpha }$ is the empirical cutoff corresponding to the top $\alpha$ fraction of posterior heights. This construction is preferable to separate marginal intervals for H and $\sigma$, since the posterior draws exhibit substantial dependence and lie on a curved ridge rather than in an approximately rectangular region.

To evaluate repeated-sampling performance, we repeat the above procedure over B simulated datasets. For each replication, we check whether the true parameter pair $(H_c,{\sigma }_0)$ lies inside the estimated joint HDR:

$$I_b=\mathbf{1}\left\{(H_c,{\sigma }_0)\in {\mathcal{R}}_{\alpha }^{\left(b\right)}\right\},b=1,\dots ,B.$$

(18)

The empirical coverage is then

$${\widehat{\mathrm{Cov}}}_{\alpha }={\displaystyle \frac{1}{B}}{\displaystyle \sum _{b=1}^B}{I}_b.$$

(19)

In addition, we summarize posterior medians, biases and posterior correlations between H and $\sigma$.

2.6. Posterior Option Pricing

To connect parameter uncertainty with option pricing, each retained posterior draw $(H^{\left(m\right)},{\sigma }^{\left(m\right)})$ is mapped through the fractional Black–Scholes pricing framework [6,16,27]. The pricing equation is used here as a sensitivity-based mechanism for propagating posterior uncertainty in $(H,\sigma )$ into option values. Let

$${\lambda }_H=2H{T}^{2H-1}.$$

(20)

Then

$$\begin{array}{cc}\hfill d_1& ={\displaystyle \frac{log(S_0/K)+\left(r+\frac{{\lambda }_H}{2}{\sigma }^2\right)T}{\sigma \sqrt{{\lambda }_{H}T}}},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill d_2& ={\displaystyle \frac{log(S_0/K)+\left(r-\frac{{\lambda }_H}{2}{\sigma }^2\right)T}{\sigma \sqrt{{\lambda }_{H}T}}},\hfill \end{array}$$

(22)

and the call price is

$$C(H,\sigma )=S_0\Phi \left(d_1\right)-K{e}^{-rT}\Phi \left(d_2\right),$$

(23)

where $\Phi (·)$ denotes the standard normal distribution function.

Applying (23) to every posterior draw produces a posterior sample of option prices,

$$C^{\left(m\right)}=C(H^{\left(m\right)},{\sigma }^{\left(m\right)}),m=1,\dots ,M,$$

(24)

from which posterior medians and credible intervals for the option price are obtained.

3. Simulation Study

This section evaluates the repeated-sampling performance of the proposed Bayesian procedure under fractional Gaussian noise. For each replication, we (i) generate an exact Gaussian fGn sample from the true covariance structure, (ii) estimate $(\mu ,\sigma ,H)$ via Bayesian MCMC under the multivariate normal model, (iii) construct an empirical joint HDR by ranking posterior draws according to posterior height, (iv) record whether the true pair $(H_c,{\sigma }_0)$ falls inside the estimated HDR, and (v) propagate posterior draws through the fractional Black–Scholes formula to obtain a posterior distribution for the call price.

Posterior inference was based on 3,000 MCMC iterations following an initial pilot run of 1,000 iterations. The first 500 iterations of the main chain were discarded as burn-in, leaving 2,500 posterior draws for inference. For each design point, 100 simulated datasets were generated. Computation times are based on a MacBook Air with an Apple M4 processor and 24GB of RAM. Approximate runtimes for the full simulation grid under the sequential implementation were 6.6 minutes for $n=50$, 16.9 minutes for $n=100$, and 4.5 hours for $n=300$, demonstrating the substantial increase in computational cost as sample size grows. The largest setting, $n=500$, required approximately 72 hours under the sequential implementation. To improve computational efficiency, a parallelized version of the simulation procedure was also implemented, reducing the runtime for $n=500$ to approximately four hours while producing qualitatively consistent results.

Across all simulation settings, estimation accuracy improved as the sample size increased, with narrower posterior regions and reduced bias observed for larger values of n. However, the overall pricing behavior and uncertainty structure remained qualitatively similar for $n\ge 300$. This suggests that moderate sample sizes already provide stable inference for the proposed Bayesian framework while avoiding the substantially higher computational cost associated with larger datasets. Similar computational trade-offs have been noted in recent studies of Bayesian Hurst exponent estimation for long-memory models [24].

Figure 1 and Figure 2 illustrate the joint posterior draws and corresponding 95% HDR clouds for the smallest and largest sample sizes considered. These figures demonstrate the substantial posterior dependence between H and $\sigma$ and show the improvement in posterior concentration as sample size increases. Corresponding results for the intermediate sample sizes $n=100$ and $n=300$ are provided in Appendix A. Figure 3a–d summarize empirical coverage and posterior dependence across sample sizes on a common scale. Coverage remains close to the nominal 0.95 level across most settings, while $\mathrm{corr}(H,\sigma )$ generally increases with $H_c$, indicating that separate marginal intervals may fail to capture important joint posterior structure.

It is well known that fractional Brownian motion models may violate the semimartingale structure underlying classical arbitrage-free pricing theory [6,30]. Consequently, the fractional Black–Scholes framework is used here primarily as a sensitivity-based tool for studying how posterior uncertainty in $(H,\sigma )$ propagates into option values, rather than as a complete market equilibrium model under all fractional specifications. Nevertheless, recent developments in rough volatility modeling further supports the usefulness of fractional dynamics as flexible representations of financial dependence structures [3,14].

To quantify the implications for option valuation, Figure 4 reports posterior medians and 95% credible intervals for the call price obtained by mapping posterior draws $(H,\sigma )$ through the fractional Black–Scholes pricing equation. The results show a clear decline in the posterior median call price as $H_c$ increases, with the decrease becoming especially pronounced when $H_c$ approaches one. This behavior indicates that, under the chosen option-pricing parameters, stronger persistence can substantially reduce the recommended posterior median call price.

The uncertainty bands also widen noticeably for larger values of $H_c$, particularly near $H_c=0.90$ and $H_c=0.95$. In this region, the posterior price distributions become more dispersed and appear left-skewed, with lower posterior medians observed for smaller sample sizes. This suggests that when the true dependence is very strong, limited sample sizes may lead to greater uncertainty and more conservative pricing recommendations.

From a practical perspective, increasing the sample size improves posterior concentration, but the gain must be balanced against computational cost. The difference between $n=300$ and $n=500$ is relatively small compared with the substantial increase in runtime required for $n=500$. In contrast, $n=100$ produces wider intervals, but these intervals often contain the corresponding intervals for larger sample sizes. Therefore, $n=100$ may be sufficient for exploratory analysis or preliminary pricing decisions, while $n=300$ provides a more reliable compromise between statistical precision and computational feasibility. The use of $n=500$ is recommended only when very high precision is required or when parallel computation is available. Table 1 provides the corresponding numerical summaries, including the 95% interval $[q_{0.025},q_{0.975}]$ and the 90% band $[q_{0.05},q_{0.95}]$, which may be interpreted as a conservative bid–ask range.

Figure 5 displays the posterior lower, median, and upper quantile curves of the option price as functions of the true Hurst parameter $H_c$ for each sample size. The lower and upper curves form a 95% credible interval. The shrinkage of the gap between $C_{0.025}$ and $C_{0.975}$ as n increases illustrates posterior concentration in the joint parameter space. The median call price decreases with $H_c$ under the chosen option parameters, while uncertainty is substantially larger for small samples and tightens as n increases. These quantile curves also facilitate decision-oriented summaries, such as one-sided credible bounds, which may be preferable when overpricing and underpricing have asymmetric consequences.

4. Empirical Application

This section illustrates the practical performance of the proposed Bayesian framework using real financial data from energy markets. We consider daily log returns of WTI crude oil and natural gas futures, which are known for their pronounced volatility and relevance in derivative pricing applications.

To investigate the effect of market conditions, the oil data are divided into two distinct periods representing different volatility regimes. The first period corresponds to a high-volatility environment, while the second reflects relatively stable market behavior. Natural gas is included as an additional example of a persistently volatile asset, allowing us to examine the robustness of the proposed approach across different settings.

For each dataset, we estimate the Hurst parameter and volatility using the Bayesian methodology described in the previous section. These estimates are then incorporated into the fractional Black–Scholes framework to obtain posterior distributions of option prices. This approach allows us to assess not only point estimates but also the uncertainty associated with pricing decisions.

The empirical analysis is designed to distinguish between volatility-driven effects and genuine long-range dependence, and to quantify their impact on option pricing. This distinction is particularly important, as traditional methods may confound volatility with persistence, leading to misleading inference.

4.1. WTI Crude Oil: Regime Comparison

We begin with WTI crude oil, which exhibits distinct regimes of market behavior. The period 2020–2022 is characterized by elevated volatility, while 2023–2025 reflects comparatively more stable market conditions. This provides a natural setting for examining how changes in market volatility influence posterior inference and option pricing uncertainty.

Table 2 shows that the earlier period exhibits substantially higher volatility, while the average returns remain close to zero in both regimes. This suggests that differences between the two periods are driven primarily by changes in market variability rather than systematic shifts in mean returns.

For pricing, we consider a European call option with $S_0=100$, $K=98$, $T=90/365$, and $r=0.05$, holding these contract inputs fixed when comparing classical Black–Scholes and fractional Black–Scholes prices.

Table 3 reports posterior summaries for the Hurst parameter and volatility. Across both periods, the estimated Hurst parameter remains close to 0.5, indicating weak long-range dependence in returns. In contrast, posterior volatility estimates differ substantially between periods, with the 2020–2022 regime exhibiting considerably higher uncertainty. The wider posterior pricing interval during this period, with a width 1.63 compared to 0.59 in the later regime, indicates that elevated volatility translates directly into greater uncertainty in option valuation.

Figure 6 further illustrates that the posterior distribution of the Hurst parameter remains relatively stable across periods, whereas volatility exhibits substantial variation between market regimes. The posterior clouds also reveal noticeable dependence between H and $\sigma$, emphasizing the importance of joint uncertainty quantification in the Bayesian framework.

The fractional Black–Scholes formula is used here as a sensitivity-based pricing framework to study how posterior uncertainty in $(H,\sigma )$ propagates into option values, rather than as a claim of complete market arbitrage-free dynamics under all fractional specifications.

Figure 7 shows that the high-volatility period produces a substantially wider posterior distribution of option prices, reflecting increased pricing uncertainty. In contrast, the more stable 2023–2025 regime yields a considerably tighter distribution. Although the posterior median prices remain reasonably close to the corresponding classical Black–Scholes prices, the range of plausible option values differs substantially across market regimes.

Overall, these findings suggest that differences in market regimes affect option pricing uncertainty primarily through changes in volatility rather than through strong long-range dependence in returns.

4.2. Natural Gas: High-Volatility Market

We next consider natural gas as an example of a persistently high-volatility commodity market. In contrast to the regime-based behavior observed in oil, natural gas exhibits consistently high volatility over the entire sample period.

Table 4. Summary statistics for the natural gas return series used in the empirical analysis.

Asset	Ticker	Period	n	Mean Return	Annualized Volatility
Natural Gas Futures	NG=F	2022–2024	711	$-0.00041$	0.832

Table 4. Summary statistics for the natural gas return series used in the empirical analysis.

Asset	Ticker	Period	n	Mean Return	Annualized Volatility
Natural Gas Futures	NG=F	2022–2024	711	$-0.00041$	0.832

For pricing, we consider a European call with $S_0=100$, $K=98$, $T=90/365$, and $r=0.05$, holding these contract inputs fixed when comparing BSM and fBSM prices.

Table 5 reports posterior summaries of the Hurst parameter and volatility. The estimated Hurst parameter remains close to 0.5, indicating weak long-range dependence in returns. In contrast, the volatility is relatively high, reflecting the pronounced variability of the natural gas market. The corresponding posterior pricing interval width of 1.62 further indicates substantial uncertainty in option valuation under highly volatile market conditions.

The joint posterior distribution in Figure 8 reveals a concentration of H around 0.5, while the volatility parameter displays considerable variability. This pattern suggests that the observed dynamics are primarily driven by volatility rather than persistence. The substantial posterior dispersion in $\sigma$ also contributes directly to the wider range of plausible option prices observed under the fractional pricing framework.

Figure 9 illustrates the posterior distribution of option prices. The posterior distribution is relatively wide, reflecting substantial uncertainty in pricing outcomes. Although the classical Black–Scholes price lies within the posterior pricing interval, it represents only a single point estimate and therefore does not capture the full uncertainty associated with the joint estimation of $(H,\sigma )$.

Overall, these results indicate that, in highly volatile markets, parameter uncertainty can translate into meaningful variation in option values, even in the absence of strong long-range dependence.

5. Discussion and Conclusions

This paper develops a Bayesian framework for estimating the Hurst parameter and volatility in the fractional Black–Scholes model and examines their impact on option pricing in energy markets. The empirical results provide several important insights.

First, across all datasets considered, the estimated Hurst parameter remains close to 0.5, suggesting weak evidence of long-range dependence in returns despite elevated market volatility. This finding is consistent across different assets and market conditions, indicating that persistence plays a limited role in describing return dynamics. Although the estimated Hurst parameter for returns remains close to 0.5, financial markets may still exhibit persistence in volatility dynamics, which is consistent with the broader literature on volatility clustering and rough volatility models. This suggests that fractional dynamics may remain useful for uncertainty quantification even when strong long-range dependence is not directly observed in return series.

However, the posterior distributions also demonstrate non-negligible uncertainty surrounding these estimates. Consequently, the analysis does not support imposing the classical Brownian motion assumption $H=0.5$ with certainty. From a risk management perspective, it is preferable to allow for the possibility that the true persistence structure may deviate from standard Brownian motion, even if the posterior center lies near 0.5. The Bayesian framework, therefore, provides a natural mechanism for incorporating parameter uncertainty directly into inference and option pricing, thereby avoiding overconfidence in a single fixed persistence assumption. In this sense, the fractional model serves not only as a tool for detecting strong long-range dependence, but also as a framework for quantifying uncertainty about persistence in financial markets.

An important contribution of the proposed framework is not necessarily the detection of strong long-range dependence itself, but rather the quantification of uncertainty surrounding persistence. Although the posterior distributions of the Hurst parameter were centered near 0.5 in the empirical applications, the results do not support treating the classical Brownian motion assumption H = 0.5 as known with certainty. Instead, the Bayesian framework provides a probabilistic assessment of persistence by propagating uncertainty in the Hurst parameter directly into option pricing inference. Consequently, the methodology remains informative even when strong long-range dependence is not conclusively supported by the data, since uncertainty regarding persistence can still influence financial decision-making and risk assessment.

Second, volatility emerges as the primary driver of market behavior. In the case of WTI crude oil, the comparison of two distinct periods reveals that differences in market regimes are characterized by substantial changes in volatility, while the Hurst parameter remains relatively stable. Similarly, natural gas exhibits consistently high volatility, leading to greater uncertainty in parameter estimates.

Third, these differences in volatility translate directly into differences in option pricing uncertainty. While the median fractional Black–Scholes prices remain close to the classical Black–Scholes values, the posterior distributions reveal that the range of possible option prices can vary considerably, particularly in high-volatility environments. This highlights that classical pricing models may underestimate the uncertainty associated with option valuation.

From a methodological perspective, the results emphasize the importance of jointly modeling dependence and volatility. The Bayesian approach provides a flexible framework that captures parameter uncertainty and avoids over-reliance on point estimates. These findings are consistent with recent studies emphasizing the advantages of Bayesian approaches for uncertainty quantification in long-memory estimation problems [23,24]. In contrast, methods based solely on moment conditions may be more sensitive to volatility-driven variation, potentially leading to misleading interpretations of persistence.

Overall, the findings suggest that, in financial markets, uncertainty in volatility plays a more significant role than long-range dependence in determining option prices. This has important implications for risk management, as it indicates that accounting for parameter uncertainty is essential for accurate pricing and hedging decisions.

This study has several limitations. The return series are modeled using a Gaussian fractional Gaussian noise likelihood with a constant volatility parameter within each sample period. Consequently, the framework does not explicitly account for stochastic volatility, jumps, heavy tails, or structural breaks that may arise in financial markets. In addition, the empirical analysis is based on daily data, which may mask microstructure effects and intraday dependence patterns.

Future research may extend this framework to incorporate stochastic volatility models, alternative prior structures, and heavy-tailed likelihoods. Additional directions include applications to high-frequency financial data and the integration of fractional dynamics with more complex market mechanisms such as jumps, leverage effects, and regime-switching behavior. Overall, the proposed Bayesian framework demonstrates that uncertainty quantification in fractional financial models can provide valuable insight into the stability and reliability of option pricing decisions under complex market conditions.