Reduced {\alpha}-stable dynamics for multiple time scale systems forced with correlated additive and multiplicative Gaussian white noise

TL;DR

This paper investigates how a slow process driven by correlated additive and multiplicative noise can be approximated as being influenced by an alpha-stable noise, extending stochastic averaging to infinite-variance forcing scenarios.

Contribution

It introduces a method to approximate a slow process driven by CAM noise with an alpha-stable process, enabling stochastic averaging in infinite-variance contexts.

Findings

Slow process behaves as driven by alpha-stable noise under large time-scale separation.

Conditions identified for CAM process to influence the slow process.

Numerical results demonstrate the approximation's validity.

Abstract

Stochastic averaging problems with Gaussian forcing have been studied thoroughly for many years, but far less attention has been paid to problems where the stochastic forcing has infinite variance, such as an {\alpha}-stable noise forcing. It has been shown that simple linear processes driven by correlated additive and multiplicative (CAM) Gaussian noise, which emerge in the context of atmosphere and ocean dynamics, have infinite variance in certain parameter regimes. In this paper, we study a stochastic averaging problem where a linear CAM noise process in a particular parameter regime is used to drive a comparatively slow process. It is shown that the slow process exhibits properties consistent with being forced by a white {\alpha}-stable noise in the case of large time-scale separation. We identify the conditions required for the fast linear CAM process to have such an influence in driving a slower process, and then derive an (effectively) equivalent fast, infinite-variance process for which an existing stochastic averaging approximation is readily applied. These results are illustrated using a set of representative numerical results.

Equations

$$\psi (k)=\exp \left[-{\sigma }^{\alpha }|k{|}^{\alpha }\Xi (k;\alpha ,\beta )\right],$$

$$\Xi (k;\alpha ,\beta )$$

$$\phi (k)$$

$$d{x}_t$$

$$d{y}_t$$

$$d{y}_t=\frac{1}{ϵ}\left(L+E^2/2\right)y_{t}dt+\frac{1}{{ϵ}^{1/2}}\left[(E{y}_t+g)d{W}_{1,t}+bd{W}_{2,t}\right].$$

$$\frac{\partial p}{\partial t}=\mathcal{A}p,\mathcal{A}p$$

$$+\frac{1}{2}\frac{{\partial }^2}{\partial y^2}\left([(Ey+g{)}^2+b^2]p\right).$$

$$p_s(y)=$$

$$\times \exp \left(\frac{2g\nu }{b}\arctan \left(\frac{Ey+g}{b}\right)\right),$$

$$\nu =-\left(\frac{L}{E^2}+\frac{1}{2}\right)$$

$$\mathcal{N}$$

$$=\frac{2\pi (2b{)}^{-(2\nu +1)}\Gamma (2\nu +1)}{E\Gamma (\nu +1-i\frac{2g\nu }{b})\Gamma (\nu +1+i\frac{2g\nu }{b})}$$

$$p_s(y)\sim \frac{h(\mathrm{sgn}\left(y\right))}{|y{|}^{2(\nu +1)}}\text{mbox}as|y|\to \infty ,$$

$${\mathrm{ACD}}_y(\tau )=\log \left[\frac{\mathbb{E}\left[\exp (i(y_{t+\tau }-y_t))\right]}{\mathbb{E}\left[\exp (i{y}_{t+\tau })\right]\mathbb{E}\left[\exp (-i{y}_t)\right]}\right].$$

$$d{z}_t$$

$$d{L}_t^{(\alpha ,\beta )}$$

$${\mathrm{ACD}}_z(\tau )$$

$$\frac{1}{n^{1/\alpha }}\sum _{j=1}^n\left(R_j-{\overline{R}}_n\right)\underset{D}{\to }{\mathcal{S}}_{{\alpha }^{\ast }}({\beta }^{\ast },{\sigma }^{\ast })\text{mbox}for\alpha ={\alpha }^{\ast }.$$

$${\alpha }^{\ast }=2\nu +1=-2L/E^2.$$

$${\beta }^{\ast }=\tanh \left(\frac{\pi g\nu }{b}\right).$$

$${\sigma }^{\ast }$$

$$={\left(\frac{2\cosh (\pi g\nu /b)}{E^{({\alpha }^{\ast }+1)}{\alpha }^{\ast }\mathcal{N}}\Gamma (1-{\alpha }^{\ast })\cos \left(\frac{\pi {\alpha }^{\ast }}{2}\right)\right)}^{1/{\alpha }^{\ast }}.$$

$$\frac{\partial P}{\partial t}$$

$${\mathcal{D}}_z^{(\alpha ,\beta )}$$

$${\psi }_z(k)$$

$${\psi }_v(m,t)$$

$$v_t={\int }_0^t{z}_{s/ϵ}ds\underset{D}{\to }{\mathcal{S}}_{\alpha }\left(\beta ,{ϵ}^{\gamma }{t}^{1/\alpha }\frac{{\sigma }_z}{\theta }\right),$$

$${\int }_0^T$$

$$Y_j$$

$$u_{Y}(r)\sim\begin{cases}q^{-}|r|^{-(\alpha^{*}+1)}&{\mbox{for $r<-a$ }}\\ q^{+}r^{-(\alpha^{*}+1)}&\mbox{for $r>a$}\end{cases},$$

$${\beta }_Y\equiv \frac{q^{+}-q^{-}}{q^{+}+q^{-}}\approx \frac{h^{+}-h^{-}}{h^{+}+h^{-}}={\beta }^{\ast }.$$

$$CD(Y_k,Y_l)$$

$$CS(Y_k,Y_l)$$

$$S_{N_Y}=\left(N_Y^{-1/{\alpha }^{\ast }}\sum _{j=1}^{N_Y}{Y}_j\right),$$

$$p_S(r)=\frac{1}{2\pi }{\int }_{\mathbb{R}}\exp \left(-ikr\right){\psi }_S\left(k\right)dk,\text{mbox}where{\psi }_S(k)\sim {\psi }_Y^{N_Y}\left(\frac{k}{N_Y^{1/{\alpha }^{\ast }}}\right).$$

Full text

In this paper, we study a stochastic averaging problem where a linear CAM noise process in a particular parameter regime is used to drive a comparatively slow process. It is shown that the slow process exhibits properties consistent with being forced by a white $\alpha$-stable noise in the case of large time-scale separation. We identify the conditions required for the fast linear CAM process to have such an influence in driving a slower process, and then derive an (effectively) equivalent fast, infinite-variance process for which an existing stochastic averaging approximation is readily applied. These results are illustrated using a set of representative numerical results.

I Introduction

Stochastic differential equation (SDE) modelling techniques are used broadly in modern quantitative research in disciplines including, but not limited to, physics Gardiner (1985), finance Cont and Tankov (2004), and biology Bressloff (2014). They allow researchers to incorporate elements into dynamical models that are impossible or impractical to model explicitly due to their unpredictable nature and/or complexity, by representating their effects with stochastic processes. The most commonly used stochastic driving process is Gaussian white noise. The ubiquity of Gaussian white noise follows from the Central Limit Theorem (CLT), which states that the sum of a sufficiently large set of independent, identically-distributed random variables with finite variance converges in distribution to a Gaussian random variable. However, the assumptions necessary for the CLT are not always satisfied. For example, if the random variables have a density with power law tails such that the variance is not finite, then we must consider the Generalized Central Limit Theorem (GCLT) Feller (1966). The GCLT states that a sum of $n$ independent, identically-distributed random variables of this type, not necessarily having finite variance, converges in distribution to an $\alpha$-stable random variable as $n\to\infty$. The distribution $\mathcal{S}_{\alpha}(\beta,\sigma)$ of an $\alpha$-stable random variable depends on three parameters: the stability index $\alpha\in(0,2]$, the skewness parameter $\beta\in[-1,1]$, and the scale parameter $\sigma\in(0,\infty)$. Such distributions do not have in general a closed-form expression for their probability density functions (PDFs), but their characteristic functions have the form

