LRD spectral analysis of multifractional functional time series on manifolds

This paper addresses the estimation of the second–order structure of a manifold cross–time random field (RF) displaying spatially varying Long Range Dependence (LRD), adopting the functional time series framework introduced in [28]. Conditions for the asymptotic unbiasedness of the integrated periodogram operator in the Hilbert–Schmidt operator norm are derived beyond structural assumptions. Weak–consistent estimation of the long–memory operator is achieved under a semiparametric functional spectral framework in the Gaussian context. The case where the projected manifold process can display Short Range Dependence (SRD) and LRD at different manifold scales is also analyzed. The performance of both estimation procedures is illustrated in the simulation study, in the context of multifractionally integrated spherical functional autoregressive–moving average (SPHARMA(p,q)) processes.


Introduction
The literature on weakly dependent functional time series has been widely developed in the last few decades, allowing the statistical analysis, and inference on stochastic processes under a Markovian framework (see, e.g., [7]; [16]).Nowadays, spectral analysis of functional time series models constitutes an open research area.Under suitable functional cumulant mixing conditions, and the summability in time of the trace norm of the elements of the covariance operator family, in [21], a weighted periodogram operator estimator of the spectral density operator is derived.Its asymptotic analysis is addressed.Particularly, the asymptotic normality of the functional discrete Fourier transform (fDFT) of the curve data is proved (see also [29]).In [22], a harmonic principal component analysis of functional time series, based on Karhunen-Loéve-like decomposition in the temporal functional spectral domain is proposed, the so-called Cramér-Karhunen-Loéve representation (see also [26]; [27]).Some recent applications in the context of functional regression are obtained in [23].Hypothesis testing for detecting modelling differences in functional time series dynamics is achieved in [30] in the functional spectral domain.
Recently, an attempt to extend spectral analysis of functional time series to the context of LRD functional sequences has been presented in [28], covering, in particular, some examples of the LRD funtional time series family analyzed in [17] in the temporal domain.The application of harmonic analysis in this more general context entails important advantages as given in [28].Particularly, under stationary in time, the temporal dependence range can be approximated from the behavior in a neighborhood of zero frequency of the spectral density operator family at different spatial resolution levels.Moreover, a more flexible modeling framework can be introduced in this setting.Particularly, the projected process can display LRD and SRD depending on the spatial scale, according to the support of the spectral measure of the LRD operator characterizing the distribution of its eigenvalues.In [17], Functional Principal Component Analysis (FPCA), based on the long-run covariance function, is applied in the consistent estimation of the dimension and the orthonormal functions spanning the dominant subspace, where the projected curve process displays the largest dependence range.Fractionally integrated functional autoregressive moving averages processes constitute an interesting example (see [17]).The multifractionally integrated version of this process family can be analyzed under the modeling framework introduced in [28].
Connected and compact two-point homogeneous spaces constitute an example of manifold, with isometrically equivalent properties to the sphere, locally resembles an Euclidean space.Here, we will denote it as M d , with d being its topological dimension.The isotropy or invariance of a kernel with respect to the group of isometries of M d allows its diagonal representation in terms of a fixed orthonormal basis given by the eigenfunctions of the Laplace-Beltrami operator on L 2 (M d , dν).Thus, the separable Hilbert space H = L 2 (M d , dν) of square integrable functions on M d is considered in our functional time series analysis.Here, dν is the normalized Riemannian measure on M d .Particularly, this Hilbert space framework has been adopted by several authors in the current literature for the special case of the sphere.That is the case of [10], where estimation and asymptotic analysis of spherical functional time series is achieved, introducing new model families (see [8]).Also, in the LRD framework, sphere cross-time random fields are analyzed in [20], investigating the asymptotic behavior, under temporal increasing domain, of the empirical measure of a excursion area at any threshold.Time-dependent RF solution to a fractional pseudodifferential equation on the sphere is introduced in [14] (see also [3]).The eigenfunctions of the Laplace Beltrami operator on L 2 (M d , dν) , and the corresponding zonal functions play a crucial role in the analysis of manifold cross-time random fields (see, e.g., [18]).For example, Sobolev regularity and Hölder continuity of Gaussian RFs on a connected and compact two-point homogeneous space are studied in [12] and [13], by exploiting asymptotic properties of the pure point spectra of invariant kernels, and projection of functions into the eigenfunctions of the Laplace Beltrami operator.Some motivating real data applications can be found in the field of Cosmic Microwave Background (CMB) radiation (see, e.g., [19] and references therein).It is well-known the interest of these RFs in climatic change analysis (see, e.g., [1]).
The spectral domain allows to characterize LRD in functional time series in terms of the unboundedness at zero frequency of the corresponding element of the family of spectral density operators.Specifically, the divergence of the eigenvalues of the elements of the spectral density operator family at a neighborhood of zero frequency leads to different levels of singularity at this frequency depending on the spatial scale (see, e.g., [28]).The case of SRD and LRD at different spatial scales can also be analyzed, in the case of non-trivial null space of the LRD operator.Thus, the projected process displays SRD in this subspace, and LRD in the eigenspaces associated with the elements of the support of the spectral measure of the LRD operator.
From a theoretical point of view, this paper contributes providing a sufficient condition for the asymptotic unbiasedness of the integrated periodogram operator, in the Hilbert-Schmidt operator norm, for the class of zero-mean, stationary and isotropic, mean-square continuous Gaussian, or elliptically contoured, spatiotemporal RFs on M d .This result provides a suitable setting for applying the LRD spectral functional time series framework introduced in [28], where the weak-consistent estimation of the second-order structure in the functional spectral domain is achieved under a Gaussian scenario.Our scenario is a bit different, since the Hilbert-Schmidt operator scenario has been considered.Note that this scenario has usually been adopted in the current literature on functional time series (see [7]; [10]).The case where SRD and LRD are displayed at different manifold scales is also addressed, combining semiparametric estimation of the spectral density operator based on minimum contrast at manifold scales where LRD is displayed, and a nonparametric estimation based on the weighted periodogram operator in the remaining manifold scales where SRD is observed.The preliminary result derived in the Supplementary Material on spatiotemporal Karhunen-Loéve expansion of the restriction to a bounded closed temporal interval of a zero-mean, stationary in time and isotropic in space, mean-square continuous Gaussian spatiotemporal RF on M d has been applied.
From a practical point of view, the simulation study undertaken in Section 5 illustrates the performance beyond the Gaussian scenario of the two proposed estimation approaches in the functional spectral domain, in the context of multifractionally integrated SPHARMA(p,q) processes.Two cases are analyzed respectively corresponding to a decreasing and increasing positive bounded eigenvalue sequence of the LRD operator.Particularly, we study the scenario where the projected process displays SRD and LRD at different spherical scales.This last case is illustrated when the eigenvalues of the LRD operator vanish at high discrete Legendre frequencies.Details on the implementation of both estimation methodologies, and some conclusions on the results obtained beyond the Gaussian scenario are also included in Sections 5 -6 (see also the Supplementary Material).
The outline of the paper is the following.Section 2 introduces notation, technical tools and preliminary elements.Asymptotic unbiasedness of the integrated periodogram operator, and weak consistent minimum contrast estimation are studied in Section 3.An extended formulation of this estimation methodology to the case where temporal LRD and SRD are displayed at different manifold scales is proposed in Section 4. A simulation study is undertaken in Section 5 to illustrate the performance of the proposed estimation methodologies.A summary of conclusions about the simulation study is given in Section 6.The Supplementary Material complements the theoretical material about the paper, and provides the results of the simulation study in the remaining scenarios analyzed that are not displayed in the paper.

