* (15) S. Monnier and G. W. Moore, “Remarks on the Green-Schwarz terms of six-dimensional supergravity theories,” Commun. Math. Phys. 372, no. 3, 963 (2019) arXiv:1808.01334 [hep-th].
  * (16) S. J. Lee and T. Weigand, “Swampland Bounds on the Abelian Gauge Sector,” Phys. Rev. D 100, no. 2, 026015 (2019) arXiv:1905.13213 [hep-th].
  * (17) C. Angelantonj, Q. Bonnefoy, C. Condeescu and E. Dudas, “String Defects, Supersymmetry and the Swampland,” JHEP 11, 125 (2020) [arXiv:2007.12722 [hep-th].
  * (18) N. Raghuram, W. Taylor and A. P. Turner, “Automatic enhancement in 6D supergravity and F-theory models,” JHEP 07, 048 (2021) arXiv:2012.01437 [hep-th].
  * (19) H. C. Tarazi and C. Vafa, “On The Finiteness of 6d Supergravity Landscape,” arXiv:2106.10839 [hep-th].
  * (20) A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis and R. Rattazzi, “Causality, analyticity and an IR obstruction to UV completion,” JHEP 0610, 014 (2006) arXiv:hep-th/0602178.
  * (21) H. Ooguri and C. Vafa, “On the geometry of the string landscape and the swampland,” Nucl. Phys. B 766, 21 (2007) arXiv:hep-th/0605264.
  * (22) N. Arkani-Hamed, L. Motl, A. Nicolis and C. Vafa, “The string landscape, black holes and gravity as the weakest force,” JHEP 0706, 060 (2007) arXiv:hep-th/0601001.
  * (23) J. Polchinski, “Monopoles, duality, and string theory,” Int. J. Mod. Phys. A 19S1, 145 (2004) arXiv:hep-th/0304042.
  * (24) T. Banks and N. Seiberg, “Symmetries and Strings in Field Theory and Gravity,” Phys. Rev. D 83 084019 (2011) arXiv:1011.5120 [hep-th].
  * (25) D. Harlow and H. Ooguri, “Symmetries in quantum field theory and quantum gravity,” arXiv:1810.05338 [hep-th].
  * (26) P. S. Aspinwall and D. R. Morrison, “Nonsimply connected gauge groups and rational points on elliptic curves,” JHEP 9807, 012 (1998) hep-th/9805206.
  * (27) P. S. Aspinwall, S. H. Katz, and D. R. Morrison, “Lie groups, Calabi-Yau threefolds, and F theory”, Adv. Theor. Math. Phys. 4, 95–126 (2000) hep-th/0002012.
  * (28) C. Mayrhofer, D. R. Morrison, O. Till, and T. Weigand, “Mordell-Weil Torsion and the Global Structure of Gauge Groups in F-theory”, JHEP 1410, 016 (2014) arXiv:1405.3656 [hep-th].
  * (29) F. Apruzzi, M. Dierigl and L. Lin, “The Fate of Discrete 1-Form Symmetries in 6d,” arXiv:2008.09117 [hep-th].
  * (30) M. Cvetič, M. Dierigl, L. Lin and H. Y. Zhang, “String Universality and Non-Simply-Connected Gauge Groups in 8d,” Phys. Rev. Lett. 125, no.21, 211602 (2020) arXiv:2008.10605 [hep-th].
  * (31) M. B. Green, J. H. Schwarz and P. C. West, “Anomaly Free Chiral Theories In Six-Dimensions,” Nucl. Phys. B 254, 327 (1985).
  * (32) A. Sagnotti, “A Note on the Green-Schwarz mechanism in open string theories,” Phys. Lett. B 294, 196 (1992) arXiv:hep-th/9210127.
  * (33) J. Erler, “Anomaly Cancellation In Six-Dimensions,” J. Math. Phys. 35, 1819 (1994) arXiv:hep-th/9304104.
  * (34) V. Kumar and W. Taylor, “A bound on 6D ${\cal N}=1$ supergravities,” JHEP 0912, 050 (2009) arXiv:0910.1586 [hep-th].
  * (35) V. Kumar, D. S. Park and W. Taylor “6D supergravity without tensor multiplets,” JHEP 1104, 080 (2011). arXiv:1011.0726 [hep-th].
  * (36) A. Grassi and D. R. Morrison, “Anomalies and the Euler characteristic of elliptic Calabi-Yau threefolds,” Commun. Num. Theor. Phys. 6, 51 (2012) arXiv:1109.0042 [hep-th].
  * (37) D. S. Park and W. Taylor, “Constraints on 6D Supergravity Theories with Abelian Gauge Symmetry,” JHEP 1201, 141 (2012) arXiv:1110.5916 [hep-th].
  * (38) D. S. Park, “Anomaly Equations and Intersection Theory,” JHEP 1201, 093 (2012) arXiv:1111.2351 [hep-th].
  * (39) S. Monnier, G. W. Moore and D. S. Park, “Quantization of anomaly coefficients in 6D $\mathcal{N}=(1,0)$ supergravity,” JHEP 1802, 020 (2018) arXiv:1711.04777 [hep-th].
  * (40) N. Seiberg and E. Witten, “Comments on String Dynamics in Six Dimensions,” Nucl. Phys. B 471, 121 (1996) arXiv:hep-th/9603003.
  * (41) J.J. Heckman, D.R. Morrison and C. Vafa, “On the Classification of 6D SCFTs and Generalized ADE Orbifolds,” JHEP 1405, 028 (2014) arXiv:1312.5746 [hep-th].
  * (42) M. Del Zotto, J.J. Heckman, A. Tomasiello and C. Vafa, “6d Conformal Matter,” JHEP 1502, 054 (2015) arXiv:1407.6359 [hep-th].
  * (43) D. R. Morrison, “TASI lectures on compactification and duality,” Strings, Branes, and Gravity, pp. 653–719 hep-th/0411120.
  * (44) F. Denef, “Les Houches Lectures on Constructing String Vacua,” Les Houches 87 (2008) pp. 483–610 arXiv:0803.1194 [hep-th].
  * (45) W. Taylor, “TASI Lectures on Supergravity and String Vacua in Various Dimensions,” arXiv:1104.2051 [hep-th].
  * (46) T. Weigand, “F-theory,” PoS TASI 2017, 016 (2018) arXiv:1806.01854 [hep-th].
  * (47) D. R. Morrison, “What is F-theory?,” to appear.
  * (48) K. Kodaira, “On compact analytic surfaces. II, III,” Ann. Math. 77 563 (1963); Ann. Math. 78 1 (1963).
  * (49) V. Sadov, “Generalized Green-Schwarz mechanism in F theory,” Phys. Lett. B 388, 45 (1996) arXiv:hep-th/9606008.
  * (50) A. Grassi, D. R. Morrison, “Group representations and the Euler characteristic of elliptically fibered Calabi-Yau threefolds”, J. Algebraic Geom. 12 (2003), 321-356 arXiv:math/0005196.
  * (51) V. Kumar, D. R. Morrison and W. Taylor, “Mapping 6D ${\cal N}=1$ supergravities to F-theory,” JHEP 1002, 099 (2010) arXiv:0911.3393 [hep-th].
  * (52) T. W. Grimm and A. Kapfer, “Anomaly Cancelation in Field Theory and F-theory on a Circle,” JHEP 1605, 102 (2016) arXiv:1502.05398 [hep-th].
  * (53) F. Bonetti and T. W. Grimm, “Six-dimensional (1,0) effective action of F-theory via M-theory on Calabi-Yau threefolds,” JHEP 1205, 019 (2012) arXiv:1112.1082 [hep-th].
  * (54) L. Lin and T. Weigand, “Towards the Standard Model in F-theory,” Fortsch. Phys. 63, no. 2, 55 (2015) arXiv:1406.6071 [hep-th].
  * (55) M. Cvetič, D. Klevers, D. K. Mayorga Peña, P. K. Oehlmann and J. Reuter, “Three-Family Particle Physics Models from Global F-theory Compactifications,” JHEP 1508, 087 (2015) arXiv:1503.02068 [hep-th].
  * (56) C. Lawrie, S. Schäfer-Nameki and J. M. Wong, “F-theory and All Things Rational: Surveying U(1) Symmetries with Rational Sections,” JHEP 1509, 144 (2015) arXiv:1504.05593 [hep-th]
  * (57) T. W. Grimm, A. Kapfer and D. Klevers, “The Arithmetic of Elliptic Fibrations in Gauge Theories on a Circle,” JHEP 1606, 112 (2016) arXiv:1510.04281 [hep-th].
  * (58) L. Lin and T. Weigand, “G4-flux and standard model vacua in F-theory,” Nucl. Phys. B 913, 209 (2016) arXiv:1604.04292 [hep-th].
  * (59) M. Cvetič and L. Lin, “The Global Gauge Group Structure of F-theory Compactification with U(1)s,” JHEP 1801, 157 (2018) arXiv:1706.08521 [hep-th].
  * (60) M. Cvetič, L. Lin, M. Liu and P. K. Oehlmann, “An F-theory Realization of the Chiral MSSM with $\mathbb{Z}_{2}$-Parity,” JHEP 1809, 089 (2018) arXiv:1807.01320 [hep-th].
  * (61) M. Cvetič, J. Halverson, L. Lin, M. Liu and J. Tian, “Quadrillion $F$-Theory Compactifications with the Exact Chiral Spectrum of the Standard Model,” Phys. Rev. Lett. 123, no. 10, 101601 (2019) arXiv:1903.00009 [hep-th].
  * (62) W. Taylor and A. P. Turner, “Generic construction of the Standard Model gauge group and matter representations in F-theory,” arXiv:1906.11092 [hep-th].
  * (63) N. Raghuram, W. Taylor and A. P. Turner, “General F-theory models with tuned $(SU(3)\times SU(2)\times U(1))/\mathbb{Z}_{6}$ symmetry,” arXiv:1912.10991 [hep-th].
  * (64) M. Bershadsky, K. A. Intriligator, S. Kachru, D. R. Morrison, V. Sadov and C. Vafa, “Geometric singularities and enhanced gauge symmetries,” Nucl. Phys. B 481, 215 (1996) hep-th/9605200.
  * (65) S. H. Katz and C. Vafa, “Matter from geometry,” Nucl. Phys. B 497, 146 (1997) arXiv:hep-th/9606086.
  * (66) D. R. Morrison and W. Taylor, “Matter and singularities,” JHEP 1201, 022 (2012) arXiv:1106.3563 [hep-th].
  * (67) D. Klevers, D. R. Morrison, N. Raghuram and W. Taylor, “Exotic matter on singular divisors in F-theory,” JHEP 1711, 124 (2017) arXiv:1706.08194 [hep-th]
  * (68) T. W. Grimm and W. Taylor, “Structure in 6D and 4d N=1 supergravity theories from F-theory,” JHEP 1210, 105 (2012) arXiv:1204.3092 [hep-th].
  * (69) E. Witten, “Phase transitions in $M$-theory and $F$-theory,” Nuclear Phys. B 471 (1996) 195–216, arXiv:hep-th/9603150.
  * (70) J. Tate, “On the conjectures of Birch and Swinnerton-Dyer and a geometric analog,” Séminaire Bourbaki, Années 1964/65-1965/66, Exp. No. 306, W. A. Benjamin, New York-Amsterdam, pp. 415–440.
  * (71) T. Shioda, “On elliptic modular surfaces,” J. Math. Soc. Japan 24 (1972) 20–59.
  * (72) R. Wazir, “Arithmetic on elliptic threefolds,” Compos. Math., 140, pp. 567-580, (2004).
  * (73) D. R. Morrison and W. Taylor, “Classifying bases for 6D F-theory models,” Central Eur. J. Phys. 10, 1072 (2012), arXiv:1201.1943 [hep-th].
  * (74) O. Zariski, “The theorem of Riemann-Roch for high multiples of an effective divisor on an algebraic surface,” Ann. of Math. (2) 76, pp. 560–615, (1962).
  * (75) T. Bauer, M. Caibăr, and G. Kennedy, “Zariski decomposition: a new (old) chapter of linear algebra,” Amer. Math. Monthly 119, pp. 25–41, (2012), arXiv:0911.4500 [math.AG].
  * (76) G. Martini and W. Taylor, “6D F-theory models and elliptically fibered Calabi-Yau threefolds over semi-toric base surfaces,” JHEP 1506, 061 (2015) arXiv:1404.6300 [hep-th].
  * (77) Y. N. Wang, “Tuned and Non-Higgsable U(1)s in F-theory,” JHEP 1703, 140 (2017) arXiv:1611.08665 [hep-th].
  * (78) D. R. Morrison, D. S. Park and W. Taylor, “Non-Higgsable abelian gauge symmetry and $\mathrm{F}$-theory on fiber products of rational elliptic surfaces,” Adv. Theor. Math. Phys. 22, 177 (2018) arXiv:1610.06929 [hep-th].
  * (79) W. Taylor and A. P. Turner, “Generic matter representations in 6D supergravity theories,” JHEP 1905, 081 (2019) arXiv:1901.02012 [hep-th].
  * (80) J. H. Silverman and J. Tate, “Rational Points on Elliptic Curves,” Springer, 1992.
  * (81) J. H. Silverman, “The Arithmetic of Elliptic Curves,”, Springer, 1986.
  * (82) D. Kubert, “Universal bounds on the torsion of elliptic curves,” Proc. London Math Soc. 33 193-237 (1976).
  * (83) N. Hajouji and P. K. Oehlmann, “Modular Curves and Mordell-Weil Torsion in F-theory,” arXiv:1910.04095 [hep-th].
  * (84) B. Mazur, “Modular curves and Eisenstein ideal,” IHES Publ. Math. 47, 33-186 (1977); “Rational isogenies of prime degree,” Invent. Math. 44 129-162 (1978).
  * (85) N. Bourbaki, Groupes et algèbres de Lie, Ch. IV, V, VI, Hermann, Paris, 1968.
  * (86) S. Katz, D. R. Morrison, S. Schäfer-Nameki and J. Sully, “Tate’s Algorithm and F-theory,” JHEP 1108, 094 (2012) arXiv:1106.3854 [hep-th].
  * (87) L. Anderson, J. Gray, N. Raghuram and W. Taylor, “Matter in Transition,” JHEP 1604, 080 (2016). arXiv:1512.05791 [hep-th]
  * (88) Y. Kimura, “Unbroken $E_{7}\times E_{7}$ nongeometric heterotic strings, stable degenerations and enhanced gauge groups in F-theory duals,” arXiv:1902.00944 [hep-th].
  * (89) M. Cvetič, J. J. Heckman and L. Lin, “Towards Exotic Matter and Discrete Non-Abelian Symmetries in F-theory,” JHEP 1811, 001 (2018) arXiv:1806.10594 [hep-th].
  * (90) Y. C. Huang and W. Taylor, “Comparing elliptic and toric hypersurface Calabi-Yau threefolds at large Hodge numbers,” JHEP 1902, 087 (2019) arXiv:1805.05907 [hep-th].
  * (91) D. R. Morrison and D. S. Park, “F-Theory and the Mordell-Weil Group of Elliptically-Fibered Calabi-Yau Threefolds,” JHEP 1210, 128 (2012) arXiv:1208.2695 [hep-th].
  * (92) D. R. Morrison and W. Taylor, “Sections, multisections, and U(1) fields in F-theory,” J. Singularities 15 pp. 126–149 (2016) arXiv:1404.1527 [hep-th].
  * (93) D. Klevers, D. K. Mayorga Pena, P. K. Oehlmann, H. Piragua and J. Reuter, “F-Theory on all Toric Hypersurface Fibrations and its Higgs Branches,” JHEP 1501, 142 (2015) arXiv:1408.4808 [hep-th]
  * (94) N. Raghuram, “Abelian F-theory Models with Charge-3 and Charge-4 Matter,” JHEP 1805, 050 (2018) arXiv:1711.03210 [hep-th].
  * (95) F. M. Cianci, D. K. Mayorga Peña and R. Valandro, “High U(1) charges in type IIB models and their F-theory lift,” JHEP 1904, 012 (2019) arXiv:1811.11777 [hep-th].
  * (96) A. Collinucci, M. Fazzi, D. R. Morrison and R. Valandro, “High electric charges in M-theory from quiver varieties,” JHEP 1911, 111 (2019) arXiv:1906.02202 [hep-th].
  * (97) N. Raghuram and W. Taylor, “Large U(1) charges in F-theory,” JHEP 1810, 182 (2018) arXiv:1809.01666 [hep-th].
  * (98) D. R. Morrison and W. Taylor, “Toric bases for 6D F-theory models,” Fortsch. Phys. 60, 1187 (2012) arXiv:1204.0283 [hep-th].
  * (99) L. B. Anderson, X. Gao, J. Gray and S. J. Lee, “Fibrations in CICY Threefolds,” JHEP 1710, 077 (2017) arXiv:1708.07907 [hep-th].
  * (100) U. Persson. “Configurations of Kodaira fibers on rational elliptic surfaces.” Mathematische Zeitschrift 205.1 (1990): 1-47.
  * (101) R. Miranda. “Persson’s list of singular fibers for a rational elliptic surface.” Mathematische Zeitschrift 205.1 (1990): 191-211.
  * (102) B. C. Hall, Lie groups, Lie algebras, and representations: An elementary introduction, Second Edition, Graduate Texts in Mathematics, vol. 222, Springer-Verlag, New York, 2015.
  * (103) T. Bröcker and T. tom Dieck, Representations of compact Lie groups, Graduate Texts in Mathematics, vol. 98, Springer-Verlag, New York, 1995, Corrected reprint of the 1985 translation.
  * (104) J. C. Baez and J. Huerta, “The algebra of grand unified theories,” Bull. Am. Math. Soc. 47 (2010) 483–552, arXiv:0904.1556 [hep-th].
  * (105) H. Georgi and S. Glashow, “Unity of All Elementary Particle Forces,” Phys. Rev. Lett. 32 (1974) 438–441.


# Pricing Multivariate European Equity Option Using Gaussians Mixture
Distributions and EVT-Based Copulas.

Hassane Abba Mallam${}^{1},$ Diakarya Barro2Yaméogo WendKouni3 and Bisso
Saley${}^{1}.$ 1UAM, Université Abdou Moumouni (Niger) 2 Université Thomas
Sankara (Burkina Faso). 3 LANIBIO, Université Joseph Ki-Zerbo (Burkina Faso)

###### Abstract

In this article, we present an approach which allows to take into account the
effect of extreme values in the modeling of financial asset returns and in the
valorisation of associeted options. Specifically, the marginal distribution of
assets returns is modeled by a mixture of two gaussiens distributions.
Moreover , we model the joint dependence structure of the returnsusing a
copula function, the extremal one, which is suitable for our financial data,
particularly the extreme value copula. Applications are made on the Atos and
Dassault Systems actions of the CAC40 index. Monte-Carlo method is used to
compute the values of some equity options such as the call on maximum, the
call on minimum, the digital option and the spreads option with the basket
(Atos, Dassault systems) as underlying.

Keywords: Options, Extremes values, Gaussians mixture; Copulas; Monte-Carlo;
European market.

2010 Mathematical Subject Classification: 60E05; 91B24; 91G20; 91G60; 91G70.

2010 JEL classification: C02; F31; G13.

##  1 Introduction

Since the pionering work of Black-Scholes [1] and Cox and al. [6]
(respectively in the continuous and discrete case), option pricing has become
a crucial topic in finance. Indeed, considering an european-type option on an
underlying asset with a price $S_{t}$; strike $K$ and expiration $T$, Black
and Scholes have made it possible to determine a formula for the price of such
options under certain assumptions, the fundamental of which are the lack of
arbitrage opportunity and that on the price $S_{t}$ of the asset underlying
($S_{t}$ follows geometric Brownian motion), i.e.

$dS_{t}=\mu S_{t}dt+\sigma S_{t}W_{t},$

where $\mu$ and $\sigma$ are constant and $W_{t}$ is a standard geometric
Brownian motion.

Thus, the formula of the relative theoretical value $\dfrac{C_{t}}{S_{t}}$ of
a call option is then given by:

$C^{all}_{t}=N(d_{1})-\kappa e^{-r(T-t)}N(d_{2}),$ (1)

and the relative theoretical value of a put is given by:

$P^{ut}_{t}=-N(-d_{1})+\kappa e^{-r(T-t)}N(-d_{2}),$ (2)

where $N(.)$ is the standard normal distribution;
$d_{1}=\dfrac{(r+\frac{\sigma^{2}}{2})(T-t)-\ln(k)}{\sigma\sqrt{T-t}}$;
$d_{2}=d_{1}-\sigma\sqrt{T-t}$ and $\kappa=\frac{K}{S_{t}}$ the relative
strike and $r$ the risk-free interest rate.

Options are essential financial products allowing to their holders to hedge
against the risk of falls in their investments. This is how we are
increasingly seeing the creation of several types of options such as exotic
options; multivariate options; etc., with the aim of providing more security.
As a result, valuation models are also evolving. Of all the multiple option
pricing models, it turns out that each one is primarily based on the dynamics
of underlying asset pricing model (for options with only one underlying) or
the asset portfolio (for options on multiple assets), when market assumptions
are known. In fact, since the assumption of no arbitrage opportunity (NAO) in
the markets is the basis of the fundamental results obtained in finance, it’s
considered by default 111There are markets on which the arbitration assumption
is considered.. The advantage under this NAO assumption is that, associated
with that of market completeness, there is a single risk-neutral probability
for which the discounted flows are martingales. In the univariate case, one of
the most interesting results obtained in this direction on valuation is that
of Breeden and al. [3]. It states that the second derivative (when it exists
and continuous) of the price of a standard option relatively to the strike
coincides with the risk-neutral density. Indeed, if $D_{t}$ is the price of a
european option of an underlying asset with price $X_{t}$ having for pay-off
$g(X_{T})$, $T$ the time to expiration and $r$ the risk-free interest rate;
then the risk-neutral density $f^{*}(X_{T})$ is linked to $D_{t}$ by:

$D_{t}=e^{-r(T-t)}\mathbb{E}^{*}\\{g(X_{T})\\}=e^{-r(T-t)}\int
g(X_{T})f^{*}(X_{T})dX_{T}.$ (3)

For valuation in the multivariate framework, this risk-neutral formula is a
simple generalization. Talponen and al. [19] recently gave the multivariate
version of the univariate result.

Multivariate options (rainbow; digital; quantos;etc.), which will be the main
subject of our study in this paper, constitute the central themes of current
research on financial risk coverage. The advantage lies in the fact that they
offer better coverage against risks. Indeed, the basic idea is that when the
option is a function of several assets, the fall in value of one asset is
compensated by the rise of another asset in the portfolio. Thus the
association or dependence between assets plays a major role in the pricing of
these types of options. To take such an aspect into account in the valuation,
the use of the copula is a good alternative.

The valuation of multivariate options by copulas is in full development. The
copula gives the advantage of joining the marginal and the dependence
structure. This is the case for many works on the valuation of options with
copulas, the emphasis is first on marginal risk-neutral densities and then on
the joint risk-neutral density (risk-neutral copula). For example, we can cite
the work of, Cherubuni and al. [4]; Cherubuni and al [5], Rosenberg and al.
[11], Salmon and al. [13]; Slavchev and al. [17]. However, all this work did
not take into account the effects of extreme values in the marginal, which is
not without effect on valuation (risk of over-valuation or under-valuation).
However, there are other copula modeling approaches based on volatility
dynamics as in Goorbergh and al. [20]. The reader can consult this last one
for full details on the literature on this approach.

In this present study, we propose a valuation method for multivariate options
allowing to take into account the effects of extreme values in the marginal
and the joint structure on the basis of the works of Edier et al [7] and the
use of extreme values copulas.
In the rest of this work, the first section we give the results obtained by
Idier et al[7]. which will be necessary and somes essentielles notions on
copula. In the second section, we expose the methodology used for leading
properly the application of the demarche. Then, the obtained results of
different estimations and simulations are presented, with their analysis and
interpretations. The last section present a conclusion and discussion.

##  2 Preliminaries

###  2.1 Results of an approach of modeling financial assets

It is prouved that the empirical distribution of financial assets returns has
thicker tails than that of the Gaussian distribution. This indicates the
presence of extreme values. This fact shows also that the normal distribution
does not make it possible to model rigorously the returns of financial assets
because it does not take into account the extremes. This is the case with the
method proposed by Black-Scholes.

To take into account the effects of extreme values, Idier et al [7] proposed,
as an alternative to the normal distribution, to model the distribution of the
rates of return of the underlying asset of a univariate option, under NAO
assmption, by a mixture of Gaussian distribution in the continuous framework
222Their method is a generalization of the method in the discrete case of
Bertholon, Monfort and Pegoraro(2006); Pegoraro(2006). They justified their
choice by the fact that a mixture of Gaussian distributions makes it possible
to approximate all the distributions usually used (Gaussian; alpha-stable;
Student; hyperbolic; etc.); also that it has certain theoretical properties
allowing easy handling in the frame a theoretical model for valuing asset
price; that it is easy to simulate and can reproduce various sets (mean,
variance, skewness and kurtosis) observed in the data.

Under the assumption that the historical distribution of the returns of the
underlying $X_{t+1}=\ln(\frac{S_{t+1}}{S_{t}})$ where $S_{t}$ is the price at
time $t$ of the underlying asset, is a mixture of $2$ Gaussian distribution.
Its density is given by:

$f(x)=\sum_{i=1}^{2}p_{i}.n(x,\mu_{i},\sigma_{i}^{2}),$ (4)

where

$n(x,\mu_{i},\sigma_{i}^{2})=\frac{1}{\sigma_{i}\sqrt{2\pi}}\exp\bigg{\\{}\dfrac{-(x-u_{i})^{2}}{2\sigma_{i}^{2}}\bigg{\\}},$

is the density of a Gaussian distribution with mean $\mu_{i}$ and standard
deviation $\sigma_{i}$; $0<p_{i}<1\;and\;\sum_{i=1}^{2}p_{i}=1.$

Moreover, the stochastic discount factor is characterized by an affine
exponential form, i.e.

$M_{t;t+1}=\exp\\{\alpha_{t}X_{t}+\beta_{t}\\}.$ (5)

They establish, under these assumptions, that the risk-neutral distribution is
also a Gaussian mixture and that its density $f^{*}$ is defined by:

$f^{*}(x)=\sum_{i=1}^{2}v_{i}.n(x,\mu_{i}+\alpha.\sigma_{i}^{2},\sigma_{i}^{2}),$
(6)

where

$v_{i}=(p_{i}\exp\\{\mu_{i}\alpha+\frac{\alpha^{2}}{2}\sigma_{i}^{2}\\})/(\sum_{i=1}^{2}p_{i}\exp\\{\mu_{i}\alpha+\frac{\alpha^{2}}{2}\sigma_{i}^{2})\;with\;0<v_{i}<1,\sum_{i=1}^{2}v_{i}=1.$

Thus, they derive the relative theoretical price of an European call with a
one-period maturity ($T=t+$ 1) and a relative strike $k$:

$c_{t}(T,k)=\sum_{i=1}^{2}v_{i}\gamma_{i}c_{bs}(\sigma^{2}_{i},\frac{k}{\gamma_{i}}),$
(7)

where $c_{bs}(.,.)$ is the Black-Scholes one period $(T=t+1)$ formula for the
relative price of a call and
$\gamma_{i}=\exp\\{\mu_{i}+\alpha.\sigma_{i}^{2}-r+\frac{\sigma^{2}_{i}}{2}\\}$,
for $J=2$.

###### Remark 1

The existence of the call-put parity relation makes it possible to simplify
the task in calculating option price. It is then sufficient to calculate the
price of the call to deduce that of the corresponding put (or vice versa) by
the relation:

$C_{t}(T,K)+Ke^{-r(T-t)}=P_{t}(T,K)+S_{t}.$ (8)

###  2.2 A survey of copulas

In this section we recall the basics notions on copulas. These are the
definitions and properties essential for our study. For more details on
copulas, see Nelsen [10].

####  2.2.1 Definitions and properties

The copula is a function allowing to capture the structure of dependence
between several random variables.

A function $C:[0,1]^{d}\rightarrow[0,1]$ is a d-copula if it satisfies the
following properties:

  * i)

For all $u\ in[0,1]$, $C(1,...1,u,1,...1)=u$;

  * ii)

