#### 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.

#### 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, of level , is given by Rice's formula [1, 2], expressed as where is the joint probability density function (j-pdf) of the response and its instantaneous time derivative . The applicability of (1.1) lies in the availability of the information on the j-pdf . This is, however, rarely available.

Exact information about the j-pdf, , 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 where and and indicate the variance and the expected value, respectively.

The probability distribution of the extreme response in a fixed period , namely, , can be conservatively estimated by means of the inequality See, for example, [3, 4]. For stationary Gaussian responses the stronger result that as 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 force , the response can be written as a sum where

In (1.5), and , respectively, denote the linear and quadratic transfer functions of the system. Here, it can be assumed that is a smooth Gaussian process, given by where is a Brownian motion while 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 as in (1.4). Statistical analysis of extremes of 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, , and a continuous and increasing function . Then one forms a non-Gaussian process, , according to the transformation . In this way, the process can have a non-Gaussian marginal distribution. Different strategies to choose the function 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 , as defined in (1.4), with assumed to be an LMA process.

The paper is organized as follows. First, in Section 2, 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 , given by , 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
where is a kernel function and is a stochastic process with independent and stationary increments having a generalized asymmetric Laplace distribution. The process is referred to as *Laplace motion* and the resulting process is called the *Laplace driven moving average* (LMA). Thus may be thought of as a convolution of with the increments of the process . A process generated in this way is stationary and ergodic. In the special case where is chosen to be a Brownian motion, 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 variable is said to have a generalized asymmetric Laplace distribution if its characteristic function is given by Here, , and , are parameters of the Laplace distribution and . If , 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 is a process that starts at zero and whose distribution at is given by 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 (), 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: where 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 are given by Here, denotes the Fourier transform. It must be noted that (2.5) is valid for any square integrable kernel. As symmetrical kernels have real Fourier transform , this implies that if 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 , we get However, after having chosen the kernel 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 and excess kurtosis (if ) of the marginal distribution of . By using the expression for the characteristic function in (2.4), these are given by

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, . In fact, a Gaussian process can be obtained from the Laplace driven MA as a limiting case as and , for example, by letting and in such a way that 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 kernel . In full generality, following [21], the asymmetric Laplace motion , with drift , asymmetry parameter , and variance can be represented as Here, is a gamma process characterized by independent and homogeneous -increments having a gamma distribution with shape parameter and scale parameter while is Brownian motion. Using this representation, a simple algorithm for simulating the Laplace driven moving average with kernel is given by the following. Pick and so that is well approximated by its values onPick in order to simulate the values of the response process at the points . Simulate identically and independently distributed (i.i.d.) random variables and store them in a vector . Simulate i.i.d. zero mean standard normal random variables and store them in a vector . Compute , where , denotes convolution and the integral is computed by some numerical method.

The above simulation strategy is based on the fact that conditional on the Gamma process increments, the increment in are independent and Gaussian, where 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, , were studied for stationary process , as defined in (1.4), when is a stationary Gaussian process. Here, we consider the more general case where is modeled as an LMA process; see (2.1). Combining (1.4) and (2.1), the response process can be rewritten as Here,

For most real life engineering applications, the kernel is symmetrical. Further, we assume that the kernels are square integrable and hence vanish at infinity. Thus, by choosing sufficiently large, we may approximate the kernels by letting and for and . Under such assumptions, the Kac-Siegert technique based on the representation of the truncated kernel through its eigenfunctions and eigenvalues , can be employed. Let the eigenfunctions and eigenvalues of the kernel be defined by

For a symmetrical kernel , the eigenfunctions corresponding to the different eigenvalues are orthogonal. By further normalization, we assume that are orthonormal with eigenvalues . Suppose that the eigenfunctions are ordered according to . Both eigenvalues and eigenfunctions are real, as , and see [22]. Further, for simplicity of presentation, we assume that and hence can be expanded in a series using the orthonormal eigenfunctions , namely, From (3.6), it follows that for any where is computed from the condition that . Let us introduce In quadratic mean, the response in (1.4) can be rewritten as where are LMA processes.

Often only a few of the eigenvalues are significantly nonzero. Assuming that the number of such eigenvalues is , (3.9) can be rewritten as For notational convenience, we drop the tilde and rewrite (3.11) as where . The instantaneous time derivative process is given by Note that (3.12) and (3.13) are functions of the vectors of LMA processes and . A procedure for estimating the upcrossing intensity for in (3.12) is discussed in the following section.

#### 4. Computing the Upcrossing Intensity