[TABLE]

where

[TABLE]

The case where $\alpha=2$ is the only case of an $\alpha$-stable random variable with finite variance and corresponds to a Gaussian random variable with mean 0 and variance $2\sigma^{2}$.

In this paper, we consider stochastic dynamical systems with multiple time scales of the form

[TABLE]

where $t\geq 0$, $0<\epsilon\ll 1$, $\rho,E,g,b$, and $L<0$ are real non-zero constants, $f_{1},f_{2}$ are functions with $f_{2}(\tilde{x})\neq 0$ for any $\tilde{x}$ in the domain of $x_{t}$, and $dW_{1,t},\,dW_{2,t}$ are independent Gaussian white noise processes with $\mathbb{E}\left[dW_{p,t}\right]=0$, $\mathbb{E}\left[dW_{p,s}\,dW_{q,t}\right]=\mathbf{1}_{p=q}\delta(t-s)dt$ for $p,q=1,2$. The process $y_{t}$ is referred to as a CAM noise process Sardeshmukh and Sura (2009); Penland and Sardeshmukh (2012), described by a SDE of standard form with linear drift term and correlated additive and multiplicative (CAM) Gaussian white noise forcing. Such processes emerge naturally in highly-truncated projections of fluid mechanical systems Sardeshmukh and Sura (2009). The variable $y_{t/\epsilon}$ is a fast version of the process $y_{t}$, specifically satisfying

[TABLE]

Thus the system (4)-(5) represents a slow process $x$ coupled to a fast process $y_{t/\epsilon}$ reduced from a larger nonlinear system. In this study, we restrict our attention to $g\neq 0$, noting that for $g=0$ the additive and multiplicative noise is not correlated. As shown in Sardeshmukh and Sura (2009), the stationary distribution for $y_{t}$ is non-Gaussian, and has infinite variance if $-L\leq E^{2}$. Example time series for $y_{t}$ are plotted in Figure 1 to illustrate the variability that this process can exhibit in this parameter regime. One feature to notice is that the distribution of fluctuations in the realization are increasingly asymmetric around zero $y_{t}$ for larger values of $g>0$. This skewness is discussed further in Section II.

Our goal is to derive a one-dimensional stochastic dynamical system for a variable $X_{t}$ that weakly approximates $x_{t}$ for parameter ranges where $y_{t}$ does not possess second-order moments. Such a stochastic process $X_{t}$ is referred to as a stochastic averaging approximation (or stochastic homogenization approximation) for $x_{t}$. Stochastic averaging techniques offer systematic ways of reducing a dynamical system with multiple time scales to one with reduced dimension and fewer time scales whose properties weakly approximate those of the slow variable(s) of the full system. The prototypical stochastic averaging problem is to determine a SDE model governing the evolution of random variable $X_{t}$, such that $X_{t}$ weakly approximates $x_{t}$ where $dx_{t}=f(x_{t},z_{t/\epsilon})\,dt$, $0<\epsilon\ll 1$, and the process $z_{t}$ evolves on a faster time scale relative to $x_{t}$. The case where $z_{t}$ is driven by Gaussian white noise and has finite moments has been well-studied, with several established results for stochastic averaging approximations Khas’minskii (1966a, b); Papanicolaou and Kohler (1974); Borodin (1977); Freidlin and Wentzell (1984); Givon et al. (2004); Pavliotis and Stuart (2007). Some studies consider a fast chaotic driving process rather than a stochastic one Kantz et al. (2004); Mitchell and Gottwald (2012). The case where the fast variable is driven by an $\alpha$-stable noise process has received less attention, with some stochastic averaging approximations obtained in Srokowski (2011); Thompson et al. .

Stochastic averaging techniques have particular importance to climate modelling problems characterized by high dimensionality, chaotic dynamics, and multiple time scales Hasselmann (1976); Majda et al. (2001); Arnold et al. (2003); Monahan and Culina (2011); Mitchell and Gottwald (2012). The multiple time scale nature of climate variability is two-fold. Climate variability involves interactions between different Earth system components with broadly separated dominant timescales (e.g. atmosphere, ocean, cryosphere, land surface). Furthermore, each of these components are high-dimensional systems with a broad range of timescales of internal variability Saltzman (2002). Linear CAM noise processes of the form (5) have been used to understand certain various features of the probability distributions of observed climate variables Sura et al. (2005); Penland and Ewald (2008); Sura and Sardeshmukh (2008); Sardeshmukh and Sura (2009); Sura (2011). Since infinite variance processes can arise in simple systems like (5) with additive and multiplicative noise, the GCLT suggests that CAM noise dynamics could be a potential source of the appearance of $\alpha$-stable noise forcing in the long-term climate record Ditlevsen (1999). However, the results in Penland and Sardeshmukh (2012) demonstrate that linear CAM noise processes are not equivalent to $\alpha$-stable noise processes, so that the recent results of Thompson et al. can not be applied directly to (4)-(5). Rather, the GCLT indicates that similarities in the behavior of the time integrals of the CAM noise and $\alpha$-stable processes can be demonstrated over sufficiently long time scales.

The main result of this paper is to identify a system composed of $x_{t}$ in (4) coupled with a Ornstein-Uhlenbeck-Lévy process (OULP) $z_{t}$ that, under a stochastic averaging approximation, gives a slow proess $X(t)$ that weakly approximates the slow dynamics of $x_{t}$ in (4)-(5). Analogous to an Ornstein-Uhlenbeck (OU) process forced by Gaussian noise, the OULP $z_{t/\epsilon}$ is a one-dimensional process with a linear drift driven by an $\alpha$-stable white noise forcing $dL_{t}^{(\alpha,\beta)}\sim\mathcal{S}_{\alpha}(\beta,dt^{1/\alpha})$. The key to this main result is a comparison of the integral of $y_{t}$ in (5) with the integral of $z_{t}$. Once such a process $z_{t}$ is identified, the (N+) approximation for systems driven with additive $\alpha$-stable noise given in Thompson et al. can be used to derive a stochastic averaging approximation for the slow dynamics.

The remainder of the paper is organized as follows. In Section II, we give some properties of the linear CAM process $y_{t}$ and discuss its relationship with $\alpha$-stable distributions through the GCLT. In Section III, we identify an OULP $z_{t/\epsilon}$ whose integral weakly approximates that of $y_{t/\epsilon}$ in (5). Replacing $y_{t/\epsilon}$ with $z_{t/\epsilon}$, in Section IV we derive the stochastic averaging approximation for the slow variable $x_{t}$ in (4) by applying the (N+) stochastic averaging approximation as in Thompson et al. . As a minimum requirement, the characteristic time scale of $y_{t/\epsilon}$ must be at least an order of magnitude faster than $x_{t}$ for $X_{t}$ to provide an accurate weak approximation to $x_{t}$. However the exact time scale separation needed depends on the parameters of the CAM noise process, as the rate of convergence to the limiting $\alpha$-stable distribution depends on the tail behaviour of the random variables being summed in the integral of $y_{t}$. In Section V, we apply our approximation to specific examples, one linear system and two nonlinear systems of the form (4)-(5), illustrating the performance of the approximation.

