# Crossings of Second-Order Response Processes Subjected to LMA Loadings

Authors: Thomas Galtier; Sayan Gupta; Igor Rychlik
Journal: Journal of Probability and Statistics
 (2010)

---

## Abstract

The focus of this paper is on the estimation of the crossing intensities of responses for second-order dynamical systems, subjected to stationary, non-Gaussian external loadings. A new model for random loadings—the Laplace driven moving average (LMA)—is used. The model is non-Gaussian, strictly stationary, can model any spectrum, and has additional flexibility to model the skewness and kurtosis of the marginal distribution. The system response can be expressed as a second-order combination of the LMA processes. A numerical technique for estimating the level crossing intensities for such processes is developed. The proposed method is a hybrid method which combines the saddle-point approximation with limited Monte Carlo simulations. The performance and the accuracy of the proposed method are illustrated through a set of numerical examples.


---

## Body

## 1. Introduction

Failures in randomly vibrating systems occur primarily in two different modes of gradual deterioration of the material properties resulting in fatigue type failure and/or due to overloading, when the structure response exceeds specified threshold levels for the first time. Quantification of the risk associated with a structural system requires probabilistic characterization of the structure response. The probability of the first passage type of failures can be estimated from the statistics of the extreme structure response. On the other hand, predicting the risk against fatigue type of failures requires the probability distribution of the amplitudes of the response cycles corresponding to various ranges. In either case, the corresponding statistic is related to the intensity of the upcrossing of levels. For smooth stationary processes, the upcrossing intensity,μ(u) of level u, is given by Rice's formula [1, 2], expressed as(1.1)μ(u)=∫0∞zfY(0),Ẏ(0)(u,z)dz,
where fY(0),Ẏ(0)(u,z) is the joint probability density function (j-pdf) of the response Y(0) and its instantaneous time derivative Ẏ(0). The applicability of (1.1) lies in the availability of the information on the j-pdf fY(0),Ẏ(0)(u,z). This is, however, rarely available.Exact information about the j-pdf,fY(0),Ẏ(0)(u,z), is available when the response is stationary and Gaussian. This is usually applicable when stationary Gaussian loads act on systems with very weak nonlinearities, enabling approximating such systems as time invariant linear systems. This simplification implies that the response is also stationary and Gaussian. The corresponding upcrossing intensity can thus be evaluated using (1.1), leading to(1.2)μ(u)=fze-(1/2)(u-E[Y(0)])2/V(Y(0)),
where fz=(1/2π)V(Ẏ(0))/V(Y(0)) and V(·) and E[·] indicate the variance and the expected value, respectively.The probability distribution of the extreme response in a fixed periodT, namely, MT=max0≤t≤TY(t), can be conservatively estimated by means of the inequality(1.3)P(MT>u)≤P(Y(0)>u)+Tμ(u).
See, for example, [3, 4]. For stationary Gaussian responses the stronger result that P(MT>u)/(Tμ(u))→1 as u tends to infinity is true, see [5]. Hence, for a long time the study of random loads has been dominated by Gaussian processes, that is, the dynamics of the system were linearized while external loads were modeled by means of Gaussian processes.However, there are situations where a simple linearization of weakly nonlinear, time invariant systems leads to approximations that are too crude. Such systems are often represented by means of Volterra functional expansion that is truncated after the second-order term. More precisely, we assume that with input forceX(t), the response Y(t) can be written as a sum(1.4)Y(t)=Y1(t)+Y2(t),
where (1.5)Y1(t)=∫-∞∞h1(s)X(t-s)ds,Y2(t)=12∬-∞∞h2(s1,s2)X(t-s1)X(t-s2)ds1ds2.In (1.5), h1(·) and h2(·,·), respectively, denote the linear and quadratic transfer functions of the system. Here, it can be assumed that X(t) is a smooth Gaussian process, given by(1.6)X(t)=∫-∞∞f(t-x)dB(x),
where B(x) is a Brownian motion while f(x) is a suitably chosen kernel. The pdf of responses and crossing properties of processes defined by (1.4), with Gaussian forcing, have been studied by many authors; see, for example, [6–9] and the more recent studies [10–14].However, many real loads, for example, ocean waves in shallow water or during heavy storms, show considerable non-Gaussian features, such as, a skewed marginal distribution with heavy tails. These waves are sometimes modeled by Volterra series expansions with Gaussian input, that is, a process of the same type asY(t) in (1.4). Statistical analysis of extremes of Y(t) when the forcing is quadratic is a difficult task. One approach would be to employ Monte Carlo simulations. However, to estimate the crossing intensities of very high levels, which in turn imply rare events, would require large number of simulation runs making Monte Carlo simulations prohibitively expensive.An alternative approach to modeling non-Gaussian forcing is to use a class of transformed Gaussian processes [15]. These processes take their starting point in a Gaussian process, Z(t), and a continuous and increasing function g(·). Then one forms a non-Gaussian process, X(t), according to the transformation X(t)=g(Z(t)). In this way, the process X(t) can have a non-Gaussian marginal distribution. Different strategies to choose the function g(·) have been proposed and studied in [16–19]. The drawback of this class of models is the inability to exactly model the spectral density function.In this paper, we consider another class of processes, the-so called Laplace moving averages (LMA), to model the forcing. These models are characterized by mean, spectrum (as in the Gaussian case), and two more parameters for skewness and kurtosis of the marginal distribution [20]. In this way, LMA processes offer an alternative to the transformed Gaussian models that is preserving the correct spectrum. Both simulating from the model and passing through linear filters are straightforward as the linear filtering does not lead outside of this class. In this paper we will study crossings of response Y(t), as defined in (1.4), with X(t) assumed to be an LMA process.The paper is organized as follows. First, in Section2, we introduce the LMA process and review some simple properties of this model. In Section 3, we define the response process, (1.4), with LMA forcing and develop the necessary equations. In Section 4, we present a method based on the saddle-point approximation to compute the crossing intensity of Y(t), given by μY(u), when the joint moment generating function of the response and its instantaneous time derivative is available. Subsequently, some numerical examples are presented in Section 5 to highlight the applicability of the developments proposed in this paper, and discussions on the accuracy of the estimates are presented. The salient features of the study carried out in this paper are highlighted in the concluding section.


## 2. The LMA Process

### 2.1. The Laplace Driven Moving Average Model

The model we propose for loads is a continuous time moving average which may be written as(2.1)X(t)=∫-∞∞f(t-x)dΛ(x),
where f(x) is a kernel function and Λ(x) is a stochastic process with independent and stationary increments having a generalized asymmetric Laplace distribution. The process Λ(x) is referred to as Laplace motion and the resulting process X(t) is called the Laplace driven moving average (LMA). Thus X(t) may be thought of as a convolution of f(·) with the increments of the process Λ(x). A process generated in this way is stationary and ergodic. In the special case where Λ(x) is chosen to be a Brownian motion, X(t) becomes a Gaussian process; otherwise, in general it is non-Gaussian.The generalized Laplace distribution is compactly defined by its characteristic function. More precisely, a random variableZ is said to have a generalized asymmetric Laplace distribution if its characteristic function is given by(2.2)ϕZ(v)=E[eivZ]=eivθ(1-iμv+σ2v2/2)1/ν.
Here, θ, μ∈ℝ and ν, σ>0 are parameters of the Laplace distribution and i=-1. If μ=0, the distribution is symmetric; otherwise it is asymmetric. An extensive overview of Laplace distributions is available in [21]. The generalized asymmetric Laplace distribution can be used to construct a process with independent and stationary increments—the previously mentioned Laplace motion. The Laplace motion Λ(x) is a process that starts at zero and whose distribution at x is given by(2.3)ϕΛ(x)(v)=E[eivΛ(x)]=eivζx(1-iμv+σ2v2/2)x/ν,
where ζ is a parameter representing the drift of the process. The Laplace motion can be extended to the whole real line by basically taking two independent copies of it and mirroring one of them in the origin. The extended process can then be used to define the moving average in (2.1). Since the increments of the Laplace motion are allowed to have an asymmetric distribution (μ≠0), it turns out that the corresponding moving average process will also have a nonsymmetric marginal distribution. In fact, the marginal distribution of the Laplace driven MA has the following characteristic function:(2.4)ϕX(t)(v)=exp(∫-∞∞iζvf(x)-1νlog(1-iμvf(x)+σ2f2(x)v22)dx),
where log(·) is the complex logarithm function.For the Laplace driven MA defined in (2.1), one can show that the mean and the two-sided spectral density S(ω) are given by(2.5)E[X(t)]=(ζ+μν)∫-∞∞f(x)dx,S(ω)=σ2+μ2ν12π|ℱf(ω)|2.
Here, ℱ denotes the Fourier transform. It must be noted that (2.5) is valid for any square integrable kernel. As symmetrical kernels f(·) have real Fourier transform F(ω), this implies that if S(ω)0.5 is integrable, a symmetrical kernel is obtained from (2.5). This means that by choosing different kernels one can, in principle, model any spectrum. By normalizing the kernels, such that ∫f(x)2dx=1, we get(2.6)V(X(t))=σ2+μ2ν.
However, after having chosen the kernel f(·) and fitting the mean and variance, there are still two free parameters, out of the four original ones. These “two degrees of freedom” can be used, for example, to fit skewness s and excess kurtosis κ (if κ>3) of the marginal distribution of Y(t). By using the expression for the characteristic function in (2.4), these are given by(2.7)s=μν1/22μ2+3σ2(μ2+σ2)3/2∫-∞∞f3(x)dx,κ=3ν(2-σ4(μ2+σ2)2)∫-∞∞f4(x)dx.This ability to fit both spectrum and the marginal skewness and kurtosis can be very useful when modeling second-order processes. Note that for a Gaussian process both skewness and excess kurtosis equal zero, that is,s=κ=0. In fact, a Gaussian process can be obtained from the Laplace driven MA as a limiting case as s=0 and κ→0, for example, by letting μ=0 and ν→0 in such a way that V(X(0)) in (2.6) is constant; see [21, page 183] for more detailed discussion. Consequently, in the following, we consider Gaussian moving averages as a special case of Laplace moving averages.