Preliminaries
The family of manifold cross-time RFs analyzed in this paper is introduced in this section.The connected and compact two-point homogeneous spaces are briefly described, providing some preliminary algebraic notions, with reference to the invariant probabilistic measure, and the Laplace Beltrami operator.
Let X = {X(x, t), x ∈ M d , t ∈ T} be a wide sense stationary in time and isotropic in space zero-mean, and mean-square continuous Gaussian, or elliptically contoured, RF on the basic probability space (Ω, A, P ).Here, T denotes the temporal domain, which usually is Z or R (see also the Supplementary Material for the case of a bounded temporal interval [0, T ]).Assume also that the map X t : (Ω, A) −→ L 2 (M d , dν), B(L 2 (M d , dν)) is measurable, with X t (x) := X(x, t) for every t ∈ T and x ∈ M d .Here, B(L 2 (M d , dν)) denotes the σ-algebra generated by all cylindrical subsets of L 2 (M d , dν).
Let d M d be the geodesic distance induced by the isometry with the unit sphere, preserving the spherical distance ρ(x, y) = cos −1 (⟨x, y⟩), x, y ∈ S d , and let ω d = M d dν(x).In the following, R , with P (α,β) n being the Jacobi polynomial of degree n with a pair of parameters (α, β) (see, e.g., [2]).For each n ∈ N 0 , {S d n,1 , . . ., S d n,δ(n,d) } is the orthonormal basis of eigenfunctions of the eigenspace H n of Laplace-Beltrami operator ∆ d on L 2 (M d , dν), associated with the eigenvalue λ n = −nε(nε ) .
Note that α = (p + q − 1)/2, β = (q − 1)/2, and ε = 2 if M d = P d (R), the projective space over the field R, and ε = 1, otherwise.Parameters q and p are the dimensions of some root spaces connected with the Lie algebras of the groups G and K, with M d ≃ G/K, being G the connected component of the group of isometries of M d , and K the stationary subgroup of a fixed point o in M d (see, e.g., Table 1 in [18]).The next lemma is applied along the paper.
Lemma 1 (See [15, Theorem 3.2.]and [2, p 455]) For every n ∈ N 0 , the following addition formula holds: The following diagonal series expansion holds under the conditions of Theorem 4 in [18]: 3 Operator based minimum contrast parameter estimation of LRD This section adopts the semiparametric spectral LRD functional time series framework introduced in [28] (see Condition C1 below), for inference on an LRD manifold cross-time RF, constructed from a zero-mean, stationary in time, and isotropic in space, mean-square continuous Gaussian, or elliptically contoured spatiotemporal RF X = {X(x, t), x ∈ M d , t ∈ T = Z} (see, e.g., [20] for sphere cross-time RFs).Note that under the conditions assumed in Section 2, X = {X(x, t), x ∈ M d , t ∈ T = Z} defines a functional time series { X t (•), t ∈ Z}.In Section 3.1 below, we will work under this scenario to implement minimum contrast parameter estimation under Condition C1 (providing the semiparametric functional spectral framework introduced in [28]).Condition C0 below establishes a sufficient condition for the asymptotic unbiasedness of the integrated periodogram operator of X in the Hilbert -Schmidt operator norm beyond structural assumptions (see Theorem 1).The weak-consistency of the minimum contrast estimator of the LRD operator is then obtained in Theorem 2, under new Condition C0 and Condition C1, adopting the L 2 (M d , dν)-valued time series framework.
The following condition will be assumed in the subsequent development.