II The CAM noise process and its relationship to $\alpha$-stable distributions

In this section, we state some basic properties of the CAM noise process $y_{t}$ (many of which are derived in Sardeshmukh and Sura (2009); Penland and Sardeshmukh (2012)), identify the parameter domains of interest, and discuss the relationship between the CAM noise distribution and the $\alpha$-stable distribution for these parameters.

II.1 Properties of the linear CAM noise process

First we note that $y_{t}$ can display a range of different behaviours, despite its relatively simple appearance. It reduces to an Ornstein-Uhlenbeck process (OUP) as $E\to 0$ and to a geometric Brownian motion when both $E/g,\,E/b\to\infty$. We focus on parameter ranges for which $y_{t}$ has a stationary probability density function without a finite variance. The relevant parameter ranges can be determined by considering the probability density function (PDF) $p(y,t)$ for (5), which satisfies the FKE

[TABLE]

The stationary PDF, $p_{s}(y)$ is obtained by solving the time-independent FKE, $\mathcal{A}p_{s}=0$. This function is given in Penland and Sardeshmukh (2012),

[TABLE]

where

[TABLE]

and $\mathcal{N}$ is the normalization constant, shown in Penland and Sardeshmukh (2012) to be given by

[TABLE]

where $\Gamma$ is the complex Gamma function. If $\nu\leq 0$, the first moment of the distribution (8) does not exist. We consider $\nu>0$ for the remainder of this paper, and hence the stationary mean of $y_{t}$ is equal to 0. The distribution is heavy-tailed Penland and Sardeshmukh (2012), since for large $|y|$, $p_{s}$ (8) decays according to a power law,

[TABLE]

where $h(s)=\frac{\exp\left(\frac{\pi g\nu}{b}s\right)}{\mathcal{N}E^{2(\nu+1)}}$. Eq. (12) indicates that $y$ with stationary density $p_{s}$ in (8) does not have a finite variance $0<\nu\leq\frac{1}{2}$. By the GCLT, this case is related to non-Gaussian $\alpha$-stable processes and so we restrict ourselves to corresponding values of $L$ and $E$ in (9), that is, $\frac{E^{2}}{2}<(-L)<E^{2}$.

To capture the serial dependence of $y_{t}$ (or, less formally, the memory of $y_{t}$) we use the autocodifference function (ACD function), the quantity for stochastic processes without a finite variance analogous to the autocovariance for processes with finite variance. The ACD for the stationary process $y_{t}$ is defined in terms of characteristic function of the process at different times Taqqu and Samarodnitsky (1994); Wyłomań ska et al. (2015):

[TABLE]

Figure 2 presents sample estimates of (13) for CAM processes corresponding to two sets of parameter values. For comparison, also shown are ACDs for the OULPs $z(t)$

[TABLE]

The parameter values of these OULPs were selected so that the asymptotic decay rates of their ACD functions match those of the corresponding CAM processes. We observe that $\log[{\rm Re}(\operatorname{ACD}_{y}(\tau))]$ of the CAM noise process is evidently a nonlinear function of the lag $\tau$, in contrast to $\log[{\rm Re}(\operatorname{ACD}_{z})]$ for the OULP which is close to being linear in $\tau$. Specifically, for shorter lag times $\tau$, $\operatorname{ACD}_{y}(\tau)$ decays more rapidly than the $\operatorname{ACD}_{z}$ for the OULP with the same long-time dependence structure, while for larger values of $\tau$, $\log[{\rm Re}(\operatorname{ACD_{y})]}$ asymptotes to approximately linear behavior. This serial dependence plays an important role in considering the asymptotic behavior of integrals of $y_{t}$ in Section III below.

II.2 Generalized central limit theorem applied to independent CAM random variables

The GCLT states that a scaled and shifted sum of independent, identically distributed (iid) random variables $R_{j}$ with density $u_{R}(r)\propto|r|^{-(\alpha+1)}$ as $|r|\rightarrow\infty$ with $\alpha\in(0,2)$, converges in distribution to an $\alpha$-stable random variable (rather than a Gaussian random variable) Feller (1966), denoted

[TABLE]

We illustrate how the parameters $\alpha^{*}$, $\beta^{*}$, and $\sigma^{*}$ are determined in the context where the independent random variables $\{R_{j}\}_{j=1}^{n}$ are drawn from the distribution $p_{s}$ given by (8). The stability index of the attracting distribution is determined from the exponent of the tail behaviour of $p_{s}$ in (12)

[TABLE]

For $0<\nu<\frac{1}{2}$, $p_{s}$ has infinite variance and mean zero as discussed in Section II.1. The skewness parameter $\beta^{*}$ and scale parameter $\sigma^{*}$ are determined following Keller and Kuske (2001). First, $\beta^{*}=(h^{+}-h^{-})/({h^{+}+h^{-}})$ where $h^{\pm}=h(\pm 1)$ with $h(s)$ in (12), which can be written simply as

[TABLE]

The scale parameter $\sigma^{*}$ is determined by comparing the arguments of the characteristic function of the normalized sum of $\{R_{j}\}_{j=1}^{n}$ (17) in the limit $n\to\infty$ to those of an $\alpha$-stable random variable (1) and is given by

[TABLE]

The analysis in Keller and Kuske (2001) of the sum (17) gives an approximate expression for its PDF in the form of the attracting $\alpha$-stable PDF plus a correction of $O(n^{1-2/\alpha^{*}})$ for $\alpha^{*}\in(1,2)$ and large $n$. Then $n^{1-2/\alpha^{*}}$ determines the rate of convergence to the $\alpha$-stable density for large $n$, with slower convergence for $\alpha^{*}$ closer to 2. If $\alpha^{*}=2$ ($\nu=1/2$), then the sum of $\{R_{j}\}_{j=1}^{n}$ converges to a Gaussian distribution with a correction that is asymptotically $O(1/\log(n))$ for large $n$ Keller and Kuske (2001). We do not treat the case $\alpha^{*}=2$ here, due to its logarithmic rate of convergence and the fact that the limiting distribution does not have power law tails.

The stochastic averaging approximation to which we now turn involves integrals of serially-dependent stochastic processes rather than sums of independent, identically-distributed random variables. Our central ansatz is that Eqns. (18)-(19) determine the stability index and skewness parameter of such integrals, so that only the scale parameter needs to be computed.

III Approximation of a linear CAM process by an OULP