### 2.2. Simulation of the Laplace Driven MA

The Laplace driven moving average can be simulated in several different ways. The simplest and most straightforward one is to first simulate the increments of the Laplace motion over an equally spaced grid and then convolve it with the kernelf(·). In full generality, following [21], the asymmetric Laplace motion Λ(x), with drift ζ, asymmetry parameter μ, and variance σ2 can be represented as(2.8)Λ(x)=ζx+μΓ(x)+σB(Γ(x)).
Here, Γ(x) is a gamma process characterized by independent and homogeneous dx-increments having a gamma distribution with shape parameter dx/ν and scale parameter 1 while B(x) is Brownian motion. Using this representation, a simple algorithm for simulating the Laplace driven moving average with kernel f(·) is given by the following. (1)
Pickm and dx so that f(·) is well approximated by its values on(2.9)-mdx<⋯<-dx<0<dx<⋯<mdx.(2)
Pickn≫2m+1 in order to simulate the k=n-2m values of the response process Y at the points 0<dx<2·dx<⋯<(k-1)·dx.(3)
Simulaten identically and independently distributed (i.i.d.) Γ(dx/ν,1) random variables and store them in a vector G=[Gj].(4)
Simulaten i.i.d. zero mean standard normal random variables and store them in a vector Z.(5)
ComputeX=ζ∫f(x)dx+μf*G+σf*(G·Z), where G·Z=[Gj·Zj], * denotes convolution and the integral ∫f(x)dx is computed by some numerical method.The above simulation strategy is based on the fact that conditional on the Gamma process increments, the increment inGZ are independent and Gaussian, where G is the standard deviation. The advantage with the above simulation procedure is that it is very fast and efficient and that it works for long simulations and for most values of the parameters. The disadvantage is that one looses some resolution where the jumps in the Gamma process occur, due to taking an equally spaced grid.


## 2.1. The Laplace Driven Moving Average Model

The model we propose for loads is a continuous time moving average which may be written as(2.1)X(t)=∫-∞∞f(t-x)dΛ(x),
where f(x) is a kernel function and Λ(x) is a stochastic process with independent and stationary increments having a generalized asymmetric Laplace distribution. The process Λ(x) is referred to as Laplace motion and the resulting process X(t) is called the Laplace driven moving average (LMA). Thus X(t) may be thought of as a convolution of f(·) with the increments of the process Λ(x). A process generated in this way is stationary and ergodic. In the special case where Λ(x) is chosen to be a Brownian motion, X(t) becomes a Gaussian process; otherwise, in general it is non-Gaussian.The generalized Laplace distribution is compactly defined by its characteristic function. More precisely, a random variableZ is said to have a generalized asymmetric Laplace distribution if its characteristic function is given by(2.2)ϕZ(v)=E[eivZ]=eivθ(1-iμv+σ2v2/2)1/ν.
Here, θ, μ∈ℝ and ν, σ>0 are parameters of the Laplace distribution and i=-1. If μ=0, the distribution is symmetric; otherwise it is asymmetric. An extensive overview of Laplace distributions is available in [21]. The generalized asymmetric Laplace distribution can be used to construct a process with independent and stationary increments—the previously mentioned Laplace motion. The Laplace motion Λ(x) is a process that starts at zero and whose distribution at x is given by(2.3)ϕΛ(x)(v)=E[eivΛ(x)]=eivζx(1-iμv+σ2v2/2)x/ν,
where ζ is a parameter representing the drift of the process. The Laplace motion can be extended to the whole real line by basically taking two independent copies of it and mirroring one of them in the origin. The extended process can then be used to define the moving average in (2.1). Since the increments of the Laplace motion are allowed to have an asymmetric distribution (μ≠0), it turns out that the corresponding moving average process will also have a nonsymmetric marginal distribution. In fact, the marginal distribution of the Laplace driven MA has the following characteristic function:(2.4)ϕX(t)(v)=exp(∫-∞∞iζvf(x)-1νlog(1-iμvf(x)+σ2f2(x)v22)dx),
where log(·) is the complex logarithm function.For the Laplace driven MA defined in (2.1), one can show that the mean and the two-sided spectral density S(ω) are given by(2.5)E[X(t)]=(ζ+μν)∫-∞∞f(x)dx,S(ω)=σ2+μ2ν12π|ℱf(ω)|2.
Here, ℱ denotes the Fourier transform. It must be noted that (2.5) is valid for any square integrable kernel. As symmetrical kernels f(·) have real Fourier transform F(ω), this implies that if S(ω)0.5 is integrable, a symmetrical kernel is obtained from (2.5). This means that by choosing different kernels one can, in principle, model any spectrum. By normalizing the kernels, such that ∫f(x)2dx=1, we get(2.6)V(X(t))=σ2+μ2ν.
However, after having chosen the kernel f(·) and fitting the mean and variance, there are still two free parameters, out of the four original ones. These “two degrees of freedom” can be used, for example, to fit skewness s and excess kurtosis κ (if κ>3) of the marginal distribution of Y(t). By using the expression for the characteristic function in (2.4), these are given by(2.7)s=μν1/22μ2+3σ2(μ2+σ2)3/2∫-∞∞f3(x)dx,κ=3ν(2-σ4(μ2+σ2)2)∫-∞∞f4(x)dx.This ability to fit both spectrum and the marginal skewness and kurtosis can be very useful when modeling second-order processes. Note that for a Gaussian process both skewness and excess kurtosis equal zero, that is,s=κ=0. In fact, a Gaussian process can be obtained from the Laplace driven MA as a limiting case as s=0 and κ→0, for example, by letting μ=0 and ν→0 in such a way that V(X(0)) in (2.6) is constant; see [21, page 183] for more detailed discussion. Consequently, in the following, we consider Gaussian moving averages as a special case of Laplace moving averages.


## 2.2. Simulation of the Laplace Driven MA

