The foundation of modern derivative valuation is the Black–Scholes–Merton (BSM) model that provides an analytical, closed-form solution assuming that the asset price follows a Geometric Brownian Motion (GBM) with constant volatility and risk-free rate [7, 8]. The failure of this assumption to hold in real markets led to the development of more sophisticated methods.

Monte Carlo (MC) simulation became essential for models without closed-form solutions, such as those involving path-dependent features or complex stochastic processes. However, standard MC still requires specifying the underlying stochastic process, often a GBM or a jump-diffusion process, which remains a main source of model risk.

A more advanced framework, Stochastic Volatility Models, try to address the persistent volatility smile by allowing the volatility parameter to evolve randomly over time [9]. While these models significantly improve pricing accuracy compared to BSM, they introduce computational complexity, via complex multi-parameter calibration and still rely on the accurate specification of the volatility dynamics [10].

Generative Adversarial Networks (GANs), which was introduced by Goodfellow et al. (2014), are a class of deep learning models designed to generate synthetic data that closely resembles a target distribution [11, 12, 13]. A GAN consists of two neural networks that are simultaneously trained: a Generator (G) and a Discriminator (D).

The generator maps random noise z ξ Pz (z) to synthetic data G(z), while the discriminator attempts to distinguish between real data x ξPdata(x) and generated data G(z).

The training involves a minimax game governed by the objective function defined as:

min max V (D, G) = Ex∼Pdata x[log ↽↼ D(x)] + Ez∼Pzz ↼ [log(1 ↽ − D(G(z)))] (1)

This adversarial setup encourages the generator to produce samples that are indistinguishable from real data, while the discriminator improves its classification accuracy. The equilibrium of this game is reached when D(x) = 0.5 for all inputs, indicating perfect confusion between real and synthetic data.

However, standard GANs often suffer from mode collapse, unstable gradients, and poor convergence, especially when applied to financial time series, which exhibit nonstationarity, volatility clustering, and heavy tails. To address these challenges, Wasserstein GANs (WGANs) were proposed [14], replacing the Jensen–Shannon divergence with the Wasserstein-1 distance (also known as Earth Mover’s Distance), which is defined as:

W(Pr, Pg) = inf E x,y ∼ ⊊ γ ↽ [∥x − y∥] (2)

γ∈⋄⇐ Pr,Pg↽

The WGAN objective functions is modified to:

min max Ex∼Pdata [D(x)] − Ez∼Pz [D(G(z))] (3)

where D represents the set of 1-Lipschitz functions, a constraint enforced via weight clipping. This formulation provides smoother gradients and better convergence, making WGANs particularly suitable for generating realistic financial trajectories. In the context of this study, WGANs are used to simulate synthetic NIFTY 50 paths that preserve empirical properties such as autocorrelation, volatility regimes, and tail behavior. These paths serve as inputs to downstream pricing models for vanilla options, enhancing robustness and realism in derivative valuation.

The application of machine learning, such as deep learning in quantitative finance has grown rapidly, spanning across areas from algorithmic trading to derivative pricing [15]. Within this domain, GANs have demonstrated promise for synthesizing realistic market data (e.g., equities, commodities, and derivatives), allowing for robust back-testing and simulation of extreme events [16]. Some of the research also shown that Time-series Generative Adversarial Networks (TimeGAN) can generate multivariate, sequence-dependent data that preserves both the cross-sectional and temporal characteristics of the original data [17]. This success forms the necessary foundation for the current study, highlighting the feasibility of using GANs to model volatility and correlation structures crucial for pricing complex derivatives.

The literature establishes that while traditional models are interpretable, they are often statistically limited, and pure machine-learning models, though powerful, can lack economically meaningful constraints. A hybrid approach is therefore desirable: one that preserves financial interpretability while leveraging the distribution-learning strengths of modern generative models [18], [19].

Classical models (Black–Scholes, Merton) assume log-normal returns and constant volatility, assumptions frequently violated in practice. Stochastic volatility models (e.g., Heston) and jump-diffusion models attempt to address these shortcomings but remain parametric and can struggle to capture complex empirical features of returns. Recent advances in machine learning, particularly GANs, have shown promise in generating realistic financial time series [20].