For all $u_{i}\ in[0,1]$, $C(u_{1},...,u_{d})=0$ if at least one of the
$u_{i}$ is zero;

  * iii)

$C$ is "grounded" and d-increasing, i.e

$\sum_{i_{1}=1}^{2}...\sum_{i_{d}=1}^{2}(-1)^{i_{1}+...+i_{d}}C(u_{1,i_{1}},...,u_{d,i_{d}}))\geq
0,$ (9)

for all $(u_{1,1},...,u_{d,1})$ and $(u_{1,2},...,u_{d,2})\in[0,1]^{d}$ with
$u_{d,1}\leq u_{d,2}$.

The fundamental result on the copula due to Sklar (1969) states that: For a
whole multivariate distribution $F$ with continuous marginal $F_{1}$, …,
$F_{d}$, then there exists a unique 333Uniqueness is not guaranteed when
marginal are not continuous. copula $C:[0,1]^{d}\rightarrow[0,1]$ such that:

$F(x_{1},..,x_{d})=C(F_{1}(x_{1}),...,F_{d}(x_{d})).$ (10)

Conversely, when $C$ is a copula and $F_{1}$, …, $F_{d}$ are marginal
distributions, the function $F$ defined by (10) is a multivariate distribution
of marginal distributions $F_{1}$, …, $F_{d}$.