The Laplace driven moving average can be simulated in several different ways. The simplest and most straightforward one is to first simulate the increments of the Laplace motion over an equally spaced grid and then convolve it with the kernelf(·). In full generality, following [21], the asymmetric Laplace motion Λ(x), with drift ζ, asymmetry parameter μ, and variance σ2 can be represented as(2.8)Λ(x)=ζx+μΓ(x)+σB(Γ(x)).
Here, Γ(x) is a gamma process characterized by independent and homogeneous dx-increments having a gamma distribution with shape parameter dx/ν and scale parameter 1 while B(x) is Brownian motion. Using this representation, a simple algorithm for simulating the Laplace driven moving average with kernel f(·) is given by the following. (1)
Pickm and dx so that f(·) is well approximated by its values on(2.9)-mdx<⋯<-dx<0<dx<⋯<mdx.(2)
Pickn≫2m+1 in order to simulate the k=n-2m values of the response process Y at the points 0<dx<2·dx<⋯<(k-1)·dx.(3)
Simulaten identically and independently distributed (i.i.d.) Γ(dx/ν,1) random variables and store them in a vector G=[Gj].(4)
Simulaten i.i.d. zero mean standard normal random variables and store them in a vector Z.(5)
ComputeX=ζ∫f(x)dx+μf*G+σf*(G·Z), where G·Z=[Gj·Zj], * denotes convolution and the integral ∫f(x)dx is computed by some numerical method.The above simulation strategy is based on the fact that conditional on the Gamma process increments, the increment inGZ are independent and Gaussian, where G is the standard deviation. The advantage with the above simulation procedure is that it is very fast and efficient and that it works for long simulations and for most values of the parameters. The disadvantage is that one looses some resolution where the jumps in the Gamma process occur, due to taking an equally spaced grid.


## 3. Quadratic Response Process with LMA Forcing

In this section, we employ a methodology developed in [6] to represent quadratic response processes with LMA forcing. The formulation closely follows the approach in [9], where asymptotical properties of the upcrossing intensity, μY(u), were studied for stationary process Y(t), as defined in (1.4), when X(t) is a stationary Gaussian process. Here, we consider the more general case where X(t) is modeled as an LMA process; see (2.1). Combining (1.4) and (2.1), the response process Y(t)=Y1(t)+Y2(t) can be rewritten as(3.1)Y1(t)=∫-∞∞q(t-x)dΛ(x),Y2(t)=12∬-∞∞Q(t-x1,t-x2)dΛ(x1)dΛ(x2).
Here, (3.2)q(t)=∫-∞∞h1(s)f(t-s)ds,Q(t,s)=∬-∞∞h2(s1,s2)f(t-s1)f(s-s2)ds1ds2.For most real life engineering applications, the kernelQ(·,·) is symmetrical. Further, we assume that the kernels q(·),Q(·,·) are square integrable and hence vanish at infinity. Thus, by choosing T sufficiently large, we may approximate the kernels by letting Q(s,t)=0 and q(s)=0 for |s|>T and |t|>T. Under such assumptions, the Kac-Siegert technique based on the representation of the truncated kernel Q(·,·) through its eigenfunctions ϕi(x) and eigenvalues λi, can be employed. Let the eigenfunctions and eigenvalues of the kernel Q(·,·) be defined by(3.3)∫-TTQ(t,s)ϕi(s)ds=λiϕi(t).For a symmetrical kernelQ(·,·), the eigenfunctions corresponding to the different eigenvalues are orthogonal. By further normalization, we assume that ϕi(·) are orthonormal with eigenvalues λi. Suppose that the eigenfunctions are ordered according to |λi|≥|λi+1|. Both eigenvalues and eigenfunctions are real, λi→0 as i→∞, and(3.4)∬-TT|Q(s1,s2)-∑i=1nλiϕi(s1)ϕi(s2)|2ds1ds2→0,asn→∞,
see [22]. Further, for simplicity of presentation, we assume that(3.5)∬-TTQ(s,t)q(s)q(t)dtds<∞,
and hence q(·) can be expanded in a series using the orthonormal eigenfunctions ϕi(·), namely,(3.6)q(s)=∑i=1∞aiϕi(s),ai=∫-TTϕi(s)q(s)ds.
From (3.6), it follows that for any n(3.7)q(s)=∑i=1n-1aiϕi(s)+∑i=n∞aiϕi(s)=∑i=1n-1aiϕi(s)+ãnϕ̃n(s),
where ãn is computed from the condition that ∫ϕ̃n2(s)ds=1. Let us introduce(3.8)W̃n(t)=∫-TTϕ̃n(t-x)dΛ(x).
In quadratic mean, the response in (1.4) can be rewritten as(3.9)Y(t)=∑i=1∞aiWi(t)+λi2Wi2(t),
where (3.10)Wi(t)=∫-TTϕi(t-x)dΛ(x)
are LMA processes.Often only a few of the eigenvaluesλi are significantly nonzero. Assuming that the number of such eigenvalues is n-1, (3.9) can be rewritten as(3.11)Y(t)=∑i=1n-1(aiWi(t)+λi2Wi2(t))+ãnW̃n(t).
For notational convenience, we drop the tilde and rewrite (3.11) as(3.12)Y(t)=∑i=1naiWi(t)+λi2Wi2(t),
where λn=0. The instantaneous time derivative process is given by(3.13)Ẏ(t)=∑i=1naiẆi(t)+λiẆi(t)Wi(t).
Note that (3.12) and (3.13) are functions of the vectors of LMA processes W(t)={Wi(t)}i=1n and Ẇ(t)={Wi̇(t)}i=1n. A procedure for estimating the upcrossing intensity for Y(t) in (3.12) is discussed in the following section.


## 4. Computing the Upcrossing Intensityμ(u)

The upcrossing intensityμ(u) of Y(t) can be computed using (1.1) if the j-pdf of Y(0) and Ẏ(0) is available. This, however, is not easy when Y(t) is defined as in (3.12). The elements in vectors W(t) and Ẇ(t) all have generalized Laplace distributions whose marginal pdfs are usually defined through their characteristic functions. Also, since the elemental processes Wi(t) and Ẇi(t) have mutual dependence, the computation of the joint characteristic function of Y(0) and Ẏ(0) is a difficult task. In the special case when Y(t) is an LMA-process, that is, Y(t)=Y1(t)=W1(t) see (1.4), and when n=1 in (3.12), it can be shown that the characteristic function can be expressed in an explicit manner; see later (5.2). In addition, as the moment generated function exists, the saddle-point method can be used for estimating the crossing intensity of the LMA processes [23]. The details of the saddle-point algorithm are available in the literature and for the sake of conciseness are not repeated here; the reader is directed to [10, 11, 14] for further details.In this paper, we extend the above method and develop a similar procedure for estimatingμ(u) for the general quadratic response process, Y(t), as defined in (1.4), but with LMA forcing. It must be noted that for general quadratic processes, Y(t), not only the characteristic functions are hard to compute but also the moment generating functions may not exist. Consequently, the application of the saddle-point method, or even methods employing characteristic functions, is not straightforward. Obviously one could use Monte Carlo (MC) approaches to simulate Y(t) or to estimate the joint density of (Y(0),Ẏ(0)) needed to compute μ(u) using (1.1). However, the MC approach is not an efficient way for computing μ(u) for high levels u, as the sample size, and in turn, the computational costs could be prohibitively large.Here, an alternative “hybrid” method is presented. The proposed method is a combination of Monte Carlo simulations and the saddle-point estimate. It uses the fact that conditionally on the Gamma process,W(t) and Ẇ(t) are normally distributed. Consequently, the computation of conditional moment generating function is straightforward and is given by(4.1)M(s,t∣γ)=E[esY(0)+tẎ(0)∣Γ(·)=γ(·)];
see Section 5.2. Now, the upcrossing intensity μ(u) can be expressed as the expectation of NY(u), that is, the number of upcrossings of level u by the process Y(t) in duration T=1. Thus, one can find the conditional upcrossing intensity of the process Y(t) when conditioned on the Gamma processes, and subsequently, the unconditional upcrossing intensity can be obtained as the first moment across the ensemble of Gamma processes. Mathematically, this can be written as(4.2)μ(u)=E[NY(u)]=E[E[NY(u)∣Γ(·)=γ(·)]].
The upcrossing intensity μ(u) can be estimated by computing the conditional moment generating function in (4.1) and using the saddle-point method to estimate E[NY(u)∣Γ(·)=γ(·)]. Subsequently, Monte Carlo simulations can be employed to estimate the unconditioned upcrossing intensity.The saddle-point algorithm is particularly efficient when the moment generating function,M(s,t), is symmetrical in t, that is, M(s,t)=M(s,-t). Note that the numerical algorithm presented in [10, 11, 14] is restricted to the symmetrical case. Unfortunately, the conditional moment generating function M(s,t∣γ) in (4.1) is not, in general, symmetrical. For the asymmetrical M(s,t), the algorithm is much slower, and further development of the method is needed before one can use it for a complex problem. As will be demonstrated in the following subsection, one can bypass this problem for time reversible processes. The sufficient condition for the time reversibility of the response process is that the kernels q(t) and Q(s,t) in (3.2) are symmetrical, which is what has been assumed in this paper.