Condition C0. The elements of the temporal coefficient sequence {B
where Note that Parseval identity leads to The convergence to zero in the Hilbert-Schmidt operator norm of the integrated bias of the periodogram operator also holds under (3.1).
where F (T ) ω denotes the mean of the periodogram operator being the functional Discrete Fourier Transform (fDFT), based on a functional sample X t , t = 1, . . ., T, of spatiotemporal RF X.Here, as before, for each t = 1, . . ., T, X t (x) := X(x, t), for every x ∈ M d .
Proof The proof follows from equation (3.3) implying Lemma 1 in [28] holds in our context.Hence, equations (4.6)-(4.8) in the proof of Theorem 1 in [28], and the remaining steps of the proof of this theorem can be obtained in a similar way.

Minimum contrast parameter estimation
For the implementation of the minimum contrast parameter estimation of the spatial-varying LRD parameter of X, in the L 2 (M d , dν)-valued time series framework introduced in [28], the following condition is assumed: , for any n ≥ 0, and θ ∈ Θ, for certain l α (θ), L α (θ) ∈ (0, 1/2).The elements of the function sequence {M n , n ∈ N 0 } are strictly positive continuous functions on [−π, π], slowly varying at zero frequency in the Zygmund's sense (see, e.g., Definition 6.6 in [6]).For ev- } defines the SRD spectral family.The associated kernel family satisfies The parameter space Θ is assumed to be compact with non null interior.For each θ ∈ Θ, {α(n, θ), n ∈ N 0 } is the system of eigenvalues of the parameterized long-memory operator A θ with kernel K A θ admitting the following diagonal series expansion: Hence, for every θ ∈ Θ, A θ has isotropic kernel and defines a strictly positive self-adjoint operator on L 2 (M d , dν), with norm in the space L(L 2 (M d , dν)) of bounded linear operators less than 1/2.Note that, since sin(ω) ∼ ω, ω → 0, where the frequency varying operator |1 − exp (−iω)| −A θ /2 is interpreted as in [11]; [24].In particular, equation (3.1) in Assumption II in [28], characterizing LRD of functional time series in the spectral domain, holds.
Remark 1 Condition C1 restricts the eigenvalues of A θ to the interval (0, 1/2), leading to a shorter range of spectral singularity at zero frequency at any spatial scale, allowing (3.3) holds.In particular, The sequence {B η n (0) ≥ 0, n ∈ N 0 } in equation (3.5) defines the eigenvalues of the trace integral autocovariance operator where Bn(0) dω, for each n ∈ N 0 .Assume that the true parameter value θ 0 lies in the non empty interior of compact set Θ, and α(•, θ 1 ) ̸ = α(•, θ 2 ), for θ 1 ̸ = θ 2 , and θ 1 , θ 2 ∈ Θ, ensuring identifiability.Let α T (n, θ) = α(n, θ T ), n ∈ N 0 be the parametric estimators of the eigenvalues of A θ by minimum contrast.The computation of the minimum contrast parametric estimator θ T requires the introduction of some operator families in the spectral domain as now briefly describe.
Operator integrals are understood here as improper operator Stieltjes integrals which strongly converge (see, e.g., Section 8.2.1 in [25]).Specifically, in the definition of our loss function, integration in the temporal spectral domain with respect to weighting operator W ω is achieved.For each ω ∈ [−π, π], the invariant kernel K Wω of W ω on M d ×M d admits the following series expansion: where the eigenvalues Under Condition C1, define the normalizing self-adjoint integral operator N θ by the kernel Hence, for each θ ∈ Θ, and ω ∈ [−π, π], ω ̸ = 0, the kernel K Υ ω,θ of the density operator Here, δ(x−y) denotes the Dirac Delta distribution.Equivalently, Given the candidate set constituted by the parametric operator families {Υ ω,θ , ω ∈ [−π, π]} , θ ∈ Θ, the theoretical loss function L(θ 0 , θ) to be minimized is defined as where the theoretical contrast operator U θ has kernel Hence, U θ ∈ L L 2 (M d , dν; C) (see also Remark 7 in [28]).
Proof The proof follows as in Theorem 2 in [28].
Global analysis.Condition (3.1) is a key condition in our approach, leading to equation (3.3) introducing the Hilbert-Schmidt operator setting, usually considered in functional time series analysis.Note that, under this setting, LRD still can be displayed (see equation (3.7)).Furthermore, equation (3.3) allows to apply the spectral analysis of LRD functional time series introduced in [28] for inference on manifold cross-time RFs.Specifically, equation (3.1) ensures Lemma 1 in [28] holds under alternative conditions for this RF family.
Then, asymptotic unbiasedness of the integrated periodogram also holds, beyond the Gaussian scenario under non structural assumptions (see Theorem 1).Furthermore, Section 3.1, and, in particular, Condition C1, provides the semiparametric functional spectral scenario introduced in [28] to be applied for minimum contrast estimation, when a functional sample of a manifold cross-time RF can be observed.Theorem 2 ensures weak-consistency under a Gaussian scenario as given in [28].

SRD-LRD estimation in the spectral domain
This section has a double, theoretical and practical, motivation.Specifically, on the one hand a wider family of spatiotemporal RF models is analyzed displaying SRD and LRD at different manifold scales, in the spirit of the LRD framework introduced in [17], but extended to the multifractional context.On the other hand, this manifold-scale varying memory behavior can be observed in some stochastic fractional or multifractional pseudodifferential in time evolution equations, defined from a local spatial differential operator, where the decay velocity of the temporal correlation function is accelerated by the decay of the spatial pure point spectrum, modifying the temporal LRD level of the model at the smallest spatial scales.The reverse situation can be observed in processes defined by evolution equations given by a local differential operator in time, and a fractional or multifractional pseudodifferential operator in space.Thus, the solution to these models displays a spatial fractal behavior, reflected in a slow decay of its spatial pure point spectrum, slowing down the decay velocity of the temporal correlation function, leading to a stronger dependence at small spatial scales.Both behaviors can be observed within the model family introduced in [3]; [4]; [5].
Let us now consider in Condition C1, l α (θ) = 0, and α(n, θ) = 0, for n ∈ D SRD ⊂ N 0 , meaning that the process projected into the eigenspaces H n , n ∈ D SRD ⊂ N 0 , of the spherical Laplace Beltrami operator displays SRD.While LRD is observed at the remaining eigenspaces.Without loss of generality let D SRD = {0, . . .n 0 }, for certain n 0 ≥ 1. (The reverse situation where the eigenvalues of the LRD operator vanish at large n has been analyzed in the simulation study in Section 5.2).The projected SRD process then admits the following expansion (see Theorem 1 in the Supplementary Material): For ω ∈ [−π, π], the fDFT projected into ⊕ n0 n=0 H n is expressed as The corresponding projected periodogram operator p then involves the tensorial product of the eigenfunctions of the Laplace Beltrami operator up to order n 0 .For ω ∈ [−π, π], the kernel estimator , with B T being the positive bandwidth parameter.Function W on R is real, and such that W is positive, even, and bounded in variation, with R W (x)dx = 1 (see [21]).Under Condition C1, the minimum contrast estimators α T (n, θ 0 ) = α(n, θ T ), n > n 0 , are computed as given in equations (3.8)-(3.16).Hence, the mixed SRD-LRD kernel estimator of the spectral density operator is given by

Simulation study
This section illustrates the results derived in the context of multifractionally integrated SPHARMA(p,q) processes (see, e.g., Example 1 in Section 3.3 in [28], and [17] in the particular case of fractionally integrated functional autoregressive moving averages processes).See also [9] and [10] for the case of SPHARMA(p,q) processes.
Considering the unit sphere S 2 in R 3 , and hence, the separable Hilbert space H = L 2 (S 2 , dν), a multifractionally integrated SPHARMA(p,q) process X t , t ∈ Z is defined by the following state space equation: (5.1) Condition C1 holds, since this model constitutes a particular case H = L 2 (S 2 , dν) of the more general formulation given in Section 3.3 in [28] in the functional time series framework.Then, as before, A θ is the LRD operator satisfying l α (θ), L α (θ) ∈ (0, 1/2).In particular, equation (3.7) holds.Here, operator (I L 2 (S2,dν) −B) A θ /2 is interpreted as in [11]; [24], and B is a difference operator satisfying E∥B j X t − X t−j ∥ 2 H = 0, for t, j ∈ Z.
Theorem 1 and Section 1.1 of the Supplementary Material have been applied in the simulation methodology adopted in the generations of some special cases within the family of multifractionally integrated SPHARMA(p,q) processes.In particular, a spherical uniform pole U = u 0 (see Figure 2) in the involved zonal functions (defined from the Legendre polynomials {P n , n ∈ N 0 }) has been independently generated of the L 2 (S 2 , dν)-valued Gaussian strongwhite noise innovation process with variance 3 of the Supplementary Material provides the specific parametric scenarios considered.Under such scenarios, Section 3 of the Supplementary Material displays generations of some special cases of multifractionally integrated SPHARMA(p,q) processes, for p = 1, 3 and q = 0, and p = 1, 3, and q = 1 in equation (5.1).The particular cases analyzed of LRD operator eigenvalue sequences can be found in Table 1 of the Supplementary Material.In the implementation of the minimum contrast estimation methodology, we consider 100 candidate systems of parametric eigenvalues for the LRD operator (see Table 2 in the Supplementary Material).These candidate sets are displayed in Figure 1 jointly with the true eigenvalue sequence.
The functional spherical values at different times of the generated multifractionally integrated SPHAR(3) process are displayed for M = 10 under decreasing eigenvalue sequence of the LRD operator in Figure 3.As commented, the remaining special cases are plotted in Section 3 of the Supplementary Material.Additionally, generations under LRD operator with eigenvalues vanishing for n ≥ n 0 = 16, i.e., the LRD-SRD case, is showed in Section 3.1 of the Supplementary Material for multifractionally integrated SPHARMA(1,1) process.
It can be observed in all generations, displayed on a spherical angular neighborhood of the selected pole (see Figure 2), that persistent in time of the local spherical behavior is stronger when a positive bounded non-decreasing sequence of eigenvalues of LRD operator is considered, since the highest positive eigenvalues are displayed at high discrete Legendre frequencies.For decreasing eigenvalue sequence of the LRD operator, the opposite LRD effect is observed through spherical scales.One can also observe a symmetry in the evolution in time of the spherical patterns, under both, decreasing and increasing LRD operator eigenvalue sequences, when spatial spherical SRD is observed.An increasing level of local linear correlation in space at intermediate times is displayed when autoregression of order 3 is considered.

Minimum contrast estimation results
We now display the minimum contrast estimation results referred to the multifractionally integrated SPHAR(3) process, when LRD operator A θ has decreasing sequence of eigenvalues as plotted at the left-hand side of Figure 1, under a truncation order M = 30.See Section 4 of the Supplementary Material, for the remaining cases of multifractionally integrated SPHAR(1), SPHAR(3), SPHARMA (1,1), SPHARMA(3,1) processes.Specifically, Figure 4 displays, for i = 1, . . ., 100 frequency nodes, operator |M ωi | 1/2 (left-hand side), and 4(sin(ω i /2)) 2 −A θ /4 |M ωi | 1/2 (right-hand side), projected into H n , n = 1, . . ., 30.These factors have been computed in the implementation of the minimum contrast estimation methodology in the functional spectral domain.Figure 5 provides the projected spectral density operator kernel at temporal Fourier frequencies ω = −π + 0.0628 (10)i, i = 1, 3, 7, 10 in the interval [−π, π].Its empirical counterpart is given in Figure 6 where the modulus of the projected fDFT, and tapered periodogram operator kernel of the generated data at temporal Fourier frequency zero are displayed.The last one over a grid of 30 × 30 Legendre frequencies.The projected empirical contrast operator U T,θ in equation (3.15) is then computed.Its minimization is performed in the bounded operator norm.  1 of the Supplementary Material), Figure 7 displays the histograms of the temporal mean of the empirical absolute errors from R = 100, 2000, 5000 independent generations of a functional sample of size T = 1000 of multifractionally integrated SPHAR(3) process.The empirical analysis performed for the remaining cases under decreasing LRD operator eigenvalue sequence, from R = 100, 2000, 5000 independent generations of each one of the functional samples of size T = 50, 500, 1000 considered, are given in Sections 4.1, 4.3, 4.5, and 4.7 of the Supplementary Material (see also Sections 4.2, 4.4, 4.6 and 4.8 of the Supplementary Material for the same analysis under increasing LRD operator eigenvalue sequence as plotted at the right-hand side of Figure 1).The results are displayed for the eigenspaces H n , n = 1, 5, 10, 15, 20, 25, 30.The empirical probabilities are computed under decreasing and increasing LRD operator eigenvalue sequence, for the multifractionally integrated SPHAR(3) model in Figure 8, considering R = 100, 2000, 5000 independent generations of each one of the functional samples of size T = 50, 500, 1000.In all scenarios displayed in Section 4 of the Supplemntary Material, the empirical probabilities (5.3) are computed from projection into the eigenspaces H n , n = 1, . . ., 30 of the Laplace Beltrami operator, and for a grid of 100 thresholds in the interval (0, 0.1).
When the number of repetitions increases, the tails of the empirical distribution of the temporal mean of the empirical absolute errors become lighter, being the shape of the empirical distribution closer to a Gaussian distribution.This effect associated with the increasing of the sample size is more pronounced at highest spatial resolution levels (i.e., at high discrete Legendre frequencies).In particular, the empirical distribution of the temporal mean of the empirical absolute errors becomes more concentrated around its mean faster at higher than at coarser resolution levels.A similar behavior is observed in the asymmetric empirical distribution of quadratic error temporal mean, displaying a larger degree of dispersion, due to the stronger effect of the functional spectral singularity at zero frequency.Regarding the empirical probability analysis, one can observe that the increasing of R has a strong effect by spherical scales on the gradually decay of the empirical probabilities to zero over the smallest threshold values in the grid in the interval (0, 0.1), while increasing parameter T enlarges the dark blue area, in the contour plots displayed, where empirical probabilities become zero over the largest threshold values in the grid analyzed.Note that the opposite effects of parameter R by spherical scales, in the increasing and decreasing LRD operator eigenvalue scenarios, are observed.All the results displayed are affected by a numerical integration error.
The performance of the SRD-LRD estimation methodology is also analyzed for the remaining multifractionally integrated spherical functional processes generated as displayed in Figure 13.Specifically, in this figure, results in terms of the empirical mean quadratic errors, associated with SRD functional spectral estimation (left hand-side), and in terms of the histograms of the temporal mean of the empirical absolute errors (right hand-side), associated with LRD spectral estimation, are respectively showed from R = 300 independent generations of a functional sample of size T = 500.Results for T = 50, T = 100, T = 500, T = 1000, and R = 100 are displayed in Section 5 of the Supplementary Material.One can observe since T = 500 the order of magnitude of the empirical mean quadratic errors is 10 −4 for the bandwidth parameter B T = 0.2 in the SRD estimation.It is well-known that the bandwidth parameter affects precision of the weighted periodogram operator estimator, and the impact of parameter T is very strong.This fact can also be observed in Section 5 of the Supplementary Material, looking at differences in the magnitude of empirical mean quadratic errors, based on 100 repetitions, associated with the weighted periodogram operator for T = 50 and T = 100, considering bandwidth parameter B T = 0.1, as well as for T = 500, and T = 1000, considering bandwidth parameter B T = 0.2, since a substantial reduction in such magnitudes occurs when increasing the functional sample size T (see, e.g., Theorem 3.6 in [21]).Similar results as in previous section are obtained in terms of the empirical distribution of the temporal mean of the empirical absolute errors in the implementation of minimum contrast estimation methodology for n = 1, . . ., 15.