This result makes it possible to deduce several properties of the copula
including invariance by any monotonic transformation. Another consequence of
Sklar’s theorem is that every copula $C$ satisfies

$max\bigg{(}\sum_{i=1}^{d}u_{i}-d+1;0\bigg{)}\leq C(u_{1},u_{2},...,u_{d})\leq
min(u_{1},u_{2},...,u_{d}).$ (11)

This relation is the variant in terms of copulas of the Frechet-Hoeffding
bounds of a multivariate distribution. The upper bound
$min(u_{1},u_{2},...,u_{d})$ is the comonotonic copula representing the
perfect positive dependence. The lower bound
$max\big{(}\sum_{i=1}^{d}u_{i}-d+1;0\big{)}$ is a copula only for $d=2$. In
this case it represents the perfect negative dependence.

###### Remark 2

If $\bar{F}$ is the multivariate survival distribution of a $F$ distribution
of marginals $F_{i};\;i=1,..,d$, then the survival copula, denoted by
$\bar{C}$, is defined by:

$\bar{F}(x_{1},...,x_{d})=\bar{C}(\bar{F}_{1}(x_{1}),...,\bar{F}_{d}(x_{d})).$

The survival copula $\bar{C}$ is related to the copula $C$, for all
$(u_{1},u_{2},...,u_{d})\in[0,1]^{d}$, by:

$\bar{C}(u_{1},u_{2},...,u_{d})=\sum_{M\subset
N}(-1)^{m}C\big{(}(1-u_{1})^{1_{1}\in M},(1-u_{2})^{1_{2}\in
M},...,(1-u_{d})^{1_{d}\in M}\big{)}$ (12)

where $N=\\{1,2,...,d\\}$, $m=\lvert M\rvert$ is the cardinality of $M$, and
$1_{i}\in M$ indicates that $i$ belongs to $M$.

It is therefore advisable not to confuse the dual copula with the survival
copula.

####  2.2.2 Sample of Copulas for finance in this study

Archimedian copulas

In the literature, there are several families of copulas and some of which are
more suited to financial modeling. Archimedean copulas family includes the
models of Clayton, Frank and Gumbel. These copulas have the advantage of
capturing the structure of positive or negative dependence between the
variables. These types of dependences are characteristics of financial
variables, which justifies the use of this copula family. In terms of option
pricing, for example, these copulas have been used in Cherubuni and Luciano
([4]) and in Slavchev and Wilkens ([17]).

Table 1: Examples of Archimediean copulas. Family | Archimedean generator | Copula $C_{\theta}(u)$, $u=(u_{1},...,u_{d})$
---|---|---
Clayton |  $\psi_{\theta}(t)=\frac{t^{-\theta}-1}{\theta}$, $\theta>0$ | $C_{\theta}(u)=\bigg{(}\sum_{i=1}^{d}(u_{i})^{-\theta}-d+1\bigg{)}^{-\frac{1}{\theta}}$
Frank |  $\psi_{\theta}(t)=\frac{-1}{\theta}\log\big{(}1-(1-e^{-\theta})e^{-t}\big{)},$ $\theta>0$ | $C_{\theta}(u)=\frac{-1}{\theta}\log\bigg{(}1-\dfrac{\prod_{i=1}^{d}(e^{-u_{i}\theta-1})}{(e^{-\theta-1})^{d-1}}\bigg{)}$
Gumbel |  $\psi_{\theta}(t)=(-\ln t)^{\theta}$,with $\theta\geq 1$ | $C_{\theta}(u)=\exp\bigg{(}-\bigg{(}\sum_{i=1}^{d}(-\log u_{i})^{\theta}\bigg{)}^{\frac{1}{\theta}}\bigg{)}$

Elliptical copulas

Other types of copulas used in finance are the normal copula and the t-copula.
They belong to the family of elliptical copulas which describe the dependence
structure of elliptical distributions. The choice of this family is justified
by the fact that elliptic distributions have long been used to model random
phenomena in many fields. Despite the demonstration of the leptokurtic
character of the returns of financial series, of which they have the weakness
to rigorously model, their continued use in finance.

The normal copula is also known as the Gaussian copula. Its expression is:

$C_{\Sigma}(u_{1},u_{2},...,u_{d})=\int_{-\infty}^{N^{-1}(u_{1})}...\int_{-\infty}^{N^{-1}(u_{d})}\dfrac{1}{(2\pi)^{d/2}\lvert\Sigma\rvert^{1/2}}\exp\big{(}-\frac{1}{2}y^{\prime}\Sigma^{-1}y\big{)}dy,$
(13)

where $N^{-1}$ is the quantile of the standard normal distribution and
$\Sigma$ is the correlation matrix.

The t-copula is the one describing the dependence of t-distributions. It has
the advantage of capturing the dependency in the distribution tails more than
the normal copula. Its expression is given by:

$C_{v,\Sigma}(u_{1},u_{2},...,u_{d})=\int_{-\infty}^{t_{v}^{-1}(u_{1})}...\int_{-\infty}^{t_{v}^{-1}(u_{d})}\dfrac{\Gamma\big{(}\frac{v+d}{2}\lvert\Sigma\rvert^{1/2}\big{)}}{\Gamma(v/2)(v\pi)^{d/2}}\big{(}1+\frac{1}{v}y^{\prime}\Sigma^{-1}y\big{)}^{\frac{v+d}{2}}dy,$
(14)

where $t_{v}^{-1}$ is the quantile of the t-distribution with $v$ degree of
freedom; $\Gamma$ is the gamma function and $\Sigma$ the correlation matrix.

####  2.2.3 Estimation-Adequacy test of a copula

The choice of the copula that rigorously describing multivariate statistical
data requires estimation and conformity testing. There are several techniques
in the literature for estimating copulas belonging to different families:
parametric; semi-parametric and non-parametric. For more details on these
methods, see Bouyé [2].

Estimation by the IFM method

The IFM method (inference functions of margins) is a two-step estimation
method of a copula. It was presented by Shih and Louis [14] in the bivariate
case and then developed in dimension greater than two by Joe and Xu [8]. It is
carried out as follows:

  * 1)

the first step consists in finding the estimators $\hat{\alpha}_{i}$ of the
parameters $\alpha_{i},\;i=1,..,d$ for marginal distributions by maximum
likelihood:

$\hat{\alpha}_{i}=argmax\sum_{j=1}^{N}\ln f_{n}(x^{j}_{i};\alpha_{i});$

  * 2)

Once the marginal have been determined, we estimate the parameter $\theta$ of
the copula that best describes these marginal by the maximum likelihood.

$\hat{\theta}=argmax\sum_{j=1}^{N}\ln
c(F_{1}(x^{j}_{1};\hat{\alpha}_{1});...;F_{d}(x^{j}_{d};\hat{\alpha}_{d});\theta).$

One of the advantages of this method is that under certain conditions of
regularity, the IFM estimator is consistent and asymptotically normal.

Also, in terms of numerical computation time, this method is better than the
"direct" maximum likelihood method since it is simpler and faster.

Fit test

To confirm whether the chosen parametric copula models the data well, it is
necessary to perform a test. The most powerful tests are based on the
processes $\sqrt{n}(\hat{C}-C_{\hat{\theta}})$, where $\hat{C}$ and
$C_{\hat{\theta}}$ are respectively the empirical copula and the parametric
copula.

The Cramer-von Mises statistic is by far the most used because it gives
satisfactory results. It is defined by:

$\int_{[0;1]^{d}}n(\hat{C}-C_{\hat{\theta}})d\hat{C}.$ (15)

Others criteria such as the Akaike Criterion (AIC) and the Bayeian Inference
Criterion(BIC) are very often used for the choice of the best copula. They are
respectively defined by:

$\displaystyle AIC=2m-log(l(\theta))$ (16) $\displaystyle
BIC=m*log(n)-log(l(\theta)).$ (17)

where $l(\theta)$ is the model likelihood for the estimated parameter
$\theta$, $m$ the number of estimated parameters and $n$ the data size.

##  3 Methodology and Application

The price of a multivariate option is a function of the density of the joint
distribution. Thus, their valuation requires the determination of the joint
risk-neutral density. To do this, it suffices to determine the marginal risk-
neutral densities and then to choose the copula that best describes their
dependence structure by using Sklar’s theorem. This perspective is possible
because the objective copula can be matched with the risk-neutral joint
copula, under certain conditions (see Rosenbergh [12]).

###  3.1 Methodology

Our approach consists firstly in determining the marginal risk-neutral
distributions by the procedure used by Idier and al. [7]. This in order to
take into account the effect of extreme values in the margins. We will also
limit ourselves to the case of a mixture of two Gaussians in this study. Then,
we will choose among the families of copulas listed in the section 2.2.2 the
one that best suits the study. And finally we will determine the prices of the
multivariate options by numerical integration (Monte-Carlo method) by using
the formulas provided below for multivariate options considered.

We will be particularly interested by rainbow options (those relating to the
maximum or the minimum of several assets, etc.). These kinds of options have
been the subject of many studies as in Stulz [18] and Jonshon [9].