### 4.1. Approximation of the Upcrossing Intensityμ(u)

Assuming thatY(t) is a stationary process, Y(t) and Y(-t) have the same expected number of upcrossings of any level u. Consequently, (4.3)Ỹ(t)=K·Y(t)+(1-K)·Y(-t),
where K is independent of the Y process and takes values 0 or 1 with probability 1/2. Additionally, Ỹ(t) has the same upcrossing intensity as the process Y(t). In the special case when Y(t) is given by (1.4) with LMA forcing, the upcrossing intensity can be expressed as(4.4)μ(u)=E[NY(u)]=E[NỸ(u)]=E[E[NỸ(u)∣Γ(·)=γ(·)]].Let the conditional crossing intensity be defined as(4.5)μ(u∣γ)=E[NỸ(u)∣Γ(·)=γ(·)].
Then, by simulating a sequence of Gamma processes, γi(·), i=1,…,N, the unconditional crossing intensity, μ(u), can be estimated by averaging μ(u∣γ), namely,(4.6)μ(u)≈1N∑i=1Nμ(u∣γi),
where N is the number of sequence of Gamma process simulated.The problem that needs to be addressed next is to develop a strategy for computing the conditional level crossing intensityμ(u∣γi). For time reversible Ỹ(t) in (4.3), since the conditional moment generating function can be expressed as(4.7)MỸ(s,t∣Γ(·)=γ(·))=12MY(s,t∣γ)+12MY(s,-t∣γ),
it is obvious that MỸ(s,t∣Γ(·)=γ(·)) is symmetrical. This enables one to use the saddle-point algorithms discussed in [10, 11, 14] to estimate the conditional upcrossing intensity μ(u∣γi).Clearly, the method to estimate the upcrossing intensityμ(u) proposed here is a hybrid method which combines Monte Carlo simulations of realizations of Gamma processes and the saddle-point approximation of upcrossing intensity. The advantage of this approach is that one can approximate crossings of extremely high levels (required when computing the extremes of responses with 100-year-return period) which is otherwise difficult if one employs Monte Carlo simulations only. The unresolved issue of the accuracy of the proposed hybrid method will be examined in the following section.


## 4.1. Approximation of the Upcrossing Intensityμ(u)

Assuming thatY(t) is a stationary process, Y(t) and Y(-t) have the same expected number of upcrossings of any level u. Consequently, (4.3)Ỹ(t)=K·Y(t)+(1-K)·Y(-t),
where K is independent of the Y process and takes values 0 or 1 with probability 1/2. Additionally, Ỹ(t) has the same upcrossing intensity as the process Y(t). In the special case when Y(t) is given by (1.4) with LMA forcing, the upcrossing intensity can be expressed as(4.4)μ(u)=E[NY(u)]=E[NỸ(u)]=E[E[NỸ(u)∣Γ(·)=γ(·)]].Let the conditional crossing intensity be defined as(4.5)μ(u∣γ)=E[NỸ(u)∣Γ(·)=γ(·)].
Then, by simulating a sequence of Gamma processes, γi(·), i=1,…,N, the unconditional crossing intensity, μ(u), can be estimated by averaging μ(u∣γ), namely,(4.6)μ(u)≈1N∑i=1Nμ(u∣γi),
where N is the number of sequence of Gamma process simulated.The problem that needs to be addressed next is to develop a strategy for computing the conditional level crossing intensityμ(u∣γi). For time reversible Ỹ(t) in (4.3), since the conditional moment generating function can be expressed as(4.7)MỸ(s,t∣Γ(·)=γ(·))=12MY(s,t∣γ)+12MY(s,-t∣γ),
it is obvious that MỸ(s,t∣Γ(·)=γ(·)) is symmetrical. This enables one to use the saddle-point algorithms discussed in [10, 11, 14] to estimate the conditional upcrossing intensity μ(u∣γi).Clearly, the method to estimate the upcrossing intensityμ(u) proposed here is a hybrid method which combines Monte Carlo simulations of realizations of Gamma processes and the saddle-point approximation of upcrossing intensity. The advantage of this approach is that one can approximate crossings of extremely high levels (required when computing the extremes of responses with 100-year-return period) which is otherwise difficult if one employs Monte Carlo simulations only. The unresolved issue of the accuracy of the proposed hybrid method will be examined in the following section.


## 5. Numerical Examples and Discussions

First, we consider an LMA process, that is, whenY(t)=Y1(t) for which the (unconditional) saddle-point method can be used. For such cases, the saddle-point method is very accurate, see [23], and the computed estimate can be used to benchmark the accuracy of the proposed method. This will allow us to study how large N in (4.6) should be in order to reach desired accuracy.Next, we study the crossings of a simple quadratic responseY(t)=Y1(t)+λY1(t)2/2. The upcrossing intensity can be computed when upcrossing intensity of the linear response Y1(t) is known. Since the intensity can be very accurately computed by means of the saddle-point method, one can now study the convergence of (4.6) with reference to the quadratic process.Finally, we consider an example ofY(t) of full complexity and estimate the upcrossing intensity. Here, 12 eigenvalues λi differ significantly from zero. The computed crossing intensity is compared with the Monte Carlo estimate. The details of these numerical examples are elaborated in the following subsections.


### 5.1. Saddle-Point Approximation of Crossing Intensity for LMA Processes