In this work we follow a hybrid, practical strategy: we use a GAN to learn realistic return-path distributions and integrate those synthetic scenarios into a Monte Carlo pricing framework for vanilla options. To maintain financial consistency we apply a small set of economically informed corrections — notably an expanding (cumulative) training window and a drift correction to align generated paths with the risk-neutral drift — rather than embedding an explicit parametric model (like GARCH [21]) into the GAN loss. This preserves interpretability (we still benchmark against Black–Scholes and Heston) while benefiting from the GAN’s flexible distributional learning.

To model realistic return dynamics for vanilla option pricing on the NIFTY 50 index, we employ a Generative Adversarial Network (GAN) trained on historical log-returns using an expanding window approach. Each yearly model is trained on all data available up to (but excluding) that year, ensuring the generation process reflects only past information—consistent with a forward-looking pricing setup.

The GAN consists of two neural networks [22]:

The generator maps a latent noise vector sampled from a standard normal distribution to a synthetic sequence of log-returns. The architecture consists of a fully connected layer with 64 neurons and ReLU activation, followed by an output layer with T neurons (no activation), where T varies by year between 83 and 87 trading days. The output is reshaped to match the temporal structure of real sequences. Mathematically, the generator function is defined as G(z; θG) = W3∆ ReLU(W2∆ z + b2) + b3 (4) where there are 32 latent dimensions, and the learnable parameters are represented by weights and biases.

The discriminator (or critic) evaluates input sequences and outputs a scalar score rather than a probability. The architecture employs 1D convolution with 32 filters (kernel size 3, same padding, ReLU activation) to capture temporal dependencies, followed by flattening and two fully connected layers (32 neurons with ReLU, then a single output neuron with no activation). The critic function is expressed as D(x; θD) = W5∆ ReLU(W4∆ Flatten(Conv1D(x)) + b4) + b5 (5) where all learnable parameters include convolutional filters, weights, and biases.

The WGAN training objective is formulated as previously described:

min max Exξ Pdata [D(x)] Γ Ezξ Pz [D(G(z))] (6)

where the real data distribution, latent noise distribution (standard normal), and the set of 1-Lipschitz functions are represented respectively. This is equivalent to minimizing the Wasserstein-1 distance:

W1(Pdata, PG) = supk Exξ Pdata [D(x)] Γ Exξ PG [D(x)] (7)

DkLκ 1

The discriminator loss is defined as

LD = Ezξ Pz [D(G(z))] Γ Exξ Pdata [D(x)] (8) which the discriminator aims to minimize, while the generator loss is:

LG = Γ Ezξ Pz [D(G(z))] (9) which the generator minimizes to fool the critic into assigning high scores to synthetic data.

For the WGAN objective to approximate the Wasserstein distance, the discriminator must lie within the space of 1-Lipschitz functions. In our implementation, this constraint is enforced via weight clipping: after each discriminator update, all parameters are clipped to the interval [Γ 0.01, 0.01]:

θDψ clip(θD, Γ 0.01, 0.01) (10) This ensures that the gradient of the discriminator with respect to its input remains bounded.

Following the WGAN protocol, the discriminator is updated twice per generator update to maintain a strong critic throughout training. For each batch, we first perform two discriminator updates by sampling real sequences and generated sequences, computing the discriminator loss, updating parameters via gradient descent with Adam optimizer (learning rate 0.0001), and clipping weights. We then perform one generator update by sampling new noise, computing the generator loss, and updating parameters with the same optimizer.

While our implementation uses weight clipping for computational simplicity, we acknowledge that gradient penalty (WGAN-GP) is generally preferred in the literature [23]. Gradient penalty offers smoother gradients, better convergence, and improved sample quality, particularly for high-dimensional distributions. However, it introduces additional computational overhead due to second-order derivatives. For our application—generating relatively short sequences (T ≈ 80–90) on CPU—weight clipping provided adequate training stability.

Both networks are optimized using the Adam optimizer. The discriminator is updated twice per generator step to ensure a strong critic during adversarial training.

The training follows a cumulative (expanding window) structure:

Training Data for Year t = {rτ | τ < t} (11)

where rτ are historical log-returns. The data are normalized before training to ensure numerical stability. Each GAN is trained for 400 epochs per year with a batch size of 32. Further details on model hyperparameters are outlined in the reproducibility appendix at the end of this document.