The upcrossing intensity of can be computed using (1.1) if the j-pdf of and is available. This, however, is not easy when is defined as in (3.12). The elements in vectors and all have generalized Laplace distributions whose marginal pdfs are usually defined through their characteristic functions. Also, since the elemental processes and have mutual dependence, the computation of the joint characteristic function of and is a difficult task. In the special case when is an LMA-process, that is, see (1.4), and when 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 for the general quadratic response process, , as defined in (1.4), but with LMA forcing. It must be noted that for general quadratic processes, , 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 or to estimate the joint density of needed to compute using (1.1). However, the MC approach is not an efficient way for computing for high levels , 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, and are normally distributed. Consequently, the computation of conditional moment generating function is straightforward and is given by see Section 5.2. Now, the upcrossing intensity can be expressed as the expectation of , that is, the number of upcrossings of level by the process in duration . Thus, one can find the conditional upcrossing intensity of the process 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 The upcrossing intensity can be estimated by computing the conditional moment generating function in (4.1) and using the saddle-point method to estimate . 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, , is symmetrical in , that is, . Note that the numerical algorithm presented in [10, 11, 14] is restricted to the symmetrical case. Unfortunately, the conditional moment generating function in (4.1) is not, in general, symmetrical. For the asymmetrical , 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 and in (3.2) are symmetrical, which is what has been assumed in this paper.

##### 4.1. Approximation of the Upcrossing Intensity

Assuming that is a stationary process, and have the same expected number of upcrossings of any level . Consequently, where is independent of the process and takes values or with probability . Additionally, has the same upcrossing intensity as the process . In the special case when is given by (1.4) with LMA forcing, the upcrossing intensity can be expressed as

Let the conditional crossing intensity be defined as Then, by simulating a sequence of Gamma processes, , , the unconditional crossing intensity, , can be estimated by averaging , namely, where 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 . For time reversible in (4.3), since the conditional moment generating function can be expressed as it is obvious that is symmetrical. This enables one to use the saddle-point algorithms discussed in [10, 11, 14] to estimate the conditional upcrossing intensity .

Clearly, the method to estimate the upcrossing intensity 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 -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, when 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 in (4.6) should be in order to reach desired accuracy.

Next, we study the crossings of a simple quadratic response . The upcrossing intensity can be computed when upcrossing intensity of the linear response 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 of of full complexity and estimate the upcrossing intensity. Here, eigenvalues 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 with symmetrical kernel . The corresponding moment generating function is given by [20] Since , 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; is the estimate of computed by means of the hybrid saddle-point method and (4.6)-(4.7), and denotes the estimate of by means of saddle-point method and defined in (5.2). Here is defined as the observed number of upcrossings of level divided by the length of the “observation” time. In all the examples, has units Hz.

We first focus on the computation of the conditional moment generated function . Let us consider two LMA processes defined by a common Laplace motion. More precisely, for two kernels and the Laplace motion , define

Here, is defined as in (2.8). Now, conditionally that , the Laplace motion can be written as and hence the conditional LMA processes, and , can be represented as respectively. Obviously for any , (here we take ), the joint pdf of and is Gaussian with means and covariances , , , given by Using (5.6), with and , leads to

*Example 5.1. *In this example, 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, , having two peaks. The kernel is computed from the spectrum 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.

**(a)**

**(b)**

We next need to estimate the parameters of the LMA process. The variance of the stress time history is obtained by integrating the spectrum, , 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 and , respectively. These values indicate that the stress process is slightly non-Gaussian.

Figure 3 illustrates the crossing intensity 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 , that is, the saddle-point method, where the moment generated function, has been defined in (5.2). The function is shown in the plot as a dashed dotted line. The agreement between and is seen to be very good, except at the highest observed values of . These discrepancies can be attributed to extremely large whipping effects, which consist of several crossings of high levels. This effect is averaged in .

In order to verify this claim, we simulated the LMA process for a much longer duration (-hour period) and computed the crossing intensities. The resulting crossing intensity, , 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 (computed by means of the saddle-point method and formulas (4.7)-(4.6)) and to predict the return values, that is, levels such that , and (b)to examine how fast converges to , which here is estimated by .

These are slightly different problems since in (a), one is interested in the horizontal distance between and , when plotted against levels , 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 , one gets relatively small errors (about %) in predictions of . However, the vertical convergence is slower and one needs about simulations of to get satisfactory distance between the two lines; see Figure 4(b), where the fractions for , and , are presented. The algorithm is relatively fast and one can use high values of to obtain satisfactory accuracy levels. We further note that in Figure 4(a), the computed crossing intensities towards the right extremes are better with samples than with 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 , and samples are, respectively, , , and , where is higher for higher values of . Thus, for a fixed value of , the standard deviations of the estimates with samples and samples are, respectively, and times the standard deviation of the estimate with samples. Figure 4(a) illustrates the computed intensities for only one set of values. If the exercise was repeated, we expect to see a much smaller spread in the computed crossing intensities for higher values of . Alternatively, a confidence band could be computed by means of the parametric bootstrap.