Consider $d$ assets whose price at maturity $T$ are denoted by
$S_{1}^{T};...S_{d}^{T}$ and denote by $X_{1}^{T};...;X_{d}^{T}$,
respectively, the returns associated to each asset at instant $T$ (with for
all $i=1,...,d$; $X_{i}^{T}=\ ln(S_{i}^{T}/S_{i}^{t})$). For a chosen strike
price $K$, we consider the following different types of rainbows: spread
option; option on the maximum; option on the minimum and digital option.

####  3.1.1 Spread Option

Having a pay-off equal to $max(S_{2}^{T}-S_{1}^{T}-K;0)$, its value is
calculated by:

$\displaystyle V_{OS}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}max(S_{2}^{T}-S_{1}^{T}-K;0)\bigg{\\}}=e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}\int_{\mathbb{R}}\mathbb{I}_{\\{K+S_{1}^{T}\leq
x\leq S_{2}^{T}\\}}dx\bigg{\\}}$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}\int_{\mathbb{R}}\mathbb{I}_{\\{K+S_{1}^{T}\leq
x\\}}-\mathbb{I}_{\\{K+S_{1}^{T}\leq x\;et\;S_{2}^{T}\leq x\\}}dx\bigg{\\}}$

which gives,

$\displaystyle V_{OS}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{-\infty}^{+\infty}\mathbb{P}^{*}{\\{K+S_{1}^{T}\leq
x\\}}-\mathbb{P}^{*}{\\{K+S_{1}^{T}\leq x\;et\;S_{2}^{T}\leq x\\}}dx$
$\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{-\infty}^{+\infty}\mathbb{P}^{*}{\\{\dfrac{S_{1}^{T}}{S_{1}^{t}}\leq\dfrac{x-K}{S_{1}^{t}}\\}}-\mathbb{P}^{*}{\\{\dfrac{S_{1}^{T}}{S_{1}^{t}}\leq\dfrac{x-K}{S_{1}^{t}}\;et\;\dfrac{S_{2}^{T}}{S_{2}^{t}}\leq\dfrac{x}{S_{2}^{t}}\\}}dx$

and finaly,

$\displaystyle V_{OS}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{-\infty}^{+\infty}\mathbb{P}^{*}{\\{X_{1}^{T}\leq\log(\dfrac{x-K}{S_{1}^{t}})\\}}-\mathbb{P}^{*}{\\{X_{1}^{T}\leq\log(\dfrac{x-K}{S_{1}^{t}})\;et\;X_{2}^{T}\leq\log(\dfrac{x}{S_{2}^{t}})\\}}dx$
(18) $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{K}^{+\infty}e^{r(T-t)}\dfrac{\partial P_{t,1}}{\partial
k_{1}}(T;\log(x_{1}-k_{1}))-$ $\displaystyle-$ $\displaystyle
C\big{(}e^{r(T-t)}\dfrac{\partial P_{t,1}}{\partial
k_{1}}(T;\log(x_{1}-k_{1}));e^{r(T-t)}\dfrac{\partial P_{t,2}}{\partial
k_{2}}(T;\log(x_{2}))\big{)}dx,$

where

$x_{i}=\dfrac{x}{S_{i}^{t}};\;k_{i}=\dfrac{K}{S_{i}^{t}},$

and $P_{t,i}$ is the put $i$ price’s, for $i=1;2.$

####  3.1.2 Call on the maximum

Its pay-off is equal to $max\\{max(S_{1}^{T};...;S_{d}^{T})-K;0\\}$. Thus, its
price at maturity is given by:

$\displaystyle V_{CMax}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\\{max\\{max(S_{1}^{T};...;S_{d}^{T})-K;0\\}\\}=e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}\int_{\mathbb{R}}\mathbb{I}_{\\{k\
\leq x\leq max(S_{1}^{T};...;S_{d}^{T})\\}}dx\bigg{\\}}$ $\displaystyle=$
$\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}\int_{K}^{+\infty}1-\mathbb{I}_{\\{max(S_{1}^{T};...;S_{d}^{T})\leq
x\\}}dx\bigg{\\}}.$

we obtain,

$\displaystyle V_{CMax}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{K}^{+\infty}1-\mathbb{P^{*}}\\{S_{1}^{T}\leq
x;...;S_{d}^{T}\leq x\\}dx$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{K}^{+\infty}1-\mathbb{P^{*}}\\{X_{1}^{T}\leq
x_{1};...;X_{d}^{T}\leq x_{d}\\}dx$

and finaly,

$\displaystyle V_{CMax}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{K}^{+\infty}1-C\bigg{(}e^{r(T-t)}\dfrac{\partial
P_{t,1}}{\partial k_{1}}(T;x_{1});...;e^{r(T-t)}\dfrac{\partial
P_{t,d}}{\partial k_{d}}(T;x_{d})\bigg{)}dx\hskip 28.45274pt$ (19)

where $x_{i}=log(\dfrac{x}{S_{i}^{t}})$ and $P_{t,i}$ is the put $i$ price’s,
for $i=1;2.$

####  3.1.3 Call on the minimum

It admits for pay-off $max\\{min(S_{1}^{T};...S_{d}^{T})-K;0\\}$ and its value
at maturity is then defined by:

$\displaystyle V_{CMin}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\\{max\\{min(S_{1}^{T};...;S_{d}^{T})-K;0\\}\\}=e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}\int_{\mathbb{R}}\mathbb{I}_{\\{k\
\leq x\leq min(S_{1}^{T};...;S_{d}^{T})\\}}dx\bigg{\\}}$

which is equal to,

$\displaystyle V_{CMin}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\bigg{\\{}\int_{K}^{+\infty}\mathbb{I}_{\\{min(S_{1}^{T};...;S_{d}^{T})\geq
x\\}}dx\bigg{\\}}=e^{-r(T-t)}\int_{K}^{+\infty}\mathbb{P^{*}}\\{S_{1}^{T}\geq
x;...;S_{d}^{T}\geq x\\}dx$

at the end, we obtain

$\displaystyle V_{CMin}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{K}^{+\infty}\mathbb{P^{*}}\\{X_{1}^{T}\geq
x_{1};...;X_{d}^{T}\geq x_{d}\\}dx$ (20) $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\int_{K}^{+\infty}\bar{C}\bigg{(}-e^{r(T-t)}\dfrac{\partial
C_{t,1}}{\partial k_{1}}(T;x_{1});...;-e^{r(T-t)}\dfrac{\partial
C_{t,d}}{\partial k_{d}}(T;x_{d})\bigg{)}dx\hskip 28.45274pt$

where $x_{i}=log(\dfrac{x}{S_{i}^{t}})$ and $C_{t,i}$ is the call $i$ price’s,
for $i=1;2.$

####  3.1.4  Digital Option

It has for pay-off $\mathbb{I}_{\\{S_{1}^{T}\geq K_{1};...;S_{d}^{T}\geq
K_{d}\\}}$. Thus, its value at maturity is given by

$\displaystyle V_{ODig}(t)$ $\displaystyle=$ $\displaystyle
e^{-r(T-t)}\mathbb{E}^{*}\big{\\{}\mathbb{I}_{\\{S_{1}^{T}\geq
K_{1};...;S_{d}^{T}\geq
K_{d}\\}}\big{\\}}=e^{-r(T-t)}\mathbb{P}^{*}\big{\\{}S_{1}^{T}\geq
K_{1};...;S_{d}^{T}\geq K_{d}\big{\\}}$

which give us,

$V_{ODig}(t)=e^{-r(T-t)}\bar{C}\bigg{(}-e^{r(T-t)}\dfrac{\partial
C_{t,1}}{\partial k_{1}}(T;k_{1});...;-e^{r(T-t)}\dfrac{\partial
C_{t,d}}{\partial k_{d}}(T;k_{d})\bigg{)}$ (21)

where $k_{i}=\dfrac{K_{i}}{S_{i}^{t}}$; $C_{t,i}$ is the call $i$ price’s, for
$i=1;...;d$ and $\bar{C}$ the survival copula.

###### Remark 3

Not to forget that the quantities $-e^{r(T-t)}\dfrac{\partial
C_{t,d}}{\partial k}(T;k)$ and $e^{r(T-t)}\dfrac{\partial P_{t,d}}{\partial
k}(T;k)$ are both equal to the risk neutral distribution wich density is given
by relation (6) for our study.

###  3.2  Applications

We will focus on the bivariate options on the pair of Atos and Dassault
Systems shares. The data were obtained from Investing.com and relate to the
components of the CAC40 index of the Paris stock exchange. The collected data
concerns the closing prices for the period from July 01, 2014 to June 30, 2020
(1534 days).

At first, for each of the two assets (Atos; Dassault Systems) the parameters
of the two Gaussian regimes constituting the Gaussian mixture are determined
(table 2) as well as the proportions of each diet.

Table 2: Estimated parameters of the Gaussian mixture and their proportion. | Regime 1 | Regime 2 | Gaussien mixture | Empirical
---|---|---|---|---
|  |  |  | distribution
| Moyenne | -0.0072328 | 0.000764489 | 0.00013704 | 0.0001370
Ecart-type | 0.0603574 | 0.013530408 | 0.02142835 | 0.0214349
Atos | Skewness | 0 | 0 | -0.6131518 | -2.823688
| Kurtosis | 3 | 3 | 15.7 | 42.81632
| Proportion | 0.07845771 | 0.921542310 | – | –
| Moyenne | 0.00110506 | -0.00101651 | 0.0007598786 | 0.0007599
Ecart-type | 0.01014017 | 0.03315738 | 0.01577694 | 0.01630201
Dassault systems | Skewness | 0 | 0 | -0.03745108 | 0.1949504
| Kurtosis | 3 | 3 | 10.00661 | 12.39682
| Proportion | 0.83729906 | 0.16270094 | – | –

The next step, we determine the parameters $(\alpha;\beta)$ of the stochastic
discount factor defined by the relation (5) thanks to the assumptions of the
model in particular that of lack of arbitration opportunity. The results
obtained for a risk-free rate $r=0.025$ are presented in the table 3.

Table 3: Stochastic discount factor parameters for data returns with a risk-free rate $r=$ 0.025. | Atos | Dassault
---|---|---
$\alpha$ | 36.1209027 | 37.70565
$\beta$ | -0.3610132 | -0.3500945

####  3.2.1  Copulas fitting results

We present the results of the copula estimates associated with our data. In
each case, we will base ourselves on the Cramer-Von Mises statistic and/or the
AIC criteria for the choice of the best copula.

In the table 4, we present the estimated parameters (and their Cramer-Von
statistics) for bivariate copulas. It then emerges that the three copulas with
the best Cramer-Von Mises statistic are Tawn’s copula, Frank’s copula,
Gumbel’s copula in that order.

The table 5 gives the AIC and BIC of the parameters estimated for the
bivariate copulas chosen. Based on these criteria, the four best candidate
copulas for our bivariate data are the normal copula, the Husler-Reiss’s
copula, the Galambos’s copula and the Gumbel’s copula.

We can notice that for our data, the performance of each fitted copula differs
according to the criterion. It is then difficult to make a particular choice
on a copula in suck situation on the basis of the two criteria (Cramer-Von
Mises statistic vs AIC) combined. Nevertheless, if there is a choice to be
made between these two criteria, it would be more judicious to base oneself on
the Cramer-Von Mises test.

Table 4: Parameters of bivariate copulas selected and their Cramer-Von Mises test statistics. | Normale | Clayton | Gumbel | Frank | Tawn | Galambos | Husler-Reiss
---|---|---|---|---|---|---|---
parameter | 0.2822 | 0.1766 | 1.344 | 2.3166 | 0.6868 | 0.5995 | 0.9798
Statistique of Test | 1.469 | 3.257 | 0.72485 | 0.6136 | 0.6136 | 0.79162 | 0.84238
p-value | 0.0005 | 0.0005 | 0.0005 | 0.0005 | 0.0005 | 0.0005 | 0.0005
Table 5: AIC and BIC of the estimated parameters (by the maximum likelihood.) for the selected bivariate copulas.  | Normale | Clayton | Gumbel | Frank | Tawn | Galambos | Husler-Reiss
---|---|---|---|---|---|---|---
$log(l(\theta))$ | 365.9 | 325.1 | 345.4 | 291.8 | 320 | 354.7 | 358
AIC | -729.8 | -648.2 | -688.8 | -581.6 | -638 | -707.4 | -714
BIC | -724.46 | -642.86 | -683.46 | -576.26 | -632.66 | -702.06 | -708.66

####  3.2.2  Options prices by Monte-Carlo approach

In this section, we give the simulation results of the prices (for one period
$T=t+1$) of all options presented in the section above based on the basket
(Atos, Dassault systemes). The tables 6 and 7 give, respectively, the prices
of the call on maximum and the call of maximum. We fixe the price of each
asset of the bivariate basket to $120$. The values of their prices are
calculated when it is out-of-the money (OTM); at-the-money (ATM) and in-the-
money (ITM). For the cases of digital option and the spread option, we fixe,
respectively, the price of the basket to $S=(120,130)$ and $S=(100,120)$. We
give the prices of these options in the table 8 and 9 obtained also for
different strikes. The prices are calculated by Monte-carlo method with
$N=10^{5}$ simulations.444The choice of the bivariate options with our
underlying is simply for academic interest. In fact, they are not exchanged on
the market.