A fundamental distinction in derivative pricing is the separation between the real-world probability measure (P) and the risk-neutral probability measure (Q) [24]. The GAN, trained on historical NIFTY 50 log-returns, inherently learns the empirical distribution under P, which reflects the actual return dynamics observed in the market. However, option pricing theory requires valuation under Q, where all assets grow at the risk-free rate when discounted appropriately.

Under the risk-neutral measure, the discounted asset price must be a martingale. For a non-dividend-paying asset with spot price S0, risk-free rate r, and maturity T, this condition requires:

EQ[ST ] = S0erT . (12)

Equivalently, for log-returns, if we denote the cumulative log-return over the period as the sum of individual returns, the risk-neutral expectation must satisfy:

" T # EQ r ˜i = rT. (13)

This ensures that the expected cumulative growth of the asset under Q equals the risk-free rate compounded over the time horizon, eliminating any risk premium present in the real-world measure.

The synthetic return paths {r ˜i}T i=1 generated by the GAN reflect P-dynamics and will generally exhibit a cumulative mean µGAN that differs from the risk-neutral drift rT. To transform these returns into the Q-measure, we apply a deterministic drift correction.

After training, 10,000 synthetic return paths are generated per year. For each year t with maturity T (in years) and risk-free rate r, we first compute the mean cumulative log-return across all generated paths:

" T # µGAN = E r ˜i , (14)

where the expectation is taken over the N = 10, 000 simulated paths.

PT To enforce the martingale condition i=1 r ˜i = rT, we apply a uniform adjustment to each time step. The drift-corrected returns are obtained by

X T rT − µGAN r ˜i = rT ⇒ r ˜i = ri + . (15) i=1 T This correction distributes the difference between the target risk-neutral drift (rT) and the empirically observed GAN drift (µGAN) uniformly across all T time steps. Each individual return is shifted by the same constant (rT µGAN)/T, ensuring that the cumulative return of every path sums to exactly rT.

This correction method satisfies several desirable properties:

PT  Martingale consistency: By construction, the adjusted cumulative returns satisfy i=1 r ˜i = rT for all paths, ensuring the martingale condition is respected.

Distribution preservation: The correction applies a uniform additive shift to all returns, preserving the relative structure of the GAN-learned distribution. Higher-order moments—including volatility, skewness, and kurtosis—remain governed by the generative model’s learned dynamics, while only the location (mean) is adjusted.

Path-wise enforcement: Unlike methods that only adjust the ensemble mean, our approach ensures each individual synthetic path exhibits the correct risk-neutral drift, providing consistency at both the path and distribution levels.

Computational efficiency: The correction is applied post-generation as a simple additive constant, requiring no retraining of the GAN.

Alternative approaches, such as directly training the GAN on risk-neutral returns or embedding the drift constraint into the loss function, were considered but would compromise the GAN’s ability to learn the rich empirical features of the return distribution; particularly volatility clustering and tail behavior. Our approach maintains the flexibility of the generative model while ensuring theoretical consistency with no-arbitrage pricing.

This diagram below illustrates the architecture and training flow of the WGAN used to generate synthetic log-returns. The generator transforms latent noise into time-series outputs via dense and reshape layers. The discriminator evaluates both real and generated returns using convolutional and dense layers to produce a critic score. The training loop computes losses, updates both networks, and applies weight clipping to enforce Lipschitz continuity. Arrows indicate the adversarial feedback between components, forming a closed optimization cycle.

The drift-corrected returns are then used to construct synthetic price paths for Monte Carlo option pricing, as described below:

St = S0 exp r ˜i (16) and fed into a Monte Carlo (MC) pricing module for vanilla European options. This hybrid GAN-MC framework allows the model to generate realistic, non-Gaussian market scenarios while maintaining consistency with financial pricing constraints. Finally, the expected value under the risk-neutral measure (Q) of the discounted payoffs is used for pricing, as described below:

C = e rT EQ [max(ST K, 0)] (17)

To benchmark the performance of the proposed GAN-based scenario generation framework, we compare it against two classical option pricing models widely used in financial engineering: the Black–Scholes model and the Heston stochastic volatility model.