**(a)**

**(b)**

##### 5.2. Computation of 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 Example 5.1. First, we need to find the conditional moment generating function which can be written by an explicit formula; see (5.11) derived below. Then one can simulate a sequence of gamma processes, , and as before approximate by means of (4.7).

Let be a diagonal matrix with the diagonal elements being denoted by , , and the rest of the elements being zero. Using matrix notation, the response process can be written as where (). As discussed earlier, conditionally on , the vectors and are normally distributed with means and and covariance matrices , , and , where for , Once the matrices and vectors and are computed, it is a straightforward task to compute , see [10], which is given by Here,

*Remark 5.2. *Example 5.1 is obtained as a special case when , with and while 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, , using the proposed hybrid method, for quadratic response in (5.9) for the special case when and , that is,
Considering the case provides certain advantages which can be exploited to benchmark the accuracy of the estimates, , using the proposed method. Using (4.3) and (5.13), it can be shown that the crossing intensity can be expressed as
As can be seen from (5.14), the accuracy of the estimate depends on the estimate of the crossing intensity . This, however, poses no problem as this can be very accurately obtained using the direct saddle-point method. Thus, replacing in (5.14) by the saddle-point estimate, , the expression in (5.14) can be used to benchmark the accuracy of the level crossing estimate, , obtained using the proposed hybrid method.

As in Example 5.1, is a stress time history of duration of 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 . For the quadratic response, we choose . This value is chosen so that the contribution of linear and quadratic parts to is similar; note that standard deviation of is about MPa. An estimate of the crossing intensity is obtained using (5.14) and is shown in Figure 5(a). The accuracy of the crossing intensities for the corresponding levels, , obtained using the proposed method are determined by comparing with these values.

In order to compute , one needs the expression for the conditional moment generating function . This is given in (5.11)-(5.12), with , , , . All parameters have the same values as in Example 5.1. A comparison of the crossing intensity estimates, , using the proposed hybrid method is illustrated in Figure 5(a). As in Example 5.1, we consider the three cases where , , and , where 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 .

**(a)**

**(b)**

*Example 5.4. *In this example, we consider a more general quadratic response process, such that the number of terms in (5.9) is more than one. We consider the response process defined in (3.1), where , and
The parameters in Laplace motion, , are chosen in such a way that the linear response, , has mean zero, variance one, skewness , and kurtosis . For the kernel , the first eigenvalues were found to be significantly nonzero. To determine the number of such eigenvalues, the first 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 % of the total sum. This led to for this example.

Based on experience from Examples 5.1 and 5.3, we expect that simulations of are needed for arriving at a reasonably accurate estimate of 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) independent samples of pairs were first simulated. (b)Subsequently, an approximation for the joint pdf was statistically determined. (c)Finally, an estimate of the upcrossing intensity is obtained by numerically integrating Rice's formula in (1.1).

Figure 6 illustrates the comparison of the level crossing estimates obtained using the proposed hybrid method, when , and those obtained using Monte Carlo simulations. The three dashed lines are independent estimates of , and we observe that the variability between them is small, confirming the assumption that assuming 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,
is also computed. Note that for the kernel , the variance of the linear response remains unchanged, that is, is equal to one, while skewness and kurtosis are, respectively, zero and . The corresponding crossing intensities are computed using the same algorithm as for the proposed hybrid method, but for , 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 for 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 years, then the -year-response, defined as the level crossing intensity approximately equal to , can be examined from Figure 6. We observe that while for the Gaussian forcing the level is approximately , the corresponding level for the skewed non-Gaussian loading is , 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 response 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, , computed for smaller number of simulated processes. In order to answer this question, the crossing intensities were estimated using the proposed hybrid method with gamma process simulations. Thirty independent estimates of were calculated and are represented as thin solid lines in Figure 7. From the figure, one can see that the variability of is quite large indicating that is probably too small a sample size for the statistical fluctuations to die down.

#### 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 of 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.

#### Acknowledgments

The research presented in this paper has been partially supported by the Gothenburg Stochastic Center and the Swedish foundation for Strategic Research through GMMC, Gothenburg Mathematical Modelling Center. The authors also express their gratefulness to DNV crew, management company, and owner for providing measurement data. The first author acknowledges the support of the European Union project SEAMOCS.