As $y_{t}$ is not an $\alpha$-stable process, we cannot directly apply the averaging results of Thompson et al. to obtain the weak approximation $X_{t}$ to $x_{t}$ in (4). However, we note that the solution $x_{t}$ of (4) for $t\in(0,T)$ involves the integral of $y_{t/\epsilon}$ over this interval. The connection between $\alpha$-stable random variables and sums of independent random variables with the stationary distribution of the CAM process $y_{t}$ described in Section II.2 suggests that for $\epsilon\ll 1$, the integral $\int_{0}^{T}y_{s/\epsilon}ds$ should have a distribution that is close to $\alpha$-stable. In this limit, we can determine the parameters of an OULP $z_{t}$ such that the distribution of $\int_{0}^{T}z_{s/\epsilon}\,ds$ approximates that of $\int_{0}^{T}y_{t/\epsilon}ds$, and then apply the stochastic averaging results from Thompson et al. to the two-timescale system $(x_{t},z_{t/\epsilon})$. Our hypothesis is that the resulting stochastic process is a weak approximation to the slow variable $x_{t}$ in the original system $(x_{t},y_{t/\epsilon})$.

We cannot simply invoke the GCLT to make this claim however, since the integral of $y_{t}$ is not a sum of independent, identically-distributed random variables. Demonstrating this result involves arguing that the integral $\int_{0}^{T}y_{s/\epsilon}ds$ can be decomposed into the sum of a large number of effectively independent and identically-distributed random variables. As well, previous stochastic averaging results for fast-slow systems such as Monahan and Culina (2011) and Thompson et al. indicate that details of the serial dependence of the fast process appear directly in the stochastic averaging approximation. Specifically, when an $\alpha$-stable process with characteristic decay time scale $\tau$ drives the fast process, then the $\alpha$-stable process in the the stochastic averaging approximation depends on $\tau$. This dependence appears in the scale parameter used in the approximation, but not in its stability index and skewness parameter. Thus, we expect that any $\alpha$-stable stochastic forcing replacing the fast CAM noise process would depend on the memory of $y_{t}$ through the scale parameter.

With these facts in mind, we proceed to investigate the hypothesis that the system (4), (5) can be reduced by relating the fast driving linear CAM process $y_{t/\epsilon}$ to an OULP, $z_{t/\epsilon}$ with approximately equally-distributed integrals over the interval $(0,T)$. Then the stochastic averaging result from Thompson et al. is used to approximate the weak properties of the fast fluctuations in the equation for $x_{t}$ when there is a separation of the slow and fast time scales for $\epsilon\ll 1$. In order to identify the appropriate parameter values in the SDE for $z_{t/\epsilon}$, we review the statistical properties of $\int_{0}^{t}z_{s/\epsilon}\,ds$ as derived in Thompson et al. . These are complemented by the conditions under which the integral $\int_{0}^{t}y_{s/\epsilon}\,ds$, can be weakly approximated by an $\alpha$-stable random variable.

III.1 The distribution for the integral of an OULP

We derive the characteristic function for the integral $v_{t}=\int_{0}^{t}\,z_{s/\epsilon}\,ds$ ($dv_{t}=z_{t/\epsilon}\,dt$) where $z(t)$ is the OULP in (15). The $(v_{t},z_{t/\epsilon})$ system has a joint time-dependent PDF $P(v,z,t)$ satisfying the FKE,

[TABLE]

where $P(v,z,0)=\delta(v)\delta(z-z_{0}),\,t>0$. The operator $\mathcal{D}_{z}^{(\alpha,\beta)}$ is the fractional differentiation operator Chaves (1998). As shown in Appendix A.1, we obtain the joint characteristic function $\psi_{v,z}(m,k,t)$ for $(v_{t},z_{t/\epsilon})$ from (59), the Fourier transform of (21). We find that $\psi_{v,z}(m,k,t)=\psi_{z}(k)\psi_{v}(m,t)$ for $t=O(1)$ and $\epsilon\ll 1$, where

[TABLE]

for $\Xi(m;\alpha,\beta)$ as given in (2). The expression for $\psi_{z}$ in (23) is the characteristic function corresponding to the stationary density of the OULP $z_{t}$. We note that the ratio $\sigma_{z}/\theta^{1/\alpha}$ appears in $\psi_{z}(k)$, analogous to the expression for the standard deviation in a stationary OU process. For $z_{0}=0$, the form of $\psi_{v}$ is the same as the characteristic function of an $\alpha$-stable random variable with mean [math] and the scale parameter $\epsilon^{1-1/\alpha}\sigma_{z}\theta^{-1}t^{1/\alpha}$. As convergence in characteristic function implies convergence in distribution, the asymptotic distribution for $v_{t}$ is $\alpha$-stable,

[TABLE]

where $\gamma=1-1/\alpha$.

III.2 The distribution of the integral of CAM noise

We begin by considering the integral $\int_{0}^{t}y_{s/\epsilon}\,ds$ as a sum of $N_{Y}$ partitions of length $\Delta$, denoted by $\{Y_{j}\}_{j=1}^{N_{Y}}$,

[TABLE]

where $T=N_{Y}\Delta$. We make the following ansatz for the tail behaviour of $Y_{j}$: for any fixed non-zero value $\Delta$ and sufficiently small $\epsilon$, $Y_{j}$ has a PDF $u_{Y}$, where

[TABLE]

for some $a>0$ and $q^{+},q^{-}>0$.. The ansatz for $u_{Y}$ follows from noting that the integral (27) for $Y_{j}$ can be approximated by a sum of random variables $y_{m}$, for $Mh=\Delta$ and sufficiently small $h$. In Appendix B we use a Euler-Maruyama approximation of $y_{m}$ in (5), with $dW_{i,t}$ approximated by sequences of independent random variables $\xi_{i,m}\sim N(0,\sqrt{h})$ for $i=1,2$. We express $u_{Y}(Y_{j})$ in terms of $p_{s}(y)$ via a linear relationship between $Y_{j}$ and $y_{1}$ (68). This linear transformation leads not only to the exponent $\alpha^{*}$ in the ansatz (28), but also to an approximation of a parameter $\beta_{Y}$ by the skewness parameter $\beta^{*}$ in (12). Specifically, if $q^{\pm}$ are constant multiples of $h^{\pm}$ in (12), then

[TABLE]

Figure 3 illustrates the ansatz (28) and the skewness parameter $\beta_{Y}\approx\beta^{*}$ for representative examples of $Y_{j}$. The results show that indeed the values of $Y_{j}$ have the tail behaviour assumed by the ansatz that we assume. The rate of decay of both the positive and negative tails is consistent with a power law with exponent $-1-\alpha^{*}$ and the relative weight of the tails is as predicted for $|r|$ sufficiently large. Thus, we have empirical evidence for the validity of our ansatz, ((28) and (29)). Note that while we can predict the value of the ratio $(q^{+}-q^{-})/(q^{+}+q^{-})$ from the value of $\beta^{*}$, we do not have expressions for the values of the numerator or denominator separately.

The GCLT implies that a sum of independent variables $Y_{j}$ distributed according to (28) converges in distribution to an $\alpha$-stable random variable with stability index $\alpha^{*}$ and skewness parameter $\beta_{Y}\approx\beta^{*}$ as $N_{Y}\to\infty$. Thus, we expect to be able to approximate the distribution of $\int_{0}^{T}y_{s/\epsilon}\,ds$ with that of an $\alpha$-stable random variable with an error that decreases with $N_{Y}$, if the following two conditions are satisfied.

• •