The Black–Scholes (BS) model assumes that the underlying asset price follows a geometric Brownian motion with constant drift and volatility [25]:

dSt = rSt dt + σSt dWt (18)

where r is the continuously compounded risk-free rate, σ is the volatility, and Wt is a standard Wiener process.

The analytical solution for the European call option price is:

C = S0N(d1) Ke rT N(d2) (19)

ln(S0/K) + (r + 1 σ2)T p d1 = p 2 , d2 = d1 σ T (20)

where N( ) is the cumulative distribution function (CDF) of the standard normal distribution.

In our implementation, the volatility parameter σ is both estimated historically (theoretical prices) and inferred from the market—that is, the implied volatility that equates the BS theoretical price to the observed market option price (calibrated prices).

The Heston model extends the BS framework by introducing stochastic volatility. The asset price and variance evolve under correlated stochastic processes:

dSt = rSt dt + vtSt dWt p v (21) dvt = κ(θ vt) dt + σv vt dWt

where vt is the instantaneous variance, κ is the mean-reversion rate, θ is the long-run mean of variance, σv is the volatility of volatility, and WtS and Wtv are correlated Wiener processes with correlation coefficient ρ.

Unlike BS, the Heston model captures volatility clustering, leverage effects, and skewness in the return distribu-tion—phenomena commonly observed in equity markets.

Model calibration involves minimizing the error between theoretical option prices (obtained via the semi-analytical Heston pricing formula) and observed market prices across multiple strikes and maturities, yielding an optimal parameter set:

fκ, θ, σv, ρ, ν0g (22)

Once calibrated, simulated Heston paths are generated and used to compute model-based distributions of returns, forming a benchmark for comparison against GAN-generated distributions. The distribution process is replicated for theoretical and calibrated BS models, though pricing is carried out via analytical solutions.

The Heston model parameters were calibrated to market option prices using a two-stage optimization procedure. First, a global search was performed using differential evolution [26] with 40 iterations and population size 10 to explore the parameter space and avoid local minima. The algorithm minimized the mean squared error between model-predicted and market-observed option prices, subject to bounds ensuring positivity constraints and correlation bounds. Soft penalties were applied for violations of the Feller condition to maintain numerical stability. The best solution from differential evolution was then refined using local L-BFGS-B optimization [27] with a maximum of 500 iterations to achieve precise convergence. This hybrid approach balances global exploration with local refinement.

This study utilizes two primary datasets:

(1) NIFTY 50 Index (2015–2025) Daily closing prices were obtained from Yahoo Finance (ticker:ˆNSEI) [28] and used to compute log returns for training the Generative Adversarial Network (GAN). The same dataset provides the underlying spot price S0 required for option pricing. Returns are normalized and structured into rolling sequences corresponding to the empirical trading lengths for each year (ranging from 83 to 87 days).

(2) NSE Option Market Data Option prices for the NIFTY 50 were extracted from the National Stock Exchange (NSE) [29]. NIFTY 50 index options traded on the NSE are European-style contracts, eliminating concerns about early exercise that would complicate pricing under standard frameworks. This characteristic ensures compatibility with the Black–Scholes, Heston, and Monte Carlo pricing approaches employed in this study.

Model performance was assessed across three complementary dimensions:

(1) Pricing Accuracy Mean Squared Error (MSE) between model-derived and market-observed option prices serves as the primary metric for evaluating pricing performance [31]:

1 X n 2 MSE = Pimodel Pimarket (24) n i=1 where Pimodel and Pimarket are theoretical and actual option prices, respectively.

(2) Distributional Fidelity Beyond price-level accuracy, the statistical similarity between model-simulated and empirically realized return distributions was analyzed. For each year, the first four moments—mean, standard deviation, skewness, and kurtosis—were computed for the synthetic distributions generated by each model (Black– Scholes, Heston, GAN) and compared against market distributions. The total absolute deviation across these moments provides a Match Score, where lower values indicate closer alignment with the true data-generating process.

(3) Computational Efficiency Training and calibration times were recorded to assess the practical feasibility of each model:

Black–Scholes: Executed instantaneously due to its closed-form nature.