We consider the crossings of a linear response process, given by(5.1)Y(t)=∫q(t-x)dΛ(x),
with symmetrical kernel q(·). The corresponding moment generating function is given by [20](5.2)M(s,t)=exp(ζ∫-∞∞sq(x)+tq̇(x)dx)·exp(-1ν∫-∞∞log(1-μ(sq(x)+tq̇(x))-σ22(sq(x)+tq̇(x))2)dx).
Since M(s,t)=M(s,-t), one can use the efficient algorithm of the saddle-point method discussed in [10, 11, 14].In order to simplify the presentation we introduce the following notations;μNs(u) is the estimate of μ(u)=E[N(u)] computed by means of the hybrid saddle-point method and (4.6)-(4.7), and μs(u) denotes the estimate of E[N(u)] by means of saddle-point method and M(s,t) defined in (5.2). Here N(u) is defined as the observed number of upcrossings of level u divided by the length of the “observation” time. In all the examples, N(u) has units Hz.We first focus on the computation of the conditional moment generated functionM(s,t∣γ). Let us consider two LMA processes defined by a common Laplace motion. More precisely, for two kernels f1(·),f2(·) and the Laplace motion Λ(x), define(5.3)X1(t)=∫f1(t-x)dΛ(x),X2(t)=∫f2(t-x)dΛ(x).Here,Λ(x) is defined as in (2.8). Now, conditionally that Γ(·)=γ(·), the Laplace motion can be written as(5.4)λ(x)=ζx+μγ(x)+σB(γ(x))
and hence the conditional LMA processes, X1(t) and X2(t), can be represented as(5.5)X1(t)=∫f1(t-x)dλ(x),X2(t)=∫f2(t-x)dλ(x),
respectively. Obviously for any t, (here we take t=0), the joint pdf of X1(0) and X2(0) is Gaussian with means and covariances mi, σij, i,j=1,2, given by(5.6)mi=ζ∫fi(x)dx+μ∫fi(x)dγ(x),σij=σ2∫fi(x)fj(x)dγ(x).
Using (5.6), with f1(x)=q(x) and f2(x)=q̇(x), leads to(5.7)M(s,t∣γ)=E[esY(0)+tẎ(0)∣Γ(·)=γ(·)]=exp(sm1+tm2+0.5s2σ11+0.5t2σ22+stσ12).Example 5.1.
In this example,30 minutes of measured stress in a ship under stationary severe sea conditions is modeled as an LMA process. A part of the stress is shown in Figure 1(a). One can clearly see the existence of high-frequency oscillations, likely due to whippings, which get superimposed with the wave-induced stress. Figure 1(b) illustrates an estimated spectrum, S(ω), having two peaks. The kernel q(x) is computed from the spectrum S(ω) by inverting (2.5). The choice of kernels for LMA processes is quite important. Though there can be many kernels which give the same spectrum, there can be only one symmetrical kernel. In this paper, we consider only symmetrical kernels. Thus, this unique symmetrical kernel is obtained by finding the inverse of (2.5). Such spectra give time reversible loads. However, these loads do not always match with observed loadings. More work is needed to determine the kernel from observed data. The symmetrical kernel obtained by imposing these conditions is shown in Figure 2.Example5.1: (a) A part of the measured stress. (b) An estimate of the spectral density for the measured stress.
(a)(b)Figure 2
Example5.1: Estimated kernel function q(x).We next need to estimate the parameters of the LMA process. The variance of the stress time history is obtained by integrating the spectrum,S(ω), with respect to ω. Additionally, we assume that stress time history is mean zero. In order to identify the remaining parameters of the LMA process, we compute the skewness and excess kurtosis which are 0.13 and 0.21, respectively. These values indicate that the stress process is slightly non-Gaussian.Figure3 illustrates the crossing intensity N(u) for the measured stress (solid irregular line). In prediction of extremes, the crossing intensity needs to be extrapolated to much higher levels. Here, the LMA model is used for the extrapolation. The crossings of LMA are estimated by means of μs(u), that is, the saddle-point method, where the moment generated function, M(s,t) has been defined in (5.2). The function μs(u) is shown in the plot as a dashed dotted line. The agreement between N(u) and μs(u) is seen to be very good, except at the highest observed values of N(u). These discrepancies can be attributed to extremely large whipping effects, which consist of several crossings of high levels. This effect is averaged in μs(u).Figure 3
Example5.1: Observed crossing intensity N(u) in the measured stress-solid line; simulated crossing intensity in 100 times longer signal than measured (50 hours)—solid line with dots; the saddle-point approximation μs(u) of E[N(u)]—dashed dotted line.In order to verify this claim, we simulated the LMA process for a much longer duration (50-hour period) and computed the crossing intensities. The resulting crossing intensity, N(u), is superimposed in Figure 3 by solid line with dots. We observe that the estimated crossing intensities follows closely those computed using the saddle-point method. This confirms the accuracy of the saddle-point method.The primary objectives of this example are(a)
to study the applicability of the approximationμNs(u) (computed by means of the saddle-point method and formulas (4.7)-(4.6)) and to predict the return values, that is, levels uT such that E[N+(uT)]=1/T, and(b)
to examine how fastμNs(u) converges to E[N+(uT)], which here is estimated by μs(u).These are slightly different problems since in (a), one is interested in the horizontal distance betweenμNs(u) and μs(u), when plotted against levels u, while in (b), one examines the vertical distance between the functions. The conclusions of these studies are illustrated by means of Figures 4(a) and 4(b). In Figure 4(a), we observe that even for as low N=102, one gets relatively small errors (about 10%) in predictions of uT. However, the vertical convergence is slower and one needs about N=104 simulations of γi to get satisfactory distance between the two lines; see Figure 4(b), where the fractions μNs(u)/μ∞s(u) for N=102, 103 and 104, are presented. The algorithm is relatively fast and one can use high values of N to obtain satisfactory accuracy levels. We further note that in Figure 4(a), the computed crossing intensities towards the right extremes are better with 102 samples than with 103 samples. This apparent contradiction can be explained by the fact that the crossing intensities determined in the procedure above are statistical estimates. The standard deviations of the estimates obtained by 102, 103 and 104 samples are, respectively, c/102, c/103, and c/104, where c is higher for higher values of u. Thus, for a fixed value of c, the standard deviations of the estimates with 103 samples and 104 samples are, respectively, 0.32 and 0.1 times the standard deviation of the estimate with 102 samples. Figure 4(a) illustrates the computed intensities for only one set of N values. If the exercise was repeated, we expect to see a much smaller spread in the computed crossing intensities for higher values of N. Alternatively, a confidence band could be computed by means of the parametric bootstrap.Example5.1: (a) Crossing intensities μNs(s,t) computed using the proposed hybrid method: Sample size for simulated gamma processes γi; N=102: dotted line; N=103: dashed dotted line; N=104: dashed line; Crossing intensity μs(u): solid line. (b) Corresponding relative errors μNs(u)/μs(u).
(a)(b)


### 5.2. Computation ofM(s,t) for the Quadratic Response