Final comments
Results displayed in Sections 5.1 and 5.2 (see also Sections 4 and 5 of the Supplementary Material) are based on the computation of the empirical distribution of the temporal mean of the absolute errors, and the empirical probabilities.The interaction between parameters n, R and T in the asymptotic analysis of the two estimation methodologies proposed is illustrated beyond the Gaussian scenario in the simulation study undertaken.Specifically, from the empirical distributions plotted, one can conclude that their rate of convergence is a function of the spherical scale n, the functional sample size T, and the number of repetitions R. Differences between empirical distributions of absolute errors by scales are more pronounced for decreasing sequence of LRD operator eigenvalues than in the case of increasing LRD operator eigenvalue sequence.As expected (see also Supplementary Material), the effect of the element of SPHARMA(p,q) process family considered is negligible for the minimum contrast estimation methodology.Namely, a slightly increasing of the concentration rate of the empirical errors when the parameter of autoregression p increases is observed.Additionally, the empirical probability analysis also reflects the interaction of these three parameters through the rate of converge to zero.Under nondecreasing eigenvalue sequence of the LRD operator, a smoother decay to zero of the empirical probabilities than in the case of decreasing eigenvalue sequence of the LRD operator is observed.A new battery of limit results will be investigated beyond the Gaussian scenario in a subsequent paper for the asymptotic analysis of the proposed estimators of the second-order structure of the LRD manifold cross-time RFs studied here.