Heston Model: Requires numerical calibration of five parameters via iterative optimization (e.g., L-BFGS-B), resulting in moderate computation time.

GAN: Involves deep neural network training across multiple epochs, leading to significantly higher computational cost.

The comparative analysis of option pricing models—Black–Scholes (BS) Theoretical, BS Calibrated, Heston Calibrated, and the Hybrid GAN-Monte Carlo (MC) Framework—yielded several critical observations regarding model performance, computational efficiency, and adherence to market distributions. Visually, all four modeling approaches showed a rough fit to the market-based option prices across the range of strikes. However, quantifiable performance using Mean Squared Error (MSE) revealed significant differences (see Figure 2, Figure 3, Figure 4, and Figure 5).

Calibrated Black–Scholes, using implied volatility, was repeatedly the best-performing model, consistently achieving the lowest MSE across most years, which highlights the powerful predictive quality of implied volatility derived from the market itself. The Heston Model offered competitive results, often placing second; it is plausible that a more finely tuned calibration process, particularly around strikes closer together, could have further improved its accuracy heston. Conversely, the BS Theoretical model, relying solely on historical volatility, performed the worst due to its fundamental inability to capture the volatility smile—the variation of implied volatility across strikes. Notably, the GAN-MC Framework demonstrated strong, comparable performance, achieving the lowest MSE in 2022, with error rates mostly commensurate with the other calibrated models, proving its viability as a data-driven alternative.

NIFTY 50 Historical Data: Yahoo Finance (ticker: ˆNSEI), daily closing prices from January 1, 2015 to August 31, 2025 [28].

Option Market Data: National Stock Exchange of India (NSE), European-style call options with approximately 3-month maturities [29].

Risk-Free Rate: Indian 3-Month Treasury Bill Yield [30].

The complete implementation code, including GAN training, risk-neutral drift correction, and Monte Carlo pricing modules, is available upon request from the corresponding author.

[20] Wiese, M., Knobloch, R., Korn, R., & Kretschmer, P. (2020). Quant GANs: Deep generation of financial time series. Quantitative fi - nance, 20(9), 1419–1440. https://doi.org/10.1080/14697688.2020.1730426 [21] Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of econometrics, 31(3), 307–327. https://doi.org/10.1016/0304-4076(86)90063-1 [22] Aggarwal, A., Mittal, M., & Battineni, G. (2021). Generative adversarial network: An overview of theory and applications. Interna-tional journal of information management data insights, 1(1), 100004. https://doi.org/10.1016/j.jjimei.2020.100004 [23] Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., & Courville, A. (2017). Improved training of wasserstein GANs. Proceedings of the 31st international conference on neural information processing systems (pp. 5769 – 5779). Curran Associates Inc. https://doi.org/10.5555/3295222.3295327 [24] Van Dijk, M. T. P., De Graaf, C. S. L., & Oosterlee, C. W. (2018). Between P and Q. Journal of risk and financial management, 11 (4), 67. https://doi.org/10.3390/jrfm11040067 [25] Hossain, A. S., Noushin, M. T., & Hasan, K. (2019). A comparison of the Black-Scholes option pricing model and its alternatives. Dhaka university journal of science, 67(2), 105-110. https://www.banglajol.info/index.php/DUJS/article/view/54581/38354 [26] Storn, R., & Price, K. (1997). Differential evolution. Journal of global optimization, 11(4), 341–359.

https://doi.org/10.1023/A:1008202821328 [27] Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM journal on scientific computing, 16(5), 1190–1208. https://doi.org/10.1137/0916069 [28] Yahoo Finance. (2025). Historical data for Nifty 50 [data set]. https://finance.yahoo.com/quote/%5ENSEI/history/?guccounter [29] National Stock Exchange of India. (2026). NIFTY 50 option market data [data set]. https://www.nseindia.com/ [31] Wang, Z., & Bovik, A. C. (2009). Mean squared error: Love it or leave it? IEEE signal processing magazine, 26(1), 98–117. https://doi.org/10.1109/MSP.2008.930649.

