The problem of the analysis of nonlinear time series models has attracted a great deal of interest in finance, business, physics, and other sciences. Granger and Anderson (1978) and Tong and Lim (1980) noticed, in many time series modeled by Box and Jenkins (1970) , that the squared residuals are significantly autocorrelated even though the residual autocorrelations are not. This indicates that the innovation of these models might be uncorrelated but not independent. The authors suggested using the autocorrelation function of the squared values of the series in order to detect the nonlinearity. In this respect, Engle (1982) showed that the classical portmanteau tests proposed by Box and Pierce (1970) and Ljung and Box (1978) , based on the autocorrelation function of the residuals, fail to detect the presence of the Autoregressive Conditional Heteroscedasticity, arch , in many financial time series models. He introduced a Lagrange multiplier statistic, based on the autocorrelations of the squared-residuals, to test for the presence of arch process. McLeod and Li (1983) proposed a portmanteau test, which was also based on the squared-residuals autocorrelations, to detect nonlinearity and arch effect. Their statistic is similar in spirit to the modified test of Box and Pierce (1970) proposed by Ljung and Box (1978) .
Since then, several authors have developed goodness-of-fit tests based on the autocorrelation function of the squared-residuals to detect nonlinear structures in time series models (e.g., Li and Mak (1994) ; Rodríguez (2002, 2006) ; Rodríguez and Ruiz (2005) ; Fisher and Gallagher (2012) ). Simulation studies demonstrate the usefulness of using these statistics to test for arch effects but they tend to lack power compared to other types of nonlinear models that do no have such effects. A possible reason for such lack of power could be the fact that none of these tests have considered the cross correlation between the residuals at different powers where the correlations at positive and negative lags might be vary. Lewis (1985, 1987) introduced the idea of using the generalized correlations 1 of the residuals to detect nonlinearity without taking into account the distribution of the parameters estimation. Recently, Psaradakis and Vávra (2019) implemented the Lewis (1985, 1987) idea and they proposed four goodness-of-fit tests, based on the generalized correlations of the residuals of linear time series model, to test for linearity of stationary time series. Two of these tests are similar in spirit to Box and Pierce (1970) where they replaced the autocorrelations of the residuals in Box and Pierce (1970) by the cross-correlations of the residuals at different powers getting two tests: one is associated with positive lags; and the other is based on negative lags 2 . Similarly, the other two tests were proposed by replacing, respectively, the autocorrelations of the residuals in Ljung and Box (1978) (or square-residuals in McLeod and Li (1983) ) tests by the cross-correlations of the residuals at different powers. Their preliminary analysis indicated that the McLeod and Li (1983) modification tests control the Type I error probability somewhat more successfully than the other tests. Hence, they restricted their simulation study focusing on comparing these two tests (where r, s ∈ {1, 2}, r = s) with McLeod and Li (1983) test (where r = s = 2) and they suggested, in most cases, that at least one of the two cross-correlation tests tends to have more 1 Generalized correlation is defined as the autocorrelation between the residuals to the power r at time t and the residuals to the power s at time t + k where r, s are positive s and k is the lag time.
2 The values of autocorrelations of the residuals at positive lags are the same values at negative lags. On the other hand the values of the cross-correlations between the residuals to the power r at lag time t and the residuals to the power s at lag time t + k when k > 0 are not the same when k < 0.
power than the test based on the autocorrelation of the squared-residuals in several time series models.
To our knowledge, none of the tests in the time series literature has simultaneously combined the autocorrelations of the residuals and the autocorrelations of their square values with the cross-correlation between them. In this article, we fill this gap by proposing four goodness-of-fit tests. The first part in each of the proposed tests is based on the autocorrelations of the residuals, which is used to test for linearity (adequacy of fitted arma model); the second is based on the autocorrelation of the squared-residuals, which can be used to test for arch , and the third one is based on the cross correlations between the residuals and their squared values, which can be used to test for other types of nonlinear models in which the residuals and their values are crosscorrelated. Similar to Psaradakis and Vávra (2019) , the cross-correlations between the residuals and their squared values propose two different tests. One is based on positive lags and other is based on negative lags.
In Section 2 we discuss the generalized correlations of residuals and review some test statistics that have been commonly used to detect nonlinearity structure in stationary time series models.
In Section 3, we propose new goodness-of-fit (auto-and-cross-correlated) tests that can be used to detect, simultaneously, linear, bilinear, and nonlinear dependency in time series models, and derive their asymptotic distribution as a chi-squared distribution. In Section 4, we provide a Monte Carlo study comparing the performance of the proposed statistics with those from the literature. We show that the empirical significance level of the proposed tests successfully controls the Type I error probability and tends to have higher power than others. An illustrative application is given in Section 5 to demonstrate the usefulness of the proposed test for a real world data. We finish this article in Section 6 by providing a concluding remarks.
2 Auto-and-cross-correlated portmanteau tests Let z 1 , z 2 , · · · , z n denote a realization of real-valued stochastic process {z t }
where 
where the polynomials Φ p (B) and Θ q (B) are assumed to have all roots outside the unit circle and to have no common roots and η = (φ 1 , · · · , φ p , θ 1 , · · · , θ q ).
In time series literature, the notion of linearity is commonly used when the process {z t } admits the moving-average, ma , representation (1) where {ε t } are independent and identically distributed, i.i.d., random variables with a zero mean and a constant variance σ 2 . This is the notion considered by McLeod and Li (1983) ; Lewis (1985, 1987) and Psaradakis and Vávra (2019) which will also be adopted in this article. On the other hand, there are several time series that do not exhibit a linear behavior. For example, when the innovations {ε t } are uncorrelated but not independent, Engle (1982) proposed the Autoregressive Conditional Heteroscedasticity, arch , model that is widely used for analyzing financial time series. This has been generalized by Bollerslev (1986) to the generalized autoregressive conditional heteroscedasticity
where {ξ t } is a sequence of i.i.d. random variables with a mean value of 0 and a variance value
has been proposed to model the dynamic behavior of conditional heteroscedasticity in real time series. e.g.; the exponential garch (egarch ) model proposed by Nelson (1991) to allow for asymmetric effects between positive and negative financial time series. Tsay (2005) and Carmona (2014) provided nice reviews of the garch models.
Another popular nonlinear model is the threshold autoregressive (tar ) model of (Tong, 1978 (Tong, , 1983 (Tong, , 1990 Tsay, 1989) , which was generalized by Chan and Tong (1986) and Terasvirta (1994) to the smooth transition autoregressive (star ) model. The two regime-switching star model of order (2; p, p) takes the form
i ), i = 0, 1, · · · , p are the autoregressive coefficients, 0 ≤ F (.) ≤ 1 is a transition continuous function that allows the dynamics of model to switch between regimes smoothly, and s t is a transition variable. A common formulations for the transition function, which was proposed by Terasvirta (1994) , is the first-order logistic function and can represented as:
where γ > 0 denotes the smoothness parameter of the transition from one regime to the other, d ≥ 1 is the delay parameter, c is a threshold variable that separates the two-regimes, and σ z is the standard deviation of z t .
An alternative choice, which was also proposed by Terasvirta (1994) , is the exponential function given by
When the value of γ increases, the transition function F (z t ; γ, c) approaches the indicator function
In this case, the model in (4) reduces to the tar model of (Tong, 1978) . Arguably the tar models is consider to be the most popular class of nonlinear time series model. Therefore, testing for tar and star models have attracted much attention (Xia et al., 2020) .
Volterra series is another form of the nonlinear stationary process that is widely used in many applications (Wiener, 1953, lecture 10) . These models have the form
where µ is the mean level of z t and {ε t , −∞ < t < ∞} is a strictly stationary process of i.i.d. random variables.
Let β = (η, µ, σ 2 ) denote the true parameter values in (1),β = (η,μ,σ 2 ) denote the estimated values, andε i , i = 1, · · · , n denote the residuals and define the correlation coefficient at lag time k betweenε r t andε s t+k (r, s = 1, 2) aŝ
whereγ rs (k) = n −1 n−k t=1 f r (ε t )f s (ε t+k ) for k ≥ 0, γ rs (k) = γ sr (−k) for k < 0, is the autocovariance (cross-covariance), at lag time k, between the residuals to the power r and the residuals to the power s for r, s = 1, 2, and f j ( (1) of the squared residuals from the linear model given in (2) have been widely used in several portmanteau statistics to test for nonlinearity. The commonly employed test, which is similar in spirit to Box and Pierce (1970) , isQ
where 0 < m < n is the maximum lag considered for a significant autocorrelation 4 . This test is asymptotically distributed as chi-square with m degrees of freedom.
McLeod and Li (1983) showed that theQ 22 test can be improved if the autocorrelation coefficients in (8) are replaced with their standardized values r rs (k) = n + 2 n − |k|r rs (k), (r = s = 1, 2), k = ±1, · · · , ±m.
McLeod and Li (1983) proposed a portmanteau test, for detecting the presence of the arch effects.
Their test statistic is given by
where the Q 22 test is asymptotically distributed as chi-square with m degrees of freedom.
Simulation studies have showed that theQ 22 and Q 22 tests respond well to arch models but tend to lack power compared to other types of nonlinear models that do not have the arch effects.
3 The case for absolute residuals is beyond the scope of this article and we focus our attention to the squaredresiduals case. 4 The value of m depends on the sample size, and as a rule of thumb, we consider m ∈ {1, 2, · · · , √ n }, where n denotes the largest integer not exceeding n.
A reasonable justification for the lack of powers could be neglecting the cross-correlations between the residuals at different exponent powers (generalized correlations betweenε r t andε s t+k , where k ≥ 0 and r + s > 2). In this regard, Lewis (1985, 1987) proposed the idea of implementing the sample generalized correlations to detect the nonlinear dependency without considering the distribution of the parameter estimation. Recently, Psaradakis and Vávra (2019) used the idea of Lewis (1985, 1987) under the assumption in (1). They proposed portmanteau tests based on the generalized correlations and demonstrated the usefulness of using the test statistics based on the cross-correlation, between the residuals at different powers, for detecting linearity in time series. Their test statistics arê
and
whereQ rs and Q rs (r = s) are asymptotically distributed as χ 2 m . Psaradakis and Vávra (2019) suggested that the tests based on the cross-correlations tend to be more powerful against many types of nonlinearity compared to other statistics based on squared-residual autocorrelations.
Motivated by the ideas of Lewis (1985, 1987) , and Psaradakis and Vávra (2019), we propose in the next section new test statistics that can be used to detect nonlinearity in time series models.
In this section we made some assumptions that are needed to derive the asymptotic distribution of the proposed statistics. In general, the limiting distribution of the portmanteau test statistic required the following assumptions:
A2. The polynomial Ψ ∞ (B) is differentiable with respect to η in an open neighborhood of the closed disc |B| ≤ 1. This guarantees that the 1/Ψ ∞ (B) series is converge and the process {z t } admits the following Autoregressive, AR, representation with order ∞
where φ j , j = 0, 1, · · · which can be found by solving
. Hereβ can be estimated by the least squares method or the maximum (or quasi-maximum) likelihood method (see e.g., Hannan (1973) ; Hosoya and Taniguchi (1982) ; Kuersteiner (2001)).
A4. ∂γ rs (k)/∂β = O p (1/ √ n) for k = 1, 2, · · · , n − 1 and r, s ∈ N such that r + s ≥ 2 and
If the model in (2) is correctly identified and the assumptions A1-A4 are held, we propose the portmanteau test statisticsĈ
and
where the asymptotic distribution ofĈ rs and C rs is chi-square with 3m − (p − q) degrees of freedom, where r = s ∈ {1, 2}.
Remark. Each of the two test statistics in (16) can be seen as a linear combination of three existent tests, Ljung and Box (1978) , McLeod and Li (1983) , and Psaradakis and Vávra (2019), modifying the corresponding tests in (15).
Theorem 1. If {z t } satisfies (1) and the assumptions A1-A4 are held, then, for any integer m < n, the asymptotic distribution of √ n(r rr (1), · · · ,r rr (m),r rs (1), · · · ,r rs (m),r ss (1), · · · ,r ss (m)) , where r, s = 1, 2, r = s, as n → ∞ is Gaussian with zero mean vector and the covariance matrix equals to I 3m .
Proof. For a fixed m < n, where n is large, the first m sample autocorrelations of the innovations are asymptotically normal with a mean of (ρ 11 (1), ρ 11 (2), · · · , ρ 11 (m)) and covariance matrix n −1 W , where the ( , k)th element of W are given by
where ρ 11 (m) is the population autocorrelation of the innovations at lag m. Under the assumption A1-A3, Bartlett (1946) showed that ρ 11 (m) = 0 for all m = 0; hence the asymptotic distribution of r 1 = √ n(r 11 (1), · · · ,r 11 (m)) is multivariate normal with mean vector zero and identity covariance matrix (Anderson and Walker, 1964) . Similarly, under the assumption A1-A4 and from
Theorem 14 of Hannan (1971, p.228 ) and the result 2c.4.12 of Rao (1973) , McLeod and Li (1983) showed that the limiting distribution r 2 = √ n(r 22 (1), · · · ,r 22 (m)) is multivariate normal with mean vector zero and identity covariance matrix. Also, under the same assumptions, Psaradakis and Vávra (2019) showed that the asymptotic distribution of r 12 = √ n(r 12 (1), · · · ,r 12 (m)) is
Gaussian with zero mean vector and identity covariance matrix. Furthermore, the previous assumptions imply that r 1 , r 2 , and r 12 are independent. Therefore, we may conclude that, as n → ∞, the distribution of √ n(r 11 (1), · · · ,r 11 (m),r 12 (1), · · · ,r 12 (m),r 22 (1), · · · ,r 22 (m)) is Gaussian with zero mean vector and covariance matrix equals to I 3m . Similarly, it is straightforward to show that √ n(r 22 (1), · · · ,r 22 (m),r 21 (1), · · · ,r 21 (m),r 11 (1), · · · ,r 11 (m)) converges weakly to the standard normal distribution on R 3m .
Using the results from Theorem 1, under the adequacy of the fitted arma (p, q) model, and from the theorem on quadratic forms given by Box (1954) , it is straightforward to conclude that C rs and C rs will be asymptotically distributed as chi-square with 3m−(p−q) degrees of freedom.
This section presents the simulation results regarding the finite-sample properties of the asymptotic results of the proposed tests. For illustrative purposes, we also consider the three statistics Q 22 and (Q 12 , Q 21 ) given by (11) and (13). The comparative study is similar to that provided in Psaradakis and Vávra (2019) where the data-generating processes (DGP) are based on the following eighteen models: (2019)).
For each of the M1-M18 models, we use R to simulate 1,000 independent trajectories, each is a series of length n + n/2 with n ∈ {200, 500, 1, 000}, but only the last n data points are used to carry out portmanteau tests to the residuals at different lags m ∈ {5, 10, 15, 20, 25, 30}
(R Core Team, 2020). Figure 1 illustrates the accuracy of the approximation of the empirical distribution ofĈ rs and C rs , r = s ∈ {1, 2} defined by (15) and (16) the sample size), and our preliminary analysis 5 indicate that the portmanteau tests based on the statistics C rs control the Type I error probability more successfully than the tests based on the statisticsĈ rs ; hence, we recommend the use of C rs , r = s ∈ {1, 2}.
In this section, we calculate the type I error probability (of nominal size 1% and 5% 6 ) based on the five test statistics, (C 12 , C 21 ), (Q 12 , Q 21 ), and Q 22 defined by (16) of fitting a true model, for each n ∈ {200, 500, 1, 000}, are averaged across the non-Gaussian linear (M1-M5) and nonlinear (M6-M14) models. In addition, the relative rejection frequencies outside the 99% significant limits are put in boldface. As shown in Table 2 , the simulation results clearly indicate insensitivity to non-Gaussianity noise and the empirical level (regardless of the sample size) agrees very well with the corresponding nominal size. 
The powers of the proposed tests C 12 and C 21 are compared with the powers of Q 12 , Q 21 , and Q 22 statistics for nominal levels α = .01, .05, and .10 7 . Similar to the simulation design in Section 4.1, we generate artificial series of lengths n ∈ {200, 500, 1, 000} from M1-M18 models with innovations having either normal or non-normal distributions. For the non-normal distribution, As seen in Figures 5-7 , the proposed tests C 12 and C 21 both, substantially, have higher rejection frequencies than the Q 12 , Q 21 , and Q 22 statistics and in all cases, the test C 21 is the more powerful test.
In this section, the aforementioned statistics are applied to detect nonlinear dependency in the daily adjusted log-returns of Amazon.com, Inc., spanning the period January 02, 2019 to Decem- estimation. The order of the best of fitted model is selected by minimizing the Bayesian Information Criterion (BIC) over the arma (p, q) models, where (p, q) ∈ {0, 1, · · · , 8(n/100) 1/4 } as suggested by Ng and Perron (2005) . The BIC suggested that no arma process (linear) is detected, in the log-returns of the Amazon.com, Inc., neither before nor after the spread of the disease. The simulation study suggested the proposed tests, C 12 and C 21 , appeared to be as good as the McLeod and Li (1983) statistic in detecting long memory processes and more powerful than Psaradakis and Vávra (2019). As seen in Table 3 , the two tests of Psaradakis and Vávra (2019) fail to detect the nonlinearity structure in both series. However, when the McLeod and Li (1983) and the proposed statistics are used, a clear indication of nonlinear structure appears.
In this article, we propose four goodness-of-fit tests to detect various types of linear and nonlinear dependency in stationary time series models. The proposed tests are based on a linear combination of three auto-and-cross-correlation components. The first and the second components are based on the autocorrelations of the residuals and their squares, respectively, whereas the third is based on the cross-correlations between the residuals and their squares. Two tests can be seen as an extended modification version to the Ljung and Box (1978) test, and the others might be considered as an extended modification version to the Box and Pierce (1970) test. Our simulation study recommends to use the former tests as they are insensitive with respect to non-Gaussianity assumption of the noise, controlling the Type I error probability more successfully than the other tests, and almost always having more power than McLeod and Li (1983) and Psaradakis and Vávra (2019) tests. The idea discussed in this article might be extended to detect seasonality in time series and to identify various types of nonlinearity dependency in multivariate time series.
When the underlying process is assumed to be uncorrelated but yet is not independent (weak assumptions), one might be able to propose a new portmanteau test. In such a case, the Monte Carlo significance test as suggested by Lin and McLeod (2006) and Mahdi and McLeod (2012) might be used to calculate the p-values based on approximating the sampling distribution of the proposed test.