Fig. 2 :
Fig. 2: The selected pole in the zonal functions

Fig. 4 : 30 Fig. 5 :Fig. 6 :
Fig. 4: The square root (s.r.) of the modulus of the projected regular factor of the spectral density operator (left-hand side), and the product of this factor with the s.r. of the modulus of the singular factor of the spectral density operator (right-hand side) at Legendre frequencies n = 1, . . ., 30

30 Fig. 7 :
Fig. 7: Histograms of the temporal mean of the empirical absolute errors from a functional sample of size T = 1000 (LRD operator decreasing eigenvalue sequence)

Fig. 13 :
Fig. 13: The empirical mean quadratic errors (E.M.Q.E.s) associated with the projected weighted periodogram operator (P.W.P.O.) estimator (eigenspaces Hn, n = 16, . . ., 30, of the Laplace Beltrami operator), are displayed at the left-hand side.The remaining plots provide the histograms of the temporal mean of the empirical absolute errors, associated with the minimum contrast parameter estimation of LRD operator defining multifractional integration (M.I), for eigenspaces Hn, n = 10, 15.All the results displayed are based on R = 300 independent generations of a functional sample of size T = 500.The bandwidth parameter B T = 0.2 has been chosen in all the cases ). E.M.Q.E. of P.W.P.O..M.Q.E. of P.W.P.O. E