Condition A) The value of $\Delta$ must be large enough such that $Y_{j}$ and $Y_{j+1}$ are (to a good approximation) independent for all $j\in\{1,2,\dots,N_{Y}-1\}$.

• •

Condition B) The number of partitions $N_{Y}=T/\Delta$ must be large enough such that the distribution of the normalized sum of $\left\{Y_{j}\right\}_{j=1}^{N_{Y}}$ is well-approximated by that of an $\alpha$-stable random variable.

If Conditions A) and B) are satisfied for $T=N_{Y}\Delta$, then $N_{Y}^{-1/\alpha^{*}}\sum_{j=1}^{N_{Y}}Y_{j}$ can be approximated by a $\alpha$-stable random variable (34). Furthermore, Conditions A) and B) taken by themselves appear to suggest that $T$ should be large; however, when considered within the larger context of approximating $x$ in (4) on the slow time scale, we see that $T$ must be of the same order of the characteristic time scale of $x$. In fact, these conditions can be satisfied by taking $\epsilon$ to be sufficiently small. The asymptotic approximations that use Conditions A) and B) are discussed below.

1. A)

The value of $\Delta$ must be large enough such that $Y_{j}$ and $Y_{j+1}$ are (effectively) independent for all $j\in\{1,2,\dots,N_{Y}-1\}$.

For the GCLT to apply to the sum $\sum_{j}Y_{j}$, the individual terms $Y_{j}$ must be independent. We intuitively expect $Y_{j}$ and $Y_{j+1}$ to be asymptotically independent (and thus also $Y_{j}$ and $Y_{j+m}$ for $m>1$) provided that $\Delta/\epsilon\gg\tau_{y}$, for $\tau_{y}$ a characteristic memory time scale of $y_{t}$ and $\Delta/\epsilon$ in the definition of $Y_{j}$ (27). For example, $\tau_{y}$ could be the $e$-folding time of the ACD function: $ACD_{y}(\tau_{y})=e^{-1}ACD_{y}(0)$. For $\operatorname{ACD}_{y}(\tau)$ decreasing with $\tau$ as shown in Figure 2, we see that $\tau_{y}=O(1)$.

Two random variables with infinite variance $Y_{k}$ and $Y_{l}$ are independent if both the codifference ($CD$) and the cosum ($CS$) vanish Gross (1994); Wyłomań ska et al. (2015), where

[TABLE]

Hence, we consider the asymptotic behaviour of $CD(Y_{j},Y_{j+1})$ and $CS(Y_{j},Y_{j+1})$, which decrease with increasing $\Delta/\epsilon$ as illustrated in Figure 4 for example CAM processes. It should be noted that the choice to normalize the values of $Y_{j}$ in the figure is done for stylistic purposes and does not affect the relative values of the cosum and codifference functions for different values of $\Delta/\epsilon$. We expect that the codifference for an $\alpha$-stable random variable with itself is equal to 2 when the scale parameter is normalized to 1 Wyłomań ska et al. (2015), as the results in Figure 4. For $CD(Y_{j},Y_{j+1})$ and $CS(Y_{j},Y_{j+1})$ approaching zero, we treat $Y_{k}$ for $k>j+1$ as independent $Y_{j}$, treating all $Y_{j}$ effectively independent. It is clear that for $|L|=O(1)$, we can expect effective independence of $Y_{j}$ and $Y_{j+1}$ provided that $\Delta/\epsilon\gtrsim 5$.

2. B)

The number of partitions $N_{Y}=T/\Delta$ must be large enough such that the distribution of the sum of $\left\{Y_{j}\right\}_{j=1}^{N_{Y}}$ is well-approximated by that of an $\alpha$-stable random variable. To characterize the error in relating $\sum_{j=1}^{N_{Y}}Y_{j}$ to an $\alpha$-stable random variable, we review the results from Keller and Kuske (2001) that gives the distribution of a sum $S_{N_{Y}}$ of heavy-tailed random variables $Y_{j}$,

[TABLE]

via its characteristic function $\psi_{S}$. The rate of convergence to this distribution depends on its tail behaviour. For $\Delta/\epsilon$ large enough for effective independence of subsequent values of $Y_{j}$, we can write the PDF $p_{S}(r)$ of $S_{N_{Y}}$ as

[TABLE]

and $\psi_{Y}$ is the characteristic function of $Y_{j}$. Then for large $N_{Y}$, using the result of Keller and Kuske (2001) for a sum of independent, identically distributed heavy-tailed random variables, the characteristic function $\psi_{S}(k)$ can be approximated from an expansion of $\psi_{Y}$ for small argument. Using the expansion (66) for $\psi_{Y}$ as derived in Appendix A.2 yields the approximation

[TABLE]

for $Q$ a constant. As $N_{Y}$ increases, (33) approaches the form of a characteristic function for an $\alpha$-stable distribution as given in (1)

[TABLE]

Then the distribution of $S_{N_{Y}}$ converges to an $\alpha$-stable distribution, implied by the pointwise convergence of their characteristic functions (as per Lévy’s continuity theorem Feller (1966)). The error term in (33) has coefficient $N_{Y}^{1-2/\alpha^{*}}$, characterizing the rate of convergence of $S_{N_{Y}}$ to an $\alpha$-stable distribution, as noted in Section II.2. This error increases substantially for $\alpha^{*}<2$ and $(2-\alpha^{*})\ll 1$.

Putting these results together, it follows from (34) that we can specify the parameters for the approximate distribution of $S_{N_{Y}}$ as

[TABLE]

and $T=N_{Y}\Delta$, and $N_{Y},\Delta$ are sufficiently large for Conditions A) and B) to hold. Note that because we do not know the value of $q^{+}+q^{-}$, Eqn. (35) does not allow us to compute the value of $\sigma_{Y}$.

Now, we demonstrate that $\sigma_{Y}$ is proportional to $\epsilon^{\gamma^{*}}\Delta^{1/\alpha^{*}}$ for $\gamma^{*}=1-1/\alpha^{*}$, which allows us to determine the dependence on $\epsilon$ and $T$ of the scale parameter of $\int_{0}^{T}y_{s/\epsilon}\,ds$. The integral (27) can be expressed by multiplying the sum $S_{N_{Y}}$ by $N_{Y}^{1/\alpha^{*}}$ which gives us an $\alpha$-stable distribution for $\int_{0}^{T}y_{s/\epsilon}\,ds$:

[TABLE]

It follows that

[TABLE]

for sufficiently large $T$. The integral $\int_{0}^{T}y_{\hat{s}}\,d\hat{s}$ has no dependence on $\epsilon$ or $\Delta$, so that

[TABLE]

This dependence of ${\sigma}_{Y}$ on $\epsilon$ and $\Delta$ is illustrated in Figure 5. Using this result, it is useful to rewrite (37) as

[TABLE]

Note that by (39), $\Sigma$ is independent of $\epsilon$ and $\Delta$ and depends only on the parameters of the CAM noise in (5). However, (35) and (39) do not provide the value of $\Sigma$, since the coefficients $q^{\pm}$ are not specified. We numerically estimate $\Sigma$ via the empirical characteristic function obtained from simulations of $Y_{j}$, as described in Appendix D.

Finally, we obtain the approximation result that on a sufficiently long time $T$,

[TABLE]

if the parameters in (15) are taken to be

[TABLE]