[32] Hasnaine, Q. R., Fouda, M. M., & Chiu, S. C. (2025, April). GPU vs. CPU: A comparative study of performance, architecture, and NVIDIA's innovations in modern computing. In 2025 IEEE 4th international conference on computing and machine intelligence (ICMI) (pp. 1-6). IEEE. https://doi.org/10.1109/ICMI65310.2025.11141315 [33] Raj, P., & Sekhar, C. (2020). Comparative study on CPU, GPU and TPU. international journal of computer science and information technology for education, 5, 31-38. https://doi.org/10.21742/IJCSITE.2020.5.1.04.

The experiments were performed on the following hardware and software configuration:

Operating System: Windows 10 (64-bit)

Processor: Intel Core i5-10350U CPU @ 1.20 GHz

RAM: 8 GB

Python Version: 3.12.12

TensorFlow Version: 2.19.0

NumPy Version: 2.0.2

Computing Environment: Google Colab (CPU runtime)

The Wasserstein GAN architecture and training procedure employed the following hyperparameters:

Latent Dimension: 32

Training Epochs: 400 per year

Batch Size: 32

Learning Rate: 0.0001 (both generator and discriminator)

Optimizer: Adam (default parameters: β1 = 0.9, β2 = 0.999, ϵ = 1 10 8)

Discriminator Updates per Generator Update: 2 (ncritic = 2)

Weight Clipping Range: [-0.01, 0.01]

Number of Synthetic Paths Generated: 10,000 per year

Input Layer: Latent vector of dimension 32

Hidden Layer: Fully connected, 64 neurons, ReLU activation

Output Layer: Fully connected, T neurons (no activation), where T 2 f 83 , 85, 86, 87g

Reshape: Output reshaped to (T, 1)

Input Layer: Sequence of shape (T, 1), Output Layer: Single neuron (no activation)

Convolutional Layer: 32 filters, kernel size 3, same padding, ReLU activation

Flatten Layer: Convert 2D feature maps to 1D

Hidden Layer: Fully connected, 32 neurons, ReLU activation

Historical NIFTY 50 log-returns were normalized using z-score standardization before training:

r µtrain rnorm = (A1) σtrain where µtrain and σtrain are computed from all available training data up to (but excluding) each target year. After generation, synthetic returns were denormalized and scaled to match the empirical volatility of the training data, followed by risk-neutral drift correction as described in Section 3.4.

The expanding window training approach resulted in the following dataset sizes per year.

Table A1. Training Dataset Sizes per Yearly Window

Year Sequence Length (T) Training Observations Training Sequences 2019 86 days 1,000 914 2020 85 days 1,250 1,165 2021 83 days 1,500 1,417 2022 87 days 1,750 1,663 2023 86 days 2,000 1,914 2024 87 days 2,250 2,163 2025 85 days 2,500 2,415

WGAN training exhibited stable convergence across all years. Both discriminator and generator losses approached zero as training progressed, indicating successful adversarial equilibrium. Representative training dynamics are shown below for selected epochs:

Table A2. WGAN Loss Convergence by Epoch

Year Epoch 0 Epoch 100 Epoch 200 Epoch 300 2019 = 0005 = 0003 = 0007 = 0007 = 0010 = 0012 2020 = 0004 = 0003 = 0006 = 0008 = 0008 = 0009 2021 = 0002 = 0002 = 0005 = 0006 = 0006 = 0004 2022 = 0002 = 0002 = 0002 = 0004 = 0004 = 0004 2023 = 0003 = 0003 = 0004 = 0004 = 0007 = 0007 2024 = 0004 = 0004 = 0006 = 0008 = 0016 = 0013 2025 = 0005 = 0004 = 0008 = 0007 = 0012 = 0012

No mode collapse or training instability was observed. Losses converged smoothly without divergence or oscillation, confirming the robustness of the WGAN framework with weight clipping.

Total training and path generation time across all seven years was approximately 8.37 minutes (CPU-only). Individual year training times ranged from approximately 1.0 to 1.3 minutes per model. This excludes subsequent option pricing calculations via Monte Carlo simulation.

No explicit random seed was set in the original experiments reported in this paper. Consequently, exact numerical replication will produce different (but statistically similar) results due to stochastic initialization of network weights and random sampling during training.

For deterministic reproducibility in future work, researchers should add the following code at the beginning of their implementation:import numpy as np import tensorflow as tf RANDOM_SEED = 42 np.random.seed(RANDOM_SEED)