For the case of the call on maximum (table 6), the price obtained by normale
is superior than the prices with all others copulas (archimedian and extreme)
in all the three situations without only the case when it is at the money with
Cayton’s copula. We notice that the prices obtained with the others are
approximatily the same (weak discrepencies). The normale copula over-estimate
the price compared to the Tawn copula with has the good fitness test.

Table 6:  Prices of the call on maximum with different strikes. | Normale | Clayton | Gumbel | Frank | Galambos | Tawn | Husler-Reiss
---|---|---|---|---|---|---|---
OTM (K=130) | 2.7126 | 2.6947 | 2.6642 | 2.6888 | 2.6593 | 2.6664 | 2.6617
ATM(K=120) | 7.328 | 7.417 | 7.0226 | 7.101 | 7.0216 | 7.0264 | 7.0245
ITM(K=110) | 17.211 | 17.057 | 16.743 | 16.794 | 16.674 | 16.660 | 16.673

In the case of the call on minimum (table7), the normale copula present a
price which is lower than that obtained by any extreme value copula when it is
at-the-money or out-of-the money. We notice the contrary when it is in-the-
money.
Particularly, the normale copula give the highest price that of Tawn’s copula
when it is in-the-money with a small discrepency. And, when it is at-the-money
or out-of-the money the Tawn’s copula produice a price superior to that of
normale copula with a fairly large gap than the first situation.

Table 7:  Prices of the call on minimum with different strikes K. | Normale | Clayton | Gumbel | Frank | Galambos | Tawn | Husler-Reiss
---|---|---|---|---|---|---|---
OTM(K=130) | 0.0386 | 0.02188 | 0.05882 | 0.03639 | 0.06085 | 0.05437 | 0.06224
ATM(K=120) | 1.181 | 0.995 | 1.361 | 1.268 | 1.356 | 1.369 | 1.353
ITM(K=110) | 10.628 | 10.183 | 10.559 | 10.451 | 10.536 | 10.523 | 10.543

For prices of the digital option (table 8), that obtained by normale copula is
inferior to the others calculated by extreme copula in the three situations of
valuation.

Table 8:  Prices of the bivariate digital option (paying one unit of the money) with different strikes $K$. | Normale | Clayton | Gumbel | Frank | Galambos | Tawn | Husler-Reiss
---|---|---|---|---|---|---|---
OTM | 0.02283713 | 0.014745475 | 0.03258418 | 0.02359763 | 0.03306132 | 0.03181797 | 0.03329281
ATM | 0.5208968 | 0.507406 | 0.5273577 | 0.5344094 | 0.5258449 | 0.53006026 | 0.5247992
ITM | 0.97519641 | 0.9751979 | 0.9751963 | 0.9751963 | 0.9751963 | 0.9751964 | 0.9751963

For the spreads option (table 9), we notice that when it is at-the-money the
normale copula give the high price and the Tawn copula has the small price.
Other wise, when its in-the-money the normale copula give the smallest price.
Finaly, when the option is out-of-the money, the normale copula give the
second greatest price. Comparatively to the price obtained with Tawn’s copula,
the price calculated with normale copula is greatest when the option is at-
the-money and out-of-the-money but smallest when it is in-the-money.

Table 9:  Prices of the bivariate spread option with different strikes K. | Normale | Clayton | Gumbel | Frank | Galambos | Tawn
---|---|---|---|---|---|---
OTM(K=30) | 0.04626 | 0.05037 | 0.01244 | 0.02615 | 0.00898 | 0.01909
ATM (K=20) | 2.1429 | 1.379 | 0.9609 | 1.0818 | 1.006 | 0.8878
ITM (K=10) | 6.5442 | 7.8945 | 7.4151 | 7.8331 | 7.5363 | 7.4094

###### Remark 4

When $X$ and $Y$ are two random variables that modeling the returns of two
shares, having a extreme value copula, one can compute the discordance
function for more information about the dependance between these variables.
For more details on this measure, see Simplice and al. [15].

##  4 Conclusion and discussions

This paper proposes an approach that allows to take in amount the effect of
extreme values in the marginal and joint distribution of the underlying for
the valorisation of multivariate options. For doing so, at first, each
marginal distribution of the returns of any underlying asset is modelled by a
mixture of gaussiens as in Idier et al.[7] and the dependence structure is
modelled by a copula. The choice of the best copula is confirmed by fitting
and goodness test fit.
An application is made on the basket (Atos, Dassault systems) of the financial
market CAC40 reveal that the Tawn’s copula is the best for modelling the
dependence structure of their returns. Thus, the prices of four type of
options are calculated by use of Monte-Carlo simulation. The simulations
results show that the normale copula over-estimate the prices for the call on
maximum and the spread option when they are at-the money. In the case of
digital and call on minimum, this copula under-estimate the prices when the
options are is at-the money.

## References

  * [1] Black, F. and Scholes, M.  The Pricing of Options and Corporate Liabilities.  The Journal of Political Economy, 1973. DOI: 10.1086/260062
  * [2] Bouyé, E.  Copulas for Finance: A reading guide and some applications.  Groupe crédit Lyonnais, 2000.
  * [3] Breeden, D. T. and Robert, H. Litzenberger.  Prices of State-contingent Claims Implicit in option prices.  The Journal of Business, Vol.51, num. 4, Oct. 1978.
  * [4] Cherubuni, U. and Luciano, E.  Multivariate Options Pricing with Copulas.  Internationnal Centre for Economic Reearch, 2000.
  * [5] Cherubuni, U. and Luciano, E.  Pricing Vulnerable Options With Copulas.  The Journal of Risk Finance, 2003. http://dx.doi.org/10.1108/eb022977
  * [6] Cox, J. and Ross, S.A. and Rubinstein, M.  Option Pricing: A Simlified Approach.  Journal of Financial Economics, 1979. https://doi.org/10.1016/0304-405X(79)90015-1
  * [7] Idier, J. and Jardet, C. and Gaëlle le Fol and Alain Monfort and Fulvio Pegoraro.  La prise en compte des événements extrêmes pour la valorisation d’options européennes.  Revue de la stabilité financière, 2008.
  * [8] Joe, H. and Xu, James J.  The estimation method of Inference Functions for Margins.  Tecnical Report, The University of British Columbia Vancouver, Canada .
  * [9] Johnson, H.  Options on the minimum or the maximum of several assets.  Journal of Financial and Quantitative Analysis 22: 277–283, 1987.
  * [10] Nelsen, R.  An introduction to Copulas.  Springer, 2006.
  * [11] Rosenberg, J.V.  Semiparametric pricing of multivariate contingent claims.  Working Paper, New York University, New York, August 1999.
  * [12] Rosenberg, J.V.  Nonparametric Pricing of Multivariate Contingent Claims.  Federal Reserve Bank of New York Staff Reports, no. 162, March 2003.
  * [13] Salmon, M. and Christoph Schleichere.  Pricing Multivariate Currency Options with Copulas.  Warwick Business School and Bank of England, 2000.
  * [14] Shih, J. H. and Thomas A. Louis.  Inferences on the Association Parameter in Copula Models for bivariate Survival Data.  International Biometrics Society , Vol.51, num. 4, Dec. 1995, pp. 1384-1399.
  * [15] Simplice Dossou-Gbete, Blaise Somé , Diakarya Barro -(2009).  Modelling the Dependence of Parametric Bivariate Extreme Value Copulas.  Asian Journal of Mathematics $\&$ Statistics,2(3), 41-54. http://dx.doi.org/10.3923/ajms.2009.41.54
  * [16] A. Sklar.  Fonctions de répartition à n dimensions et leurs marges.  Publ. Inst. Statist. Univ. Paris 8, 229-231, 1959.
  * [17] Slavchev, I. and Wilkens, S.  The valuation of multivariate equity options by means of copulas: Theory and application to the European derivatives market.  Journal of Derivatives and Hedge Funds Vol. 16, 4, 303–318, 2009. doi:10.1057/jdhf.2010.5
  * [18] Stulz, R.M.  Options on the minimum or the maximum of two risky assets: analysis and applications.  J. Finan. Econ. 10, 161-185, 1982.
  * [19] Talponen, J. and Viitasaari, L.  Note on Multidimensional Breeden-Litzenberger Representation for State Price Densities.  Springer, math. Finan. Econ., Janvier 2014. DOI 10.1007/s11579-014-0113-5
  * [20] van den Goorbergh, R. W. J. and Genest, C. and Werker, B. J. M.  Multivariate Option Pricing Using Dynamic Copula Models.  , 2004.


# Quartic solitons of a mode-locked laser distributed model

D. Malheiro  Departamento de Física, Universidade de Aveiro, Campus
Universitário de Santiago, 3810-193 Aveiro, Portugal  M. Facão  Departamento
de Física, Universidade de Aveiro, Campus Universitário de Santiago, 3810-193
Aveiro, Portugal  I3N-Aveiro,Campus Universitário de Santiago, 3810-193
Aveiro, Portugal<EMAIL_ADDRESS>M. I. Carvalho  DEEC/FEUP and INESC TEC,
Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal

###### Abstract

Dissipative quartic solitons have gained interest in the field of mode-locked
lasers for their energy-width scaling which allows the generation of
ultrashort pulses with high energies. Pursuing the characterization of such
pulses, here we found soliton solutions of a distributed model for mode-locked
lasers in the presence of either positive or negative fourth order dispersion
(4OD). We studied the impact the laser parameters may have on the profiles,
range of existence and energy-width relation of the output pulses. The most
energetic and narrowest solutions occur for negative 4OD, with the energy
having an inverse cubic dependence with the width in most cases. Our
simulations showed that the spectral filtering has the biggest contribution in
the generation of short (widhts as low as 39 fs) and very energetic (392 nJ)
optical pulses.

††journal: opticajournal

Quartic solitons (QS) are shape-preserving pulses that occur in media with
Kerr nonlinearity and fourth order dispersion (4OD), in the presence or in the
absence of second order dispersion (group velocity dispersion - GVD). The
conservative quartic solitons exist for negative 4OD and have been first
studied in the 1990s [1, 2, 3] However, the interest on these solutions has
increased more recently. Quartic solitons in silico-based slot waveguides were
explored in [4], soliton-like solutions were observed in photonic crystal
waveguides with negative 4OD and negligible GVD (normal or anomalous) [5]. Tam
et al. [6, 7] have compiled the previous results of conservative models,
including results in [8], by studying the nonlinear Schrödinger equation (NLS)
plus negative 4OD, showing regions of existence of solitons with exponentially
decaying tails, purely exponential or with oscillations. In 2017-19, quartic
solitons were studied in dissipative Kerr optical cavities [9, 10] and in
2020, a mode-locked laser using a cavity dispersion dominated by 4OD was
demonstrated [11], thus opening the chapter of dissipative quartic solitons.
Contrary to the conservative QS, the dissipative QSs have been shown to exist
for both negative [11, 12] and positive [13] 4OD.

One of the advantages that has been pointed for the QS in ultrashort laser
applications is its energy, $E$, inversely proportional to the third power of
temporal width, $\tau$. This would yield very energetic short pulses. However,
this energy-width law should be only universally valid for the conservative
QSs. In fact, recent works on dissipative QSs have reported two different
energy-width relations for the quartic soliton laser, Runge et al. [11] found
$E\propto\tau^{-3}$ for negative 4OD and Qian et al. [13] found
$E\propto\tau^{3}$ for positive 4OD, both by varying the gain saturation
energy. These contradictory results could be attributed to the sign of 4OD,
but further studies are needed to establish if the law $E\propto\tau^{-3}$ is
always valid for dissipative QSs in the presence of negative 4OD. Indeed, it
is well known that the characteristics of dissipative solitons are fixed by
the equation parameters, an argument that was already referred for quartic
solitons in mode-locked lasers in [13].

The works that have studied the QS laser, referred above, have used a lumped
model, which has the advantage of a more accurate modelling of the real laser.
The cubic quintic complex Ginzburg-Landau equation (CGLE) has also been
extensively studied in connection with passively mode-locked lasers [14, 15,
16, 17], with the advantage of having a lower number of parameters and
allowing for semi-analytical approaches. Between those two approaches there is
a distributed model, proposed in the works of Zaviyalov et al., which does not
approximate the saturable absorber (SA) term that is present in the lumped
model [15, 16]. Here, we exploit the QS solutions of this distributed model
including positive and negative 4OD, showing parameter ranges of existence,
presenting their amplitude profiles, the dependence of the energy on temporal
width and searching for effective ways to maximize pulse energy while
simultaneously minimizing the width.

Let us consider an equation similar to equation (21) of [16], representing a
distributed model for a mode-locked laser