Here $\alpha^{*}$ and $\beta^{*}$ are given by Eqns (18)-(19) and $\Sigma$ is estimated empirically as described in Appendix D. Note that for the averaging approximation the choice of $\theta$ is arbitrary, as only $\Sigma$ enters the expression for the scale parameter of the integral of $y_{t/\epsilon}$. A reasonable choice is to use a characteristic inverse timescale of Eqn. (5) for $y_{t}$ and set $\theta=L+E^{2}/2$.

IV The stochastic averaging approximation

We now use the results from the previous section to find the SDE for a slow process $X(t)$ that weakly approximates $x(t)$. To accomplish this we find a weak approximation for the fast fluctuations that involve $y_{t/\epsilon}$ in the equation for $x(t)$. Following Thompson et al. , we make a change of variable $\eta_{t}=\mathcal{U}(x_{t})$ for $\mathcal{U}^{\prime}(x)=\frac{1}{f_{2}(x)}$, so that $y_{t/\epsilon}$ enters as an additive term in the equation for $\eta(t)$:

[TABLE]

Notice that the transformation $\mathcal{U}$ is invertible, since $f_{2}(x)\neq 0$ for any $x$ in the domain. We write the dynamics (44) in integral form,

[TABLE]

and consider the integral of $y_{t/\epsilon}$ over a time interval of length $T$, $\int_{0}^{T}y_{s/\epsilon}\,ds$.

Then we approximate $\int_{0}^{T}y_{s/\epsilon}\,ds$ in the SDE for $\eta(t)$ (45) using $\int_{0}^{T}z_{s/\epsilon}\,ds$ where $z_{t}$ satisfies (15) with $\theta$ and $\sigma_{z}$ in (42),

[TABLE]

While $T$ must be sufficiently large for (41) to hold, we also require that the approximation is appropriate on the characteristic timescale $\tau_{x}$ of $x_{t}$, consistent with the overarching goal to find an approximate equation for the dynamics on the slow variable $x_{t}$.

The system (46) and (15) is a slow-fast system for $\eta_{t}$ and $z_{t}$ for which we seek a stochastic averaging approximation, as an intermediate step to weakly approximating $x_{t}$. For completeness, we give the key components of the (N+) approximation based on Thompson et al. for the general slow-fast system,

[TABLE]

where $t\geq 0,0<\epsilon\ll 1$, $\rho$ is a constant, $f_{2}({u})\neq 0$ and $g_{2}({u})>0$ for any ${u}$ in the domain of $u_{t}$. Note that (15) is a special case of (48). The process $U(t)$ that weakly approximates the slow dynamics of $u(t)$ satisfies

[TABLE]

where $\overline{f}(U)=f_{1}(U)+[f_{2}(U)g_{1}(U)/g_{2}(U)]$. The ‘$\diamond$’ symbol indicates the Marcus interpretation of the stochastic differential terms Marcus (1978), which is analogous to the Stratonovich interpretation for Gaussian white noise driven systems. If $\sigma f_{2}/g_{2}$ is a constant, then (49) reduces to the Itō interpretation. If $\alpha=2$, then (49) is interpreted in the sense of Stratonovich. Details of theory and simulations are given in the Section V.3 and the references Thompson et al. ; Chechkin and Pavlyukevich (2014); Cont and Tankov (2004).

Applying the (N+) stochastic averaging approximation (49) to the $(\eta,z)$ system defined in (46) and (15), we obtain the SDE for the weak approximation $\tilde{\eta}_{t}$

[TABLE]

Taking the inverse of the transformation $\mathcal{U}$ as in Thompson et al. , we obtain the SDE for $X_{t}=\mathcal{U}^{-1}(\tilde{\eta}_{t})$ that weakly approximates $x_{t}$

[TABLE]

The value for $\Sigma$ used in this approximation is determined from the relationship (40) and is approximated numerically as discussed in Appendix D.

V Sample systems

In this section, we simulate one linear and two nonlinear systems of the form (4), (5) and compare the stationary PDFs and ACDs of the numerically simulated trajectories of $x_{t}$ to those of the corresponding stochastic averaging approximation, $X_{t}$ (51). We take $\rho=\gamma^{*}$ for convenience, in which case the power of $\epsilon$ in the noise coefficient of (51) vanishes.

V.1 Linear system

Let the dynamics of $x_{t},\,t\geq 0$ be given by

[TABLE]

where $\mu>0$, $\zeta$ is a constant, $x_{0}$ is known and $y_{t/\epsilon}$ is a fast linear CAM noise process with dynamics (5). Our stochastic averaging result (51) gives the reduced system,

[TABLE]

where $\sigma^{*}$ is given in (20). Not surprisingly, given the linear form for (52), $X_{t}$ is an OULP. We simulate $x_{t}$ numerically and compare its estimated PDF and ACD functions with the known PDF and ACD for the corresponding process $X_{t}$. The results are shown in Figure 6 for small $g$ (and thus small skewness parameter $\beta^{*}$) and in Figure 7 for larger $g\neq 0$. It is evident in both cases that the stationary distribution and ACD of the reduced model converge to that of the full system as $\epsilon$ decreases.

V.2 Nonlinear system 1: nonlinear potential

We illustrate our reduction method in a system with cubic nonlinearity in the slow equation and additive linear CAM noise driving:

[TABLE]

where $y_{t}$ is given by (5). This system is close to linear for small $x_{t}$, but experiences a stronger nonlinear drift for large $x_{t}$. For $\mu>0$, the system is globally attracted to $x=0$ in the noise-free limit, and for $\mu<0$, the origin is an unstable (repelling) equilibrium and the system is locally attracted to one of two stable equilibria at $x=\pm\sqrt{-\mu}$. According to our reduction results, we expect that the dynamics for $x_{t}$ can be weakly approximated by

[TABLE]

where $\Sigma$ is numerically estimated in the same fashion as described in Section III.2. We use a predictor-corrector method described in Thompson et al. to simulate both the full and reduced systems, avoiding numerical instabilities due to the cubic nonlinearity for large values of $x_{t}$ or $X_{t}$. The CAM noise process is integrated as described in Appendix C. The numerically estimated PDFs for $\mu>0$ and $\mu<0$ are compared in Figures 8 and 9, respectively. In both cases, the stochastic averaging approximation works well for sufficiently small $\epsilon$. The PDF tails from both the full and reduced systems show modest fluctuations due to sampling variability.

V.3 Nonlinear system 2: bilinear interaction

The second nonlinear system that we consider has a bilinear term in which the CAM noise process appears multiplicatively in the slow equation:

[TABLE]

It is clear from considering the limit as $x\to 0$ that the dynamics of $x_{t}$ are restricted to the positive real line, due to the inclusion of the state-independent drift term $c>0$. Applying our stochastic averaging result, we can weakly approximate the dynamics of (56) by $X_{t}$ where

[TABLE]

where ‘$\diamond$’ denotes the Marcus interpretation of the stochastic driving term. Details and references on the calculation of Marcus stochastic integrals used in simulating these SDE’s can be found in Thompson et al. ; Chechkin and Pavlyukevich (2014); Cont and Tankov (2004). To summarize, the Marcus stochastic term, $\zeta\Sigma X_{t}\diamond dL_{t}^{(\alpha^{*},\beta^{*})}$, can be written as $\lambda(1;\Delta L,X_{t-})-X_{t-}$ where $\lambda(s;\Delta L,X_{t-})=\lambda(s)$ satisfies