The general quadratic response is only notationaly more complex and we will proceed in a similar way as for the LMA process discussed in Example5.1. First, we need to find the conditional moment generating function(5.8)M(s,t∣γ)=E[esY(0)+tẎ(0)∣Γ(·)=γ(·)],
which can be written by an explicit formula; see (5.11) derived below. Then one can simulate a sequence of gamma processes, γi(·), i=1,…,N and as before approximate M(s,t) by means of (4.7).LetΛ be a diagonal matrix with the diagonal elements being denoted by λi, i=1,…,n, and the rest of the elements being zero. Using matrix notation, the response process can be written as(5.9)Y(t)=aW(t)T+12W(t)ΛW(t)T=∑j=1najWj(t)+12∑j=1nλjWj2(t),
where a=(a1,…,an) (an=1,λn=0). As discussed earlier, conditionally on Γ(·)=γ(·), the vectors W=W(0) and Ẇ=Ẇ(0) are normally distributed with means m and ṁ and covariance matrices Σ11, Σ12, and Σ22, where for 1≤i,j≤n, (5.10)σ11(i,j)=σ2∫ϕi(x)ϕj(x)dγ(x),σ12(i,j)=σ2∫ϕi(x)ϕ̇j(x)dγ(x),σ22(i,j)=σ2∫ϕ̇i(x)ϕ̇j(x)dγ(x),m(i)=ζ∫ϕi(x)dx+μ∫ϕi(x)dγ(x),ṁ(i)=ζ∫ϕ̇i(x)dx+μ∫ϕ̇i(x)dγ(x).
Once the matrices Σij and vectors m and ṁ are computed, it is a straightforward task to compute M(s,t∣γ), see [10], which is given by(5.11)M(s,t∣γ)=1det(Σ)exp(ms+ṁt+12t2m̃Vm̃+12tTΣ-1t).
Here, (5.12)Σ=Σ11-1-sΛ-t(ΛΣ21Σ11-1+Σ11-1Σ12Λ),t=sm̃+tṁΛ+tm̃Σ11-1Σ12+t2m̃VΛ,m̃=a+mΛ,m=amT+12mΛmT,ṁ=aṁT+ṁΛmT.Remark 5.2.
Example5.1 is obtained as a special case when n=1, with ϕ1(s)=q(s) and λi=0 while a1=1 in (5.9). Under these conditions, using simple algebraic manipulations, it can be shown that the conditional moment generated function is equal to the expression in (5.7).Example 5.3.
In this example, we focus on checking the accuracy of the estimates of the level crossing intensity,μNs(u), using the proposed hybrid method, for quadratic response Y(t) in (5.9) for the special case when n=2 and ϕn=0, that is,
(5.13)Y(t)=Y1(t)+λY12(t)2=λ(Y1(t)+1/λ)22-12λ.
Considering the case n=2 provides certain advantages which can be exploited to benchmark the accuracy of the estimates, μNs(u), using the proposed method. Using (4.3) and (5.13), it can be shown that the crossing intensity μY(u)=E[NY(u)] can be expressed as
(5.14)μY(u)=μY1(-1λ+2uλ+1λ2)+μY1(-1λ-2uλ+1λ2).
As can be seen from (5.14), the accuracy of the estimate μY(u) depends on the estimate of the crossing intensity μY1. This, however, poses no problem as this can be very accurately obtained using the direct saddle-point method. Thus, replacing μY1 in (5.14) by the saddle-point estimate, μY1s, the expression in (5.14) can be used to benchmark the accuracy of the level crossing estimate, μNs(u), obtained using the proposed hybrid method.
As in Example5.1, Y1(t) is a stress time history of duration of 30 minutes measured in a particular location of a ship impinged by ocean waves during the course of its journey; see Figure 1(a). We use the LMA process described in Example 5.1 to model Y1(t). For the quadratic response, we choose λ=0.01. This value is chosen so that the contribution of linear and quadratic parts to Y1 is similar; note that standard deviation of Y1(t) is about 47 MPa. An estimate of the crossing intensity μY(u) is obtained using  (5.14) and is shown in Figure 5(a). The accuracy of the crossing intensities for the corresponding levels, μNs(u), obtained using the proposed method are determined by comparing with these values.
In order to computeμNs(u), one needs the expression for the conditional moment generating function M(s,t∣γ). This is given in (5.11)-(5.12), with Σ11=σ11, a=1, Σ22=σ22, Σ12=σ12. All parameters have the same values as in Example 5.1. A comparison of the crossing intensity estimates, μNs, using the proposed hybrid method is illustrated in Figure 5(a). As in Example 5.1, we consider the three cases where N=102, 103, and 104, where N is the number of gamma process simulations in the proposed hybrid method. A comparison of the relative errors is shown in Figure 5(b). As in Example 5.1, we observe that the estimates are in fairly good agreement with the accuracy expectedly improving for larger values of N.Example5.3: (a) Crossing intensities μNs(u) computed using the proposed hybrid method: sample size for simulated gamma processes γi; N=102: dotted line; N=103: dashed dotted line; N=104: dashed line; μs(u): solid line. (b) Corresponding relative errors μNs(u)/μs(u).
(a)(b)Example 5.4.
In this example, we consider a more general quadratic response process, such that the number of termsn in  (5.9) is more than one. We consider the response process Y(t)=Y1(t)+Y2(t) defined in (3.1), where q(s)=exp(-s2/50)/25π, -25≤s≤25 and
(5.15)Q(t,s)=0.01exp(-(s-t)250).
The parameters in Laplace motion, Λ(x), are chosen in such a way that the linear response, Y1(t)=∫-TTq(t-s)dΛ(s), has mean zero, variance one, skewness 0.5, and kurtosis 4.5. For the kernel Q(t,s), the first 12 eigenvalues were found to be significantly nonzero. To determine the number of such eigenvalues, the first 100 eigenvalues were found and ordered according to their absolute values, and their corresponding ratios with respect to their total summation were calculated. It was assumed that the series could be truncated when the sum of the absolute value of the eigenvalues exceeded 99.9% of the total sum. This led to n=12 for this example.
Based on experience from Examples5.1 and 5.3, we expect that N=1000 simulations of γi are needed for arriving at a reasonably accurate estimate of μNs(u) using the proposed hybrid method. In the absence of any closed-form analytical solutions for the crossing intensities of the quadratic response process, we compare the estimates obtained using the proposed hybrid method with those obtained from Monte Carlo simulations. For Monte Carlo simulations, simulating a large number of response processes and checking for their crossing intensities would be computationally very expensive and time consuming. Instead, we adopt the following MC procedure. (a)
1×107 independent samples of pairs (Y(0),Ẏ(0)) were first simulated.(b)
Subsequently, an approximation for the joint pdffY(0),Ẏ(0) was statistically determined.(c)
Finally, an estimate of the upcrossing intensity is obtained by numerically integrating Rice's formula in  (1.1).
Figure6 illustrates the comparison of the level crossing estimates obtained using the proposed hybrid method, when N=1000, and those obtained using Monte Carlo simulations. The three dashed lines are independent estimates of μNs(u), and we observe that the variability between them is small, confirming the assumption that assuming N=1000 leads to estimates that are reasonably free from statistical fluctuations. The irregular solid line is obtained from Monte Carlo simulations and a fairly good agreement between the crossing intensities is observed. Though the required computation time in the Monte Carlo method is of the same order as in the proposed hybrid method, it is clear from Figure 6 that the estimates from the proposed method are more accurate for higher levels.
We next focus on examining the errors induced in estimating upcrossing intensities for high levels when the non-Gaussian features of the response processes are neglected. Consequently, the upcrossing intensity of the response with Gaussian loading, namely,(5.16)YG(t)=∫-TTq(t-s)dB(s)+∬-TTQ(t-s1,t-s2)dB(s1)dB(s2),
is also computed. Note that for the kernel q(·), the variance of the linear response remains unchanged, that is, is equal to one, while skewness and kurtosis are, respectively, zero and 3. The corresponding crossing intensities are computed using the same algorithm as for the proposed hybrid method, but for N=1, the response process is unconditionally Gaussian and no simulation of gamma processes is required. The results are illustrated in the same plot, see Figure 6, as the thicker solid line. For completeness, the corresponding level crossing intensities were also computed using the Monte Carlo technique used in this example. These estimates are shown in Figure 6 as the irregular thick line. Based on these observations, one can conclude the following.
(i) One can see that the extremal responses forYG(t) are much smaller than the ones under LMA forcing, even though in both cases mean and variance are equal. For example, if one assumed that the two forcing are stationary and last for 100 years, then the 100-year-response, defined as the level crossing intensity approximately equal to 3×10-10, can be examined from Figure 6. We observe that while for the Gaussian forcing the level is approximately 10, the corresponding level for the skewed non-Gaussian loading is 23, a difference of more than 100%. It is quite obvious that neglecting the non-Gaussian features of the response leads to an underestimation of the level crossing intensities. This highlights the importance of modeling the non-Gaussian features of the response, especially in the context of risk analysis against high levels (rare events).
(ii) The close agreement between the level crossing estimates for the responseYG(t) using the saddle-point method (whose performance has already been examined in detail in other studies) and the Monte Carlo simulation approach used in this example provides confidence on the accuracy of the level crossing estimates obtained using the proposed MC approach.
Finally, one may ask about the accuracy of the estimates,μNs(u), computed for smaller number N of simulated γi processes. In order to answer this question, the crossing intensities were estimated using the proposed hybrid method with N=100 gamma process simulations. Thirty independent estimates of μNs(u) were calculated and are represented as thin solid lines in Figure 7. From the figure, one can see that the variability of μ100s(u) is quite large indicating that N=100 is probably too small a sample size for the statistical fluctuations to die down.Figure 6
Example5.4: The dashed lines (3 in number) indicate the crossing intensities μNs(u) using the proposed hybrid method, with sample size for simulated gamma processes γi being N=1000; the corresponding irregular solid line is MC estimate; the thicker dashed dotted line is the saddle-point estimate μs(u) with Gaussian forcing; the corresponding irregular solid line is the corresponding MC estimate.Figure 7
Example5.4: The thin solid lines (30 in number) are the crossing intensities μNs(u) using the proposed hybrid method, with sample size for simulated gamma processes γi being N=100; the corresponding irregular solid line is MC estimate; the thicker dashed dotted line is the saddle-point estimate μs(u) with Gaussian forcing; the corresponding irregular solid line is the corresponding MC estimate.


## 5.1. Saddle-Point Approximation of Crossing Intensity for LMA Processes