$i\frac{\partial W}{\partial
z}-\frac{1}{2}\left(\beta_{2}+ig_{0}T_{2}^{2}\right)\frac{\partial^{2}W}{\partial
t^{2}}+\frac{\beta_{4}}{24}\frac{\partial^{4}W}{\partial t^{4}}=\\\
=\left(\frac{g_{0}-k_{\text{OC}}/L}{2}\right)W-\frac{i}{2}\frac{\delta_{0}L_{\text{SA}}/L}{1+|W|^{2}/\bar{P}_{\text{sat}}}W-\bar{\gamma}|W|^{2}W,$
(1)

where $t$ is the retarded time, $z$ is the propagation distance, $W(z,t)$ is
the slowly varying pulse envelope, $\beta_{2}$ and $\beta_{4}$ are the second
and fourth order dispersion parameters (the latter is introduced here),
$g_{0}$ is the small signal gain, $T_{2}$ is the inverse linewidth of the
parabolic gain and spectral filtering, $k_{\text{OC}}$ represents losses of
the output coupler, $L$ is the cavity length and $\delta_{0}L_{\text{SA}}$ is
the modulation depth of the SA, with $\delta_{0}$ representing the small
signal loss and $L_{\text{SA}}$ the length of the SA. The parameters
$\bar{\gamma}$ and $\bar{P}_{\text{sat}}$ are effective parameters of the
distributed model associated with the nonlinear parameter $\gamma$ and
saturation power $P_{\text{sat}}$ of the SA, respectively, and given by
$\bar{\gamma}=\gamma(\exp(g_{0}L)-1)/g_{0}L$,
$\bar{P}_{\text{sat}}=P_{\text{sat}}\exp(-g_{0}L)$. In this work we focus on
the role played by $g_{0}$, $P_{sat}$, $T_{2}$ and $\beta_{4}$. All the
following results were obtained with fixed $\delta_{0}L_{\text{SA}}=0.3$,
$\gamma=0.005$ W-1m-1, $L=1$ m and $k_{\text{OC}}=-\text{ln}(0.3)$.

We searched for stationary solutions of the form
$W(z,t)=F(t)e^{i\theta(t)+i\sigma z}$, where both $F$ and $\theta$ are real
and $\sigma$ is a propagation constant. Introducing this ansatz into (1)
yields a system of nonlinear ordinary differential equations, which was solved
using the Newton conjugate gradient method [18, 19] with hyperbolic secant
inputs as well as soliton solutions from previous simulations. To verify that
the solutions were stable, they were taken as inputs in (1) which was
subsequently integrated through a pseudo-spectral method [20] evolving in $z$
with a fourth order Runge-Kutta method, up to distances where the pulse was
stationary. This approach is not very efficient in the identification of
parameter regions where stable solutions exist and thus, we also calculated
the eigenvalues of the linear stability operator for stationary solutions
under (1).

Stationary quartic soliton solutions were found for both positive and negative
4OD, with different signs of $\beta_{4}$ leading to different pulse shapes. As
an example, in Fig. 1 the pulse profiles for $\beta_{4}=\pm 0.012$ and
$0\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$ when
$\beta_{2}=-0.024\leavevmode\nobreak\ \text{ps}^{2}\text{m}^{-1}$,
$g_{0}=1.45\leavevmode\nobreak\ \text{m}^{-1}$, $T_{2}=100\leavevmode\nobreak\
\text{fs}$ and $P_{\text{sat}}=80\leavevmode\nobreak\ \text{W}$ are
represented. The solution in the absence of 4OD takes the approximate shape of
a hyperbolic secant, as is typical in GVD only solitons. For negative 4OD, the
pulses are gaussian with symmetrically oscillating tails, which are
particularly noticeable in logarithmic scale (Fig. 1). For positive 4OD, the
pulses come in the shape of a hyperbolic secant with broader tails. In
logarithmic scale, these tails take the shape of straight lines, showing that
they are, in fact, exponentially decaying.

Figure 1: Pulse profiles in 1 linear and 1 logarithmic scale with negative,
null and positive 4OD. The inset shows more detail in the pulse shape for
$\beta_{4}\leq 0$. Pulse energy 1 and energy-width relation 1 obtained from
sweeping $\beta_{4}$ from $-0.12$ to $0.12\leavevmode\nobreak\
\text{ps}^{4}\text{m}^{-1}$. All solutions were obtained using
$\beta_{2}=-0.024\leavevmode\nobreak\ \text{ps}^{2}\text{m}^{-1}$,
$g_{0}=1.45\leavevmode\nobreak\ \text{m}^{-1}$, $T_{2}=100\leavevmode\nobreak\
\text{fs}$ and $P_{\text{sat}}=80\leavevmode\nobreak\ \text{W}$.

The pulse energy is also greatly influenced by the value of the 4OD
coefficient. Fig. 1 shows that the energy increases with
$\left|\beta_{4}\right|$, with the growth rate being greater when
$\beta_{4}<0$. Fig. 1 contains the energy-width relation of the pulses
obtained by sweeping $\beta_{4}$ and keeping all other parameters constant.
The pulses obtained with $\beta_{4}<0$ are, in general, the most energetic as
well as the narrowest, with pulses as short as 300 fs having energies close to
1.2 nJ. For negative 4OD there is also a very sharp decrease of energy with an
increase of the pulse width. In opposition, the energy increases with the
width for positive $\beta_{4}$ but at a very slow rate. Therefore, the
exploitation of negative 4OD in mode-locked laser cavities is much more
favorable to achieve higher energies at lower pulse widths.

Figure 2: Pulse shapes for different $T_{2}$ for 2, 2
$\beta_{4}=-0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$,
$P_{\text{sat}}=3\leavevmode\nobreak\ \text{W}$ and
$g_{0}=1.36\leavevmode\nobreak\ \text{m}^{-1}$, and 2, 2
$\beta_{4}=0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$,
$P_{\text{sat}}=80\leavevmode\nobreak\ \text{W}$ and
$g_{0}=1.48\leavevmode\nobreak\ \text{m}^{-1}$, in linear and logarithmic
scale, respectively. In both cases, $\beta_{2}=0$.

The sign of 4OD is not the only factor defining the shape of the pulses. All
of the parameters under investigation here were found to have some effect on
the pulse shape, amplitude and width. Fig. 2 shows the impact of growing
$T_{2}$ in the shape of pulses obtained with
$\beta_{4}=-0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$. For
$T_{2}=100\text{ fs}$ (blue solid curve), the pulse takes the same shape as
the one with negative 4OD in Fig. 1. As $T_{2}$ increases to 150 fs, and
further so to 200 fs (red dashed and green dash-dot curves respectively), the
oscillations in the tails begin to flatten. For $T_{2}=250\leavevmode\nobreak\
\text{fs}$ the pulse becomes hyperbolic secant shaped, maintaining this shape
for all greater $T_{2}$ values where stationary solutions exist. The
flattening of the tail oscillations is further evidenced in logarithmic scale
(Fig. 2), with the oscillations in the tails becoming less and less noticeable
until the tails become approximately straight lines. A similar behaviour
exists in the presence of positive 4OD. In Fig. 2, considering
$\beta_{4}=0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$,
$g_{0}=1.48\leavevmode\nobreak\ \text{m}^{-1}$,
$P_{\text{sat}}=80\leavevmode\nobreak\ \text{W}$ and
$T_{2}=115\leavevmode\nobreak\ \text{fs}$, a pulse with exponentially decaying
tails is obtained, similar to the one in Fig. 1. Changing $T_{2}$ to 200 fs
results in a hyperbolic secant shaped pulse. Fig. 2 shows the different pulse
shapes in logarithmic scale, highlighting a much more rapid decay in the tails
for $T_{2}=200\leavevmode\nobreak\ \text{fs}$.

Figure 3: Energy dependence with $g_{0}$ (top row), $P_{\text{sat}}$ (middle
row) and $T_{2}$ (bottom row) for $\beta_{4}=-0.08\leavevmode\nobreak\
\text{ps}^{4}\text{m}^{-1}$ (left) and $\beta_{4}=0.08\leavevmode\nobreak\
\text{ps}^{4}\text{m}^{-1}$ (right). When constant,
$P_{\text{sat}}=80\leavevmode\nobreak\ \text{W}$ and
$T_{2}=100\leavevmode\nobreak\ \text{fs}$ were considered in all cases, while
$g_{0}$ took values of $1.36\leavevmode\nobreak\ \text{m}^{-1}$ for
$\beta_{4}<0$ and of $1.45\leavevmode\nobreak\ \text{m}^{-1}$ for
$\beta_{4}>0$.

To design a mode-locked fiber laser capable of emitting high energy pulses, it
is important to find how the different laser parameters impact the pulse
energy. In Fig. 3, the energy dependence with $g_{0}$, $P_{\text{sat}}$ and
$T_{2}$ is plotted for $\beta_{4}<0$ (Figs. 3, 3, 3, respectively) and
$\beta_{4}>0$ (Figs. 3, 3, 3, respectively). It immediately becomes apparent,
reinforcing the results from Fig. 1, that the pulses with highest energies
occur in the negative 4OD regime. In fact, the highest energy obtained in a
$g_{0}$ sweep for negative 4OD is about 9 times greater than the one found for
positive 4OD. More striking is the case of the $T_{2}$ sweep, where the
highest energy with $\beta_{4}<0$ is approximately 3 orders of magnitude
larger than the one for $\beta_{4}>0$. It is also clear that the highest
energies occur for higher $g_{0}$ values as well as low $P_{\text{sat}}$ and
$T_{2}$ values. It is important to note however, that the choice of just one
of these parameters is enough to completely alter the range where the other
parameters allow stationary solutions. As an example, in Fig. 3, with a
saturation power of 3 W, stationary solutions were found from $T_{2}=50$ fs to
$T_{2}>1$ ps. However, by simply increasing $P_{\text{sat}}$ to 80 W, the
maximum $T_{2}$ value reached was of approximately 140 fs, and became of $\sim
110$ fs when $P_{\text{sat}}=160\leavevmode\nobreak\ \text{W}$ was considered
(see inset). The use of positive 4OD also leads to a much shorter valid
$g_{0}$, $P_{\text{sat}}$ and $T_{2}$ ranges for the existence of stationary
soliton solutions of (1). Finally, Fig. 3 also shows that the different
parameters also have very different impacts in the output pulse energy. For
instance, with $\beta_{4}<0$, very small changes for $T_{2}\lesssim 150$ fs,
lead to drastic variations on the pulse energy, while small increases in
$P_{\text{sat}}$ do not necessarily lead to very significant energy decreases.
In fact, in Fig. 3, the change in $P_{\text{sat}}$ altered the soliton
stability range, but changes in energy were not very significant.

As noted above, quartic soliton lasers have gained interest due to their
ability to generate highly energetic pulses with short widths. In the
distributed-model, it was found, in most cases, that the energy does scale
inversely with the width cubed, in agreement with previous studies on the
subject using other models [6, 11]. In Fig. 4 the energy-width relation of the
soliton solutions of (1) are presented, with Fig. 4 containing the results for
the sweeps of $g_{0}$, $P_{\text{sat}}$ and $T_{2}$ with negative 4OD
respectively, while in Fig. 4 the same curves are plotted for positive 4OD.
Interestingly, most of the energy-width points for negative 4OD fall
approximately in the same curve that, in fact, follows the relation
$E\propto\tau^{-3}$ (see gray-dashed line). The exception is for one of the
branches occurring for $P_{\text{sat}}=3$ W while $T_{2}$ is swept above 150
fs, for which the energy increases with width. For positive 4OD however, this
relation is not verified, nor is the cubic dependence predicted in [13]. In
fact, our curve fittings to data for $\beta_{4}>0$ yielded different relations
as may be observed by the dashed lines of Fig. 4 that were obtained by fitting
the data to the function $a\tau^{b}$ (where $a$ and $b$ are fitting
coeficients). The pulses for $\beta_{4}>0$ not only have lower energies but
they are also much wider. In this case, the shortest and most energetic pulses
have widths close to 1.6 ps and energies close to 100 pJ, while the pulses
with $\beta_{4}<0$ reach energies in the nJ range, and widths close to 300 fs.

Figure 4: Energy-width curves using logarithm base 10 scale and associated
fittings, obtained by sweeping $g_{0}$, $P_{\text{sat}}$ and $T_{2}$ for
$\beta_{4}=-0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$ 4 and
$\beta_{4}=0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$ 4. When
constant, values of $T_{2}=100\leavevmode\nobreak\ \text{fs}$ were considered
, while $g_{0}$ took values of $1.36\leavevmode\nobreak\ \text{m}^{-1}$ for
$\beta_{4}<0$ and of $1.45\leavevmode\nobreak\ \text{m}^{-1}$ for
$\beta_{4}>0$ and $P_{\text{sat}}$ is 80 W in all cases not referred in the
legends.