[TABLE]

with $\Delta L=dL_{t}^{(\alpha^{*},\beta^{*})}$ is the jump of the $\alpha$-stable process at time $t$ and $X_{t-}=\lim_{s\to t^{-}}X_{s}$. In this case, (58) can be solved for $\lambda(1)$, giving $\lambda(1;\Delta L,X_{t-})=\exp(\zeta\Sigma\,\Delta L)X_{t-}$. Having an analytic expression for the Marcus stochastic increment facilitates a straightforward and accurate simulation of (57). The estimates of the stationary PDFs of $x_{t}$ and the reduced process $X_{t}$ compare well (Figure 10).

V.4 Discussion of numerical results

From the results displayed in Figures 6 - 10, we see that the stationary behaviour of $x_{t}$ in each system under consideration can be well-approximated by a system influenced by $\alpha$-stable noise for sufficiently small values of $\epsilon$. Both the PDFs and the ACD functions of the full and reduced systems match well.

As expected for stochastic averaging approximations for systems with multiple time scales, as $\epsilon$ decreases, we see improvement of the approximation of the weak properties of $x_{t}$ in the full system by $X(t)$ in the reduced system. However, the difference between the PDFs and ACDs of the reduced systems and the full ones depends not only on $\epsilon$, but also on $\alpha^{*}$. Specifically, for $\alpha^{*}$ close to 2 the value of $\epsilon$ needs to be noticeably smaller, as we can see by reviewing Conditions A) and B) from Section III.2 for $\int_{0}^{T}y_{s/\epsilon}ds$ to be approximately distributed as $\alpha$-stable. A sum of independent heavy-tailed random variables converges to an $\alpha$-stable variable more slowly for values of $\alpha$ closer to 2 Keller and Kuske (2001). As such, the number of approximately independent terms in the discretization of the integral $\int_{0}^{T}y_{s/\epsilon}ds$ for fixed $T$ must be larger for $\alpha$ closer to two in order for this integral to have the desired properties - requiring a smaller value of $\epsilon$. The quality of the weak approximation of $x_{t}$ by $X_{t}$ is determined not only by the timescale separation, but also by the tail behaviour of the limiting $\alpha$-stable process.

VI Conclusions

This study demonstrates that $\alpha$-stable forcing can appear in the asymptotic slow dynamics of fast-slow systems where the fast process is linear and forced with a combination of additive and multiplicative noise such that its distribution has infinite variance. We studied the linear CAM process in this analysis due to its use in applied research problems and the fact that many analytical results are available Sardeshmukh and Sura (2009); Penland and Sardeshmukh (2012). A particularly valuable aspect of the linear CAM process is that it provides a simple dynamical form in which the interaction of multiplicative terms and Gaussian white noise forcing results in a process possessing a distribution with power law tails and infinite variance. Due to the infinite variance stationary behaviour of the linear CAM noise process and the explicit predictions of the power law tail behaviour, we hypothesized that the linear CAM noise process would appear to be equivalent to an $\alpha$-stable forcing term when used to drive a slower process. The derivation of the approximation required us to focus our attention on the properties of the integral of the linear CAM noise process and determining an OULP having integral behaviour with similar statistics. Formulas for the stability index and skewness parameter of the corresponding OULP in terms of the CAM process parameters are straightforward. In contrast, a coefficient in the OULP drift term needs to be estimated numerically as it depends on the serial dependence properties of the linear CAM noise process for which analytic results are not available. When we apply our stochastic averaging approximation to linear systems as well as systems with nonlinearities in the slow variable, we observe good agreement between simulations of the full and reduced dynamics, which improve as $\epsilon$ decreases. For values of the parameters such that the corresponding OULP has a stability index close to 2, the ratio of time scales between the fast and slow processes needs to be extremely large to observe the distributional convergence of the slow variable to the predicted distribution.

Besides presenting a method of approximating fast-slow systems that are forced with correlated additive and multiplicative noise processes, this analysis also suggests a possible mechanism through which $\alpha$-stable forcing can emerge in the modelling of physical problems. For example, as mentioned above, CAM noise processes emerge from considerations of the dynamics of quadratically nonlinear systems like those describing atmospheric motion Majda et al. (2001); Sardeshmukh and Sura (2009) and could result on longer time scale processes experiencing forcing terms distributed according to an approximately $\alpha$-stable law. This result offers a possible explanation for the observation of $\alpha$-stable noise in various climatic and fluid dynamical time series Huber et al. (2001); Seo and Bowman (2000); Ditlevsen (1999); Del-Castillo-Negrete (1998).

There are various extensions to this research that are worth exploring. All of the analysis given is for the case where both the fast and slow subsystems are univariate. Extensions of the results of this paper to higher-dimensional fast slow systems may yield unexpected challenges, but are required before this stochastic averaging method could be applied to more general research problems. Also, while there is broad utility in studying models where the fast linear CAM noise process perturbs the slow variable, it would be worth exploring the implications of nonlinear fast perturbations to the slow process or where the fast, infinite-variance process is something other than the linear CAM noise process. Another interesting situation would be to study the situation where the fast linear CAM process is conditionally dependent of the slow process, or in other words the slow variable influences the parameters of the linear CAM noise process. This would be a particularly useful scenario to consider in the context of climate modelling.

Acknowledgements.

AHM and RK acknowledge partial support from the Natural Sciences and Engineering Research Council (NSERC) Discovery Grant program. RK also was partially supported by a grant from the Simons Foundation for work carried out in part at the Isaac Newton Institute. WFT acknowledges support from the NSERC Alexander Graham Bell Canadian Graduate Scholarship program and the UBC Faculty of Science Graduate Award program.

Appendix A Characteristic functions

A.1 Calculations for $\psi$

We provide some details of the calculations for the characteristic functions arising in Section III.1.

Taking the Fourier transform of (21), gives a quasi-linear partial differential equation whose solution is the joint characteristic function $\psi_{v,z}(m,k,t)=\mathcal{F}[P](m,k,t)=\iint_{\mathbb{R}^{2}}\exp(ikz+imv)P(v,z,t)\,dvdz$.

[TABLE]

where $\Xi(\Lambda;\alpha,\beta)=1-i\beta\operatorname{sgn}\left(\Lambda\right)\tan(\pi\alpha/2)$ as given in (2). Solving (59) via the approach in Thompson et al. , the method of characteristics gives the solution for $\psi_{v,z}$

[TABLE]

where $\Lambda(r)=\frac{\epsilon m}{\theta}+\left(k-\frac{\epsilon m}{\theta}\right)e^{-\theta r/\epsilon}$. Following Thompson et al. , the integral term in (61) has the asymptotic behavior for $t=O(1)$ and $\epsilon\ll 1$,

[TABLE]

Note that this expression can be factored into separate functions of $k$ and $m$ for $\theta t/\epsilon\gg 1$, resulting in the expressions for $\psi_{v}$ and $\psi_{z}$ in (23 - 24).

A.2 Asymptotic behaviour of the characteristic function $\psi_{Y}(k)$ for small $k$

The characteristic function of $Y_{j}$ is given by the Fourier transform, $\psi_{Y}(k)$, of the PDF $u_{Y}(r)$ in (28)