We consider the crossings of a linear response process, given by(5.1)Y(t)=∫q(t-x)dΛ(x),
with symmetrical kernel q(·). The corresponding moment generating function is given by [20](5.2)M(s,t)=exp(ζ∫-∞∞sq(x)+tq̇(x)dx)·exp(-1ν∫-∞∞log(1-μ(sq(x)+tq̇(x))-σ22(sq(x)+tq̇(x))2)dx).
Since M(s,t)=M(s,-t), one can use the efficient algorithm of the saddle-point method discussed in [10, 11, 14].In order to simplify the presentation we introduce the following notations;μNs(u) is the estimate of μ(u)=E[N(u)] computed by means of the hybrid saddle-point method and (4.6)-(4.7), and μs(u) denotes the estimate of E[N(u)] by means of saddle-point method and M(s,t) defined in (5.2). Here N(u) is defined as the observed number of upcrossings of level u divided by the length of the “observation” time. In all the examples, N(u) has units Hz.We first focus on the computation of the conditional moment generated functionM(s,t∣γ). Let us consider two LMA processes defined by a common Laplace motion. More precisely, for two kernels f1(·),f2(·) and the Laplace motion Λ(x), define(5.3)X1(t)=∫f1(t-x)dΛ(x),X2(t)=∫f2(t-x)dΛ(x).Here,Λ(x) is defined as in (2.8). Now, conditionally that Γ(·)=γ(·), the Laplace motion can be written as(5.4)λ(x)=ζx+μγ(x)+σB(γ(x))
and hence the conditional LMA processes, X1(t) and X2(t), can be represented as(5.5)X1(t)=∫f1(t-x)dλ(x),X2(t)=∫f2(t-x)dλ(x),
respectively. Obviously for any t, (here we take t=0), the joint pdf of X1(0) and X2(0) is Gaussian with means and covariances mi, σij, i,j=1,2, given by(5.6)mi=ζ∫fi(x)dx+μ∫fi(x)dγ(x),σij=σ2∫fi(x)fj(x)dγ(x).
Using (5.6), with f1(x)=q(x) and f2(x)=q̇(x), leads to(5.7)M(s,t∣γ)=E[esY(0)+tẎ(0)∣Γ(·)=γ(·)]=exp(sm1+tm2+0.5s2σ11+0.5t2σ22+stσ12).Example 5.1.
In this example,30 minutes of measured stress in a ship under stationary severe sea conditions is modeled as an LMA process. A part of the stress is shown in Figure 1(a). One can clearly see the existence of high-frequency oscillations, likely due to whippings, which get superimposed with the wave-induced stress. Figure 1(b) illustrates an estimated spectrum, S(ω), having two peaks. The kernel q(x) is computed from the spectrum S(ω) by inverting (2.5). The choice of kernels for LMA processes is quite important. Though there can be many kernels which give the same spectrum, there can be only one symmetrical kernel. In this paper, we consider only symmetrical kernels. Thus, this unique symmetrical kernel is obtained by finding the inverse of (2.5). Such spectra give time reversible loads. However, these loads do not always match with observed loadings. More work is needed to determine the kernel from observed data. The symmetrical kernel obtained by imposing these conditions is shown in Figure 2.Example5.1: (a) A part of the measured stress. (b) An estimate of the spectral density for the measured stress.
(a)(b)Figure 2
Example5.1: Estimated kernel function q(x).We next need to estimate the parameters of the LMA process. The variance of the stress time history is obtained by integrating the spectrum,S(ω), with respect to ω. Additionally, we assume that stress time history is mean zero. In order to identify the remaining parameters of the LMA process, we compute the skewness and excess kurtosis which are 0.13 and 0.21, respectively. These values indicate that the stress process is slightly non-Gaussian.Figure3 illustrates the crossing intensity N(u) for the measured stress (solid irregular line). In prediction of extremes, the crossing intensity needs to be extrapolated to much higher levels. Here, the LMA model is used for the extrapolation. The crossings of LMA are estimated by means of μs(u), that is, the saddle-point method, where the moment generated function, M(s,t) has been defined in (5.2). The function μs(u) is shown in the plot as a dashed dotted line. The agreement between N(u) and μs(u) is seen to be very good, except at the highest observed values of N(u). These discrepancies can be attributed to extremely large whipping effects, which consist of several crossings of high levels. This effect is averaged in μs(u).Figure 3
Example5.1: Observed crossing intensity N(u) in the measured stress-solid line; simulated crossing intensity in 100 times longer signal than measured (50 hours)—solid line with dots; the saddle-point approximation μs(u) of E[N(u)]—dashed dotted line.In order to verify this claim, we simulated the LMA process for a much longer duration (50-hour period) and computed the crossing intensities. The resulting crossing intensity, N(u), is superimposed in Figure 3 by solid line with dots. We observe that the estimated crossing intensities follows closely those computed using the saddle-point method. This confirms the accuracy of the saddle-point method.The primary objectives of this example are(a)
to study the applicability of the approximationμNs(u) (computed by means of the saddle-point method and formulas (4.7)-(4.6)) and to predict the return values, that is, levels uT such that E[N+(uT)]=1/T, and(b)
to examine how fastμNs(u) converges to E[N+(uT)], which here is estimated by μs(u).These are slightly different problems since in (a), one is interested in the horizontal distance betweenμNs(u) and μs(u), when plotted against levels u, while in (b), one examines the vertical distance between the functions. The conclusions of these studies are illustrated by means of Figures 4(a) and 4(b). In Figure 4(a), we observe that even for as low N=102, one gets relatively small errors (about 10%) in predictions of uT. However, the vertical convergence is slower and one needs about N=104 simulations of γi to get satisfactory distance between the two lines; see Figure 4(b), where the fractions μNs(u)/μ∞s(u) for N=102, 103 and 104, are presented. The algorithm is relatively fast and one can use high values of N to obtain satisfactory accuracy levels. We further note that in Figure 4(a), the computed crossing intensities towards the right extremes are better with 102 samples than with 103 samples. This apparent contradiction can be explained by the fact that the crossing intensities determined in the procedure above are statistical estimates. The standard deviations of the estimates obtained by 102, 103 and 104 samples are, respectively, c/102, c/103, and c/104, where c is higher for higher values of u. Thus, for a fixed value of c, the standard deviations of the estimates with 103 samples and 104 samples are, respectively, 0.32 and 0.1 times the standard deviation of the estimate with 102 samples. Figure 4(a) illustrates the computed intensities for only one set of N values. If the exercise was repeated, we expect to see a much smaller spread in the computed crossing intensities for higher values of N. Alternatively, a confidence band could be computed by means of the parametric bootstrap.Example5.1: (a) Crossing intensities μNs(s,t) computed using the proposed hybrid method: Sample size for simulated gamma processes γi; N=102: dotted line; N=103: dashed dotted line; N=104: dashed line; Crossing intensity μs(u): solid line. (b) Corresponding relative errors μNs(u)/μs(u).
(a)(b)


## 5.2. Computation ofM(s,t) for the Quadratic Response