The results presented in Fig. 3 indicate how a single parameter, $g_{0}$,
$P_{sat}$ or $T_{2}$, could be manipulated to obtain pulses with higher
energy. Nevertheless, the maximization of the pulse energy must take into
account the simultaneous variation of all these parameters. Not only is the
dependence of the pulse energy on these parameters quite different but, more
importantly, stable solutions only exist within a certain parameter region. In
effect, an increase of $P_{\text{sat}}$ implies that the SA saturates for
higher power and thus, the nonlinear gain decreases. For a higher value of
$P_{\text{sat}}$, the nonlinear gain may not be sufficient to counter balance
other losses and thus, soliton solutions are not found. This same effect also
occurs with lower $g_{0}$ values, which lead to higher effective saturation
powers, $\bar{P}_{\text{sat}}$. For higher $g_{0}$ values, the linear losses
($-g_{0}/2+k_{\text{OC}}/2L+\delta_{0}L_{\text{SA}}/2L$) tend to zero and
linear gain would produce unstable background. On the other hand, lower values
of $P_{\text{sat}}$ lead to a sharp increase in the nonlinear gain which would
be unbalanced by the other parameters, also leading to unstable solitons.

In order to have a better understanding of the combinations of these
parameters that lead to highly energetic pulses with short widths, we searched
for solutions and calculated stability eigenvalues on a $g_{0}-T_{2}$ region
for $\beta_{2}=0\leavevmode\nobreak\ \text{ps}^{2}\text{m}^{-1}$ and
$\beta_{4}=-0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$. Fig. 5 shows
the energy (in logarithmic scale) and width (fs) values as contour plots for
the stable solutions when $P_{\text{sat}}$ is $80\leavevmode\nobreak\
\text{W}$ and $160\leavevmode\nobreak\ \text{W}$. Note that, despite the
energy increasing with the increase of $g_{0}$ for fixed $T_{2}$ and
$P_{\text{sat}}$, which is in agreement with Fig. 3, lower values of $g_{0}$
allow solutions to exist and evolve steadily for lower values of $T_{2}$, and
solutions in that region have the highest energy. Similarly, even though for
fixed $T_{2}$ and $g_{0}$, a decrease of $P_{\text{sat}}$ was favourable for
higher energy, if $T_{2}$ and $g_{0}$ are allowed to vary we may achieve
higher energies for higher values of $P_{\text{sat}}$. When it comes to the
widths, represented in Fig. 5, it is shown that the pulse width decreases as
$T_{2}$ decreases. $T_{2}$ being the parameter of the spectral filtering is
influencing the pulse shape and, as expected, the pulse width is always
changing in the same direction as $T_{2}$. Therefore, the most energetic
pulses are also the shortest ones. In fact, using $T_{2}=7.5$ fs,
$g_{0}=1.3\leavevmode\nobreak\ \text{m}^{-1}$ and
$\beta_{4}=-0.08\leavevmode\nobreak\ \text{ps}^{4}\text{m}^{-1}$ yielded
quartic solitons with energy of 392 nJ and 39 fs width. These energy values
are much greater than the ones presented in previous works [6, 10, 11].

Figure 5: Stability border solitons in the $g_{0}-T_{2}$ plane. Above the
curve, 5 and 5 the base 10 logarithm of the energy in nJ, in 5 and 5 the pulse
width in fs are represented in a contourplot. Equation parameters:
$\beta_{2}=0$, $\beta_{4}=-0.08\leavevmode\nobreak\
\text{ps}^{4}\text{m}^{-1}$ and $P_{\text{sat}}=80\leavevmode\nobreak\
\text{W}$ (top row) and $P_{\text{sat}}=160\leavevmode\nobreak\ \text{W}$
(bottom row).

In conclusion, we found quartic soliton solutions of a mode-locked laser
distributed model in the presence of positive and negative and 4OD, as well as
in the absence or in the presence of negative GVD. The selection of laser
parameters is of paramount importance for the existence and stability of
solitons, as well for defining their shape. Most of profiles for negative 4OD
have oscillations in their tails and for positive 4OD are hyperbolic secant in
the center but upon a larger pedestal. Nevertheless, there exist profiles very
similar to sech for both negative and positive 4OD if larger spectral
filtering $T_{2}$ is considered. We studied the energy-width relation of the
different solutions, finding that the energy dependence with the width varies
with laser parameters. With positive 4OD, our results did not show a single
universal trend, but with negative 4OD, most of our results followed an
$\text{energy}\propto\text{width}^{-3}$ trend. The exceptions were for the
sweep in $\beta_{4}$ as well as the sweep in $T_{2}$ for low $P_{\text{sat}}$
, where a non-monotonous trend was found. The highest energies and shortest
widths were found for negative 4OD and for low $T_{2}$ values. To obtain
solutions for such low $T_{2}$, $P_{\text{sat}}$ should be increased and
$g_{0}$ decreased. In effect, our simulations of the distributed model show
that the solution for surpassing current energy limits may be through quartic
solitons by carefully optimizing gain bandwidth and spectral filtering within
the available gain and saturation power levels.

Funding Content in the funding section will be generated entirely from details
submitted to Prism. Authors may add placeholder text in the manuscript to
assess length, but any text added to this section in the manuscript will be
replaced during production.

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

## References

  * [1] M. Karlsson and A. Höök, Optics Communications 104, 303 (1994).
  * [2] N. Akhmediev, A. Buryak, and M. Karlsson, Optics Communications 110, 540 (1994).
  * [3] A. Buryak and N. Akhmediev, Physical Review. E 51 4, 3572 (1995).
  * [4] S. Roy and F. Biancalana, Phys. Rev. A 87, 025801 (2013).
  * [5] A. Blanco-Redondo, C. de Sterke, J. Sipe, T. Krauss, B. Eggleton, and C. Husko, Nature Communications 7, 10427 (2016).
  * [6] K. K. K. Tam, T. J. Alexander, A. Blanco-Redondo, and C. M. de Sterke, Opt. Lett. 44, 3306 (2019).
  * [7] K. K. K. Tam, T. J. Alexander, A. Blanco-Redondo, and C. M. de Sterke, Phys. Rev. A 101, 043822 (2020).
  * [8] V. I. Kruglov and J. D. Harvey, Phys. Rev. A 98, 063811 (2018).
  * [9] C. Bao, H. Taheri, L. Zhang, A. Matsko, Y. Yan, P. Liao, L. Maleki, and A. E. Willner, Journal of the Optical Society of America B 34, 715 (2017).
  * [10] H. Taheri and A. B. Matsko, Optics Letters 44, 3086 (2019).
  * [11] J. Runge, D. Hudson, K. Tam, C. Sterke, and A. Blanco-Redondo, Nature Photonics 14, 492–497 (2020).
  * [12] Z.-X. Zhang, M. Luo, J.-X. Chen, L.-H. Chen, M. Liu, A.-P. Luo, W.-C. Xu, and Z.-C. Luo, Opt. Lett. 47, 1750 (2022).
  * [13] Z.-C. Qian, M. Liu, A.-P. Luo, Z.-C. Luo, and W.-C. Xu, Opt. Express 30, 22066 (2022).
  * [14] H. Haus, J. Fujimoto, and E. Ippen, IEEE Journal of Quantum Electronics 28, 2086 (1992).
  * [15] A. Zavyalov, R. Iliew, O. Egorov, and F. Lederer, Phys. Rev. A 79, 053841 (2009).
  * [16] A. Zaviyalov, R. Iliew, O. Egorov, and F. Lederer, J. Opt. Soc. Am. B 27, 2313 (2010).
  * [17] P. Grelu and N. Akhmediev, Nature Photon 6, 84–92 (2012).
  * [18] J. Yang, Journal of Computational Physics 228, 7007 (2009).
  * [19] J. Yang, Studies in Applied Mathematics 134, 420 (2015).
  * [20] B. Fornberg and G. B. Whitham, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 289, 373 (1978).

solitoes


# Development of an interactive code for quick data analyses between STOR-M
tokamak experimental plasma discharges.

Masaru Nakajimaa,SS, Debjyoti Basua,b,§, A.V. Melnikovc,d, David McColla,
Chijin Xiaoa aPlasma Physics Laboratory, University of Saskatchewan,
Saskatoon, Canada  b Present address: University of Maryland, Baltimore
County, MD-21250, USA  cNRC ‘Kurchatov Institute’, 123182, Moscow, Russia.  d
National Research Nuclear University MEPhI, 115409, Moscow, Russia.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Saskatchewan Torus-Modified (STOR-M) is a small tokamak, well known for
various fusion related basic experimental studies such as edge turbulent
heating, different instabilities, AC (alternating current) tokamak operation,
Ohmic H-mode triggering by the electrode biasing, fueling and momentum
injection by Compact Torus(CT) injection, and effects of Resonance Magnetic
Perturbations (RMP), among others. Some of those experiments require real time
visualization of magnetic surface reconstructions either through EFIT or quick
analyses and visualization of experimental data during experiments. Recently
experimental studies of Geodesic Acoustic Mode(GAM) and zonal flows had been
performed in STOR-M tokamak. The GAM experiments strongly require collection
of fluctuations data from different Langmuir probes installed at different
poloidal locations, but on the same magnetic surfaces. This is need of the
adjustment of radial locations between discharges. It is therefore important
to analyze and visualize the features of all probe data quickly during
discharges. For this purpose, a Python code has been developed and used for
quick analyze of data. This article will describe the development of the code
using Python and its use in detail.

###### pacs:

Valid PACS appear here

##  I Introduction

Discovery of the concept of zonal flows has great scientific value in the
field of basic turbulent plasma physics and magnetically confined plasma
systems aiming to nuclear fusion in a controlled manner. A basic requirement
to achieve self-sustaining burning plasmas in a magnetically confined fusion
reactor is to control anomalous transportation of particle and energy. It is
believed that anomalous transport happens due to small scale drift turbulent
fluctuations in a magnetically confined toroidal system like tokamaks or
stellarators. The anomalous transport can be controlled by radially sheared
poloidal flow which is generated by either radially sheared mean or
oscillating electric field. The poloidal mean flow may be generated by
external means but in general the oscillating poloidal flows, confined within
a small radial width, is a self-organized process which is known as ’zonal
flows’. Theoretically and experimentally it has been shown that the drift wave
turbulence is a sole cause of the origination of zonal flows. Two types of
zonal flows have been found in a magnetically confined toroidal plasma system
which are i)near zero or very low frequency Zonal Flows(ZF) and ii)Geodesic
Acoustic Mode(GAM) which has been found only in toroidal geometry. The
conventional GAM is featured by i) radially localized, ii) toroidally
axisymmetric for potential and density fluctuations $(n=0;m=0)$ and
$k_{\phi}=0$, iii) poloidally symmetric potential fluctuations, so, $(n=0)$,
$k_{\theta}=0$ for potential fluctuations but asymmetric density fluctuation
$(m=1)$ and $k_{\theta}=\pi$ for density fluctuations kn:Diamond , kn:Fujisawa
.

Since, the GAM mode is originated in a particular radius and normally confined
within a narrow radial width, the experimental detection of GAM mode demands
that diagnostics should collect data from the same magnetic surfaces where the
mode is confined. The main issue of the tokamaks which are not facilitated by
good feedback systems is the reproducibility of particular shot in a same
magnetic surface with same characterizations. Similarly, the radial positions
of GAM formation of a particular Tokamak plasma discharge may change shot by
shot. In that case, the detection of radial locations of GAM mode should be
done during real time operations through quick visualization of experimental
data of the density and potential fluctuations whether they have been
collected from the same magnetic surface by the help of highly advanced
software like Equilibrium Fit(EFIT)kn:lao . A high performance Python code has
been developed for quick analysis the data taken by diagnostic system and
derive mode structure between consecutive plasma discharges in STOR-M due to
absence of EFIT.

In this article, detailed description of code development, preliminary
experimental data analysis through the code during STOR-M discharges and brief
description of mode detection will be discussed.

##  II Experimental Setup

The STOR-M with circular plasma cross-section has major and minor radii 46cm
and 12cm respectively. The plasma discharge duration is around 30 ms. The
plasma boundary is defined by an elliptical stainless steel poloidal limiter,
placed in a way that its internal major $(13cm)$ and minor $(12cm)$ radii have
aligned along horizontally and vertically. Three Langmuir probe systems(top,
radial-2, bottom) were installed through the top, equatorial and bottom ports
of the vessel in the same poloidal cross-section. A Fourth one(radial-1) is
placed equatorially, offset toroidally by $90^{0}$ in clockwise direction (top
view) with respect to the other three sets. Detailed description of
diagnostics arrangement can be found from the article of ‘GAM detection in
STOR-M tokamak’kn:DB . Now, the necessary and very important requirement is to
determine that all Langmuir probes will collect data from the same magnetic
field surface for the GAM detection experimental studies. The main purpose of
this code development is to quick determination through fast analysis whether
those Langmuir probes data are collected from the same magnetic surface. The
detailed descriptions of the code development are given in following sections.

##  III Code Development