[TABLE]

Here $\Xi$ is given in (2), $\beta_{Y}=\frac{q^{+}-q^{-}}{q^{+}+q^{-}}$, and we have used the facts $\int u_{Y}(r)dr=1$ and $E[Y_{j}]=0$ to rewrite the integrals as in Keller and Kuske (2001) in order to facilitate an expansion for small $k$. Using the result from Sato (1999)

[TABLE]

in an expansion for small $k$ as in Keller and Kuske (2001)) yields

[TABLE]

Appendix B The probability density $u_{Y}$

For each of the $y_{m}$ in $Y_{j}\approx\sum_{m}y_{m}h$, $m=1,2,\ldots M$, we use the Euler-Maruyama approximation for (5) to express $y_{m}$ in terms of $y_{1}$ and the random variables $\xi_{k,m}\sim N(0,\sqrt{h})$, $k=1,2$ in the approximation. Following this straightforward but tedious calculation, $Y_{j}$ takes the form

[TABLE]

where ${\cal C}_{1}$ is linear in $\Delta$, ${\cal{K}}_{i}$ for $i=1,2$ are $O(1)$ constants dependent on $E$. We use (68) in the expression for the density of $Y_{j}$ obtained via conditioning on $y_{1}$ and $\xi_{i},m$

[TABLE]

assuming $y_{1}$ is taken from the stationary distribution for $y$, $\rho_{i}$ is the density for $\xi_{i,m}$, and $\langle\cdot\rangle_{\xi}$ indicates expected value with respect to $\xi$. Using that ${\cal G}$ is linear in $y_{1}$, and that ${\cal C}_{i}$ for $i=2,3$ are sums of products of independent random variables $\xi_{i,m}$, we conclude that $u_{Y}$ has the behavior given in (28).

Appendix C Simulating the CAM noise process

We use the weak order 2.0 explicit method Kloeden and Platen (1992) to simulate the CAM noise process in (5). This weak numerical approximation takes the form for $\hat{y}_{n}=y_{n\delta{t}}$,

[TABLE]

where

[TABLE]

Here $\tilde{L}=L+E^{2}/2$ and $\delta{t}$ is the size of the discrete time step. The terms $\delta{W_{1,n}},\,\delta{W_{2,n}},\,n=0,1,2,\dots$ are independent Gaussian random variables with mean 0 and variance $\delta{t}$.

C.1 Consistency of simulations

Since $y_{t/\epsilon}$ evolves on the fast time scale, $\delta t$ must be chosen smaller than $\epsilon$ in order to resolve the fast dynamics. We compared the stationary PDF $p_{s}$ for $y_{t/\epsilon}$ given by (8) to numerical approximations to $p_{s}$ based on simulations of (5). As expected, $\delta t$ must be an order of magnitude smaller than $\epsilon$ in order to obtain relative errors of the numerical simulations to $p_{s}$ that are $O(10^{-2})$ or smaller. We found that the relative error of two cases, $\delta{t}=\epsilon/10$ and $\epsilon/100$ were of the same order of magnitude. Furthermore, we also considered whether the value of $\delta t$ has an effect on the estimates of $Y_{j}=\int_{(j-1)\Delta}^{j\Delta}y_{s/\epsilon}\,ds$. We compared the behaviour of the density of $Y_{j}$ for decreasing values of $\epsilon$. As in the case of the CAM noise process, we do not see noticeable changes in the density for $Y_{j}$ when we choose stepsize values smaller than $\delta t\approx\epsilon/10$ for simulating $y_{t}$. Therefore we simulate the fast CAM noise process with time discretization $\delta{t}\leq\epsilon/10$ throughout the paper.

Appendix D Estimating $\Sigma$

We estimate $\Sigma$ based on simulations of the integral of $\int_{0}^{T}\,y_{s/\epsilon}\,ds$ with a value of $T=O(1)$ and a trapezoidal method, simulating $y_{t}$ in (5) with a step size $\delta{t}\leq\epsilon/10$. To ensure that the estimate for $\Sigma$ avoids the potential sensitivities described in Section V.4 we use a value of $\epsilon$ smaller than the values used in the numerical examples of Section ($\epsilon=10^{-5}$). By choosing a small value of $\epsilon$, we seek an approximation for $\Sigma$ based on large $N_{Y}$. As highlighted in Section V.4, $N_{Y}$ increases with decreasing $\epsilon$ for $T$ fixed and $\Delta$ satisfying Conditions $A$ and $B$ in Section III. Then the distribution of $\int_{0}^{T}y_{s/\epsilon}\,ds$ is close to the $\alpha$-stable distribution (41), and the approximation for $\Sigma$ is obtained from (41) for given values of $\epsilon$ and $T$.

The estimate of $\Sigma$ is based on a least squares fit of the characteristic function of $\int_{0}^{T}\,y_{s/\epsilon}\,ds=\sum_{j}Y_{j}$, $\psi_{S}(l;\sigma_{Y})=\mathbb{E}\left[\exp\left(il\int_{0}^{T}\,y_{s/\epsilon}\,ds\right)\right]$, to the empirical characteristic function

[TABLE]

based on realizations of $S_{j}^{(T)},\,j=1,2,\dots,N_{S}$. The characteristic function $\psi_{S}$ is given in (66),

[TABLE]

The least squares fit of $\hat{\psi}_{S}(l)$ to $\psi_{S}(l;\sigma_{S})$ is obtained by the minimization

[TABLE]

From the estimate of $\sigma_{S}$, we obtain our estimate for $\sigma_{Y}=\sigma_{S}/N_{Y}^{1/\alpha^{*}}$ where $\alpha^{*}$ is defined in (18). We note that the sum is proportional to integral estimates of the square of the difference between $\psi_{S}$ and $\hat{\psi}_{S}$ on the interval $[l_{\rm min},l_{\rm max}]$. This numerical method is similar to the method of Koutrouvelis which estimates the parameters for $\alpha$-stable random variables via a least squares method and the characteristic function Koutrouvelis (1980). In our estimator, we use uniformly spaced points in the domain of $l$, but other studies have shown that a non-uniform spacing of points is optimal for such an estimator Besbeas and Morgan (2008). While we opted for an unweighted estimator, there are also weighted versions of the Koutrouvelis method that may be worth considering. Nonetheless, our estimator is sufficient for obtaining accurate estimates of $\sigma_{Y}$ under the assumption that $S_{j}$ are distributed according to an $\alpha$-stable law. We use this estimator to verify the scaling relationships for $\Sigma$ with respect to $\Delta$ and $\epsilon$ (i.e. (39), as per Figure 5), and estimate a value of $\theta$ via the relationship $\Sigma=\sigma^{*}/\theta$ (42) where $\sigma_{z}=\sigma^{*}$ to determine an equivalent OULp for the $y_{t}$ (as per Section III).

We use our estimator based on the empirical characteristic function $\hat{\psi}_{S}(l)$, rather than packages that can simultaneously estimate all parameters for a stable distribution such as the STBL package for MATLAB for values, which is based on the Koutrouvelis method Veillette (2012). Tests with STBL showed it to be unreliable in determining an estimate of the stability index $\alpha$ (which we know exactly) for values of $\alpha$ near 2, and consequently estimates of the corresponding scale parameter are not reliable.