The general quadratic response is only notationaly more complex and we will proceed in a similar way as for the LMA process discussed in Example5.1. First, we need to find the conditional moment generating function(5.8)M(s,t∣γ)=E[esY(0)+tẎ(0)∣Γ(·)=γ(·)],
which can be written by an explicit formula; see (5.11) derived below. Then one can simulate a sequence of gamma processes, γi(·), i=1,…,N and as before approximate M(s,t) by means of (4.7).LetΛ be a diagonal matrix with the diagonal elements being denoted by λi, i=1,…,n, and the rest of the elements being zero. Using matrix notation, the response process can be written as(5.9)Y(t)=aW(t)T+12W(t)ΛW(t)T=∑j=1najWj(t)+12∑j=1nλjWj2(t),
where a=(a1,…,an) (an=1,λn=0). As discussed earlier, conditionally on Γ(·)=γ(·), the vectors W=W(0) and Ẇ=Ẇ(0) are normally distributed with means m and ṁ and covariance matrices Σ11, Σ12, and Σ22, where for 1≤i,j≤n, (5.10)σ11(i,j)=σ2∫ϕi(x)ϕj(x)dγ(x),σ12(i,j)=σ2∫ϕi(x)ϕ̇j(x)dγ(x),σ22(i,j)=σ2∫ϕ̇i(x)ϕ̇j(x)dγ(x),m(i)=ζ∫ϕi(x)dx+μ∫ϕi(x)dγ(x),ṁ(i)=ζ∫ϕ̇i(x)dx+μ∫ϕ̇i(x)dγ(x).
Once the matrices Σij and vectors m and ṁ are computed, it is a straightforward task to compute M(s,t∣γ), see [10], which is given by(5.11)M(s,t∣γ)=1det(Σ)exp(ms+ṁt+12t2m̃Vm̃+12tTΣ-1t).
Here, (5.12)Σ=Σ11-1-sΛ-t(ΛΣ21Σ11-1+Σ11-1Σ12Λ),t=sm̃+tṁΛ+tm̃Σ11-1Σ12+t2m̃VΛ,m̃=a+mΛ,m=amT+12mΛmT,ṁ=aṁT+ṁΛmT.Remark 5.2.
Example5.1 is obtained as a special case when n=1, with ϕ1(s)=q(s) and λi=0 while a1=1 in (5.9). Under these conditions, using simple algebraic manipulations, it can be shown that the conditional moment generated function is equal to the expression in (5.7).Example 5.3.
In this example, we focus on checking the accuracy of the estimates of the level crossing intensity,μNs(u), using the proposed hybrid method, for quadratic response Y(t) in (5.9) for the special case when n=2 and ϕn=0, that is,
(5.13)Y(t)=Y1(t)+λY12(t)2=λ(Y1(t)+1/λ)22-12λ.
Considering the case n=2 provides certain advantages which can be exploited to benchmark the accuracy of the estimates, μNs(u), using the proposed method. Using (4.3) and (5.13), it can be shown that the crossing intensity μY(u)=E[NY(u)] can be expressed as
(5.14)μY(u)=μY1(-1λ+2uλ+1λ2)+μY1(-1λ-2uλ+1λ2).
As can be seen from (5.14), the accuracy of the estimate μY(u) depends on the estimate of the crossing intensity μY1. This, however, poses no problem as this can be very accurately obtained using the direct saddle-point method. Thus, replacing μY1 in (5.14) by the saddle-point estimate, μY1s, the expression in (5.14) can be used to benchmark the accuracy of the level crossing estimate, μNs(u), obtained using the proposed hybrid method.
As in Example5.1, Y1(t) is a stress time history of duration of 30 minutes measured in a particular location of a ship impinged by ocean waves during the course of its journey; see Figure 1(a). We use the LMA process described in Example 5.1 to model Y1(t). For the quadratic response, we choose λ=0.01. This value is chosen so that the contribution of linear and quadratic parts to Y1 is similar; note that standard deviation of Y1(t) is about 47 MPa. An estimate of the crossing intensity μY(u) is obtained using  (5.14) and is shown in Figure 5(a). The accuracy of the crossing intensities for the corresponding levels, μNs(u), obtained using the proposed method are determined by comparing with these values.
In order to computeμNs(u), one needs the expression for the conditional moment generating function M(s,t∣γ). This is given in (5.11)-(5.12), with Σ11=σ11, a=1, Σ22=σ22, Σ12=σ12. All parameters have the same values as in Example 5.1. A comparison of the crossing intensity estimates, μNs, using the proposed hybrid method is illustrated in Figure 5(a). As in Example 5.1, we consider the three cases where N=102, 103, and 104, where N is the number of gamma process simulations in the proposed hybrid method. A comparison of the relative errors is shown in Figure 5(b). As in Example 5.1, we observe that the estimates are in fairly good agreement with the accuracy expectedly improving for larger values of N.Example5.3: (a) Crossing intensities μNs(u) computed using the proposed hybrid method: sample size for simulated gamma processes γi; N=102: dotted line; N=103: dashed dotted line; N=104: dashed line; μs(u): solid line. (b) Corresponding relative errors μNs(u)/μs(u).
(a)(b)Example 5.4.
In this example, we consider a more general quadratic response process, such that the number of termsn in  (5.9) is more than one. We consider the response process Y(t)=Y1(t)+Y2(t) defined in (3.1), where q(s)=exp(-s2/50)/25π, -25≤s≤25 and
(5.15)Q(t,s)=0.01exp(-(s-t)250).
The parameters in Laplace motion, Λ(x), are chosen in such a way that the linear response, Y1(t)=∫-TTq(t-s)dΛ(s), has mean zero, variance one, skewness 0.5, and kurtosis 4.5. For the kernel Q(t,s), the first 12 eigenvalues were found to be significantly nonzero. To determine the number of such eigenvalues, the first 100 eigenvalues were found and ordered according to their absolute values, and their corresponding ratios with respect to their total summation were calculated. It was assumed that the series could be truncated when the sum of the absolute value of the eigenvalues exceeded 99.9% of the total sum. This led to n=12 for this example.
Based on experience from Examples5.1 and 5.3, we expect that N=1000 simulations of γi are needed for arriving at a reasonably accurate estimate of μNs(u) using the proposed hybrid method. In the absence of any closed-form analytical solutions for the crossing intensities of the quadratic response process, we compare the estimates obtained using the proposed hybrid method with those obtained from Monte Carlo simulations. For Monte Carlo simulations, simulating a large number of response processes and checking for their crossing intensities would be computationally very expensive and time consuming. Instead, we adopt the following MC procedure. (a)
1×107 independent samples of pairs (Y(0),Ẏ(0)) were first simulated.(b)
Subsequently, an approximation for the joint pdffY(0),Ẏ(0) was statistically determined.(c)
Finally, an estimate of the upcrossing intensity is obtained by numerically integrating Rice's formula in  (1.1).
Figure6 illustrates the comparison of the level crossing estimates obtained using the proposed hybrid method, when N=1000, and those obtained using Monte Carlo simulations. The three dashed lines are independent estimates of μNs(u), and we observe that the variability between them is small, confirming the assumption that assuming N=1000 leads to estimates that are reasonably free from statistical fluctuations. The irregular solid line is obtained from Monte Carlo simulations and a fairly good agreement between the crossing intensities is observed. Though the required computation time in the Monte Carlo method is of the same order as in the proposed hybrid method, it is clear from Figure 6 that the estimates from the proposed method are more accurate for higher levels.
We next focus on examining the errors induced in estimating upcrossing intensities for high levels when the non-Gaussian features of the response processes are neglected. Consequently, the upcrossing intensity of the response with Gaussian loading, namely,(5.16)YG(t)=∫-TTq(t-s)dB(s)+∬-TTQ(t-s1,t-s2)dB(s1)dB(s2),
is also computed. Note that for the kernel q(·), the variance of the linear response remains unchanged, that is, is equal to one, while skewness and kurtosis are, respectively, zero and 3. The corresponding crossing intensities are computed using the same algorithm as for the proposed hybrid method, but for N=1, the response process is unconditionally Gaussian and no simulation of gamma processes is required. The results are illustrated in the same plot, see Figure 6, as the thicker solid line. For completeness, the corresponding level crossing intensities were also computed using the Monte Carlo technique used in this example. These estimates are shown in Figure 6 as the irregular thick line. Based on these observations, one can conclude the following.
(i) One can see that the extremal responses forYG(t) are much smaller than the ones under LMA forcing, even though in both cases mean and variance are equal. For example, if one assumed that the two forcing are stationary and last for 100 years, then the 100-year-response, defined as the level crossing intensity approximately equal to 3×10-10, can be examined from Figure 6. We observe that while for the Gaussian forcing the level is approximately 10, the corresponding level for the skewed non-Gaussian loading is 23, a difference of more than 100%. It is quite obvious that neglecting the non-Gaussian features of the response leads to an underestimation of the level crossing intensities. This highlights the importance of modeling the non-Gaussian features of the response, especially in the context of risk analysis against high levels (rare events).
(ii) The close agreement between the level crossing estimates for the responseYG(t) using the saddle-point method (whose performance has already been examined in detail in other studies) and the Monte Carlo simulation approach used in this example provides confidence on the accuracy of the level crossing estimates obtained using the proposed MC approach.
Finally, one may ask about the accuracy of the estimates,μNs(u), computed for smaller number N of simulated γi processes. In order to answer this question, the crossing intensities were estimated using the proposed hybrid method with N=100 gamma process simulations. Thirty independent estimates of μNs(u) were calculated and are represented as thin solid lines in Figure 7. From the figure, one can see that the variability of μ100s(u) is quite large indicating that N=100 is probably too small a sample size for the statistical fluctuations to die down.Figure 6
Example5.4: The dashed lines (3 in number) indicate the crossing intensities μNs(u) using the proposed hybrid method, with sample size for simulated gamma processes γi being N=1000; the corresponding irregular solid line is MC estimate; the thicker dashed dotted line is the saddle-point estimate μs(u) with Gaussian forcing; the corresponding irregular solid line is the corresponding MC estimate.Figure 7
Example5.4: The thin solid lines (30 in number) are the crossing intensities μNs(u) using the proposed hybrid method, with sample size for simulated gamma processes γi being N=100; the corresponding irregular solid line is MC estimate; the thicker dashed dotted line is the saddle-point estimate μs(u) with Gaussian forcing; the corresponding irregular solid line is the corresponding MC estimate.


## 6. Concluding Remarks

The problem of estimating the crossing intensities of the response process of second-order dynamical systems, subjected to non-Gaussian loadings, has been studied. The loads are assumed to be strictly stationary and are modeled as LMA processes. This enables retaining the non-Gaussian features, such as skewness and kurtosis, of the marginal distributions. For second-order dynamical systems, the response is expressed as a quadratic combination of the LMA processes is non-Gaussian. Direct application of Rice's formula is not possible as the joint pdf of the response and its instantaneous time derivative are not available. A numerical method is developed so that approximations for the crossing intensities can be computed with fairly reasonable accuracy. Three numerical examples have been presented to illustrate the proposed method. The salient features emerging from this study are as follows.(1)
The proposed method is a hybrid method that combines the analytical saddle-point approximation and the Monte Carlo approach. Consequently, the proposed method is much faster than Monte Carlo simulations.(2)
The accuracy levels of the proposed hybrid method depend on the number of samples of Gamma process simulations and are expectedly better for larger sample size. For the examples considered in this paper, a sample size of1000 is found to lead to estimates of fairly good accuracies.(3)
Neglecting the non-Gaussian effects of the loading can severely underestimate the crossing intensities of the response, particularly for high levels. This, in turn, implies overestimating the safety and reliability of a system subjected to rare loadings, leading to unsafe designs.(4)
The proposed method is applicable for systems with symmetric second-order kernels. Fortunately, most physical second-order dynamical systems ensure symmetric second-order kernels. Therefore, this is not a severe restriction.



---

*Source: 752452-2010-04-29.xml*