This code utilize the Graphical User Interface(GUI) feature in Python
language, specifically for the STOR-M Tokamak experiments where probes data
are recorded and stored in National Instruments(NI) digitize modules plugged
in to a ‘PXI’ crate and directly to the PSI slots in computers. There are
around 120 channels in total. LabView software is used to manage all hardware
and to transfer onboard data in digitizers to computers and save them in to
the data files in the hard drives of computers after each shot. In
experimental settings, the duty cycle time period between two discharges can
be set to 4-10 minutes depending on particular experiments whereas the time to
transfer and save data complete within a few seconds. It is imperative to
quickly visualize and analyze the experimental data after every recorded
STOR-M discharge for the complex experiment like GAM detection to determine
whether satisfactory data were obtained, or the configurations of the
diagnostics set up needs to be modified. The present development of code with
interactive data analysis tools meet such requirements. In this section the
logical frame of the application will be discussed. The resultant GUI window
is presented by figure 1. The work flow chart of this application is shown in
figure 2(a).

Figure 1: It presents the image of GUI window of the code. Figure 2: It
presents(a) work flow chart, (b)sequential organization of data store.

It describes how the programs run when the large amount of diagnostics data
related to the experiment are store periodically after each STOR-M discharge
in the pre-specified file in directories. So, the each step of the operations
and their inter connections from $BrowseData\Rightarrow LoadData\Rightarrow
VisualizeData$, the next level of analytical operations, including basic
plots, ‘Power Spectral Density’, ‘Cross Spectral Density’, and ‘Bi-coherence’
are controlled by dedicated background programs which were solely developed.
Others available applications which were used in a logical manner are
graphical interface through the PyQt package, a wrapper package for a GUI tool
kit Qt kn:pyqt ,the plotting and visualization of data through PyQtGraph
kn:pyqtgraph , a python toolkit for interactive data visualization and the
package SciPy for spectral analysis. Figure 2(b) represent the sequence of
data arrangement consecutively.

Each STOR-M plasma discharge is labeled by a shot number as a common
experimental procedure for data storage. The data files for National
Instrument digitizers are stored in the text files format with unique number
synchronized with shot number corresponding to each set of digitizers within
it. It is very important to develop a logical working procedure in the code
which includes its methods of reading, interacting and operating on the stored
experimental data. So, it requires a guiding soft file which can lead the code
for its working flows. Here, the Excel book is the pilot file from which the
code collects the instructions for loading data files. Each Excel pilot files
were named in such a way that code can understand the range of shots number
stored on the particular folder. The pilot files contain the descriptions for
each measurement channel names, conversion factors, and the final unit of the
measurement. Since multiple shots share a same measurement setup, the file
names of the DAQ information files indicate the range of shot numbers for
which the files are applicable. The data files are organized as shown in the
Figure 2(b) where it shows the pilot excel books and the data files within a
folder.

Figure 3: It presents flow charts for (a)browsing data, (b)loading data, (c)
generating basic plot.

The code starts to load experimental data in the GUI window shown in figure 1
through navigation of the proper directory where the experimental data of
interests were stored. This process is done through the logical flows of
‘browse data’, shown in figure 3(a). Figure 3(b) explains the logical flow
chart of the ‘load data’ in GUI window through interaction with pilot Excel
book which is one of the key process of the code. Once the shot numbers are
available, the shots of interest need to be ”loaded” to the program through
GUI window. The program finds all the data files associated with the selected
shot numbers. All the files of the selected shots number of interests will be
opened and stored in a table, particularly in the object pandas.DataFrame
which is called the data table, where each time series data is given a unique
data ID. After that the program opens the appropriate Excel book and extracts
information about the data such as unit, measurement name, and calibration
factor for each data time series. The parameters related to each data ID are
stored in another table which is called the parameter table. The program
completes the loading process by listing the shot numbers and measurement
names or channels. After proper loading of the data files under analysis,
figure 3(c) shows the work flow for the basic plots and visualization and the
nature of experimental data.

It is very important to analyze the fundamental signals of various channels to
check its proper nature. The program retrieve the information of relevant data
and parameters from the tables of data and parameter. The optional analysis,
such as smoothing factor, are also possible through graphical interface. An
interactive plot is then generated using a plot item in pyqtgraph package.
Based on the parameters, the program updates the unit, axis labels and the
legends.

Here, the most important feature of the code is successful using the in built
ROI tool of pyqtgraph which allows to indicate region of interest (ROI) and
corresponding analysis of the focussed region of data. The ROI can be changed
in an interactive fashion where it can be changed either by modifying the
width of ROI or moving the ROI without changing its width. The corresponding
new plot is then updated automatically. This mapping is done through a
function to a signal emitted by the ROI object. In this particular case, the
connected function simply plots the data within the selected region. In our
application, the ROI tool is a very important feature to carry out spectral
analysis which allows quick and detailed analysis of modes present in the
signals. The ROI tool have been used to efficiently identify GAM mode within a
time series.

To explain GAM features, Power spectral density (PSD) and Cross-power Spectral
analysis such as Cross Spectral Density(CSD), Cross Coherence and Cross Phase
are essential analyzing parameters which are adopted in the code through GUI
taking from Python library. Figure 4(a) and 4(b) present the logic flow chart
of the PSD and Cross spectral analysis.

Figure 4: Flow charts of (a)Power Spectral Density(PSD) and (b)Cross Spectral
Analysis.

The analysis options for PSD such as segment length and overlap percentage are
read from the graphical interface. A default time window is set and visualized
as a colored ROI object on the plot of the time series. The PSD is computed
using a function in scipy package and plotted in a separate window within a
plot. The time window, an object in pyqtgraph package, can be changed in an
interactive manner. Once the time window is changed, the PSD for the new time
range is computed and updated for the corresponding plot.

Cross spectral analysis is very useful for identifying common modes in two
time series signals in GAM detection experiment. For this process, two time
series signals, collected from two different independent channels are
compared. The program also reads the analysis options from the graphical
interface like ‘PSD’ plot. Initially, the two time series data are plotted. A
default time window is set and displayed as a segment item on the time series
data, similar process for ‘PSD’ plot. The cross spectral density, cross-
coherence, and cross-phase are computed for the time window and plotted in
three separate windows within the plot. When the time window is changed, the
above three quantities are computed again for the new time range and the
corresponding plots are updated.

##  IV Results and discussions:

The previous section discussed about the how code works. Figure 5 represents
one of the typical example about the basic plots and the corresponding plots
of ‘PSD’, ‘CSD’, ‘Cross-coherence’, Cross-phase’ at ROI of the basic data
plots. Figure 5(a) presents the ‘PSD’ of the data regime with ROI of the basic
plot of ion saturation current of Langmuir probe, as shown in figure 5(b).
Similarly, figure 5(c), 5(d) and 5(e) present the ‘CSD’, ‘Cross-coherence’,
Cross-phase’ at ROI data regime of the two different ion saturation current
signals collected from two different Langmuire Probes. Here, the ROI is marked
by the gray rectangular window on the basic plots of figure 5. So, it is clear
from the typical example of figure 5 that the mode can be easily determine in
a time series by changing the ROI positions on the basic data plots and it is
very powerful tools for GAM detections.

Figure 5: It presents (a)Power Spectral Density(PSD) of ROI of $LP-I_{Sat}$
signal, (b)corresponding fundamental $LP-I_{sat}$ signal, (c)CSD, (d) cross-
coherence, (e) cross-phase of ROI of two different $LP-I_{Sat}$ signals,
(f)corresponding two $LP-I_{sat}$ signals.

It is well known that GAM detection experiment is a very sophisticated due to
its complexity and on the other hand it is very important to understand
inherent turbulent physics as well as L-H transition in Tokamaks or
Stellarators. It is basically a high frequency zonal flows kn:Diamond ,
kn:Diamond , kn:Fujisawa , kn:DB , kn:Melnikov1  and normally confine within a
small radial width. Also, recently it is found that it can be a global
eigenmode as well which exhibits long-range spatial correlationkn:Melnikov13
but not in all Tokamaks or Stellarators. So, it is necessary to detect
experimentally the GAM mode collecting the information from approximately same
magnetic field surfaces. The same or nearly same magnitude of the $I_{sat}$
signals profiles, collected from top and bottom Langmuir Probes during the
same plasma discharge serve as a criterion to indicate that all measured
density fluctuations or floating potentials fluctuations are accumulated from
same or nearly same magnetic field surface.

Figure 6: An indication of GAM mode with frequency width 31-37MHz appears
within ROI having time window 22ms as indicated by (a)Power Spectral
Density(PSD) of ROI of $LP-I_{Sat}$ signal, (b)corresponding two $LP-I_{Sat}$
signals, (c)CSD, (d) cross-coherence, (e) cross-phase of ROI of two different
$LP-I_{Sat}$ signals, (f) corresponding two $LP-I_{Sat}$ signals.

At first, ROI has been taken over whole the flat top region of $I_{sat}$
signals, collected by top and bottom Langmuir probes of a typical STOR-M
plasma discharge with shot number 294348 to get an overall scenario. Figure 6
clearly indicates the presence of GAM mode at flat top where the ROI time
window is 22 ms. It shows the GAM-like mode having frequency width 30-40MHz in
the plasma discharge which is observed from ‘PSD’, ‘CSD’, ‘Cross-coherence’
and ‘Cross-phase’ which are presented in figure 6(a), 6(c), 6(d), 6(e)
respectively. Figure 6(b) and figure 6(f) presents the similar $I_{sat}$
signals collected from ‘top’ and ‘bottom’ Langmuir Probes.

Figure 7: A typical GAM mode with frequency width 31-37MHz appears within ROI
as indicated by (a)Power Spectral Density(PSD) of ROI of $LP-I_{Sat}$ signal,
(b)corresponding two $LP-I_{Sat}$ signals, (c)CSD, (d) cross-coherence, (e)
cross-phase of ROI of two different $LP-I_{Sat}$ signals, (f) corresponding
two $LP-I_{Sat}$ signals. Figure 8: GAM mode with frequency width 31-37MHz
disappears at another ROI indicated by (a)Power Spectral Density(PSD) of ROI
of $LP-I_{Sat}$ signal, (b)corresponding $LP-I_{Sat}$ signal, (c)CSD, (d)
cross-coherence, (e) cross-phase of ROI between two different $LP-I_{Sat}$
signals, (f) corresponding two $LP-I_{Sat}$ signals.

The ROI window has been set and started to move in different parts of those
two top and bottom Langmuir probes $I_{sat}$ signals to characterize and to
understand more precisely at which time GAM mode appears as well as its time
duration. Figure 7 indicates that a clear GAM mode of frequency width 31-37MHz
emerges around 18.8 ms within a time window of of width at around 2.9 ms which
is clearly marked by a gray rectangular area of ROI. Here, figure 7(b) and
figure 7(f) show similar $I_{sat}$ signals collected from ‘top’ and ‘bottom’
Langmuir probes which have been smoothed to get clear envelope of those
signals. The evolution of the envelopes of signals suggests that they are
located approximately on the same magnetic surfaces. The sharp peaks of ‘PSD’,
‘CSD’, ‘Cross-coherence’ and ‘Cross-phase’, presented in figure 7(a), 7(c),
7(d), 7(e) respectively, show clear indication of GAM-like mode when the time
range of ROI at around 18.8ms within a time window of at around 2.9 ms.

When the ROI window is moved in different locations for the same two signals,
it has been noticed that the GAM mode disappeared as shown in figure 8 for ROI
window within 22-24ms on the $I_{sat}$ signals, due to the lack of the clear
peaks of GAM-like mode in the ‘PSD’, ‘CSD’, ‘Cross-coherence’ and ‘Cross-
phase’, presented by the figure 8(a), 8(c), 8(d), 8(e) respectively.

So, it is cleared that the moving of ROI on the signals can easily detect the
GAM mode very precisely in terms of its time of appearance in the discharge as
well as its existing time duration through the quick and efficient
characterizations of the GAM mode.

##  V Conclusion and future work:

So, it is cleared from previous discussions that this code is highly efficient
to analyze and detect the sophisticated GAM-like mode. This code is user
friendly through its graphical interface. It has also options to add more
analytical options or components on its GUI window. Recently, the Bi-coherence
analysis has been added into the codes with its flexible parameters options
such that it can easily calculate higher order coherence. This addition will
be tested and validated soon. In future, the more analytical options related
to time series signal processing will be added and will be used for more
elaborate and larger scenario.

###### Acknowledgements.

We would like to acknowledge NSERC for supporting this work. We also would
like to acknowledge the kind support and inspiration provided by Prof. Akira
Hirose. The work of A.V. Melnikov was supported in part by the Competitiveness
Programme of the National Research Nuclear University “MEPhI”.

## References:

## References

  * [1] P.H. Diamond, et.al. 2005 Plasma Phys. Control Fusion 47 R35.
  * [2] A. Fujisawa, et.al. 2007 Nucl. Fusion 49 013001.
  * [3] L.L. Lao, et.al. 1985 Nucl. Fusion 25 1611.
  * [4] Debjyoti Basu, et.al. 2018 Nucl. Fusion 58 024001.
  * [5] PyQT, 2012, http://www.riverbankcomputing.com/static/Docs/PyQt4/html/index.html.
  * [6] PyQtGraph, 2021, www.pyqtgraph.org.
  * [7] A. V. Melnikov 2016 Nature Phys. 12 386.
  * [8] A. V. Melnikov et.al. 2015 Journal of Phys. Conf. Series 591 012003.