The nonlinear diffusion equations have found vast applications in various fields in physics [6], neuroscience [8, 34, 35] , psychology [7] , economy and marketing [22] [23] [24] [26] [27] [28] , biology and biophysics [5, 6, [8] [9] [10] , and population dynamics [13] . The examples of physical systems that are described by the nonlinear diffusion equation are the plasma systems [1, 2], surface physics [3, 4] , astrophysics [5, 6] , the polymers [14, 15] , fluids and particle beams [10] , liquid surfaces [16] , nonlinear hydrodynamics [25] , pattern formation [12] and laser physics [11] , and most important, in financial analysis [24] . Many aspects of the nonlinear diffusion equation have been analyzed and found, like stationary [9] and H-theorem [17, 18] , autocorrelations [19] , path integral formulation [20] , nonextensive maximum entropy approach [21, 31] , the distributed approximating functional method [29] , the associated entropy [30] , and anomalous diffusion [32] , for a good review see [33] . Due to the broadness of the problems involving anomalous diffusion, one needs to apply different kinds of theoretical approaches,such as the porous media equation (PME). The PME, as a subclass of the nonlinear anomalous diffusion equation, has been subjected to numerous and vast studies due to its possible applications for the porous media systems comprising of three essential equations: power-law equation of state, conservation of mass, and Darcy's law [36] . Many analytical [37] [38] [39] [40] [41] [42] [43] [44] [45] and numerical [46] [47] [48] [49] methods have been developed to study the properties of this model, which is not restricted to porous media systems, but also to the stock markets [24, 56] . The fractional PMEs (FPMEs) has been studied in many papers [51] [52] [53] 56] , aiming to study anomalous diffusion in porous media and other problems related to PME, each of which with a particular (local or non-local) "fractionalization" scheme. Finding * f.gharari@uma.ac.ir solutions for nonlinear anomalous diffusion equations is a challenge since, besides its difficulty to get exact analytical solutions, the principle of superposition is not applicable as in the linear case, so that the Fourier analysis can not be done. Despite this huge interest and theoretical studies on the problem, very limited information is available concerning the possible solutions of these equations and their properties, especially the dependence of the solutions on the fractionalization parameters.
In this paper, we aim to get to this issue by fractionalizing the PME with local and non-local fractional derivative operators. Intuitively a non-local operator is defined as the operator that needs the information in a finite interval upon its operation on a function, contrary to local operators that need only the information at one point in its close vicinity (see [75] , and Appendix A). We consider both local and non-local FPMEs focusing on three distinct cases: (LL) referring to the case where both time and space derivatives are local, (LN) or (NL) where one of them is local, and the other is non-local, and the (NN) referring to the case both derivatives of time and space are non-local. The main finding of the present paper is that the local and non-local cases admit generalized q-Gaussian functions as their Green function solutions. The difference between them is the number and the form of the fitting parameters. As an application, the local and non local generalized q-Gaussian distributions are used to describe the regimes observed during the time evolution of the probability density function (PDF) of the S&P 500 index. After addressing the q-Gaussian distribution function as a self-similar solution of the PME in the next section, we present the solution of the PME for the local fractional derivative (LL) in Sec. III. Sections IV and V, contain the analysis of the PME with (LN) and (NN) fractionalization. In Sec. VI, we present an application of the generalized q-Gaussian distribution to describe the price return of S&P 500 from the past 24 years, and we compare the results with previous solutions. [104] [105] [106] [107] [108] [109] Table I. Definitions of the most common fractional derivatives.
The PME is one of the simplest examples of a nonlinear diffusion widely used to describe processes that involve fluid flows, gas flows, and heat transfer [54] . The classical PME is:
A solution of this partial differential equation (PDE) is the Barenblatt function for q > 1 and t > 0 [83] ,
(2) where C is an integration constant. The Eq.(1) has been generalized to analyze several physical situations that present anomalous diffusion [57] . The present paper proposes a generalized form of the PME that admits a broader range of results:
The fractional derivative D of orders ξ and γ is a function of three variables: The limits (a and b), the argu-ments (t and x), and the degree order of the arguments (n and α). This last type of variable allows us to have a derivative with respect to a function when n, α = 1. A particular case of the Eq. (3), when n = α = 1 and a = b = 0, is the nonlinear fractional diffusion equation
No general solution of Eq.(4) is known. In the present paper, we aim to show that a particular solution to this PDE is the Green function. This function is obtained from the boundary condition P (x, t) = 0, when x → ±∞, and the initial condition P (x, 0) = δ(x), where δ(x) is the Dirac delta function. The Green functions can be expressed in terms of well-known distributions. Some cases are the Gaussian, the Levy-Stable [57] , and the q-Gaussian distributions. We show that the Eq.(4) admits exact solutions that vary depending on the definition of the fractional derivatives applied. The definitions of the most commonly used fractional derivatives are contained in Table I . The Eq.(4) allows space and time to scale differently, and as a consequence, different solutions can be obtained.
In searching for the solutions of Eq.(4), we exploit the fact that they follow the self-similarity law.
Green function Fractional derivatives Authors 1 Diffusion (*)
Bachelier [58] , A.Einstein [59, 60] 2 Anomalous superdiffusion
Riesz P.Levy [61, 63] , W.Feller [62] [63] [64] [65] 3 Classical PME (*)
Integer derivatives Barenblatt [66, 67] C.Tsallis [52, [68] [69] [70] 4 Space-FPME
C.Tsallis [32, 56] 5 Time-FPME (*)
Eq. (29) 6 Time-Space-FPME (*)
Eq. (12) 7 Time-Space-FPME
(time) and RL (space)
Eq. (21) 8 Particular case of the Generalized PME
ξ, γ > 0, ξ = 1 + λ + 1/α,
Eq. 
The Eq.(5) has often been used to model the price return in stock markets.
This price return obeys φ(t) ∼ t H , being H the characteristic exponent of the PDF. Additionally, for stock markets, f fits well to a q-Gaussian distribution [24] . Then, the equations presented in Table II can be used to model the detrended price return if they obey the power law and if its solution is a q-Gaussian.
A summary of these solutions are contained in Table  II . The first four equations of Table II have been solved previously by applying integer or fractional derivatives [55, 56, 68] . The last four equations have been solved in this manuscript after applying local and non-local derivatives. The local derivatives were used to solve equations N.5 and N.6 in Table II Table II . The equation N.8 has been solved by applying non-local fractional definitions on time and space. The solution of Eq. (3) is called second generalized q-Gaussian function and it presents a self-similar form, P (x, t) ∼ t −H F (xt −H ). By replacing λ = 1 1−q in the second generalized q-Gaussian, the first generalized q-Gaussian can be recovered. In the present paper, we aim to show that both local and non-local FPME admit generalized q-Gaussian solutions. For local derivation, the Katugamapola definition is applied; for non-local derivation, the Riemann-Liouville and Caputo derivatives are applied. For the non-local derivation, we consider a fractional derivative with respect to another function, in the sense of Caputo derivative.
The generalized forms of PME are obtained by replacing the first time derivative or second space derivative by fractional orders derivatives in the classical PME. These generalized PMEs may model more efficiently certain real-world phenomena, especially when the dynamics are affected by constraints inherent to the system. Typically, fractional derivatives are defined with an integral representation. Consequently, they are non-local in character. There exists several definitions for fractional derivatives and fractional integrals like the Riemann-Liouville, Caputo, Hadamard, Riesz, Grünwald-Letnikov. However, some usual properties of these fractional derivatives are different from ordinary derivatives, such as the Leibniz rule, the chain rule, and the semigroup property. Consequently, these fractional derivatives can not be applied for local scaling or differentiability properties. For further details, we refer the reader to [95, 97] and Appendix C.
The concept of local fractional derivatives keep some of the properties of ordinary derivatives. Nevertheless they lose the memory condition of fractional order derivatives [96] . There exist several definitions for local fractional derivatives like the Kolwankar, Chen, Conformable, Katugampola, see [75] for details. Recently, these local derivatives have been used to model phenomena of turbulence [84] , and anomalous diffusion [85] .
Local definitions for fractional derivatives were applied by Lenzi et al. [78] , his work contains some classes of solutions of a general fractional nonlinear diffusion equation with some observations. A similar study was made by Assis et al. [79] . More information about the use of local fractional derivatives can be found in [80] , and [81] . These examples are related to diffusion equations with nonlinear terms and fractional time derivatives which are quite few [76, 77] .
One popular type of local fractional derivatives is the "conformable operator" [86] . The "conformable operator" has been used in a wide range of applications. Some applications of the "conformable operator" are in Newtonian mechanics [87] , diffusion equation [88] , and nonlinear diffusion equation [89] . However, this local operator cannot be applied with zero as an order of the derivative. Recently, the Katugampola [91] operator has been used as a limit based fractional derivative that allows zero as a possible order of the derivative. The Katugampola operator maintains many of the familiar properties of standard derivatives such as the product, quotient, and chain rules. Throughout this section, we consider the Katugampola derivative (Katugampola operator), to solve the generalized PME. By applying this local derivative, our solution will be a generalized q-Gaussian distribution. Information about the Katugampola's definition and its properties can be found in Table I, Appendix A  and Table IV .
In this section, we solved a time-space FPME (TS-FPME) with (LL) fractionalization. The Katugampola fractional definition was applied for the time and space fractional operators. Such FPME can be written as:
where 0 < ξ ≤ 1 < γ ≤ 2, |ν| < 1 are a set of three free parameters, and D is the diffusion coefficient. To solve Eq. (6), we express the function P (x, t) in its self-similar form:
where φ(t) is a function to be identified. The Eq. (7) is consistent with a symmetric probability distribution. By considering, z = x φ(t) , and inserting Eq. (7) into Eq. (6), we have the following two equations:
In the above equation, the properties of the Katugampola derivative were used (see Appendix A and Table  IV) . Then, we arrange everything in such a way that all quantities in one side are only a function of z, and in the other side are a sole function of t. This procedure leads us to obtain the two following independent equations:
The solution of the first equation is φ ∝ t ξ ν+γ−1 , and for the second one we have:
, which can be rewritten as:
By applying the property: D µ (f ) = D µ−1 D 1 (f ) for 1 < µ 2, and taking local fractional integral with respect to z in Eq. (8), the following expression is obtained:
In the right hand side, the integration by parts is used,
ν−1 to obtain a special solution (where c 1 and c 2 are constants), the following expression is obtained:
After incorporate the previous expression into the Eq. (9), we obtain:
.
Therefore, the general solution is:
where c 1 is removed after apply the normalization condition in Eq. (10) . Then, by defining ν = 2 − q, considering α = 1−q+γ ξ , and
the following equation is reached:
where A L q is a normalization factor. In most of physical systems cases, the P (x, t) is symmetric with respect to x. This point leads us to make x → |x| (i.e. its absolute value), or we can consider some values of γ that satisfy this property. Then, the normalization factor was identified as follows:
where
For γ = 2, these parameters become α = 3−q ξ ,
.
We call Eq.(12) the local q-Gaussian (Lq-Gaussian) distribution, which is the Green function of Eq.(6) obtained from a TS-FPME with the Katugampola fractional derivative (local fractional definition). The Lq-Gaussian has been defined as, g γ q (x), Equation N.6 in Table II .
In Figure 1 , we show the Lq-Gaussians for different q values as indicated in the plot for t = 1, ξ = 1, γ = 2. In subfigure 1a, by increasing q, the peak of curves increase and the distribution becomes narrower (the tails become heavier). In the case of subfigure 1b something similar occurs, where we denote the PDFs of Lq-Gaussian for different ξ values as indicated in the plot for t = 1, q = 1.5 and γ = 2. The reverse occurs for Figure 1c , by increasing γ (considering t = 1, q = 1.5, ξ = 1), the peak of the PDFs decrease. Also, the time evolution of the Green function of Eq. 6 is shown in subfigure 1d.
In the recent subsection, we obtained Lq-Gaussian distributions by solving the TS-FPME. In fact, these Lq-Gaussians are a generalized q-Gaussians by considering |x| γ/2 , q > 1, i.e., a q−exponential in the variable |x| γ . From the definition of the q-exponential, it follows that f ∼ C f |x| −γ/(q−1) , C f > 0, as |x| → ∞. Analogously, for any q-Gaussian, g ∼ C g |x| −2/(q−1) , C g > 0, as |x| → ∞. By comparing the order of the power law of the asymptotes, we verify that for a fixed 1 < γ < 2, and for any 1 < q < 2 there exists a proportionality from Lq-Gaussians to q-Gaussian. For further details, see [99] .
Let us denote the class of random variables with
has a symmetric density f (x) with asymptotes f ∼ C|x| −(1+γ)/(1+γ(q−1)) , |x| → ∞, where 1 ≤ q < 2, 1 < γ < 2, and C is a positive constant. On the other hand, any Lq-Gaussian behaves asymptotically when C 1 /|x| γ/q−1 . Especially any Lq γ -Gaussian behaves asymptotically when C 2 /|x| γ/(qγ −1) . Hence, we obtain the following relationship :
Solving this equation with respect to q γ , we have
Three parameters were linked: γ, the parameter of the γ-stable Levy distributions, q, the parameters of correlations, and q γ , the parameters of attractors in terms of Lq γ -Gaussians. Then under Eq. (15) the density of X ∈ L q [γ] is asymptotically equivalent to Lq γ -Gaussian. The Lq-Gaussians have an interesting property. Its successive derivatives, and integrations with respect to |x| γ correspond to q γ,n −exponentials in the same variable |x| γ , where q γ,n = γq−n(q−1) γ−n(q−1) [99] . By considering G qγ,n [γ] as a set of functions {be −β|ξ| γ qγ,n , b > 0, β > 0}, and F qγ,n be the q γ,n −Fourier transform, the following expression is obtained:
This is similar to the q−exponential with the variable |ξ| γ , i.e. e −β|ξ| γ q , β > 0, which is the q−Fourier transform of (q, γ)−stable distributions [99] .
The drift is often an inevitable part of stochastic systems , that should be analyzed in details for every case study to control its effects. Although it is suggested to define the equations for the general drift term. For the case where it depends only on time (as the case for many physical systems of interest), the situation becomes easier. In this case, the governing equation is:
By change of variable τ = t ξ and the definition of the Katugampola derivative, we have:
where a (τ ) = 1 ξ a(t(τ )). By using the change of variable (s, y) ≡ (τ, x−x 0 −f (τ )), where f (τ ) = τ 0 a (τ )dτ , and using the fact that ∂y ∂τ = −a (τ ) and ∂ τ + a (τ )∂ x = ∂ s , one finds that the governing equation P (y, τ ) is:
for which the solution is (x 0 ≡ 0),
Let us equate the P (y, τ ):
P (x, t) = ∂y ∂x P (y, τ (t)).
Then, we obtain that:
where f (t) = f (τ (t)), A is a normalization factor, and c 1 is a constant. The Eq. (17) is a Lq-Gaussian solution with a drift.
To be self-contained, we consider the case where one fractional derivative is local, and the other is non-local. Therefore, in this section we solve a TS-FPME with (LN) fractionalization ,
where ∂ ξ ∂t ξ and ∂ γ ∂x γ denote the Katugampola and Riemann-Liouville fractional derivatives, respectively (see Table I ). We consider the decomposition of (Eq. 7). Then, by using the property of d γ dx γ F (ax) = a γ d γ dz γ F (z), and some properties of Katugampola derivative we obtain:
To continue, similar strategies applied in Sec. III B were used. We transform the previous equation into two independent equations:
where the solution of the first equation is φ = t ξ ν+γ−1 . For the second, the solution is : 
By integrating with respect to z, we obtain :
where c is a constant,we have set c = 0. By considering F (z) = z µ (c 1 + c 2 z) λ , and using the following property for the RL operators,
we find that
We put this result into Eq. (19) , and the following expressions are obtained:
The recent expressions reveal that the master equation admits the following solution:
where α = 1+γ 1−2γ , σ = 1−γ(1−γ) β(1+γ) , c 2 = 1 and
.
The solution is not a generalized q-Gaussian distribution. This occurs due to the factor, x αγ , which is located in front of the solution assuring that this class of PDFs tend to zero when x → 0. This makes a great difference with the previous function, the Lq-Gaussian Eq. (12). Figure 2 shows the LNLq-Gaussian distribution functions. All the graphs tend to zero as x → 0 for all times, as stated above, and don't have symmetric form around their peaks in terms of x. They manifested an abrupt increase before the peak, and long-range tails beyond it, with peaks depending on q and ξ. The peak values increase by increasing both γ and ξ.
In this section we solve one particular case of the G-PME displayed in Eq.(3), when n = −1, a = 0, and b = 1. The solution is obtained by considering a hybrid case, where the time derivative is RL operator, and the space derivative is Caputo operator (Appendix B). The other cases (RL-RL, Caputo-Caputo, and Caputo-RL derivatives) are straightforward to be processed following the same lines as this study. The resulting FPME equations is
where D is the RL operator for time derivative and C D is the Caputo operator acting on "space" coordinate x.
To construct our solution, we need to restrict ourselves to the case ν = 1+ξ 1−ξ leaving two parameters free for fitting, γ and α.
We again search for the solutions of the form of Eq. (7), where the parameters were defined in the Section III . In the following we show that the above equation admits the solution of the form F (z) = (1 + bz α ) λ , where b = 0, and α = 1. When stablished, this solution serves as another variant of the generalized q-Gaussian solution. By inserting this form in the right side of Eq. (22) we obtain:
To obtain the above equation, we have used the following property of the fractional derivatives, that is valid for all the fractional differential operators considered here,
Additionally, the following property of the fractional derivative of a function with respect to another function [93] was applied :
where n = α, if α ∈ N and n = [α] + 1, if α / ∈ N. By inserting Eq. (7) in the left side of Eq. (22) and then applying Eq. (20) for the RL operators with δ = α+β +1, we get
where b = bx α . By equating Eq. (23) and Eq. (25), one find that they match each other, yielding to:
Therefore, we see that the solution is:
where the pre-factor A N L q is a normalization constant, and η N L q (α, λ, D) is a constant depending on A N L q . If we again suppose that the distribution is symmetric with respect to x. Therefore, x → |x| (i.e. its absolute value), then the normalization constant is:
where η(α, λ, D) = η N L q (α, λ, D) λ . From this expression one can calculate the final expression of η N L q (α, λ, D) Table II . The plots for NLq-Gaussian solutions are shown in Figure 3 , for various α and q values. With respect to the local case, the behavior for the non-local case is more complicated. As is seen from the subfigure 3a, the case where α is kept constant and q increases is evaluated. For q < 1.3, the peak rises, and for q ≥ 1.3, it decreases. In subfigure 3b, for a constant q, however, the peak increases when α increases. The subfigure 3c shows the time evolution of the Green function of Eq. 22, which solution is the NLq-Gaussian distribution. The example was made for the values of q = 1.5 and α = 2. We have shown that the distribution widens as time goes on.
The price return in stock markets exhibits remarkable characteristic features. The most largely observed feature in recent studies is the self similarity law , where the PDF obeys:
in which f is a normalized distribution that is usually fit to a q-Gaussian [82] . In earlier work, the function f was assumed as a Levy-stable distribution function, L γ , which has the drawback that it presents infinite standard deviation and it does not obey the empirical power law tails [100, 101] . For long time returns it will be proved that f is a Gaussian distribution function, where the price return behaves like independent and identically distributed random variables but still following the self similar principle. Particular cases of the generalized PME are presented in Table III . The solutions of each of these partial differential equations obey the self-similar law presented in Eq. (28) and are related to the q-Gaussian distribution function.
In this part, we analyze the S&P500 stock market data during the 24-year period from January 1996 to August 2020 with a frequency of one minute. The detrended price return is defined as,
where I * (t o ) is the detrended stock market index at time t o , and I * (t) is the detrended stock market index for any time t > t o . The PDF of the detrended price return has Classical (C-PME)
Time-Space fractional (TS-FPME) Particular case of the Generalized PME (G-PME ) Forms of PME
Parameters after comparing with General PME q = 2 − ν, γ = 2, Table III . Summary of particular forms of the generalized PME Eq. 3 to model the time evolution of the PDFs of price return. The C-PME, T-FPME, TS-FPME and G-PME are obtained after set a specific value of the parameters in the general form of the PME. The fittings were perform by setting the parameters as shown in the table. Figure 4 . Three different zones were determined based on an abrupt slope change of the fitting parameters α and q. The contour plot represents the PDF of the detrended price return. The black circles represent the end points of the strong superdiffusion regime (zone A) from t = 0 to t = 35 minutes. These points are the ends of the bump obtained from the two points at the PDF with an abrupt change of slope (Figure 1-c in [24] ). The remaining area during the first 35 minutes corresponds to a weak super-diffusion regime (Zone B1). From 80 min to 10 days, the zone corresponds to a weak super-diffusion regime (zone B2). A normal diffusion process is reached after 30 days. The gray dashed lines represent the transitions between each zone.
been fitted using the q-Gaussian distribution [82] , where different zones had been captured from strong to super diffusion regime previously [82, 94] . Figure 4 shows the evolution of the PDF and its regimes in the (x, t) space. Initially, the PDF has a pronounced bump in the center that fully disappears close to 80 min. During the first 35 min there is a power law relation of time against the end points of the bumps at the PDFs (See Subfigure 1-e in [24] ). This zone is defined as the strong super diffusion regime (Zone A). The remaining area during that time and the following next area close to 10 days corresponds to the weak super diffusion regime (zone B). Finally, the last regime corresponds to a normal diffusion process (Zone C) and is reached after 30 days approximately.
We have reconstructed the time evolution of the PDFs of the detrended price return and we collapsed them after applying the corresponding re-scaling factor. Four equations of the Table II have been used to model this behaviour. The Classical PME (C-PME), the time FPME (T-FPME), the time-space FPME (TS-FPME) and a particular case of the generalized PME (G-PME). These solutions obey Eq. (28) , and are presented in Table III with more detail.
The TS-FPME and the particular solution of the G-PME proposed in this manuscript are new options to model the time evolution of the detrended price return. The Lq-Gaussian and NLq-Gaussian, which are the solutions of the TS-FPME and the particular case of the G-PME respectively, fit the collapse of the PDFs of the detrended price return well. Figure 5 shows the result of these fittings.
The first best option to fit the price return was ob- Table I : no-fractional(green), time-fractional(red), time-space fractional (black), and a particular form of the G-PME (purple). For the subfigure (d) the collapsed data is fitted by a Gaussian distribution function, which is a concurrent solution for the four previous PDEs when q = 1 and γ = 1.
tained by replacing f as the Lq-Gaussian (g γ q ) in Eq. (28) , this equation can be written as:
where C γ q is the normalization constant. B is related with the diffusion term, both of them are detailed in Table II .
The second best option was obtained by replacing f as the NLq-Gaussian (g α,λ q ) in Eq. (28) , this second equation can be written as:
The C λ,α q is the normalization constant, and B is related with the diffusion term. The definitions of these parameters are expressed in Table II Table III . The best fitting for the four cases is the NLq-Gaussian. However, the four solutions constitutes an acceptable solution for the correspondent collapses of the PDFs. A converge to a Gaussian normal distribution is observed for long time returns.
We provided different solutions for the generalized form of the FPME. To this end, we had considered the generalized PME (G-PME) as the master equation. We introduced the anomalous PME as a nonlinear fractional diffusion equation, which is a particular case of G-PME. The solutions were built by considering the local and non-local fractional derivatives assuming a Dirac's delta function as the initial condition. Our analysis proved that the solution are given by a generalized q-Gaussian, which obey a self similar law. The fractional derivatives were classified as local and non-local, where the Katugampola's is the local fractional derivative and the Riemman-Lioville, Caputo and Riesz are non-local fractional derivatives. First, we considered the case where the derivatives are local. The resulting solution is Lq-Gaussian, which is the first generalized q-Gaussian function. This solution fits the PDF of the detrended price return well (Sec. VI). The second analyzed class of G-PME is the one in which the time and space derivatives are given by the non-local fractional generalizations: Riemann-Lioville and Caputo, both of them based on the Laplace transform. For this second class, the fractional derivatives are evaluated with respect to another function, and proved that they admit the second generalized q-Gaussian solution. This second solution is symmetric about its mean (peak). The NLq-Gaussians hold a different self-similar law than the Lq-Gaussians. The main difference is that the self-similarity of the Lq-Gaussians depends on α and γ only, for the NLq-Gaussians the self-similarity depends on α, γ and one extra global exponent λ, where the first generalized q-Gaussian is recover for λ = 1 1−q . A hybrid equation has also been considered, where the time dimension is local fractional derivative, and the spatial dimension is non-local. The solution that we found is again proportional to the generalized q-Gaussian (which we named LNLq-Gaussian), but they obey a power in x that causes the PDF to vanish in the limit x → 0. The shape of these PDFs are very different from the q-Gaussian distribution and they are not symmetrical distribution.
The Lq-Gaussian (first generalized q-Gaussian) and NLq-Gaussian (second generalized q-Gaussian) have been used to model the detrended price return of S&P500. Both distribution functions describe well the fitting of the detrended price return. The NLq-Gaussian is the best model to fit the probability of the detrended price return. For future work the generalized form of the PME will be solved by applying the q-Fourier analysis. The ordinary Fourier analysis only applies for linear operators. The generalized PME contains nonlinear operators, preventing us from using the ordinary Fourier analysis.
Appendix A: Properties of Katugampola derivative In here, we give a brief summary of the definition of the Katugampola fractional operator and some of its properties. This local fractional operator is used to construct the TS-FPME in sec. III. If 0 ≤ α < 1, the Katugampola operator generalizes the classical calculus properties of polynomials [110] . Furthermore, if α = 1, the definition is equivalent to the classical definition of the first order derivative of the function f . The Katugampola derivative is defined as:
for t > 0 and α ∈ (0, 1]. When α ∈ (n, n + 1] (for some n ∈ N, and f is an n−differentiable at t > 0,), the above definition generalizes to
and if f is (n + 1)−differentiable at t > 0, then
In continue, we review some properties of the Katugampola derivative in Table IV . If f is α− differentiable in some (0, a), a > 0, and f (α) (0) = lim t→0 + D α f (t) exists, the following properties hold for Katugampola derivative. For f, g, be α− differentiable at a point t > 0.
One can define the inverse of the D α operator as a fractional integral,
where the (.) symbol is serving as place holder for the function to be operated upon. One verifies that
where f vanishes at the lower limit. Then: This section contains definitions of non-local fractional operators that are used in this paper to construct the TS-FPMEs. The Riemann-Liouville fractional derivative is a fractional operator that is used in Sections IV and V as a non-local fractional operator to construct the TS-FPMEs with (LN) and (NN) fractionalizations. The integral representation for this operator is:
where n − 1 < α ≤ n [95] . A recent variation of RL operator is the Caputo derivative [98] , defined as:
where C stands for Caputo and f (n) is the nth derivative of f . The main advantage of the Caputo derivative is that the derivate of a constant is zero, which is not the case of the RL operator. Substantially, this kind of fractional derivative is a formal generalization of the integer derivative under Laplace transform [56] .
A generalized fractional operator that we used to construct the TS-FPME with (NN) fractionalization in Sec. V is the Caputo fractional derivative of a function with respect to another function [93] , and defined as:
...
Note that the recent integral representation in the special case ψ(x) = x is reduced to the integral representation of the Caputo derivative. We solve a particular case of the generalized PME, Eq. (3), described by a fractional derivative of a function with respect to another function. This innovative approach will be useful to solve other physical problems that present a self-similar pattern and can be modelled by a q-Gaussian.
In this section, we give a short review of the properties of the fractional derivatives: Katugampola, Riemann-Lioville, and Caputo. The Katugampola is one definition for the local fractional derivative. The Riemann-Lioville and Caputo are definitions of non-local derivatives. A comparison between each property of these fractional derivatives are presented in Table IV . Katugampola's corresponds to the ordinary derivative when α = 0 and α = 1. The Riemann-Lioville and Caputo are an analytical continuation of the ordinary derivatives. The main difference between them is that the Caputo derivative of a constant is zero, a property that does not hold for Riemann-Lioville derivative. This desirable property is satisfied by the Katugampola local derivative, too.
[1] Escobedo, M., Herrero, M. and Velazquez, J., 1998.
A 
Linearity D α (af (t) + g(t)) = aD α f (t) + D α g(t)
D α (f (t)g(t)) = f (t)D α g(t) + g(t)D α f (t)
RL D α (f (t)g(t)) = ∞ k=0 α k ( RL D α−k f (t))g k (t) C D α (f (t)g(t)) = RL D α (f (t)g(t)).... 