MIXTURE INNOVATIONS’ BASED AUTO-REGRESSIVE PROCESSES WITH APPLICATION TO SEA LEVEL RISE DATA By Fatimah Alshahrani A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Statistics - Doctor of Philosophy 2020 ABSTRACT MIXTURE INNOVATIONS’ BASED AUTO-REGRESSIVE PROCESSES WITH APPLICATION TO SEA LEVEL RISE DATA By Fatimah Alshahrani Sea Level Rise (SLR) is one of the major challenges the world is facing in times of climate change. Using different approaches in obtaining, studying and modeling SLR data has often resulted in different conclusions for the same data set. This led to a debate between researchers some of whom reported accelerations whereas others advocated decelerations or linear trend. In this project, we model SLR as a linear mean reverting time series, also known as order one autoregressive processes (AR(1)), together with an in-depth study of the innovations associated with the AR process. In this direction, we define and evaluate moment ratio test to distinguish between innovations which come either from a Gaussian distribution, a Gumbel distribution or a mixture of the two. We were motivated by a claim from members of the Actuarial Climate Index (ACI) collaborative regarding the distribution of sea-level rise in the US Atlantic states. They claim that the monthly data, recorded over decades, have no inherent time structure and can be modeled as independent and identically observations from Gumbel distribution. Based on an exploratory data analysis, fitting a stationary, Markovian, mean-revering process to the data was essential, to understand its structure. Since the Ornstein–Uhlenbeck process (OU) is the high-frequency liming process for AR(1) processes, and our data has a fixed frequency of observations, we investigated fitting an AR(1) process where the innovations are either Gumbel noise, Gaussian noise, or a mixture of the two distributions. The result of these fittings, based on our moment-ratio tests and other diagnostic checks, demonstrates the preference for an AR(1) process with innovations from a mixture of Gaussian and Gumbel distributions. ACKNOWLEDGMENTS I am very grateful to my advisor and my committee Co-Chair Prof. Frederi Viens for his endless support, encouragement and patience. We had so many enlightened meetings and discussions. Prof. Viens introduced me like others to Malliavin Calculus, which is a big help in proving many theorems in Mathematical Statistics. He encouraged me to participate, attend and give talks in conferences and seminars. Prof. Viens was very supportive during my three years studying at MSU. I also would like to thank my committee Co-Chair Prof. Shrijita Bhattacharya for her sup- port and time. She and I had exchanged emails and held meetings together for two years while working on this project. Throughout our meetings, I learnt many techniques in Time Series Analysis and how to implement them. I also would like to extend my gratitude to my committee members Prof. Alexandra Kravchenko and Prof. Lyudmila Sakhanenko for accepting to serve on my committee and their valuable discussion. In the beginning of this project we had had couple of meetings with Mr. Steve Kolk from The Actuarial Climate Index, and he had explained to us the data and some of their work. So thank you so much Mr. Kolk; the discussions with you made it easier for us to find the data and analyze it. All of my study here at MSU was financially supported by Princess Nora bint Abdul Rahman University (PNU) through The Saudi Arabian Cultural Mission (SACM), so thank you so much PNU and SACM. To my family, my mom and dad, to my husband and kids, I am very thankful for your support, love and believing in me. iv TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The hypothesis tests that are related to the distribution of the innovations in the autoregressive models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Studying and modeling the sea level data . . . . . . . . . . . . . . . . . . . . . . . The moments of the Gumbel random variable and the theoretical moments of the AR(1) model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The Moments of the Gumbel Random Variable The theoretical moments of the AR(1) model . . . . . . . . . . . . . . . . . . . . The asymptotic variances and covariances of the empirical moments of the AR(1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . model The asymptotic distributions of (cid:98)m2, (cid:98)m3 and (cid:98)m4 and the non-central chi- The asymptotic normality of (cid:98)m2, (cid:98)m3 and (cid:98)m4 square S and K tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The asymptotic normality of the empirical skewness and kurtosis . . . . . . . . . The non-central chi-square tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . The Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The distributions of ST and KT . . . . . . . . . . . . . . . . . . . . . . . . . . . . The non-central Chi-square test S and its empirical power . . . . . . . . . . . . . . . . . . . . . . . . . . The non-central chi-square test K and its empirical power 1 5 14 21 21 29 37 48 49 66 70 77 77 80 87 93 The analysis of the sea level rise data . . . . . . . . . . . . . . . . . . . . . . . . 94 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 v LIST OF TABLES Table 1: Some non-Gaussian time series . . . . . . . . . . . . . . . . . . . . . . . . 11 Table 2: The distribution of ST using empirical moments and the moments of Xt when simulating from Gaussian AR(1) model . . . . . . . . . . . . . . . . Table 3: The distribution of ST using empirical moments and using the moments of Xt when simulating from AR(1) model with Gumbel innovations . . . . . Table 4: The distribution of KT using empirical moments and the moments of Xt when simulating from Gaussian AR(1) model . . . . . . . . . . . . . . . . Table 5: The distribution of KT using empirical moments and the moments of Xt when simulating from AR(1) model with Gumbel innovations . . . . . . . Table 7: Table 6: Applying the non-central chi-square S test to reject the Gaussian innova- tions in favor of the Gumbel. . . . . . . . . . . . . . . . . . . . . . . . . . The empirical power of the non-central chi-square test S where H0 : ξt ∼ Gaussian vs Ha : ξt ∼ Gumbel . . . . . . . . . . . . . . . . . . . . . . . . Table 8: Applying the non-central chi-square S test to the hypothesis H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian . . . . . . . . . . . . . . . . . . . . . . . . The empirical power of the non-central chi-square test S where H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian . . . . . . . . . . . . . . . . . . . . . . . . Table 9: Table 10: Applying the non-central Chisquare test S when simulating from an AR(1) with mixture innovation to reject the Gumbel and the Gaussian innovations . . . . . . . . . . . . . . . . . . . . . . . . . . . . in favor of the mixture. Table 11: Applying the non-central Chisquare test S when simulating from an AR(1) with mixture innovation to reject the Gumbel and the Gaussian innovations . . . . . . . . . . . . . . . in favor of the mixture innovations continuing. Table 12: The empirical power of the non-central chi-square test S where H0 : ξt ∼ Gumbel or Gaussian vs Ha : ξt ∼ mixture . . . . . . . . . . . . . . . . . . Table 13: Applying the non-central chi-square K test to reject the Gumbel innova- . . . . . . . . . . . . . . . . . . . . . . . . tions in favor of the Gaussian. vi 78 78 79 79 81 82 83 84 85 86 87 89 Table 14: The empirical power of the non-central chi-square test K where H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian . . . . . . . . . . . . . . . . . . . . . . . . tions in favor of the Gumbel. Gaussian vs Ha : ξt ∼ Gumbel Table 15: Applying the non-central chi-square K test to reject the Gaussian innova- . . . . . . . . . . . . . . . . . . . . . . . . . Table 16: The empirical power of the non-central chi-square test K where H0 : ξt ∼ . . . . . . . . . . . . . . . . . . . . . . . . Table 17: The estimators of (cid:98)φ, (cid:98)σ and (cid:98)p when fitting AR(1) model with different innovations, to the SLR data, . . . . . . . . . . . . . . . . . . . . . . . . . Table 18: The theoretical moments of AR(1) model when fitted with different inno- . . . . . . . . . . . . . . . . . . . . . . . . . . . vations to the SLR data. Table 19: The asymptotic distribution of ST and KT when fitting the AR(1) model with different innovations to the SLR data . . . . . . . . . . . . . . . . . . 3/2 Table 20: The confidence intervals of m3/m 2 and m4/m2 2 applied to the SLR data 90 91 92 98 98 99 99 Table 21: The results of applying the non-central chisquare test T to the SLR data under different null hypothesis scenarios. . . . . . . . . . . . . . . . . . . . 100 Table 22: The p values of the S test applied to the SLR data with different scenarios. 101 Table 23: The p value of KS test under the null hypothesis applied to the residuals of the SLR data with different scenarios. . . . . . . . . . . . . . . . . . . . 105 vii LIST OF FIGURES Figure 1: Correlogram of the sea level data at the CEA region. . . . . . . . . . . . Figure 2: The path of the de-trended data from 1990 to 2000. . . . . . . . . . . . . 2 3 Figure 3: Map shows the 12 different regions according to ACI . . . . . . . . . . . 94 Figure 4: The QQ-plot of the residuals of fitting a Gaussian AR(1) model to the STR data vs the Gaussian distribution . . . . . . . . . . . . . . . . . . . 103 Figure 5: The QQ-plot of the residuals of fitting an AR(1) model with Gumbel innovations to the STR data vs the Gumbel distribution . . . . . . . . . 104 Figure 6: The QQ-plot of the residuals of fitting an AR(1) model with mixture innovations to the STR data vs the mixture distribution . . . . . . . . . 104 viii Introduction Sea Level Rise (SLR) is one of the major consequences the world is facing in times of climate change. Since the result of SLR involves social-economic response, regional and global SLR has been widely investigated and studied. Using different approaches in obtaining, studying and modeling SLR data resulted in different conclusions.”This has led to an intensive debate on the existence of significant accelerations in regional and global sea level in recent years. Some researchers report accelerations where others find linear trends or even decelerations, sometimes based on the same data.” [87] Sea level data can be obtained either by climate models or measuring sea level, using methods such as tide gauge stations which is known as the historical observations or through satellite altimetry. According to sea level research group at university of Colorado, tide gauge data ”offers the only source of historical, precise, long-term sea level data.” We obtained tide gauge data from The Actuarial Climate Index (ACI), and they claim that the monthly data, recorded over decades, of the Central East Atlantic (CEA) region follow the Gumbel distribution in the absence of any time structure. [5]. This claim is investigated here, considering the data as a random sample, and fitting the Gumbel dis- tribution and applying the Kolmogorov–Smirnov test (KS test), we fail to reject the null hypothesis that the data as a random sample follows the Gumbel distribution. Looking at the auto-correlation function plot of the data in figure 1, we conclude that the assumption 1 that the data is a ”random sample” is not accurate. Figure 1: Correlogram of the sea level data at the CEA region. Figure 2shows the path of the data from 1990 to 2000, which looks like the path of the Ornstein–Uhlenbeck (OU) process. The OU process is a stationary, Gaussian, Markovian and mean reverting process. Since the Ornstein–Uhlenbeck process (OU) is the high-frequency limiting process for AR(1) processes, and our data has a fixed frequency of observations [30], we investigate fitting an AR(1) process where the innovations are Gumbel noise, Gaussian noise, or a mixture of both distributions. Having three different models with the earlier debate in mind we defined a non-central chi square S and K tests based on the ratio of the empirical moments of the AR(1) model, such that these tests can reject the Gumbel innovations in favor of the Gaussian, the Gaussian innovations in favor of the Gumbel and the Gumbel or Gaussian innovations in favor of the mixture of both. 2 Employing the moments and the ratio of moments started in the 1930s, under the as- sumption of independent identically random variables. Many researchers were interested in checking the skewness, kurtosis and normality [23]. Figure 2: The path of the de-trended data from 1990 to 2000. The result of fitting the three models to the SLR data, based on our non-central chi-square tests and other diagnostic checks that are based on the residuals, shows the preference for an AR(1) process where the innovations follow a mixture of Gaussian and Gumbel distributions. This thesis is organized as follows: the related work in both parts the goodness of fit tests of the AR(1) model and modeling the sea level data in chapter 1. In chapter 2, the moments of the Gumbel random variable and the theoretical moments of the AR(1) model were calculated. The asymptotic distributions of the empirical second, third and fouth moments are provided in chapter 3, also the main results which are the defined non-central 3 chi-square S and K tests are in chapter 3. The simulations are in chapter 4, while the case study ”the real data” in chapter 5. 4 The Literature review The hypothesis tests that are related to the distribution of the innovations in the autoregressive models The autoregressive process is one of the most successful statistical models for time series data due to its wide range of applications. An autoregressive process of order p (AR(p)), is written as: Xt = φ1Xt−1 + φ2Xt−2 + ... + φpXt−p + t (1) where t is white noise (WN(. , .)), i.e E(t) = 0, E(2 t ) = σ2 and t is uncorrelated with s for each s < t. Autoregressive Processes are a sub-class of a more general linear class called Autoregressive Moving Average Process, which for orders p and q, (ARM A(p, q)), has the form: Xt − φ1Xt−1 − ... − φpXt−p = t + θ1t−1 + ... + θqt−q (2) where t ∼ W N (0, σ2) and the polynomial (1−φ1z−...−φpzp) and (1+θ1z+...+θqzq) have no common factors for identifiability. It is further assumed that Xt is stationary and invertible. Let B be the backward shift operator, (BjXt = Xt−j, Bjt = t−j, j = 0,±1,±2, ...), then 2 can be written as φ(B)Xt = θ(B)t (3) 5 where φ(z) = 1− φ1z − ...− φpzp and θ(z) = 1 + θ1z + ... + θqzq. Xt is an AR(p) if θ(z) ≡ 1, and a moving average process of order q, (M A(q)) if φ(z) ≡ 1. It is further assumed that Xt is a stationary and invertible. The necessary and sufficient condition for stationarity is that for the polynomial φ(z) (cid:54)= 0 for all complex z with |z| = 1. Invertibility is equivalent to the condition θ(z) (cid:54)= 0 for all |z| ≤ 1. In most applications t is assumed to be Gaussian, hence if the length of realization is short, exact likelihood estimation can be used to estimate the parameters φi, i=1,2,...,p; θj, j=1,2,...,q. Alternatively, they can be asymptotically efficiently estimated by maximizing the conditional log-likelihood. Statistical analysis of time series started in the 1920s. In 1927, Udny Yule fitted an AR(2) to Sunspots data [69, 84], making this the first application of autoregressive models to real data. Yule introduces the concept of corrolgram ”serial correlation coefficient,” though he did not invent the name according to [66]. Anderson [3] stated ”Both Yule [94] and Bartlett [6] have shown that the ordinary test of significance are invalidated if successive observations are not independent of one another. The serial correlation coefficient has been introduced as a measure of the relationship between successive values of a variable ordered in time or space”. This statement indicates a fundamental and standard approach in diagnostic checks in time series models. In the early 1970s, with the computer progression and the so-called ”Box–Jenkins method”, which is a detailed strategy for time series forecasting, The ARMA process was studied and investigated widely [58]. The Box-Jenkins method is a method for analyzing discrete time se- ries data, incorporating an iterative cycle of identification, estimation and verification [12, 2]. In the first step, based on the information on the sample path, a preliminary ARMA model is 6 suggested. In the second step the estimation of a stationary ARMA model can be done using the exact or approximation methods. The third step is model diagnostic checking, which involves techniques like overfitting, residual plots and a group of tests that are preformed on the residuals. A good model in time series data should be able to capture the dependence structure of the data appropriately. That is, checking the whiteness of the residuals, since a good time series model should produce residuals that are approximately white noise. Box and Jenkins [12] employ Cramer statistics [24] g1 and g2 for testing the normality of the residuals. In this frame, these are the residuals skewness and kurtosis defined respectively as: (cid:16) 1 n(cid:88) n t=1 K3 = n(cid:88) n t=1 t (cid:17)(cid:46)(cid:16) 1 (cid:98)3 (cid:17)(cid:46)(cid:16) 1 n(cid:88) t=1 (cid:98)2 t n t (cid:17)3/2 (cid:98)2 (cid:17)2 − 3. and (cid:16) 1 n n(cid:88) t=1 K4 = (cid:98)4 t Under the assumption of normality and if the model is correct, K3 ∼ AN (0, 6/n) and K4 ∼ AN (0, 24/n) [58, 39]. As mentioned by earlier statisticians, one important measurement of dependence is the autocorrelation function. This led to using the residual autocorrelation function ˆrk as a test statistic, defined as follow: (cid:98)rk = n(cid:88) t=k+1 (cid:98)t(cid:98)t−k (cid:46) n(cid:88) t=1 (cid:98)2 t , (4) k = 1, ..., m. If the model is appropriate and m >> n, then approximately(cid:98)rk ∼= 0, for all k = 1, ..., m. 7 An informal analysis of these quantities is the core of the third step of the ”Box–Jenkins first to obtain the distribution under the autoregressive model”. method”. For the ARMA model, Box and Pierce [14] derived the asymptotic distribution of (cid:98)r, while for the AR model, Li [58] stated ”As noted by Hosking [42], Walker [88] was the In addition to the asymptotic distribution of(cid:98)rk, Box and Pierce [14] defined a several of statistics based on rk and(cid:98)rk, where n(cid:88) (cid:46) n(cid:88) rk = tt−k 2 t . t=k+1 t=1 Using the results of [3, 4], the limiting distribution of r = (r1, ..., rm)(cid:48) is multivariate with mean vector zero, var(rk) = (n−k)(cid:14)n(n+2) ≈ 1/n and cov(rk, rl) = 0 if k (cid:54)= l, they showed that, if the model is appropriate and the parameters are known then: (cid:101)Q(r) = n(n + 2) (n − k)−1r2 k, (5) m(cid:88) k=1 m(cid:88) would asymptotically follow χ2 m. Furthermore, they suggested the following approximation Q(r) = n k ∼ χ2 r2 m (6) based on the approximated variance. If the parameters are estimated and the model is k=1 appropriate, Box and Pierce [14] showed that, Q((cid:98)r) = n 8 m(cid:88) k=1 (cid:98)r2 k (7) would be asymptotically χ2 m−p−q yielding an approximate test for lack of fit. The statistics (7) is called the Box-Pierce statistic. Davies, Triggs and Newbold [25] and Prothero and Wallis [73] doubted the distribution of Q((cid:98)r) and its validity in practice for a moderated n. This led to the modified statistic: (cid:101)Q((cid:98)r) = n(n + 2) m(cid:88) (n − k)−1(cid:98)r2 k. (8) which was recommended by Ljung and Box [61]. The statistic (cid:101)Q((cid:98)r) is called the Ljung-Box k=1 statistic, and it has a finite sample distribution that is closer to that of χ2 m−p−q. In addition, they used the result by Anderson and Walker [4] to justify their belief in insensitivity of these tests to departures from normality of the t’s. Another modification on Box-Pierce statistic was suggested by Li and Mcleod [60] as follows: If the model is adequate then Q∗ k2m(m + 1) m = Q((cid:98)r) + Q∗ m asymptotically χ2 with k2(m − p − q) degree of freedom. (9) . 2n All the previous statistics to test weather or not a residual series is white noise, under the assumption of normality of the innovations, are based on the residuals’ or innovations’ autocorrelation functions. There are other approaches to check the whiteness of the residuals. For example, using nonparametric methods, Ducharme and Micheaux [29] employ a data- driven smooth test approach to define a goodness of fit test of normality for the innovations of an ARMA(p, q) model. Since the sum and the conditional distribution of Gaussian random variables are Gaus- 9 sian, these facts encourage widespread interest in and investigation of the Gaussian AR(1) model. But in many areas of science, the observations of linear time series models are of non-Gaussian nature and are very common. For instance, one of the earliest applications of this type of model was in hydrology. O’Connell and Jones [71] consider linear time series models driven by lognormal white noise. The ARMA models driven by non-Gaussian inno- vations have been studied and investigated by many scholars, for example, Li and McLeod [71]. In particular, the simplicity of AR(1) structure, its usefulness and interpretation in a wide range of context, led to the building and studying of more than 30 non-Gaussian AR(1) models, according to [36]. One of the approaches in building the non-Gaussian AR(1) model is to specify the dis- tribution of the innovations. In some cases the marginal distribution will be obtainable, for example the Cauchy AR(1) model by Brockwell and Davis [15]. Let Xt be an AR(1) model with t ∼ C(λ, 1) and λ ∈ R, then the marginal distribution of Xt is Cauchy with mean λ 1−φ2 . An example of non-Gaussian AR(1) with unknown marginal distribution is the Uniform AR(1) defined by Bell and Smith [7]. Let Xt be an AR(1) model with t ∼ U (0, 2λ), the Xt is a Uniform AR(1) model with unknown marginal distribution. 1−φ and variance 1 Another approach in constructing the non-Gaussian AR(1) model is starting with a specific marginal distribution of Xt, and then define a distribution for the innovations, even if it is not a known distribution. An example of this approach is the Logistic AR(1) [76]. Let Xt be an AR(1) model with t has a distribution of the form −∞ < z < ∞, then Xt ∼ Lg(µ, 1), where µ = λ sin(φπ) 2φπ[cosh(z−λ)+cos(φπ)], where 1−φ . Table (1) contains some examples of non-Gaussian AR(1), note that the parameter µ is a function of the parameter λ as in the 10 Logistic AR(1). Name Reference GENTS Lye and Martin [62] Laplace Dewald and Lewis [27] Bell and Smith [7] And˘el [1] EAR Gaver and Lewis [35] (cid:40) (cid:40) Innovation distribution GT (µ, σ2) w.p.φ2 λ La(λ, 1) w.p.1 − φ2 Marginal distribution Unknown La(µ, 1) Exp(1/λ) Unknown 0 Exp(1/µ) w.p.1 − φ w.p.φ Exp(1/µ) Table 1: Some non-Gaussian time series Jacobs and Lewis [45, 46] describe steps toward modeling and estimating non-Gaussian time series, using the predesignated marginal distribution approach. When O’Connell and Jones [71] consider fitting an autoregressive model with log-normal innovations, the Yule- Walker equation is used to estimate the autoregressive parameters. Li and McLeod [59] study the asymptotic normality of the maximum likelihood estimator for non-Gaussian ARMA models under some general conditions. Other alternative approaches to modeling non-Gaussian time series have been considered including the Bayesian forecasting models of West and Migon [90], or state space models approach of Kashiwagi and Yanagimoto [52], and through the transformation as in Forecasting non-normal time series by Swift and Janacek [82]. As in the Gaussian ARMA model, Li and McLeod [59] prove the asymptotic normality of the residual autocorrelation function for the non-Gaussian ARMA model. Under the 11 assumption in [59] it can be shown that: n(cid:88) (cid:98)rk = ((cid:98)t − ¯)((cid:98)t−k − ¯) (cid:46) n(cid:88) t=1 ((cid:98)t − ¯)2, (10) (cid:98)t n and where(cid:98)t are residuals from the fitted non-Gaussian ARMA(p,q) model, ¯ = (cid:80)n t=k+1 t=1 k = 1, ..., m , have an asymptotic multivariate normal distribution similar to that of (4) with a different information matrix. The theoretical skewness and kurtosis of ARMA models with non-Gaussian innovations have been studied by Davies, Spedding, and Watson [26, 89]. The purpose of this study is ”to deduce the appropriate skewness and kurtosis that need to be imposed on the noise in order to generate a series, Xt, with given skewness and kurtosis” [26]. In their paper, Davies et al., [26] obtained approximations of skewness and kurtosis of ARMA models with non-Gaussian innovations in terms of skewness and kurtosis of their innovations by making use of the Pearson family of distributions. Let Xt be a series generated by ARMA(1,1) process, which according to (2), can be written as: Xt − φXt−1 = t + θt−1 Davies et al., [26] claim that, (cid:112)β1(X) = {1 + θ(3φ2 + 3φθ + θ2)} {1 + θ(2φ + θ)}3/2 (1 − φ2)3/2 (1 − φ3) (cid:112)β1() (11) 12 where is the skewness of Xt and (cid:112)β1(X) = (cid:112)β1() = E[X3 t ] t ])3/2 (E[X2 E[3 t ] t ])3/2 (E[2 is the skewness of the generating series {t}. And where {(1 − φ4) + (φ + θ)4} {1 − φ2 + (φ + θ)2}2 (cid:104) (φ + θ)2(1 − φ)4 + φ2(φ + θ)4 (1 − φ2) (1 + φ)2 β2() (cid:105) β2(X) = + 6 {(1 − φ2) + (φ + θ)2}2(1 + φ2) (12) is the kurtosis of Xt and β2(X) = E[X4 t ] (E[X2 t ])2 β2() = E[4 t ] (E[2 t ])2 is the kurtosis of {t}. They showed ” it is possible to simulate and match machined metal surface profiles having significantly non normal skewness and kurtosis” by using these tech- niques. As in the diagnostic techniques of the Gaussian ARMA model, other approaches have been investigated to check the whiteness of the residuals of the non-Gaussian ARMA mod- els. Using nonparametric approach, Fan and Zhang [34] ”consider generalised likelihood ratio 13 tests of whether or not the spectral density function of a stationary time series admits certain parametric forms,” which can be applied to test the residuals’ whiteness. Taking advantage of techniques from the empirical processes, many scholars propose goodness of fit tests for the ARMA models. For AR(1) when the innovations have unknown symmetric distribution, Koul [54] gives a class of minimum L2-distance estimators of the autoregression parameter, and based on that discussed a goodness of fit test for symmetry innovations distribution. In addition, Koul and Stute [56] investigate a class of goodness of fit test for autoregressive model. ”These tests are based on a class of empirical processes marked by certain residu- als”, according to the authors. For heavy tailed innovations, Jureckova et al. [51] construct a class of nonparametric tests on the tail index of the innovation distribution in the linear autoregressive model. Studying and modeling the sea level data Studying and investigating sea level rise (SLR) has received attention since the late nine- teenth century. According to Byrne [16] ”Sea-level rise has been tracked since the late 1800s by tide level gauges, and since 1993 changes have been recorded with high precision from altimeter satellites.” Traditionally, sea level change has been estimated from tide gauge measurements of the ’still water level’ collected over the last century.The ’still water level’ (without waves) is averaged over some period of time, such as a month or a year yields Mean Sea Level (MSL) or ”Sea Level” as in some other literature. This MSL is measured relative to fixed marks on 14 the land known as ’benchmarks’. According to Sea Level Research Group at University of Colorado, ”Although the global network of tide gauges comprises of a poorly distributed sea level measurement system, it offers the only source of historical, precise, long-term sea level data”. The Permanent Service for Mean Sea Level (PSMSL) has been responsible for the collection, publication, analysis and interpretation of sea level data from the global network of tide gauges since its establishment in 1933.Alternatively, MSL in the context of satellite altimetry, are the measurements made in a geocentric reference frame (relative to the center of the Earth). Due to the heavy density of population and infrastructure that are located at the coast, sea level rise presents an urgent topic in climate research, with direct socio-economic con- sequences [17, 77, 70]. ”Impact and risk assessment, adaptation policies and long-term decision-making in coastal areas are crucially informed by projections of coastal mean sea level and extreme water level events [47].” Projection sea level rise focuses on the physical climate system – atmosphere, land surface, ocean and sea ice- which serves the Fourth Asssessment Report (AR4) of the Intergovern- mental Panel on Climate Change (IPCC). The Program for Climate Model Diagnosis and Intercomprison (PCMDI) collects model output contributed by leading modeling centers around the world. This collection of recent model output is officially known as the “WCRP CMIP3 multi-model dataset,” where WCRP is ”World Climate Research Programme” and CMIP3 stands for phase 3 of the Coupled Model Intercomparison Project ”archived data of climate models output” [63]. Such organizations and activities facilitate studying and inves- tigating projection of sea level changes widely through different approaches. According to 15 PCMDI ”over 250 journal articles, based at least in part on the dataset, have been published or have been accepted for peer-reviewed publication”. Three different approaches are commonly considered to calculate the projections of fu- ture global sea level change based on climate models. The first approach is process-based projections by estimating the contribution of the primary processes that are causing sea- level change [8], ”using specified forcing of the climate system” [9]. As stated by Jevrejeva et al. [47] ”In a process-based approach, such as that used by Church et al. [20], various sources of information, including numerical simulations, are combined to model the evolu- tion of different sea level components.” Global mean sea level (GMSL) changes over time (t) are represented as: (cid:88) i GICi(t) + GrIS(t) + W AIS(t) + EAIS(t) GM SL(t) = T (t) + (cid:88) + LW Sj(t), (13) j where T is thermostatic expansion, GICi ice mass loss from different glaciers, GrIS ice mass loss from the Greenland Ice Sheet, W AIS and EAIS and ice mass loss from the Antarctic Ice Sheet (West and East) and the contribution from land water storage LW Sj changes. Church et al. [20] estimate the median and 90% confidence interval (CI) of sea level change under specific scenarios (representative concentration pathway) as in [68], by the sum of (14) components with their uncertainties. Regional sea level (RSL) projections using process-based approach combine multiple fac- 16 tors. The first is the dynamic sea level contribution (DSL), the contribution of sea level components, scaled by spatial “fingerprints” associated with gravitational, deformation and rotational effects (as in [67]). The next is the background (Bkgd) processes, that ”unassoci- ated with contemporary climate change that contribute to local sea level changes” [47], like localized processes such as sediment compaction. RSL projections at locations x can be represented as: (cid:88) i (cid:88) j RSL(x, t) = T (t) + DSL(x, t) + FGICi (x)GICi(t) + FGrIS(x)GrIS(t) + FW AIS(x)W AIS(t) + FEAIS(x)EAIS(t) + FLW Sj (x)LW Sj(t) + Bkgd(x, t), (14) Where FGICi , FGrIS, FW AIS, FEAIS and FLW Sj are the normalized gravitational, rota- tional and deformational fingerprints associated, with the different components. As illustrated by Meehl et al. [64], this approach is both used to construct the past global sea level changes, and also to project the future. For projecting regional changes, Kopp et al. [53], Slangen et al. [78] and Spada et al. [79] build their studies based on model data from CMIP3, while the result of Perrette et al. [72] are based on model data from the latest suite of coupled climate model results to date CMIP5 [83]. An example of studies based on model data from combining both CMIP3 and CMIP3 was proposed by Slangen et al. [77]. Yin [93] applies GFDL CM2.1 climate model to model projections of rapid sea-level rise on the northeast coast of the United States. As stated earlier there is a huge amount of studies 17 in this regard, the interested reader is advised to consult [47] for more references. The second approach to project sea level changes is called the “semi-empirical” approach, this approach is based on the observed relationship between sea-level change and global temperature change [8]. Bolin et al. [9] ”construct a time series regression model to predict global sea level from global temperature”. Other investigators use the current observations to extrapolate the future behavior of the cryosphere, and use this information in projecting sea level [8, 65]. As mentioned earlier, climate warming has contributed to global sea levels rise through the past century, and they are projected to rise throughout the 21st century [22]. Instead of using climate models or environmental factors as inputs, some authors have used historical sea level records to show sea level changes in mm/yr, that is, the estimation of long- term trends. ”Trends in sea level have traditionally been calculated from long term data sets at a few locations” [44]. Several papers discussing evidence for acceleration from tide gauge records include Woodworth [91], Douglas [28], Holgate and Woodworth [41], Church and White [21], Jevrejeva et al. [48] and Woodworth et al. [92]. According to Visser [87], in general, there are three applications of trend estimation. Trends may be viewed in terms of prediction, highlight nonstationary behavior in a series, and trends are used for ”calibration of recent attempts to model historic thermal expansion as well as glacier and ice sheet melting”. There are five main methods for estimating trends in sea level records, see Chandler and Scott [18] for details on these methods, (i)Exploratory data analysis (or smoothed time series), (ii)Parametric trend estimation, (iii)Nonparametric 18 trend estimation,(iv)Stochastic trend models (examples are ARIMA models and Structural Time series Models) and (v)Miscellaneous models. Applying trend models to sea level data has been done by many investigators. For example Vaziri [86] employs artificial neural network (ANN) and ARIMA models to predict the monthly Caspian Sea level using tide gauge records, whereas Imani et al. [43] applied the same technique to satellite altimetry records. In ”Modeling regional sea level rise using local tide gauge data”, Iz et al. [44] use the empirical trigonometric model to combine various tide gauge data, regardless of their time span, to estimate the representative of regional mean sea level trends. For a review of trend models applied to sea level data, the interested reader is referred to [87]. Using different methods to estimate sea level trends, several recent studies show that the rates of sea level rise (SLR) have been accelerating along the coastal mid-Atlantic region (the so-called “hotspot” of accelerated SLR, as referred to by Sallenger [74]) [31]. Boon et al. [11] use the least-square linear curve-fit method ”simple linear regression” and quadratic regression to calculate the trend after filtering out seasonal and decadal variability, of sea level measurements at ten Chesapeake Bay water level stations (1975-2009). Sallenger et al. [74] ”present evidence of recently accelerated SLR in a unique 1,000- km-long hotspot on the highly populated North American Atlantic coast north of Cape Hatteras” by analyzing tide-gauge records along the North American Atlantic coast. They suggest that changes in ocean dynamic contribute significantly to sea level rise. Monthly mean sea level records from 8 tide gauge stations in the Chesapeake Bay were used as an application of a new method of trend estimation by Ezer et al. [33]. The new method is based 19 on Empirical Mode Decomposition (EMD) and Hilbert-Huang Transform (HHT). The same method (EMD/HHT) is applied by Ezer et al. [31] to study the relationship between sea level variations in Chesapeake Bay and the mid-Atlantic coast and the the Gulf Stream, using both monthly mean sea level records from 10 tide gauge stations and altimeter data from both regions. In another study of measurements of sea level at tide stations along the Northeast Atlantic coastal region of the United States and Canada from 1969 to 2011, Boon [10] obtains results that are in agreement with the observational study by Sallenger et al. [74] and the pattern of modeled projections (climate model) by Yin et al. [93]. According to [10] ”The observational data support the claim that sea level rise is now accelerating at locations from Norfolk, Virginia, to Halifax, Nova Scotia, at rates as high as 0.30mm/y2”. In their paper, ”Modeling of Coastal Inundation, Storm Surge, and Relative SeaLevel Rise at Naval Station Norfolk, Norfolk, Virginia, U.S.A.”,Li et al. [57] applied the coastal modeling system to calculate the local water-surface elevation and storm-surge inundation in Norfolk, Virginia, for varying future relative SLR scenarios. They describe the consequences of potential relative SLR estimates and storm surges on inundation. In the context of low probability with high impact events ”extreme event”, Sweet and Park [81] studies the sea level data from different stations in the East and West coasts of the Unite States, including Norfolk (Sewells Point gauge). Among their results, data from Norfolk shows exponential growth when fitting growth model. 20 The moments of the Gumbel random variable and the theoretical moments of the AR(1) model Since the defined hypothesis tests are based on the empirical skewness and kurtosis of the AR(1), then the theoretical moments for the AR(1) model and for the innovations are needed. In this chapter, we started by proving a general formula for the higher moments of Gumbel random variable, then obtained the first six moments. Then, the second, third, fourth, fifth, sixth and eighth moments for AR(1) model were calculated in section 2. In section 3 we calculate the variance of second, third and fourth empirical moments, and covariance of second and third empirical moments, and second and fourth empirical moments. The Moments of the Gumbel Random Variable In this section we start by proving a general formula for the higher moments of Gumbel random variable. Theorem 0.0.1. Let X be a Gumbel random variable with location parameter µ and scale parameter β, ”X ∼ G(µ, β)”. Then the higher moments around the origin τn, E(Xn), are n ∈ N (15) given by: n(cid:88) k=0 τn = k (cid:18)n (cid:19) Γ(k)(1)(−β)kµn−k, 21 where Γ(k)(1) is the k-th derivative of Gamma function evaluated at 1, and Γ(0)(1) = Γ(1). Proof. We will use mathematical induction to prove : τn = m(n)(0) = m(t)(cid:12)(cid:12)(cid:12)t=0 dtn Γ(1 − βt)eµt(cid:12)(cid:12)(cid:12)t=0 (cid:18)n (cid:19) n(cid:88) dn = = k k=0 Γ(k)(1)(−β)kµn−k When n = 1, then τ1 = m(1)(0) Γ(1 − βt)eµt(cid:12)(cid:12)(cid:12)t=0 d dt = = Γ(1 − βt)eµtµ − Γ(cid:48)(1 − βt)eµtβ = Γ(1)µ − Γ(cid:48)(1)β 1(cid:88) (cid:19) Γ(k)(1)(−β)kµ1−k. (cid:18)1 = k k=0 (cid:12)(cid:12)(cid:12)t=0 22 Now, assume the statement is correct for n = i, s.t. µi−1βeµtΓ(cid:48)(1 − βt) + ...+ τi = m(i)(0) = = + di 0 (cid:19) 1 (cid:18)i µieµtΓ(1 − βt) − dti Γ(1 − βt)eµt(cid:12)(cid:12)(cid:12)t=0 (cid:18)i (cid:19) (cid:19) (cid:18) i µβi−1eµtΓ(i−1)(1 − βt) − (cid:19) µi−1βΓ(cid:48)(1) + ... + (cid:19) (cid:18) i Γ(k)(1)(−β)kµi−k (cid:18)i 1 i(cid:88) i − 1 = µiΓ(1) − = k k=0 (cid:18)i (cid:19) (cid:18) i i i − 1 (cid:19) (cid:12)(cid:12)(cid:12)t=0 βieµtΓ(i)(1 − βt) µβi−1Γ(i−1)(1) − βiΓ(i)(1), 23 since(cid:0)i (cid:1) =(cid:0)i (cid:1) = 1 we get the second last line. And to prove it for n = i + 1, note that 0 i µieµtΓ(1 − βt) − µi−1βeµtΓ(cid:48)(1 − βt) + ...+ (cid:18)i (cid:19) 1 (cid:18)i (cid:19) µβi−1eµtΓ(i−1)(1 − βt) − (cid:19) µi+1eµtΓ(1 − βt) − (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)t=0 i βieµtΓ(i)(1 − βt) (cid:18)i (cid:19) (cid:18)i 0 µi−1β2eµtΓ(cid:48)(cid:48)(1 − βt) + ...+ (cid:18)i (cid:19) µieµtβΓ(cid:48)(1 − βt) (cid:19) (cid:18) i βi+1eµtΓ(i+1)(1 − βt) i − 1 1 i (cid:19) µiβΓ(cid:48)(1) − (cid:18)i (cid:18) i 1 (cid:19) µiβΓ(cid:48)(1) + (cid:19) 1 µβiΓ(i)(1) − i − 1 (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)t=0 (cid:18)i (cid:19) µi−1β2Γ(cid:48)(cid:48)(1) (cid:19) (cid:18)i µβiΓ(i)(1) µ2βi−1eµtΓ(i−1)(1 − βt) − µβieµtΓ(i)(1 − βt) τi+1 = (cid:12)(cid:12)(cid:12)t=0 (cid:12)(cid:12)(cid:12)t=0 1 i 0 0 = + + = = − i − 1 i − 1 di+1 dti+1 m(t) d (mi(t)) dt d dt (cid:34)(cid:18)i (cid:19) (cid:19) (cid:18) i (cid:34)(cid:18)i (cid:19) (cid:19) (cid:18)i µiβeµtΓ(cid:48)(1 − βt) + (cid:18) i (cid:19) (cid:18)i (cid:19) (cid:34)(cid:18)i (cid:18)i (cid:19) (cid:1) = 1 implies: (cid:1) =(cid:0)i µi+1Γ(1) −(cid:104)(cid:18)i (cid:18)i (cid:19) (cid:35) (cid:19) (cid:18) i (cid:18)i (cid:35) µi+1Γ(1) − βi+1Γ(i+1)(1) i − 1 + ... + (cid:19) − = + i + 1 0 0 0 i µβieµtΓ(i)(1 − βt) + notice that(cid:0)i (cid:34) 0 τi+1 = µ2βi−1Γ(i−1)(1) − i + (cid:19)(cid:105) µiβΓ(cid:48)(1) + ... + ... −(cid:104)(cid:18) i (cid:19) i − 1 (cid:18)i (cid:19)(cid:105) i µβiΓ(i)(1) + β(i+1)Γ(i+1)(1) 24 Now applying the identity(cid:0)n (cid:18)i + 1 sion yields to: k µi+1Γ(1) − τi+1 = k k−1 (cid:1) +(cid:0)n−1 (cid:1) =(cid:0)n−1 (cid:19) µiβΓ(cid:48)(1) + ... + ... − (cid:1) to the two terms inside the previous expres- (cid:35) (cid:18)i + 1 (cid:19) µβiΓ(i)(1) + β(i+1)Γ(i+1)(1) i (cid:34) i+1(cid:88) (cid:18)i + 1 1 (cid:19) Γ(k)(1)(−β)kµ(i+1)−k = k k=0 Which is the required. Since the formula contains the higher derivatives of Gamma function evaluated at 1, the following lemma states these values without prove. (see [19] for details) Lemma 0.0.2. Let Γ(z) be the Gamma function, defined on the complex plane as: (cid:90) ∞ 0 Γ(z) = e−ttz−1dt , where R(z) > 0 and let ζ(s) be the Riemann Zeta function, defined as (cid:90) ∞ ζ(s) = 1 Γ(s) 0 xs−1 ex − 1 dx Euler evaluated ζ(2n), n ∈ N , that is ζ(2) = π2 6 , ζ(4) = π4 90 , and ζ(6) = π6 945, then: 1. Γ(1) = 1. 2. Γ(cid:48)(1) = −γ. 25 3. The second derivative 4. The third derivative 5. The forth derivative 6. The fifth derivative Γ(cid:48)(cid:48)(1) = γ2 + ζ(2) π2 6 = γ2 + Γ(3)(1) = −γ3 − 3γζ(2) − 2ζ(3) = −γ3 − γ − 2ζ(3) π2 2 Γ(4)(1) = γ4 + 6γ2ζ(2) + 8γζ(3) + = γ4 + γ2π2 + 8γζ(3) + 3 20 ζ(4) 27 4 π4 Γ(5)(1) = −γ5 − 10γ3ζ(2) − 20γ2ζ(3) − 135 2 ζ(4) − 20ζ(2)ζ(3) − 24ζ(5) = −γ5 − 5 3 γ3π2 − 20γ2ζ(3) − 3 4 γπ4 − 10 3 π2ζ(3) − 24ζ(5) 26 7. The sixth derivative Γ(6)(1) = γ6 + 15γ4ζ(2) + 40γ3ζ(3) + 404 2 (ζ(4))2 + 120γζ(2)ζ(3) + 144γζ(5) + 40(ζ(3))2 + = γ6 + 5 2 γ4π2 + 40γ3ζ(3) + 45 20 γ2π4 + 20γπ2ζ(3) + 144γζ(5) + 40(ζ(3))2 + 2745 8 ζ(6) 61 168 π6 Using these values with Theorem0.0.1, the first six moments for Gumbel random variable are obtained in the following Lemma. Lemma 0.0.3. Let X be a Gumbel random variable with location parameter µ and scale parameter β, ”X ∼ G(µ, β)”, and define E(Xn) = κn then: 1. τ1 = µ + βγ. 2. τ2 = µ2 + β2γ2 + π2 6 β2 + 2βµγ 3. τ3 = µ3 + 3µ2βγ + 3µβ2(γ2 + π2 6 ) + β3(γ3 + γ π2 2 + 2ζ(3)). 4. τ4 = µ4 + 4µ3βγ + 6β2µ2(γ2 + ) + 4β3µ(γ3 + γ π2 2 + 2ζ(3)) π2 6 π4) 3 20 + β4(γ4 + γ2π2 + 8γζ(3) + 27 5. 6. τ5 = µ5 + 5µ4βγ + 10µ3β2(γ2 + + 5µβ4(γ4 + γ2π2 + 8γζ(3) + π2 6 3 20 3 4 ) + 10µ2β3(γ3 + 3γ π2 2 + 2ζ(3)) π4) γπ4 + 10 3 π2ζ(3) + 24ζ(5)) + β5(γ5 + 5 3 γ3π2 + 20γ2ζ(3) + τ6 = µ6 + 6µ5βγ + 15µ4β2(γ2 + π2 6 + 15µ2β4(γ4 + γ2π2 + 8γζ(3) + ) + 20µ3β3(γ3 + γ 3 20 π4) + 6µβ5(γ5 + π2 2 5 3 + 2ζ(3)) γ3π2 + 20γ2ζ(3) + + 3 4 45 20 γπ4 + 10 3 π2ζ(3) + 24ζ(5)) + β6(γ6 + γ4π2 + 40γ3ζ(3) 5 2 γ2π4 + 20γπ2ζ(3) + 144γζ(5) + 40(ζ(3))2 + 61 168 π6) Proof. By applying Theorem 0.0.1 and Lemma0.0.2. Now the moments of a centered ”E(X) = 0 and Var(X) = σ2” Gumbel random variable, is stated in the following corollary. Corollary 0.0.4. Suppose X ∼ G(− √ 6 π γσ, √ 6 π σ), where σ > 0 and γ is Euler’s constant, then 1. τ1 = 0 2. τ2 = σ2 √ 3. τ3 = 2·6 π3 σ3ζ(3) 6 28 4. τ4 = 3·9 5 σ4 √ π5 σ5[ 10 5. τ5 = 36 6 3 π2ζ(3) + 24ζ(5)] 6. τ6 = 36·6 π6 σ6[40(ζ(3))2 + 61 168 π6] Where ζ(.) is the Riemann Zeta function defined in Lemma 0.0.2. The theoretical moments of the AR(1) model In this section we calculate the theoretical moments of AR(1) model. Let ξt = IID(0, σ2) and which implies: Xt+1 = φXt + ξt+1 ∞(cid:88) i=0 Xt = φiξt−i, Let E(Xr t ) = mr and E(ξr t ) = τr, then Lemma 0.0.5. m2 = τ2 1 − φ2 29 Proof. m2 = E[( φiξt−i)2] ∞(cid:88) ∞(cid:88) i=0 (cid:104) ∞(cid:88) i=0 = = E ( φiξt−i)( i=0 j=0 φ2iE[ξ2 t−i] = τ2 1 − φ2 ∞(cid:88) (cid:105) φjξt−j) Since E(ξt) = 0 we keep only the case where i = j. And the third moment is: Lemma 0.0.6. Proof. m3 = τ3 1 − φ3 ∞(cid:88) ∞(cid:88) i=0 (cid:104) ∞(cid:88) i=0 = m3 = E[( φiξt−i)3] = E ( φiξt−i)( ∞(cid:88) ∞(cid:88) (cid:105) φkξt−k) φjξt−j)( i=0 j=0 k=0 φ3iE[ξ3 t−i] = τ3 1 − φ3 Since E(ξt) = 0 we keep only the case where i = j = k. The fourth moment is: 30 Lemma 0.0.7. Proof. m4 = τ4 − 3τ 2 (1 − φ4) 2 + 3τ 2 2 (1 − φ2)2 ∞(cid:88) ∞(cid:88) i=0 i=0 (cid:104) ∞(cid:88) i=0 m4 = E[( φiξt−i)4] = E ( φiξt−i)( φjξt−j)( ∞(cid:88) ∞(cid:88) φkξt−k)( (cid:105) φlξt−l) k=0 l=0 ∞(cid:88) ∞(cid:88) j=0 ∞(cid:88) φ4iE[ξ4 t−i] + 3 φ2(i+j)E[ξ2 t−i]E[ξ2 t−j] = = = (cid:16) 1 i=0 j=0,j(cid:54)=i (1 − φ2)2 − 3τ 2 2 (1 − φ2)2 2 τ4 1 − φ4 + 3τ 2 τ4 − 3τ 2 (1 − φ4) + 2 (cid:17) 1 (1 − φ4) in the second line we keep the cases where all the indices are equal and the indices when two are equal and the other two are equal, and this is the coefficient of the second term in the (cid:1)(cid:0)2 second line, from(cid:0)4 (cid:1)(cid:14)2!. And the summation 2 2 ∞(cid:88) ∞(cid:88) j=0,j(cid:54)=i i=0 φ2(i+j)E[ξ2 t−i]E[ξ2 t−j] is done using the fact that (cid:88) i(cid:54)=k (cid:88) −(cid:88) i,k i=k = 31 as follows ∞(cid:88) ∞(cid:88) j=0,j(cid:54)=i i=0 φ2(i+j)E[ξ2 t−i]E[ξ2 t−j] = τ 2 2 = τ 2 2 = τ 2 2 And the fifth moment is: i=0 ∞(cid:88) ∞(cid:88) (cid:104) i=0 ∞(cid:88) φ2j φ2i j=0,j(cid:54)=i 1 φ2i(cid:104) (1 − φ2)2 − 1 (1 − φ2) − φ2i(cid:105) (cid:105) 1 (1 − φ4) Lemma 0.0.8. Proof. m5 = τ5 − 10τ2τ3 (1 − φ5) + 10τ2τ3 (1 − φ3)(1 − φ2) ∞(cid:88) ∞(cid:88) i=0 i=0 (cid:104) ∞(cid:88) i=0 = = = m5 = E[( φiξt−i)5] = E ( φiξt−i)( φjξt−j)( ∞(cid:88) φkξt−k)( ∞(cid:88) φlξt−l)( l=0 m=0 (cid:105) φmξt−m) j=0 ∞(cid:88) (cid:18)5 (cid:19) ∞(cid:88) (cid:16) i=0 3 ∞(cid:88) ∞(cid:88) k=0 j=0,j(cid:54)=i 1 φ5iE[ξ5 t−i] + φ3i+2jE[ξ3 t−i]E[ξ2 t−j] (cid:17) (1 − φ3)(1 − φ2) − 1 1 − φ5 τ5 1 − φ5 + 10τ2τ3 τ5 − 10τ2τ3 (1 − φ5) + 10τ2τ3 (1 − φ3)(1 − φ2) 32 in the second line we keep the cases where all the indices are equal and the indices when two are equal and the other three are equal, and this is the coefficient of the second term in the (cid:1)(cid:0)2 second line, from(cid:0)5 (cid:1). Again the summation in the scond line is done using the fact that 3 2 (cid:88) = (cid:88) −(cid:88) i,k i=k i(cid:54)=k . Where as the sixth moment is: Lemma 0.0.9. m6 = τ6 − 15τ2τ4 − 10τ 2 (1 − φ6) 3 + 30τ 3 2 15τ2τ4 − 45τ 3 (1 − φ4)(1 − φ2) 2 + + 10τ 2 3 (1 − φ3)2 + 15τ 3 2 (1 − φ2)3 33 Proof. m6 = E[( φiξt−i)6] i=0 ∞(cid:88) ∞(cid:88) (cid:0)6 (cid:1) ∞(cid:88) ∞(cid:88) (cid:1) (cid:0)6 (cid:1)(cid:0)4 ∞(cid:88) ∞(cid:88) φ6iE[ξ6 i=0 i=0 j=0,j(cid:54)=i 3 2 2 2 t−i] + ∞(cid:88) (cid:16) = + + = = (cid:18)6 (cid:19) ∞(cid:88) ∞(cid:88) 2 j=0,j(cid:54)=i i=0 φ4i+2jE[ξ4 t−i]E[ξ2 t−j] φ3i+3jE[ξ3 t−i]E[ξ3 t−j] 6 i=0 j=0 k=0,k(cid:54)=j(cid:54)=i φ2i+2j+2kE[ξ2 t−i]E[ξ2 t−j]E[ξ2 t−k] (cid:17) (cid:16) − 1 1 − φ6 2 + 10τ 2 3 (cid:17) + 15τ 3 2 τ6 1 (cid:16) 1 − φ6 + 15τ2τ4 (1 − φ2)3 − τ6 − 15τ2τ4 − 10τ 2 (1 − φ6) (1 − φ4)(1 − φ2) 1 3 + (1 − φ4)(1 − φ2) 3 + 30τ 3 2 (1 − φ6) 15τ2τ4 − 45τ 3 (1 − φ4)(1 − φ2) + 2 + 10τ 2 3 (1 − φ3)2 + 15τ 3 2 (1 − φ2)3 (cid:17) 1 (1 − φ3)2 − 1 1 − φ6 Again, by repeating the same logic to count the indices, and using the fact i(cid:54)=k . Lastly the eighth moment is: (cid:88) (cid:88) −(cid:88) = i,k i=k 34 Lemma 0.0.10. m8 = = + + 28τ2τ6 − 420τ 2 35τ 2 4 − 210τ 2 (1 − φ4)2 τ8 − 28τ2τ6 − 56τ3τ5 − 35τ 2 2 τ4 − 280τ2τ 2 (1 − φ6)(1 − φ2) 2 τ4 + 315τ 4 2 + 4 + 420τ 2 (1 − φ8) 3 + 840τ 4 2 2 τ4 + 560τ2τ 2 2 3 − 630τ 4 56τ3τ5 − 560τ2τ 2 (1 − φ5)(1 − φ3) 3 + 2 τ4 − 630τ 4 210τ 2 (1 − φ4)(1 − φ2)2 + 2 280τ2τ 2 3 (1 − φ3)2(1 − φ2) + 105τ 4 2 (1 − φ2)4 35 φ5i+3jE[ξ5 t−i]E[ξ3 t−j] + φ4i+4jE[ξ4 t−i]E[ξ4 t−j] φ6i+2jE[ξ6 t−i]E[ξ2 t−j] (cid:0)8 (cid:1) 4 2 ∞(cid:88) ∞(cid:88) j=0,j(cid:54)=i i=0 φ4i+2j+2kE[ξ4 t−i]E[ξ2 t−j]E[ξ2 t−k] φ3i+3j+2kE[ξ3 t−i]E[ξ3 t−j]E[ξ2 t−k] Proof. m8 = E[( (cid:18)8 (cid:19) ∞(cid:88) ∞(cid:88) 6 j=0,j(cid:54)=i i=0 i=0 i=0 t−i] + φ8iE[ξ8 φiξt−i)8] ∞(cid:88) ∞(cid:88) (cid:18)8 (cid:19) ∞(cid:88) ∞(cid:88) (cid:1)(cid:0)4 (cid:0)8 (cid:1) ∞(cid:88) ∞(cid:88) (cid:1)(cid:0)5 (cid:0)8 (cid:1) ∞(cid:88) ∞(cid:88) (cid:0)8 (cid:1)(cid:0)6 (cid:1)(cid:0)4 (cid:1) ∞(cid:88) ∞(cid:88) j=0,j(cid:54)=i j=0 j=0 i=0 i=0 i=0 2 5 2 4 2 3 3 2 2 2 24 τ8 ∞(cid:88) ∞(cid:88) ∞(cid:88) (cid:16) k=0 k=0,k(cid:54)=j(cid:54)=i k=0,k(cid:54)=j(cid:54)=i ∞(cid:88) (1 − φ)8 + 28τ2τ6 1 (cid:16) + 56τ3τ5 + 210τ 2 2 τ4 + 280τ2τ 2 3 + 105τ 4 2 1 (cid:16) (cid:16) − (1 − φ5)(1 − φ3) (1 − φ4)(1 − φ2)2 − − (1 − φ3)2(1 − φ2) (1 − φ2)4 − 1 6 1 (cid:16) τ8 − 28τ2τ6 − 56τ3τ5 − 35τ 2 28τ2τ6 − 420τ 2 2 τ4 − 280τ2τ 2 (1 − φ6)(1 − φ2) 35τ 2 3 56τ3τ5 − 560τ2τ 2 (1 − φ5)(1 − φ3) 2 τ4 − 630τ 4 210τ 2 (1 − φ4)(1 − φ2)2 + + 2 = + + + + = = + + + i=0 j=0 l=0,i(cid:54)=j(cid:54)=k(cid:54)=l φ2i+2j+2k+2lE[ξ2 t−i]E[ξ2 t−j]E[ξ2 t−k]E[ξ2 t−l] 1 (1 − φ6)(1 − φ2) 1 (1 − φ8) − (cid:17) 1 (1 − φ8) (cid:16) + 35τ 2 4 2 (1 − φ2)(1 − φ6) 1 (1 − φ6)(1 − φ2) − (cid:17) (cid:17) 1 (1 − φ4)2 − − 1 1 (1 − φ8) 2 (1 − φ4)2 + (1 − φ8) (cid:17) (cid:17) 2 (1 − φ3)(1 − φ5) 3 8 + 2 (1 − φ8) 6 (1 − φ4)2 − (1 − φ8) (cid:17) (1 − φ4)(1 − φ2)2 + (1 − φ6)(1 − φ2) 3 − 630τ 4 2 τ4 + 560τ2τ 2 2 + 4 + 420τ 2 (1 − φ8) 3 + 840τ 4 2 4 − 210τ 2 2 τ4 + 315τ 4 2 (1 − φ4)2 280τ2τ 2 3 (1 − φ3)2(1 − φ2) + 105τ 4 2 (1 − φ2)4 36 End of the proof. The asymptotic variances and covariances of the empir- ical moments of the AR(1) model In this section we calculate the variances and covariances of (cid:98)mi = (cid:80)T Xt =(cid:80)∞ i=0 φiξt−i, and ξt = IID(0, σ2) with E(ξr) = τr and τr < ∞ for r = 2, 3, 4, 5, 6, 8. t=1 Xi t /T where We start by the asymptotic variance of the second empirical moment. Theorem 0.0.11. Proof. T Var[(cid:98)m2] = (m4 − m2 2) (cid:18) 2 (1 − φ2) (cid:19) − 1 + o(1) Var[(cid:98)m2] = (m4 − m2 2) T + 1 T 2 T(cid:88) T(cid:88) i=1 j=1 Cov(X2 i , X2 j ) Now verify there are exactly T − r terms whose covariance values are same and equal to Cov(X2 1 , X2 1+r). Thus, Var[(cid:98)m2] = (m4 − m2 2) T T−1(cid:88) r=1 + 2 T 2 (cid:124) 37 (T − r)Cov(X2 1 , X2 1+r) (cid:125) (cid:123)(cid:122) CT,r Thus, the problem boils down to the computation of Cov(X2 1 , X2 1+r). Now note that X1+r = φrX1 + where Rr is independent of X1. Thus, r(cid:88) (cid:124) i=1 φr−iξ1+i (cid:125) (cid:123)(cid:122) Rr Cov(X2 1 , X2 1+r) = Cov(X2 1 , φ2rX2 1 + 2φrX1Rr + R2 r) = φ2rVar(X2 = φ2r(m4 − m2 1 ) + 2φrCov(X2 1 , X1Rr) 3) + 2φr(m3E[Rr] − m2m1E[Rr]) First note that m1 = 0 and E[Rr] = 0 which allows us to simplify a bit. Cov(X2 1 , X2 1+r) = φ2r(m4 − m2 3) This implies: CT,r = = = r=1 2 T 2 T−1(cid:88) T−1(cid:88) 2 T 2 2(m4 − m2 3) r=0 T (T − r)φ2r(m4 − m2 3) 3) − 2(m4 − m2 (T − r)φ2r(m4 − m2 (cid:17) (cid:16) 1 (cid:17) (cid:16) T 1 3) − 1 + o (1 − φ2) T 38 The following theorem is for the asymptotic variance of the third empirical moment. Theorem 0.0.12. (cid:104) T Var[(cid:98)m3] = 2 (m6 − m2 3) (cid:18) 1 (1 − φ3) (cid:19) − 1/2 (cid:18) + 3m4τ2 1 (1 − φ)(1 − φ2) − 1 (1 − φ3)(1 − φ2) (cid:19)(cid:105) + o(1) Proof. Var[(cid:98)m3] = (m6 − m2 3) T + 1 T 2 T(cid:88) T(cid:88) i=1 j=1 Cov(X3 i , X3 j ) Now verify there are exactly T − r terms whose covariance values are same and equal to Cov(X3 1 , X3 1+r). Thus, Var[(cid:98)m3] = (m6 − m2 3) T T−1(cid:88) r=1 + 2 T 2 (cid:124) (T − r)Cov(X3 1 , X3 1+r) (cid:125) (cid:123)(cid:122) CT,r Thus, the problem boils down to the computation of Cov(X3 1 , X3 1+r). Now note that r(cid:88) (cid:124) i=1 φr−iξ1+i (cid:125) (cid:123)(cid:122) Rr X1+r = φrX1 + 39 where Rr is independent of X1. Thus, Cov(X3 1 , X3 1+r) = Cov(X3 1 , φ3rX3 1 + 3φ2rX2 1 Rr + 3φrX1R2 r + R3 r) = φ3rVar(X3 = φ3r(m6 − m2 1 , X2 1 ) + 3φ2rCov(X3 1 Rr) + 3φrCov(X3 3) + 3φ2r(m5E[Rr] − m3m2E[Rr]) r] − m3m1E[R2 r]) + 3φr(m4E[R2 1 , X1R2 r) First note that m1 = 0 and E[Rr] = 0 which allows us to simplify a bit. Cov(X3 1 , X3 1+r) = φ3r(m6 − m2 3) + 3φrm4E[R2 r] We simplify further by noting E[R2 r] = τ2(1 − φ2r) 1 − φ2 r=1 T−1(cid:88) (T − r)φ3r(m6 − m2 T−1(cid:88) (cid:16) m6 − m2 (m6 − m2 (T − r) 3m4τ2 (cid:16) r=0 3 2 T 2 2 T 2 − 2 T (1 − φ3) (1 − φ3)(1 − φ2) CT,r = = = 3) + 3φrm4 3)φ3r + 3m4τ2 τ2(1 − φ2r) 1 − φ2 φr(1 − φ2r) 1 − φ2 (cid:17) − 2(m6 − m2 (cid:17) − 2(m6 − m2 T 3) 3) T (cid:16) 1 (cid:17) T + o + 3m4τ2 (1 − φ)(1 − φ2) And the following theorem calculates the covariance between (cid:98)m2 and (cid:98)m3 40 Theorem 0.0.13. T Cov[(cid:98)m3,(cid:98)m2] = (m5 − m3m2) (cid:18) + 3τ2m3 (1 − φ)(1 − φ2) (cid:18) 1 (cid:19) 1 − φ3 − 1 1 − φ2 + − 1 (1 − φ3)(1 − φ2) 1 1 (cid:19) + o(1) Proof. Cov[(cid:98)m3,(cid:98)m2] = 1 , X2 Cov(X3 1 ) T + T−1(cid:88) r=1 1 T 2 (cid:124) (T − r)Cov(X3 1 , X2 1+r) (cid:125) T−1(cid:88) r=1 + 1 T 2 (cid:124) (T − r)Cov(X2 1 , X3 1+r) (cid:123)(cid:122) C2 T,r (cid:123)(cid:122) C1 T,r (cid:125) By repeating logic of the previous proof: Cov(X3 1 , X2 1+r) = Cov(X3 1 , φ2rX2 1 + 2φrX1Rr + R2 r) 1 , X2 1 ) + 2φrCov(X3 = φ2rCov(X3 = φ2r(m5 − m3m2) + 2φr(m4E[Rr] − m3m1E[Rr]) = φ2r(m5 − m3m2) 1 , X1Rr) 41 since E[Rr] = 0. Also, Cov(X2 1 , X3 1+r) = Cov(X2 1 , φ3rX3 1 + 3φ2rX2 1 Rr + 3φrX1R2 r + R3 r) 1 , X3 1 ) + 3φ2rCov(X2 = φ3rCov(X2 = φ3r(m5 − m3m2) + 3φ2r(m4E[Rr] − m2 1 , X2 2E[Rr]) 1 Rr) + 3φrCov(X2 1 , X1R2 r) r] − m2m1E[R2 + 3φr(m3E[R2 r]) φr(1 − φ2r) = φ3r(m5 − m3m2) + 3τ2m3 1 − φ2 since m1 = 0, E[Rr] = 0 and E[R2 C1 T,r = (m5 − m3m2) T 2 r] = τ2(1 − φ2r)/(1 − φ2). Substituting we get, (cid:16) 1 (m5 − m3m2) (cid:16) (cid:17) 1 − φ2 − 1 (T − r)φ2r = + o T 1 (cid:17) T−1(cid:88) r=1 T T−1(cid:88) r=1 (T − r) 1 φr(1 − φ2r) 1 − φ2 − (1 − φ)(1 − φ2) (1 − φ3)(1 − φ2) 1 (cid:17) T−1(cid:88) (cid:16) r=1 (T − r)φ3r + 3τ2m3 (cid:16) (cid:17) 1 − φ3 − 1 + 3τ2m3 1 C2 T,r = (m5 − m3m2) T 2 = (m5 − m3m2) (cid:17) (cid:16) 1 T + o T Here is the asymptotic variance of (cid:98)m4 42 Theorem 0.0.14. T Var[(cid:98)m4] = 2 (cid:18) (cid:104) (m8 − m2 4) (cid:18) + 6τ2(m6 − m4m2) (cid:19) 1 (cid:18) (1 − φ4) 1 − 1/2 (1 − φ2)2 − (cid:19) 1 (1 − φ2)(1 − φ4) (cid:19)(cid:105) + 4m5τ3 1 (1 − φ)(1 − φ3) − 1 (1 − φ3)(1 − φ4) + o(1) Proof. Var[(cid:98)m4] = (m8 − m2 4) T + 1 T 2 T(cid:88) T(cid:88) i=1 j=1 Cov(X4 i , X4 j ) Now verify there are exactly T − r terms whose covariance values are same and equal to Cov(X4 1 , X4 1+r). Thus, Var[(cid:98)m4] = (m8 − m2 4) T T−1(cid:88) r=1 + 2 T 2 (cid:124) (T − r)Cov(X4 1 , X4 1+r) (cid:125) (cid:123)(cid:122) CT,r Thus, the problem boils down to the computation of Cov(X4 1 , X4 1+r). Now note that r(cid:88) (cid:124) i=1 φr−iξ1+i (cid:125) (cid:123)(cid:122) Rr X1+r = φrX1 + 43 where Rr is independent of X1. Thus, Cov(X4 1 , X4 1+r) = Cov(X4 1 , φ4rX4 1 + 4φ3rX3 1 Rr + 6φ2rX2 1 R2 r + 4φrX1R3 r + R4 r) = φ4rVar(X4 1 ) + 4φ3rCov(X4 1 , X3 1 Rr) + 6φ2rCov(X4 1 , X2 1 R2 r) + 4φrCov(X4 = φ4r(m8 − m2 + 6φ2r(m6E[R2 1 , X1R3 r) 4) + 4φ3r(m7E[Rr] − m4m3E[Rr]) r] − m4m2E[R2 r]) + 4φr(m5E[R3 r] − m4m1E[R3 r]) First note that m1 = 0 and E[Rr] = 0 which allows us to simplify a bit. Cov(X4 1 , X4 1+r) = φ4r(m8 − m2 4) + 6φ2r(m6E[R2 r] − m4m2E[R2 r]) + 4φrm5E[R3 r] Now note that E[R2 r] = E[R3 r] = τ2(1 − φ2r) 1 − φ2 τ3(1 − φ3r) 1 − φ3 44 Substituting we get, CT,r = (T − r) (T − r) r=1 2 T 2 T−1(cid:88) T−1(cid:88) − 2(m8 − m2 (cid:16) m8 − m2 4) 2 T 2 r=0 = T 4 = 2 T (1 − φ4) 8m5τ3 + T (1 − φ)(1 − φ3) (cid:16) (cid:16) φ4r(m8 − m2 4) + φ4r(m8 − m2 4) + 6τ2(1 − φ2r) 1 − φ2 6τ2(1 − φ2r) 1 − φ2 φ2r(m6 − m4m2) + φ2r(m6 − m4m2) + 4φm5τ3(1 − φ3r) 1 − φ3 4φm5τ3(1 − φ3r) 1 − φ3 (cid:17) (cid:17) − 6τ2(m6 − m4m2) (1 − φ2)(1 − φ4) − 2(m8 − m2 4) T (cid:17) 6τ2(m6 − m4m2) (1 − φ2)2 + 4m5τ3 (1 − φ3)(1 − φ4) (cid:16) 1 (cid:17) T − + o where o(1/T ) involves terms like T αT . And lastly the asymptotic covariance between (cid:98)m2 and (cid:98)m4 Theorem 0.0.15. Cov[(cid:98)m4,(cid:98)m2] = (m6 − m4m2) + 6τ2(m4 − m2 2) (cid:18) (cid:18) 1 (cid:18) 1 − φ2 + (1 − φ2)2 − − 1 1 1 (cid:19) 1 − φ4 − 1 (1 − φ2)(1 − φ4) 1 (cid:19) (cid:19) 1 + 4m3τ3 (1 − φ)(1 − φ3) (1 − φ3)(1 − φ4) + o(1) 45 Proof. Cov[(cid:98)m4,(cid:98)m2] = 1 , X2 Cov(X4 1 ) T + T−1(cid:88) r=1 1 T 2 (cid:124) (T − r)Cov(X4 1 , X2 1+r) (cid:125) T−1(cid:88) r=1 + 1 T 2 (cid:124) (T − r)Cov(X2 1 , X4 1+r) (cid:123)(cid:122) C2 T,r (cid:123)(cid:122) C1 T,r (cid:125) By repeating logic of the previous proof: Cov(X4 1 , X2 1+r) = Cov(X4 1 , φ2rX2 1 + 2φ2rX1Rr + R2 r) 1 , X2 1 ) + 2φrCov(X4 = φ2rCov(X4 = φ2r(m6 − m4m2) + 2φr(m5E[Rr] − m1E[Rr]) = φ2r(m6 − m4m2) 1 , X1Rr) since E[Rr] = 0. Also, Cov(X2 1 , X4 1+r) = Cov(X2 1 , φ4rX4 1 + 4φ3rX3 1 Rr + 6φ2rX2 1 R2 r + 4φX1R3 r + R4 r) = φ4rCov(X2 1 , X4 1 ) + 4φ3rCov(X2 1 , X3 1 Rr) + 6φ2rCov(X2 1 , X2 1 R2 r) 1 , X1R3 r) + 4φrCov(X2 = φ4r(m6 − m4m2) + 4φ3r(m5E[Rr] − m2m3E[Rr]) r] − m2 + 6φ2r(m4E[R2 = φ4r(m6 − m4m2) + 6τ2(m4 − m2 2) 2E[R2 r]) + 4φr(m3E[R3 r] − m2m1E[R3 r]) (1 − φ3r)φr (1 − φ2r)φ2r 1 − φ2 + 4τ3m3 1 − φ3 46 since E[Rr] = 0, E[R2 r] = τ2(1−φ2r) 1−φ2 , E[R3 r] = T−1(cid:88) (T − r)(m6 − m4m2)φ2r = τ3(1−φ3r) 1−φ3 . Substituting we get, (m6 − m4m2) 1 − φ2 1 T − (m6 − m4m2) T (cid:16) 1 (cid:17) T + o C1 T,r = 1 T 2 r=1 (cid:16) m6 − m4m2 − 6τ2(m4 − m2 1 − φ2 C2 T,r = + 1 (cid:16)6τ2(m4 − m2 T (1 − φ4) 1 T (1 − φ2)4 + 2) 2) (cid:17) (cid:17) − (m6 − m4m2) − 4m3τ3 1 − φ3 T (cid:16) 1 (cid:17) T + o 4m3τ3 (1 − φ)(1 − φ3) 47 The asymptotic distributions of (cid:98)m2, (cid:98)m3 and (cid:98)m4 and the non-central chi-square S and K tests Asymptotic theory plays a very important role in time series, since the exact distributions of statistics of time series can be impossible and complicated to derive, even for Gaussian pro- cesses unless in limited known very special cases. Many classical theorems were first proved under the independent and identically distribution assumptions, and extended to the depen- dent case.One way to do it, is through m-dependence approximation, which first introduced by HOEFFDING AND ROBBINS (1948) [40]. In this chapter we proved the asymptotic normality of the empirical second, third and fourth moments of a class of stationary linear process, in the first section. Then apply the Delta method a conjecture to prove the asymp- totic normality of the empirical skewness and kurtosis of the same class of processes, in section 2. In section 3, we define the non-central chisquare distribution goodness of fit tests S and K, to reject the null hypotheses of Gaussian innovations in favor of Gumbel, Gumbel innovations and Gussian in favor of mixture of both and Gumbel innovations in favor of the Gaussian. 48 The asymptotic normality of (cid:98)m2, (cid:98)m3 and (cid:98)m4 The sequence X1, X2, ... is m-dependent for some m > 0 if the vector (X1, X2, ..., Xi) is independent of (Xi+j, Xi+j+1, ...) whenever j > m. The central limit theorem (CLT) for stationary m-dependent sequence is used to prove the asymptotic normality of the second, third and fourth empirical moments of stationary linear processes. Let Xt be a linear process of the form: Xt = ∞(cid:88) i=−∞ where ξt = IID(0, σ2), and the coefficients satisfy(cid:80)∞ ical moment is: (cid:98)m2 = T(cid:88) t=1 X2 t /T, φiξt−i (16) i=−∞|φi| < ∞. then the second empir- and the third and fourth empirical moments respectively are: (cid:98)m3 = T(cid:88) t=1 X3 t /T and (cid:98)m3 = T(cid:88) t=1 X3 t /T. The CLT for stationary m-dependent sequence and ”an approximation” theorem are the main tools used in the proof. Define Xm t = m(cid:88) i=−m φi mξt−i (17) where φi m = φi if |i| ≤ m and 0 otherwise, then Xm t is 2m-dependent and strictly stationary. 49 The idea of the proof is using truncation to construct an approximating sequence ynm of random variables, that is similar to Xn, but their asymptotic distribution is obtainable. These theorems are stated as follow respectively without proofs, where the proofs can be found in [85, 75]. Theorem 0.0.16. If Xt is a strictly stationary m-dependent sequence of random variables with mean zero and autocovariance function γ(·) and if u=m(cid:88) u=−m Vm = γ(u) where Vm (cid:54)= 0, then . ¯Xn ∼ AN (0, Vm/n) Theorem 0.0.17. Let xn for n = 1, 2, .... and ymn for m = 1, 2, ... be random K × 1 vector such that: (i) ymn d−→ ym as n → ∞ for each m; (ii) ym d−→ y as n → ∞; (iii) Pr{|xn − ymn| > } = 0 for every  > 0. By Tchebycheff inequality, this m→∞ lim sup lim n→0 condition is equivalent to E{|xn − ymn|2} → 0, which is easier to verify. Then, xn d−→ y We begin by providing LLN for Stationary and Ergodic Time Series (Hamilton [38]). 50 Let Xt be a stationary process with mean E[Xt] = µ and autocovariance γj = cov(Xt, Xt−j). If ∞(cid:88) j=o |γj| < ∞ then Xt is ergodic for the mean. That is, ¯X p→ E[Xt] = µ. We state Ergodic Theorem as follows: Theorem 0.0.18. Let Xt be stationary and ergodic with E[Xi] = µ. Then T(cid:88) t=1 ¯X = 1 T p→ E[Xt] = µ. Xt Remark 0.0.19. The ergodic theorem says that for a stationary and ergodic sequence {Xt} the time average converges to the ensemble average as the sample size gets large. Remark 0.0.20. Any transformation g(·) of a stationary and ergodic process {Xt} is also stationary and ergodic. That is,{g(Xt)} is stationary and ergodic. Therefore, if E[g(Xt)] exists then the ergodic theorem gives ¯g = 1 T T(cid:88) g(Xt) p→ E[g(Xt)]. t=1 as an application by take, g(x) = xi, where i = 2, 3, 4 we get (cid:99)mi (cid:80)T Next we prove formulas for the variance and the asymptotic variance of (cid:98)m2. t=1 Xi and mi = E[Xi]. 1 T p→ mi, where (cid:99)mi = 51 Theorem 0.0.21. Let Xt be a linear process of the form: i=−∞|φi| < ∞, and E(ξ4 t ) < ∞, then ∞(cid:88) i=−∞ Xt = φiξt−i where ξt = IID(0, σ2), and the coefficients satisfy(cid:80)∞ ∞(cid:88) T Var((cid:98)m2) → γ2(r) r=−∞ where γ2(r) = Cov(X2 t , X2 t+r) = φ2r(m4 − m2 2) m4 = E(X4 t ) and m2 = E(X2 t ) Proof. Var((cid:99)m2) = = 1 T 2 1 T 2 n(cid:88) T(cid:88) T−1(cid:88) j=1 i=1 r=−(T−1) Cov(X2 i , X2 j ), (T − r)Cov(X2 1 , X2 1+r) 52 Now write Xt as: X1+r = φrX1 + r(cid:88) (cid:124) i=1 φr−iξ1+i (cid:123)(cid:122) (cid:125) Rr where Rr is independent of X1, with E[Rr] = 0 and Var[Rr] = σ2(1−φ2r) 1−φ2 . Then the cross- covariance function look like the covariance function of X1. Thus, Cov(X2 1 , X2 1+r) = Cov(X2 1 , φ2rX2 1 + 2φrX1Rr + R2 r) 1 ) + 2φrCov(X2 1 , X1Rr) 2) + 2φr(m3E[Rr] − m2m1E[Rr]) = φ2rVar(X2 = φ2r(m4 − m2 = φ2r(m4 − m2 2) Let γ2(r) = Cov(X2 t , X2 t+r), then T Var((cid:99)m2) = |φi| < ∞, this implies(cid:80) r T−1(cid:88) (cid:0)1 − |r| (cid:1)γ2(r), T r=−(T−1) i T |γ2(r)| < ∞. Note that: (cid:12)(cid:12)(cid:0)1 − |r| (cid:1)γ2(r) → γ2(r), and take fn(r) =(cid:0)1− |r| (cid:1)γ2(r)(cid:12)(cid:12) < |γ2(r)| and (cid:1)γ2(r), for |r| < n and zero for |r| > n. Let and since(cid:80) (cid:0)1− |r| Theorem we have T Var((cid:98)m2) →(cid:80)∞ The above result for the variance enables proving the following CLT for the statistic (cid:98)m2. µ(r) = 1 for r = ±1,±2, ... be the counting measure then by the Dominated Convergence r=−∞ γ2(r). n T Theorem 0.0.22. If Xt is a stationary linear process and under the assumption of Theorem 53 0.0.21, then: t ) and V =(cid:80)∞ Where m2 = E(X2 r=−∞ γ2(r). (cid:98)m2 ∼ AN (m2, V /T ) Proof. First define the stationary 2m-dependent process ym t = (Xm t )2, where Xm t is the approximation defined in (17) to Xt, then: ¯ymT = T−1 ym t T(cid:88) = T−1(cid:104) t=1 T )2(cid:105) (Xm 1 )2 + ... + (Xm also E(ym t ) = E(Xm t )2 is the second theoretical moment of the 2m dependent process Xm t . Then consider as an approximation to ymT = T 1/2 [¯ymT − E(¯ymT )] yT = T 1/2((cid:98)m2 − m2), where E(¯ymT ) is the same as E(ym t ) given above. To obtain an asymptotic distribution for yT , we apply Theorem 0.0.17, using ymT as our approximation.Now verfying (i),(ii) and (iii) of Theorem 0.0.17. (i) Applying Theorem 0.0.16, (CLT for m-dependent sequence), to the 2m dependent series 54 ymT we obtain: as T → ∞ where Vm =(cid:80)m ymT = T 1/2 [¯ymT − E(¯ymT )] d−→ ym ∼ N (0, Vm) r=−m γ2(r). (ii) Note that as m → ∞, Vm → V , using the Domi- nated Convergence Theorem, where V is defined in Theorem 0.0.22. Hence the characteristic function of ym, ϕm(λ) = exp{−1/2λ2Vm} → exp{−1/2λ2V } as m → ∞, which is the characteristic function of a random variable y ∼ N (0, V ) (character- izing convergence of distribution of Fn(.) in terms of convergence of sequence of characteristic functions ϕn(.), i.e., ϕn(λ) → ϕ(λ) ⇔ Fn(x) d−→ F (x)). (iii)To verify the last condition, note that: E(yT − ymn)2 = T E (cid:104) (cid:105)2 ((cid:98)m2 − E(X2 t )) − (¯ymT − E(¯ymT )) = T E(cid:0)(cid:98)m2 − E(X2 t )(cid:1)2 + T E(cid:0)¯ymT − E(¯ymT )(cid:1)2 − 2T E(cid:0)((cid:98)m2 − E(X2 t ))(¯ymT − E(¯ymT ))(cid:1) = T Var((cid:98)m2) + T Var(¯ymT ) − 2T Cov((cid:98)m2, ¯ymT ) using the calculations that led to the result in Theorem 0.0.21, and letting m → ∞, we get hence T Var((cid:98)m2) + T Var((cid:98)m2) − 2T Var((cid:98)m2), (V + V − 2V ) = 0 55 as T → ∞. In order to prove the limiting distribution of (cid:98)m3 we start by providing formulas for the variance and asymptotic variance of (cid:98)m3. Theorem 0.0.23. Let Xt be a linear process of the form: ∞(cid:88) i=−∞ Xt = φiξt−i where ξt = IID(0, σ2), and the coefficients satisfy(cid:80)∞ ∞(cid:88) T Var((cid:98)m3) → γ3(r) r=−∞ i=−∞|φi| < ∞, and E(ξ6 t ) < ∞, then where γ3(r) = Cov(X3 t , X3 t+r) = φ3r(m6 − m2 3) + 3φrm4 τ 2(1 − φ2r) 1 − φ2 , m6 = E(X6 t ), m3 = E(X3 t ) and m4 = E(X4 t ). 56 Proof. Var((cid:99)m3) = = 1 T 2 1 T 2 T(cid:88) T(cid:88) T−1(cid:88) j=1 i=1 r=−(T−1) Cov(X3 i , X3 j ), (T − r)Cov(X3 1 , X3 1+r) By repeating the same argument in proof of Theorem 0.0.21 we get: Cov(X3 1 , X3 1+r) = Cov(X3 1 , φ3rX3 1 + 3φ2rX2 1 Rr + 3φrX1R2 r + R3 r) = φ3rVar(X3 = φ3r(m6 − m2 1 , X2 1 ) + 3φ2rCov(X3 1 Rr) + 3φrCov(X3 3) + 3φ2r(m5E[Rr] − m3m2E[Rr]) r] − m3m1E[R2 r]) 1 , X1R2 r) + 3φr(m4E[R2 = φ3r(m6 − m2 = φ3r(m6 − m2 3) + 3φrm4E[R2 r] 3) + 3φrm4 τ 2(1 − φ2r) 1 − φ2 Let γ3(r) = Cov(X3 t , X3 t+r), then T Var((cid:99)m3) = |φi| < ∞, this implies(cid:80) (cid:1)γ3(r) → γ3(r), and take fn(r) =(cid:0)1− |r| r n and since(cid:80) (cid:0)1− |r| i T T−1(cid:88) (cid:0)1 − |r| (cid:1)γ3(r), T r=−(T−1) |γ3(r)| < ∞. Note that: (cid:12)(cid:12)(cid:0)1 − |r| (cid:1)γ3(r)(cid:12)(cid:12) < |γ3(r)| and (cid:1)γ3(r), for |r| < n and zero for |r| > n. Let T 57 µ(r) = 1 for r = ±1,±2, ... be the counting measure then by the Dominated Convergence Theorem we have T Var((cid:98)m3) →(cid:80)∞ Now the following theorem proves the asymptotic distribution for (cid:98)m3. Note that the r=−∞ γ3(r). proof will be very similar to the proof of Theorem0.0.22, where the same method and theo- rems are used, but suitable approximations and estimators. Theorem 0.0.24. If Xt is a stationary linear process and under the assumption of Theorem 0.0.23, then: t ) and V =(cid:80)∞ Where m3 = E(x3 r=−∞ γ3(r). (cid:98)m3 ∼ AN (m3, V /T ) Proof. First define the stationary 2m-dependent process ym t = (Xm t )3, where Xm t is the approximation defined in (17) to Xt, then: ¯ymT = T−1 ym t T(cid:88) = T−1(cid:104) t=1 n )3(cid:105) (Xm 1 )3 + ... + (Xm also E(ym t ) = E(Xm t )3 is the third theoretical moment of the 2m dependent process Xm t . 58 Then consider as an approximation to ymT = T 1/2 [¯ymT − E(¯ymT )] yT = T 1/2((cid:98)m3 − m3), where E(¯ymT ) is the same as E(ym t ) given above. To obtain an asymptotic distribution for yT , we apply Theorem 0.0.17, using ymT as our approximation.Now verifying (i),(ii) and (iii) of Theorem 0.0.17. (i) Applying Theorem 0.0.16, (CLT for m-dependent sequence), to the 2m dependent series ymT we obtain: ymT = T 1/2 [¯ymT − E(¯ymT )] d−→ ym ∼ N (0, Vm) as T → ∞ where Vm =(cid:80)m r=−m γ3(r). (ii) Note that as m → ∞, Vm → V , using the Domi- nated Convergence Theorem, where V is defined in Theorem 0.0.24. Hence the characteristic function of ym, ϕm(λ) = exp{−1/2λ2Vm} → exp{−1/2λ2V } as m → ∞, which is the characteristic function of a random variable y ∼ N (0, V ) (character- izing convergence of distribution of Fn(.) in terms of convergence of sequence of characteristic functions ϕn(.), i.e., ϕn(λ) → ϕ(λ) ⇔ Fn(x) d−→ F (x)). (iii)To verify the last condition, note 59 that: E(yT − ymT )2 = T E (cid:105)2 (cid:104) ((cid:98)m3 − E(X3 t )) − (¯ymT − E(¯ymT )) t )(cid:1)2 + T E(cid:0)¯ymT − E(¯ymT )(cid:1)2 = T E(cid:0)(cid:98)m3 − E(X3 − 2T E(cid:0)((cid:98)m3 − E(X3 t ))(¯ymT − E(¯ymn))(cid:1) = T Var((cid:98)m3) + T Var(¯ymT ) − 2T Cov((cid:98)m3, ¯ymT ) using the calculations that led to the result in Theorem 0.0.23, and letting m → ∞, we get T Var((cid:98)m3) + T Var((cid:98)m3) − 2T Var((cid:98)m3), (V + V − 2V ) = 0 hence as T → ∞. The same technique above applied to find the limiting distribution of (cid:98)m4. First we start by providing formulas for the variance and asymptotic variance of (cid:98)m4. Theorem 0.0.25. Let Xt be a linear process of the form: Xt = φiξt−i ∞(cid:88) i=−∞ 60 where ξt = IID(0, σ2), the coefficients satisfy(cid:80)∞ i=−∞|φi| < ∞ and E(ξ8 ∞(cid:88) T Var((cid:98)m4) → γ4(r) r=−∞ t ) < ∞, then where γ4(r) = Cov(X4 t , X4 t+r) = φ4r(m8 − m2 4) + 6τ2(1 − φ2r) 1 − φ2 φ2r(m6 − m4m2) + 4φm5τ3(1 − φ3r) 1 − φ3 mi = E(Xi t ) for i = 2, 4, 5, 6, 8, and τi = E(ξi) for i = 2, 3, 4, 5, 6, 8. Proof. Var((cid:99)m4) = = 1 T 2 1 T 2 T(cid:88) T(cid:88) T−1(cid:88) j=1 i=1 r=−(T−1) Cov(X4 i , X4 j ), (T − r)Cov(X4 1 , X4 1+r) 61 By repeating the same argument in proof of Theorem 0.0.21 we get: Cov(X4 1 , X4 1+r) = Cov(X4 1 , φ4rX4 1 + 4φ3rX3 1 Rr + 6φ2rX2 1 R2 r + 4φrX1R3 r + R4 r) = φ4rVar(X4 1 ) + 4φ3rCov(X4 1 , X3 1 Rr) + 6φ2rCov(X4 1 , X2 1 R2 r) + 4φrCov(X4 = φ4r(m8 − m2 + 6φ2r(m6E[R2 1 , X1R3 r) 4) + 4φ3r(m7E[Rr] − m4m3E[Rr]) r] − m4m2E[R2 r]) + 4φr(m5E[R3 r] − m4m1E[R3 r]) Let γ4(r) = Cov(X4 t , X4 t+r), then T Var((cid:99)m4) = |φi| < ∞, this implies(cid:80) T−1(cid:88) (cid:0)1 − |r| (cid:1)γ4(r), T r=−(T−1) i r T |γ4(r)| < ∞. Note that: (cid:12)(cid:12)(cid:0)1 − |r| (cid:1)γ4(r) → γ4(r), and take fn(r) =(cid:0)1− |r| (cid:1)γ4(r)(cid:12)(cid:12) < |γ4(r)| and (cid:1)γ4(r), for |r| < n and zero for |r| > n. Let and since(cid:80) (cid:0)1− |r| Theorem we have T Var((cid:98)m4) →(cid:80)∞ Now the following theorem proves the asymptotic distribution for (cid:98)m4. Note that the µ(r) = 1 for r = ±1,±2, ... be the counting measure then by the Dominated Convergence r=−∞ γ4(r). n T proof will be very similar to the proof of Theorem0.0.22, where the same method and theo- rems are used, but suitable approximations and estimators. Theorem 0.0.26. If Xt is a stationary linear process and under the assumption of Theorem 62 0.0.25, then: t ) and V =(cid:80)∞ Where m4 = E(x4 r=−∞ γ4(r). (cid:98)m4 ∼ AN (m4, V /T ) Proof. First define the stationary 2m-dependent process ym t = (Xm t )4, where Xm t is the approximation defined in (17) to Xt, then: ¯ymT = T−1 ym t T(cid:88) = T−1(cid:104) t=1 n )4(cid:105) (Xm 1 )4 + ... + (Xm also E(ym t ) = E(Xm t )4 is the fourth theoretical moment of the 2m dependent process Xm t . Then consider as an approximation to ymT = T 1/2 [¯ymT − E(¯ymT )] yT = T 1/2((cid:98)m4 − m4), where E(¯ymT ) is the same as E(ym t ) given above. To obtain an asymptotic distribution for yT , we apply Theorem 0.0.17, using ymT as our approximation.Now verifying (i),(ii) and (iii) of Theorem 0.0.17. (i) Applying Theorem 0.0.16, (CLT for m-dependent sequence), to the 2m dependent series 63 ymT we obtain: as T → ∞ where Vm =(cid:80)m ymT = T 1/2 [¯ymT − E(¯ymT )] d−→ ym ∼ N (0, Vm) r=−m γ4(r). (ii) Note that as m → ∞, Vm → V , using the Domi- nated Convergence Theorem, where V is defined in Theorem 0.0.25. Hence the characteristic function of ym, ϕm(λ) = exp{−1/2λ2Vm} → exp{−1/2λ2V } as m → ∞, which is the characteristic function of a random variable y ∼ N (0, V ) (character- izing convergence of distribution of Fn(.) in terms of convergence of sequence of characteristic functions ϕn(.), i.e., ϕn(λ) → ϕ(λ) ⇔ Fn(x) d−→ F (x)). (iii)To verify the last condition, note that: E(yT − ymT )2 = T E (cid:104) (cid:105)2 ((cid:98)m4 − E(X4 t )) − (¯ymT − E(¯ymT )) = T E(cid:0)(cid:98)m4 − E(X4 t )(cid:1)2 + T E(cid:0)¯ymT − E(¯ymT )(cid:1)2 − 2T E(cid:0)((cid:98)m4 − E(X4 t ))(¯ymT − E(¯ymn))(cid:1) = T Var((cid:98)m4) + T Var(¯ymT ) − 2T Cov((cid:98)m4, ¯ymT ) using the calculations that led to the result in Theorem 0.0.23, and letting m → ∞, we get hence T Var((cid:98)m4) + T Var((cid:98)m4) − 2T Var((cid:98)m4), (V + V − 2V ) = 0 64 as T → ∞. We state the following lemmas which contain the asymptotic covariances of ((cid:98)m2,(cid:98)m3) and ((cid:98)m2,(cid:98)m4), without proof, since the idea of the proofs is the same as of the previous theorems 0.0.21 , 0.0.23 and 0.0.25. Lemma 0.0.27. Let Xt be a linear process of the form: Xt = ∞(cid:88) where ξt = IID(0, σ2), the coefficients satisfy(cid:80)∞ ∞(cid:88) T Cov[(cid:98)m3,(cid:98)m2] → i=−∞ φiξt−i i=−∞|φi| < ∞, and E(ξ5 t ) < ∞, then ∞(cid:88) r=−∞ γ23(r) + γ32(r) r=−∞ where and γ23(r) = Cov(X2 t , X3 t+r) = φ3r(m5 − m3m2) + 3τ2m3 φr(1 − φ2r) 1 − φ2 γ32(r) = Cov(X3 t , X2 t+r) = φ2r(m5 − m3m2) 65 Lemma 0.0.28. Let Xt be a linear process of the form: Xt = ∞(cid:88) where ξt = IID(0, σ2), the coefficients satisfy(cid:80)∞ ∞(cid:88) T Cov[(cid:98)m4,(cid:98)m2] → i=−∞ φiξt−i i=−∞|φi| < ∞, and E(ξ6 t ) < ∞, then ∞(cid:88) r=−∞ γ24(r) + γ42(r) r=−∞ where and γ24(r) = Cov(X2 t , X4 t+r) = φ4r(m6 − m4m2) + 6τ2(m4 − m2 2) (1 − φ2r)φ2r 1 − φ2 + 4τ3m3 (1 − φ3r)φr 1 − φ3 γ42(r) = Cov(X4 t , X2 t+r) = φ2r(m6 − m4m2) The asymptotic normality of the empirical skewness and kurtosis In this section, the delta method was applied to get the asymptotic distribution of the em- pirical skewness and kurtosis of a stationary linear model to the asymptotic multivariate 66 distribution of the vectors MT = ((cid:98)m2,(cid:98)m3) and MT = ((cid:98)m4,(cid:98)m2) . We start by a conjecture for the asymptotic multivariate distribution of the vector MT = ((cid:98)m2,(cid:98)m3). 0.0.23, then: For MT = ((cid:98)m2,(cid:98)m3), we have Theorem 0.0.29. If Xt is a stationary linear process and under the assumptions of Theorem √ T (MT − µ) d =⇒ N (0, Σ) (18) where µ = [m2, m3](cid:62) r=−∞ r=−∞ ∞(cid:88) r=−∞ γ2(r) γ23(r) + γ3(r) γ32(r) ∞(cid:88) ∞(cid:88) ∞(cid:88) r=−∞ Σ11 = Σ12 = Σ22 = Now applying the delta method we get the asymptotic distribution for the empirical skewness of the linear process Xt in the following theorem. 67 Theorem 0.0.30. For ST = (cid:98)m3/(cid:98)m (cid:16) √ 3/2 2 , T ST − m3 3/2 m 2 (cid:17) d =⇒ N (0, σ2) where σ2 is given by σ2 = = (cid:104) − 3m3 9m2 5/2 2m 2 3Σ11 4m5 2 (cid:105) , 1 3/2 m 2 (cid:104) − 3m3 Σ 5/2 2m 2 (cid:105)(cid:62) , 1 3/2 m 2 − 6m3Σ12 2m4 2 + Σ22 m3 2 Proof. Applying theorem 0.0.29 and the multivariate delta method, where the function g is of the form g(s1, s2) = , we get the required. s2 3/2 s 1 The same theory of the empirical skewness is applied to the empirical kurtosis. We start by the following conjecture for the asymptotic normality of the vector ((cid:98)m2,(cid:98)m4). 0.0.25, then: For MT = ((cid:98)m2,(cid:98)m4), we have Theorem 0.0.31. If Xt is a stationary linear process and under the assumptions of Theorem √ T (MT − µ) d =⇒ N (0, Σ) (19) 68 where µ = [m2, m4](cid:62) ∞(cid:88) ∞(cid:88) ∞(cid:88) r=−∞ r=−∞ r=−∞ Σ11 = Σ12 = Σ22 = γ2(r) γ24(r) + ∞(cid:88) r=−∞ γ42(r) γ4(r) And lastly the distribution of KT Theorem 0.0.32. For KT = (cid:98)m4/(cid:98)m2 (cid:16) 2, √ T KT − m4 m2 2 (cid:17) d =⇒ N (0, σ2) where σ2 is given by (cid:104) − 2m4 m3 2 4m2 4Σ11 m6 2 σ2 = = (cid:105)(cid:62) (cid:105) (cid:104) − 2m4 , Σ 1 m2 2 − 4m4Σ12 m5 2 , m3 2 Σ22 m4 2 + 1 m2 2 Proof. Applying theorem 0.0.31 and the multivariate delta method, where the function g is of the form g(s1, s2) = s2 s2 1 , we get the required. 69 The non-central chi-square tests In this section, the empirical skewness and kurtosis are used to define a non-central chi- square distributions S and K, to reject the null hypotheses of Gaussian innovations in favor of Gumbel, Gumbel innovations and Gussian in favor of mixture of both and Gumbel innovations in favor of the Gaussian. Theorem 0.0.33. Let Xt be an AR(1) process of the form: Xt+1 = φXt + ξt+1 where ξt ∼ W N (0, σ2), with E(ξt)8 < ∞ and let S = (cid:98)m3/(cid:98)m 3/2 2 and K = (cid:98)m4/(cid:98)m2 2, then for large T , ˆ the hypothesis test: has a rejection region H0 : {ξt} follow Gaussian distribution Ha : {ξt} follow Gumbel distribution {xt : (S + bS)2 ≤ c} if (1/σ2 Sa − 1/σ2 S0 ) > 0 and {xt : (S + bS)2 > c} if (1/σ2 Sa − 1/σ2 S0 ) < 0 (20) (21) where σ2 S0 and σ2 Sa are the variances of S under the null and alternative hypotheses, 70 respectively, bS = (cid:16) µS0 (cid:17) σ2 Sa σ2 S0 −µSa σ2 S0 −σ2 Sa and c is a constant. Moreover, under the null hypothesis assumption, the test statistic (22) (23) (24) (S + bS)2/σ2 S0 ∼ χ2 1( µS0 + bS σS0 )2. H0 : {ξt} follow Gumbel distribution Ha : {ξt} follow Gaussian distribution ˆ the hypothesis test: has a rejection region {xt : (K + bK )2 ≤ c} if (1/σ2 Ka − 1/σ2 K0 ) > 0 and where σ2 K0 {xt : (K + bK )2 > c} (cid:17) (cid:16) µK0 −µKaσ2 K0 −σ2 Ka σ2 Ka σ2 K0 if (1/σ2 Ka − 1/σ2 K0 ) < 0 and σ2 Ka are the variances of K under the null and alternative hypotheses, respectively, bK = and c is a constant. Moreover, under the null hypothesis assumption, the test statistic (K + bK )2/σ2 K0 ∼ χ2 1( µK0 +bK σK0 )2. Proof. ˆ Let f0(s) and fa(s) be the density functions of the statistic S under the null and alternative hypotheses, respectively. Applying the Neyman-Pearson Lemma to S, 71 we reject H0 if ≤ c(cid:48), f0(s) fa(s) where c(cid:48) is a constant determined by α (the probability of type I error). Recall that, S ∼ AN (µS, σ2 tions follow Gaussian distribution) implies S ∼ AN (µS0 S), which under the null hypothesis assumption (the innova- ), and S ∼ AN (µSa, σ2 ) Sa , σ2 S0 under the alternative (the innovations follow Gumbel distribution). This implies(cid:32) 1(cid:113) 2πσ2 S0 −(s−µS0 e )2/2σ2 S0 (cid:44) 1(cid:113) 2πσ2 Sa (cid:33) ≤ c(cid:48), −(s−µSa )2/2σ2 Sa e hence, (cid:16) µS0 σ2 S0 − µSa σ2 Sa (cid:17) ≤ c(cid:48)(cid:48) (25) + 2s where c(cid:48)(cid:48) is a constant. If > 0 then (cid:17) (cid:17) σ2 Sa σ2 Sa − 1 σ2 S0 − 1 σ2 S0 s2(cid:16) 1 (cid:16) 1 (cid:16) s +(cid:0) µS0 (cid:17) , then (cid:1)(cid:17)2 ≤ c, − µSaσ2 S0 − σ2 Sa σ2 Sa σ2 S0 P(cid:16)(cid:0)s + bS (cid:1)2 ≤ c|H0 (cid:17) ≤ α 72 and c is a constant. (cid:16) µS0 σ2 Sa σ2 S0 −µSaσ2 S0 −σ2 Sa Let bS = where under the null hypothesis (cid:16) 1 σ2 Sa (cid:17) − 1 σ2 S0 And If which implies (cid:0)S + bS (cid:1)2|H0 ∼ χ2 1 (cid:0) µS0 + bS σS0 (cid:1). < 0, in (25), then the rejection region will be (cid:16) σ2 Sa σ2 S0 s +(cid:0) µS0 P(cid:16)(cid:0)s + bS (cid:1)(cid:17)2 > c, − µSaσ2 S0 − σ2 Sa (cid:1)2 > c|H0 (cid:17) ≤ α. ˆ Repeating the same logic in the previous part applied on the statistic K. Let f0(k) and fa(k) be the density functions of the statistic K under the null and alternative hypotheses, respectively. Applying the Neyman-Pearson Lemma to K, we reject H0 if ≤ c(cid:48), f0(k) fa(k) where c(cid:48) is a constant determined by α (the probability of type I error). Recall that, K ∼ AN (µK , σ2 vations follow Gumbel distribution) implies K ∼ AN (µK0 K ), which under the null hypothesis assumption (the inno- ), and K ∼ AN (µKa, σ2 , σ2 K0 ) Ka under the alternative (the innovations follow Gaussian distribution). 73 This implies(cid:32) 1(cid:113) 2πσ2 K0 (cid:44) −(k−µK0 e )2/2σ2 K0 1(cid:113) 2πσ2 Ka (cid:33) ≤ c(cid:48), −(k−µKa)2/2σ2 Ka e hence, where c(cid:48)(cid:48) is a constant. If (cid:17) + 2k (cid:17) σ2 Ka σ2 Ka − 1 σ2 K0 − 1 σ2 K0 k2(cid:16) 1 (cid:16) 1 (cid:16) k +(cid:0) µS0 (cid:17) P(cid:16)(cid:0)k + bK σ2 Ka σ2 K0 , then σ2 Ka σ2 K0 −µKa σ2 K0 −σ2 Ka and c is a constant. (cid:16) µK0 Let bK = (cid:16)µK0 σ2 K0 − µKa σ2 Ka (cid:17) ≤ c(cid:48)(cid:48) (26) > 0 then − µKaσ2 K0 − σ2 Ka (cid:1)(cid:17)2 ≤ c, (cid:1)2 ≤ c|H0 (cid:17) ≤ α where under the null hypothesis (cid:0)K + bK (cid:1)2|H0 ∼ χ2 1 (cid:0) µK0 + bK σK0 (cid:1). 74 (cid:16) 1 σ2 Ka (cid:17) − 1 σ2 K0 And If which implies < 0, in (26), then the rejection region will be (cid:16) σ2 Ka σ2 K0 k +(cid:0) µK0 P(cid:16)(cid:0)k + bK (cid:1)(cid:17)2 > c, − µKaσ2 K0 − σ2 Ka (cid:1)2 > c|H0 (cid:17) ≤ α. Here are some remarks regarding the non-central chisquare S and the non-central chisquare K tests: Remark 0.0.34. The non-central chi-square S test can be used to reject the Gumbel innova- tions in favor of the Gaussian innovations, and it requires less data. Remark 0.0.35. The non-central chi-square K test can be used to reject the Gaussian innova- tions, in favor of the Gumbel innovations, but it needs larger sample size than the non-central chisquare S test. Remark 0.0.36. For the mixture of the Gaussian and Gumbel innovations, the non-central chisqure test S can be used to reject the Gumbel or the Gaussian innovations in favor of the mixture of both. Note that the S test can not be used to reject the mixture of both in favor of the Gaussian or in favor of the Gumbel, since any Gaussian AR(1) model or AR(1) model with Gumbel innovations can be written as a mixture with appropriate wights. While the non-central chisquare K test can not be used with the mixture. Remark 0.0.37. In practice, starting with non-central chisquare S would be advised, that is 75 checking the normality first, and then checking the mixture assumption. That is because if there is a moderate portion of the observations are Gumbel (like the mixture), then applying the non-central chi-square S or K tests may not reject the Gumbel assumption in favor for the Gaussian since it is heavier tail than both of them. Remark 0.0.38. The empirical power for these tests are computed, for more details see the simulation studies in chapter 4. Remark 0.0.39. A very useful well known test is, the usual z-test, for example, by finding the confidence interval and use it to conclude the rejection of any of the null hypotheses the innovations are Gaussian, Gumbel or a mixture of both. 76 The Simulation This chapter contains simulations for the main results in chapter 3. We start by simulations for the asymptotic distribution of the ST and KT . And then the non-central chi-square S and K tests were illustrated by simulations for all the different null hypothesis scenarios. In addition, their empirical power was computed using simulations. The distributions of ST and KT In this section, we compare the asymptotic distributions of ST and Kt using our calculations and using ”the empirical moments”. In table 2 we calculate the asymptotic distribution of ST where we simulate from a Gaussian AR(1) model. T is the length of the realization, ”the sample size”, and φ is the AR(1) parameter. In the third column, we simulate an AR(1) model with the same length as T and φ number of iterations (say 1000 times for example). each time we calculate the empirical moments, and then take their average over the number of iterations, and then we get the variances and covariances of the empirical second and third moments, and finally apply the delta method to the result which appear in column 3. In the last column, we first simulate one realization of the Gaussian AR(1), with the same T and φ and then obtain the asymptotic distributions of (cid:98)m2 and (cid:98)m3 using the results as in theorem0.0.29 and then applying the delta method. 77 T φ Distribution of ST using empirical moments Distribution of ST using the moments of Xt 100 200 500 800 0.35 0.59 0.599 0.864 1600 0.850 N (0, 0.502) N (0, 0.472) N (0, 0.302) N (0, 0.442) N (0, 0.292) N (0, 0.512) N (0, 0.472) N (0, 0.302) N (0, 0.432) N (0, 0.292) Table 2: The distribution of ST using empirical moments and the moments of Xt when simulating from Gaussian AR(1) model As we can see from the table the distributions are almost the same. When the innovations are Gumbel, applying the same techniques above yield to the results in table 3. The results are satisfying in both cases. n φ Distribution of ST using empirical moments Distribution of ST using the moments of Xt 100 200 800 0.31 0.662 0.856 1000 0.863 1600 0.857 N (0.966, 0.692) N (0.667, 0.592) N (0.419, 0.432) N (0.408, 0.402) N (0.420, 0.312) N (1.00, 0.692) N (0.674, 0.602) N (0.4216, 0.442) N (0.410, 0.402) N (0.420, 0.312) Table 3: The distribution of ST using empirical moments and using the moments of Xt when simulating from AR(1) model with Gumbel innovations The same method in comparing the asymptotic distribution of ST using the results in chapter 3 with the empirical moments, are also applied to KT . In this case the empirical second and fourth moments are used, and theorem 0.0.31. Table 4 shows the distribution 78 of KT where the innovations are Gaussian, while for the Gumbel innovations the results are illustrated in table 5. All of these distributions are in agreement with each other. T φ Distribution of KT using empirical moments Distribution of KT using the moments of Xt 100 0.149 200 0.47 500 0.599 800 0.864 300 0.786 100 0.54 N (3.00, 0.482) N (3.00, 0.3632) N (2.99, 0.252) N (2.99, 0.2722) N (3.00, 0.4122) N (3.00, 0.5322) N (2.99, 0.492) N (2.99, 0.3652) N (2.99, 0.252) N (2.99, 0.2722) N (2.99, 0.4222) N (2.99, 0.5362) Table 4: The distribution of KT using empirical moments and the moments of Xt when simulating from Gaussian AR(1) model T φ Distribution of KT using empirical moments Distribution of KT using the moments of Xt 100 0.263 200 0.47 500 0.597 800 0.864 300 0.811 100 0.522 N (5.06, 3.682) N (4.53, 2.122) N (4.13, 1.122) N (3.367, 0.6792) N (3.494, 1.1492) N (4.33, 2.612) N (5.10, 3.742) N (4.47, 2.102) N (4.14, 1.172) N (3.36, 0.7352) N (3.496, 1.222) N (4.37, 2.82) Table 5: The distribution of KT using empirical moments and the moments of Xt when simulating from AR(1) model with Gumbel innovations 79 The non-central Chi-square test S and its empirical power In this section we show the results of applying the non-central chi-square test S to reject the Gaussian innovations in favor of the Gumbel and the mixture of both innovations of AR(1) model, and to reject the Gumbel in favor of the Gaussian and the mixture of both mixture innovations of the Ar(1) model. Also we calculate the empirical power in all hypotheses. We start with the hypothesis H0 : ξt ∼ Gaussian vs Ha : ξt ∼ Gumbel. We simulate AR(1) with Gumbel innovations where T and φ as in the first and second column, and then we fit an AR(1) model with Gaussian innovations and obtain the estimators which are used to get the distribution of ST as in theorem0.0.30. And then on the same realization, we fit a Gaussian AR(1) model and get the estimators and get the distribution of ST under the Gaussian assumption. Now the non-central chi square test S is applied as in theorem 0.0.33. As we can see the test works perfectly, we reject the null hypothesis in all of the simulations. Then when we simulate from Gaussian AR(1) with the same T and φ, and apply the same procedure above we fail to reject the null hypothesis, that is the innovations are Gaussian. 80 simulation of AR(1) simulation of AR(1) with Gumbel innova- with Gaussian innova- tions H0 : ξt ∼ Gaussian Ha : ξt ∼ Gumbel tions H0 : ξt ∼ Gaussian Ha : ξt ∼ Gumbel T 250 250 300 300 400 500 700 900 900 φ The result The result 0.168 Reject H0 Fail to reject H0 0.256 Reject H0 Fail to reject H0 0.351 Reject H0 Fail to reject H0 0.455 Reject H0 Fail to reject H0 0.519 Reject H0 Fail to reject H0 0.58 Reject H0 Fail to reject H0 0.659 Reject H0 Fail to reject H0 0.720 Reject H0 Fail to reject H0 0.783 Reject H0 Fail to reject H0 1000 0.801 Reject H0 Fail to reject H0 Table 6: Applying the non-central chi-square S test to reject the Gaussian innovations in favor of the Gumbel. Table 7 shows the empirical power of the non-central chi-square S test of the previous hypothesis, that is the hypothesis H0 : ξt ∼ Gaussian vs Ha : ξt ∼ Gumbel. As we can see the power change T and φ. The bigger the |φ| the larger T is needed to get a good power. 81 φ The empirical power T 200 200 500 200 500 300 500 800 0.17 0.30 0,30 0.50 0.50 0.60 0.60 0.60 899/1000 790/1000 982/1000 557/1000 894/1000 550/1000 743/1000 887/1000 925/1000 846/1000 1000 0.60 1200 0.70 Table 7: The empirical power of the non-central chi-square test S where H0 : ξt ∼ Gaussian vs Ha : ξt ∼ Gumbel The second scenario the hypothesis H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian. We simulated from AR(1) with Gumbel innovations and fit first Gumbel innovations and get the parameters using EM algorithm and get the distribution of ST . WE repeated the same techniques after fitting a Gaussian AR(1) model to the sample path. We applied the non- central chi-square S test to the data. 82 And we get the results in table8, with the empirical power in table 9. simulation of AR(1) with Gumbel innova- tions H0 : ξt ∼ Gaussian Ha : ξt ∼ Gumbel The result simulation of AR(1) with Gumbel innova- tions H0 : ξt ∼ Gumbel Ha : ξt ∼ Gaussian The result φ 0.17 Reject H0 Fail to reject H0 0.0.291 Reject H0 Fail to reject H0 0.529 Reject H0 Fail to reject H0 T 500 500 700 1200 0.728 Reject H0 Fail to reject H0 1800 0.798 Reject H0 Fail to reject H0 Table 8: Applying the non-central chi-square S test to the hypothesis H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian 83 And this is the empirical power. The last case where we apply the non-central chi square φ The empirical power T 500 500 700 700 0.17 0.30 0.50 0.60 999/1000 992/1000 952/1000 840/1000 789/1000 752/1000 842/1000 1200 0.70 2500 3000 0.8 0.8 Table 9: The empirical power of the non-central chi-square test S where H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian test S is to reject the Gumbel or Gaussian innovations in favor of the mixture. That is H0 : ξt ∼ Gaussian or H0 : ξt ∼ Gumbel vs Ha : ξt ∼ mixture. The same procedure above was applied here. We first simulate from AR(1) with mixture innovations, and then fit an AR(1) model with mixture, then find the estimators which enable us to calculate the distribution of ST as in theorem 0.0.30. On the same realization we got from the simulation with mixture innovations, we fit an AR(1) with Gumbel innovations and then find the estimators which we use to find the distribution of ST as in theorem 0.0.30 under the Gumbel assumption. Doing the same procedure on the same realization but with Gaussian innovations. 84 φm p (cid:98)φg (cid:98)φn H0 : ξt ∼ Gaussian Ha : ξt ∼ M ixture The result T 0.109 0.075 0.096 0.711 Reject H0 700 0.150 0.137 0.133 0.657 Fail to reject H0 700 0.138 0.142 0.134 0.501 Reject H0 700 0.186 0.207 0.192 0.685 Fail to reject H0 700 0.209 0.193 0.197 0.685 Reject H0 700 0.138 0.159 0.174 0.557 Reject H0 700 700 0.235 0.229 0.229 0.260 Reject H0 2000 0.114 0.137 0.131 0.324 Reject H0 2000 0.185 0.175 0.176 0.255 Reject H0 2000 0.158 0.173 0.173 0.843 Fail to reject H0 2000 0.186 0.145 0.161 0.815 Reject H0 2000 0.271 0.268 0.267 0.364 Reject H0 0.209 0.197 0.205 0.546 Reject H0 H0 : ξt ∼ Gumbel Ha : ξt ∼ M ixture The result Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Fail to reject H0 Reject H0 Reject H0 2000 2000 0.276 0.266 0.273 0.754 Fail to reject H0 Reject H0 Table 10: Applying the non-central Chisquare test S when simulating from an AR(1) with mixture innovation to reject the Gumbel and the Gaussian innovations in favor of the mix- ture. After obtaining the distribution of ST under the three assumptions where we simulate from a mixture we applied the non-central chi square test S as in theorem 0.0.33 and the results are in the following tables 10 and 11. 85 The table is continued here: φm p (cid:98)φg (cid:98)φn H0 : ξt ∼ Gaussian Ha : ξt ∼ M ixture The result T 3000 0.380 0.393 0.383 0.299 Reject H0 3000 0.347 0.340 0.342 0.424 Reject H0 3000 0.383 0.379 0.382 0.411 Reject H0 3000 0.308 0.321 0.327 0.592 Reject H0 3000 0.333 0.338 0.338 0.703 Reject H0 4000 0.445 0.439 0.442 0.254 Reject H0 4000 0.442 0.444 0.448 0.544 Reject H0 4000 0.540 0.542 0.539 0.508 Reject H0 4000 0.531 0.533 0.538 0.745 Fail to reject H0 9000 0.652 0.652 0.651 0.353 Reject H0 9000 0.656 0.654 0.655 0.512 Reject H0 9000 0.740 0.737 0.739 0.690 Fail to reject H0 9000 0.663 0.662 0.665 0.769 Reject H0 H0 : ξt ∼ Gumbel Ha : ξt ∼ M ixture The result Fail to reject H0 Fail to reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Reject H0 Fail to reject H0 Reject H0 Reject H0 Reject H0 Table 11: Applying the non-central Chisquare test S when simulating from an AR(1) with mixture innovation to reject the Gumbel and the Gaussian innovations in favor of the mixture innovations continuing. 86 As we see from table 10, the performance of the test depends on many factors, for example the φ, T and the p. Most of the cases where the test fails are when p is large or small and also related to the φ. But in general the test preforms better when we increase the sample size. As we see from the empirical power in table 12. T φ p The empirical power 700 0.17 0.53 1400 0.17 0.53 1800 0.17 0.53 2200 0.17 0.53 3000 0.17 0.53 1200 0.30 0.53 2500 0.30 0.53 3500 0.30 0.53 3500 0.50 0.53 4500 0.50 0.53 6000 0.70 0.53 61/100 76/100 80/100 86/100 91/100 65/100 80/100 88/100 71/100 85/100 63/100 Table 12: The empirical power of the non-central chi-square test S where H0 : ξt ∼ Gumbel or Gaussian vs Ha : ξt ∼ mixture The non-central chi-square test K and its empirical power In this section we show the results of applying the non-central chi-square K test. There are two cases where this test can be applied, to reject the Gumbel innovations in favor of Gaussian and the other way, to reject the Gaussian innovations in favor of the Gumbel. 87 We start by the hypothesis H0; ξt ∼ Gumbel vs H0; ξt ∼ Gaussian, where we simulate from a Gaussian AR(1), and then after obtaining the sample we fit a Gaussian AR(1) and use the estimators to get the distribution of Kt by applying theorem 0.0.32 . Then next, we fit an AR(1) with Gumbel innovations on the same sample and use the estimators to get the distribution of KT as in theorem 0.0.32. After having these distributions we apply theorem 0.0.33 to obtain the results of K test. 88 simulation of AR(1) with Gaussian innova- tions H0 : ξt ∼ Gumbel Ha : ξt ∼ Gaussian The result simulation of AR(1) with Gumbel innova- tions H0 : ξt ∼ Gumbel Ha : ξt ∼ Gaussian The result φ 0.171 Reject H0 Fail to reject H0 0.266 Reject H0 Fail to reject H0 0.33 Reject H0 Fail to reject H0 T 700 850 900 1000 0.483 Reject H0 Fail to reject H0 1500 0.56 Reject H0 Fail to reject H0 2000 0.508 Reject H0 Fail to reject H0 2500 0.60 Reject H0 Fail to reject H0 2500 0.762 Reject H0 Fail to reject H0 3000 0.77 Reject H0 Fail to reject H0 3000 0.84 Reject H0 Fail to reject H0 Table 13: Applying the non-central chi-square K test to reject the Gumbel innovations in favor of the Gaussian. Lastly, we simulate from AR(1) with Gumbel innovations and do the same procedure above, and apply the same test under the same assumptions. 89 As we can see in table 13 the test preforms well, where we reject the null hypothesis when we generate from Gaussian and we fail to reject the null hypothesis when the innovation is generated from the Gumbel. The following table 14 shows the empirical power of the test. φ The empirical power T 500 700 700 0.17 0.17 0,30 770/1000 908/1000 853/1000 973/1000 812/1000 926/1000 769/1000 911/1000 627/1000 817/1000 897/1000 1000 0.30 1000 0.50 1300 0.50 1300 0.60 1800 0.60 1800 0.70 2500 0.70 3000 0.70 Table 14: The empirical power of the non-central chi-square test K where H0 : ξt ∼ Gumbel vs Ha : ξt ∼ Gaussian 90 simulation of AR(1) with Gumbel innova- tions H0 : ξt ∼ Gaussian Ha : ξt ∼ Gumbel The result simulation of AR(1) with Gaussian innova- tions H0 : ξt ∼ Gaussian Ha : ξt ∼ Gumbel The result φ 0.171 Reject H0 Fail to reject H0 0.266 Reject H0 Reject H0 0.33 Reject H0 Fail to reject H0 T 700 850 900 1000 0.483 Reject H0 Fail to reject H0 1500 0.56 Reject H0 Fail to reject H0 2000 0.508 Reject H0 Fail to reject H0 2500 0.60 Reject H0 Fail to reject H0 2500 0.762 Reject H0 Fail to reject H0 3000 0.77 Reject H0 Fail to reject H0 3000 0.84 Reject H0 Fail to reject H0 Table 15: Applying the non-central chi-square K test to reject the Gaussian innovations in favor of the Gumbel. Lastly, we applied the non-central chi-square test K under the hypothesis H0 : ξ ∼ Gaussian vs Ha : ξt ∼ Gumbel to the same simulations above and we obtained the results in table15. 91 And as we see the results of applying this test are satisfactory. And the empirical power of the of test K under the previous null hypothesis is shown in table 16. φ The empirical power 870/1000 917/1000 836/1000 904/1000 757/1000 898/1000 920/1000 800/1000 885/1000 908/1000 725/1000 851/1000 930/1000 T 300 350 300 400 400 600 700 700 900 0.17 0.17 0,30 0.30 0.50 0.50 0.5 0.60 0.60 1000 0.60 1000 0.70 1400 0.70 1900 0.70 Table 16: The empirical power of the non-central chi-square test K where H0 : ξt ∼ Gaussian vs Ha : ξt ∼ Gumbel 92 The analysis of the sea level rise data It has been predicted that the global sea-level rise (SLR) will reach 0.5 to 2.00m by the end of this century (by 2100). The relative SLR (RSLR) can be greater than the global change because of local effects [57]. As mentioned in chapter one, several papers provide evidence of SLR using tide gauge data ”historical observations”, along the coastal mid-Atlantic region (a so-called ”hotspot”) [31, 10, 74]. This chapter analyzes sea level measurements obtained from the Actuarial climate index (ACI). These measurements are made available via tide gauge, which ”measure sea level rel- ative to the land below, but because the land may be moving, the ACI sea level component measures the combined effect on coastal shorelines of the generally rising seas and the rising or falling land,” According to ACI. The ACI sea level measurements are obtained from tide gague stations, which are located along 10 of the 12 regions as shown in Figure 3. These are for all regions except the CAR ”Central Arctic” and this MID ”Midwest.” This study is concerned with sea level measure- ments of the CEA ”Central East Atlantic” region, which includes the following states and districts: CT, DC, DE, MA, MD, ME, NH, NJ, NY, PA, RI, VT and WV. The population, infrastructure and economic activities in this region are dense. All this has led many au- thors to study this area as mentioned before, and our analysis shows new evidence of a rise 93 in ”hotspots” at tide stations in the CEA. There is more about the data in section one, and section two includes the methods, followed by the results and discussion in section three and four respectively. Figure 3: Map shows the 12 different regions according to ACI Data Monthly sea level standardized anomaly records from tide gauge stations in the Central East Atlantic ”CEA” were obtained from the ACI, recorded over a length of 56 years (1961-2017). 94 The ACI acquired the raw data from the Permanent Service for Mean Sea Level. Let yi be the monthly sea level measurements from 1961-2017, and let 1961-1990 be the reference period, with average µref (yi) and standard deviation σref (yi) these measurements then convert to standardized anomaly according to the following equation: yistd = ∆yi/σref (yi) where ∆yi = yi − µref (yi). Method The data (standardized anomaly records) were decomposed after converting them to time series. The path of the series {Xt} (the random part after decomposition) suggests that an autoregressive model would be a good choice. Three centered AR(1) models with different innovations distributions were fitted to the data, Gaussian, Gumbel and a mixture of both. Recall AR(1) model is written in the form: Xt = φXt−1 + ξt where ξt ∼ W N (0, σ). For the centered Gaussian AR(1), ξ ∼ N (0.σ), σ > 0, the parameters φ and σ are estimated using the Expectation Maximization, where ξt is considered to be latent variable, and the 95 log-likelihood of the Gaussian random variable ξt is: log L(0, σ2|ξt) = −T 2 log(2πσ2) − 1 2σ2 T(cid:88) t=1 ξ2 t . And for the centered AR(1) model with Gumbel innovations the ξ ∼ G(α, β) where α = −γβ, β > 0 and γ is the Euler constant. As in the Gaussian case, φ and β are estimated using the Expectation Maximization where ξt is considered a latent variable, and the log-likelihood of the Gumbel random variable ξt is: log L(α, β|ξt) = −T log(α) − T(cid:88) t=1 − T(cid:88) t=1 ξt − β α e−( ξt−β α ). And for the centered AR(1) model with mixture innovations, the innovations are the sum of the previous weighted innovations, (ξt ∼ p· N (0, σ2) + (1− p)· G(−γβ, β)), where 0 ≤ p ≤ 1. The Expectation Maximization algorithm is used to estimate the parameters φ, σ2, p, and β , and the log likelihood function for the mixture innovations is: (cid:88) log L(σ, β, p, ξt) = log[(1 − p) · f (ξt) + p · g(ξt)], t such that f (ξt) and g(ξt) are the densities for the above-centered Gaussian and Gumbel r.vs. respectively. After obtaining all estimators, the theoretical AR(1) moments, E[Xr t ] where r = 2, 3, 4, 5, 6, 8 were calculated according to formulas Lemma 0.0.5, Lemma 0.0.6, Lemma 0.0.7, Lemma 96 0.0.8, Lemma 0.0.9 and Lemma 0.0.10, where under the Gaussian assumption, the estimated φ from fitting AR(1) with Gaussian innovations used in addition to moments of Gaussian and the sum of these weighted moments is the mixture moments. random variable. The same applied under the Gumbel assumption, and for the mixture assumption, in addition to the estimated φ, (cid:98)p, is used to calculate the weighted moments, Using these moments with (cid:98)φ, the means, variances and covariances of the empirical moments (cid:98)mr =(cid:80)T t /T , where r = 2, 3, 4 were obtained according to Theorem 0.0.18, t=1 Xr Theorem 0.0.11, Theorem0.0.12, Theorem0.0.14, Theorem0.0.13 and Theorem0.0.15. Using theoremTheorem 0.0.29 and Theorem 0.0.31 imply the asymptotic multivariate normality of the vectors ((cid:98)m2,(cid:98)m3) and ((cid:98)m2,(cid:98)m4) of these statistics, and then the delta method in Theorem skewness, ST = (cid:98)m3/(cid:98)m and empirical kurtosis KT = (cid:98)m4/(cid:98)m2 0.0.30 and Theorem0.0.32 are applied to calculate the means and variances of the empirical 2, under the three different 3/2 2 assumptions on the innovations. Having these distributions enables us to apply the different tests (as in chapter three) to determine the best model. The confidence intervals were obtained, and then the z-test was applied to the means of the two statistics ST and KT under the three assumptions on ξt. To apply the non-central chi square S test as in Theorem 0.0.33, the value of the test statistic as in (22) was calculated as well as the coefficients (20) and (21) to determine the rejection regions. In addition, the p value of the non-central chisquare test S was calculated under different hypotheses (e.g. P(value of test statistic)). 97 Results In this section, the results of applying the above method to the sea level rise are summarized. Table 18 shows the estimators of φ, σ, β and p when fitting AR(1) model with Gaussian, Gumbel and mixture innovations. ξ’s distribution Gaussian Gumbel Mixture (cid:98)φ 0.172 0.190 0.170 (cid:98)σ 0.863 - 0.8321 (cid:98)β - 0.7638 0.7063 (cid:98)p - - 0.534 Table 17: The estimators of (cid:98)φ,(cid:98)σ and(cid:98)p when fitting AR(1) model with different innovations, to the SLR data, The proceeding estimators are used to calculate the AR(1) theoretical moments, which are summarized in Table 18. These moments are essential in calculating the asymptotic distribution of the empirical skewness ST and kurtosis KT . ξ’s distribution m2 Gaussian Gumbel Mixture m3 0.00 1.139 m4 3.198 5.617 m5 m6 m8 0.0025 16.484 119.321 19.146 95.670 3218.199 1.032 1.034 1.0300 0.5293 4.2999 8.8560 52.9747 1548.285 Table 18: The theoretical moments of AR(1) model when fitted with different innovations to the SLR data. Using these moments with (cid:98)φ and the moments of the Gaussian, Gumbel and of the 98 mixture of both, we calculate the variances under the three different assumptions. Table 19 contains the asymptotic distribution of the empirical skewness (ST ) and kurtosis (KT ). ξ’s distribution The distribution of Gaussian Gumbel Mixture ST N (0.00, 0.1672) N (1.08, 0.262) N (0.50, 0.252) The distribution of KT N (2.99, 0.192) N (5.24, 1.522) N (4.05, 1.122) Table 19: The asymptotic distribution of ST and KT when fitting the AR(1) model with different innovations to the SLR data Since the distributions of KT and ST under the three assumptions are found, the confi- dence intervals for the means are calculated as shown in Table 20. ξ’s distribution 95% C.I. of E(ST ) 64% C.I. of E(KT ) Gaussian Gumbel Mixture (0.173, 0.831) (−0.019, 1.024) (0.004, 1.00) (3.55, 3.90) (2.30, 5.15) (2.67, 4.78) 3/2 Table 20: The confidence intervals of m3/m 2 and m4/m2 2 applied to the SLR data To apply the hypothesis tests as in Theorem 0.0.33, the values of test statistics are obtained under different assumptions. These values are summarized in Table 21 along with (1/σ2 Sa − 1/σ2 S0 ) values and the C’s values, where α is at a level 5%. In table 22 we calculate the p value of the non-central chisquare test S under different null hypothesis scenarios. 99 The hypothesis test The Test Statistic Value 445.5 28.49 491.059 H0 : ξt ∼ Gumbel Ha : ξt ∼ mixture H0 : ξt ∼ Gaussian Ha : ξt ∼ mixture H0 : ξt ∼ mixture Ha : ξt ∼ Gumbel H0 : ξt ∼ mixture Ha : ξt ∼ Gaussian Value of − 1/σ2 S0 (1/σ2 Sa 1.438 −19.94 −1.438 The C Value ) 468.1 15.932 567.47 12.46 19.94 3.620 Table 21: The results of applying the non-central chisquare test T to the SLR data under different null hypothesis scenarios. Discussion In this section we discuss the results of fitting an AR(1) model with three different inno- vations (Gaussian, Gumbel and mixture of both) to the SLR data, and we show the best model based on the non-central chi-square S test and compare the results to other diagnostic checks based on the residuals. Table 19 contains the distributions of the empirical skewness ST and the empirical kur- tosis KT , and their confidence intervals are given in table 20. When applying the z-test to E(ST ) under the three different assumptions (ξ follow Gaussian, Gumbel and a mixture of both), using 95% confidence intervals, we reject the null hypothesis that ξt follow the Gumbel distribution, and the null hypothesis that ξt follow the Gaussian. This is because the means in both cases are not included in the confidence intervals. And since the mean of 100 The hypothesis test The P value The result H0 : ξt ∼ Gumbel Ha : ξt ∼ mixture H0 : ξt ∼ Gaussian Ha : ξt ∼ mixture H0 : ξt ∼ mixture Ha : ξt ∼ Gumbel H0 : ξt ∼ mixture Ha : ξt ∼ Gaussian 0.01489 Reject H0 0.0013 0.5067 0.493 Reject H0 Fail to reject H0 Fail to reject H0 Table 22: The p values of the S test applied to the SLR data with different scenarios. ST under the assumption that ξt follows a mixture of both the Gaussian and the Gumbel is included in the 95% confidence interval, we fail to reject the null hypothesis that the inno- vation follows a mixture of the Gaussian and Gumbel. We get the same conclusions when applying the z-test to the mean of KT under the previous three assumptions. Note that the confidence level for the E(KT ) is 65% because the sample size is small. In order to narrow the confidence interval, either increase the sample size, which leads to a decrease in the bias, or reduce the confidence level. Table 21 includes the results of applying the S test to the SLR data under different hypotheses. The purpose of this test is to reject the innovations’ hypotheses, whether Gaus- sian or Gumbel, in favor of the mixture. And it fails to reject the null hypothesis that ξt follows mixture in favor of the Gaussian or Gumbel innovations. As we see the value of (1/σ2 Sa − 1/σ2 S0 ) under the assumption H0 : ξt follows Gumbel vs Ha : ξt following the 101 mixture, which is positive. This implies that we reject the null hypothesis since the value of the test statistic falls in the rejection region, (the value of the test statistic = 445.5 < the value of C= 468.1). For the hypothesis H0 : ξt follow Gaussian vs Ha : ξt follow mixture, the value of (1/σ2 Sa − 1/σ2 S0 ) is negative, this implies we reject the null hypothesis since the value of the test statistic is greater than the value of C (28.49 > of C= 15.932), that is, the value of the test statistic falls in the rejection region. Repeating the same logic above, we fail to reject the null hypothesis under the following assumption H0 : ξt follows mixture vs Ha : ξt follows Gumbel or Gaussian. That is because for the Gumbel alternative the value of the test statistic is 491.05 which is less than the value of C, (the negative value of (1/σ2 Sa − 1/σ2 S0 ) implies the rejection region has to be greater than C). And for the Gaussian alternative, the value of the test statistic is greater than the value of C, hence it falls outside the rejection region. The p value for the non-central chi square S test under the previous null hypothesis sce- narios is included in table 22. From the table we see the values of p under the null hypotheses H0 : ξt follows Gumbel or H0 : ξt follows Gaussian are less than 5% (the p values are sig- nificant), this implies we reject these null hypotheses. Whereas we fail to reject the null hypothesis that H0 : ξt follows a mixture of both, in favor of the Gaussian or the Gumbel. That is because the p values are greater than 5%. All these results are in agreement with each others, that is the AR(1) model with mixture of both Gaussian and Gumbel fits the data more appropriately than the AR(1) with the 102 Gumbel or with the Gaussian innovations. Figure 4: The QQ-plot of the residuals of fitting a Gaussian AR(1) model to the STR data vs the Gaussian distribution These results are also in agreement with other common diagnostic check based on the residuals (cid:98)ξt. For example, the QQ plot of the residuals of fitting a Gaussian AR(1) model to the SLR data is shown in figure 4. The plot suggests a heaver tail than the Gaussian. And the QQ plot of the residuals of AR(1) with Gumbel innovations is shown in figure 5, which suggests a lighter tail than the Gumbel. Whereas the QQ plot of the residuals of the mixture A(1) as shown in 6 shows a good fit. 103 −2−1012−3−2−10123Q−Q plot for Sea Level Rise data as AR(1) with Gaussian Innovations llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll The QQ-plots of the Gumbel and mixture residuals. Figure 5: The QQ-plot of the residuals of fitting an AR(1) model with Gumbel innovations to the STR data vs the Gumbel distribution Figure 6: The QQ-plot of the residuals of fitting an AR(1) model with mixture innovations to the STR data vs the mixture distribution 104 −2−1012345−20246Q−Q plot for Sea Level Rise data as AR(1) with Gumbel Innovations llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll−2−101234−20246Q−Q plot for Sea Level Rise data as AR(1) with Mixture Innovations llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll Another diagnostic check based on the residuals is the Kolmogorov–Smirnov test (KS test). Table 23 shows the p value of the KS test applied to the residuals of the SLR data under different null hypothesis scenarios. For the Gaussian and Gumbel residuals’, H0 : (cid:98)ξt follows the Gaussian or Gumbel vs H0 : (cid:98)ξt does not follow the Gaussian or the Gumbel. As we see from the table, we reject the null hypotheses that the residuals follow Gaussian or the Gumbel. For the mixture residuals, we simulate a random sample yt from a mixture with the same parameters we obtained from fitting a mixture AR(1), and apply the KS test to the hy- pothesis H0 : (cid:98)ξt and yt come from the same distribution vs Ha the two samples are from different distributions. The p-value of the KS statistic indicates that we fail to reject the null hypothesis. The hypothesis test The p value H0 : (cid:98)ξt ∼ Gaussian H0 : (cid:98)ξt ∼ Gumbel H0 : (cid:98)ξt and yt fol- 0.010 2.2e − 16 0.760 low the same distri- bution Table 23: The p value of KS test under the null hypothesis applied to the residuals of the SLR data with different scenarios. All the results based on our methodology and based on the residuals diagnostic checks are in agreement. The advantage of our methodology is that the non-central chisqure S test 105 enables us to reject the Gaussian and to reject the Gumbel in favor of the mixture. The limitation of this test is that it needs large sample. 106 Conclusion A summary of this work as follows: we defined non-central chi-square S and K tests to reject the hypothesis that the innovations of AR(1) model follow Gumbel in favor of the Gaussian innovations, to reject the Gaussian innovations in favor of the Gumbel and to reject the Gumbel or Gaussian innovations in favor of the mixture of both of them. we modeled the SLR data as an AR(1) process where the innovations are Gumbel, Gaus- sian, or a mixture of both. Based on our non-central chi square test, we reject the Gumbel and Gaussian innovations. This work opens further problems in both theory and applications. In the theory part, proving the asymptotic multivariate normality of the vector ((cid:98)m2,(cid:98)m3) by the m-dependency or any other tools. Also in the theory, looking at the general formula of the higher moments of the Gumbel distribution, it is a function of the ”Euler Mascheroni integral” which is the derivative (and the higher derivative) of the Gamma function evaluated at 1. And note that the general formula for the moments in the Gaussian case is a function of the Hermite polynomials. The question now does the ”Euler Mascheroni integral” play other role for the Gumbel distribution? For the applied part, estimating the parameters of the centered AR(1) model with mixture innovations, without assuming the means of the mixture distributions equal to zero, this might be studied by using Bayesian analysis, in particular employing 107 RUNJAGS. 108 BIBLIOGRAPHY 109 BIBLIOGRAPHY [1] J Andel. On ar (1) processes with exponential white noise. Communications in Statistics-Theory and Methods, 17(5):1481–1495, 1988. [2] OD Anderson. The box-jenkins approach to time series analysis. RAIRO-Operations Research, 11(1):3–29, 1977. [3] Ronald L Anderson. Distribution of the serial correlation coefficient. The Annals of Mathematical Statistics, 13(1):1–13, 1942. [4] TW Anderson and AM Walker. On the asymptotic distribution of the autocorrelations of a sample from a linear stochastic process. The Annals of Mathematical Statistics, 35(3):1296–1303, 1964. [5] ANNMARIE GEDDES BARIBEAU. The slr factor: As sea levels rise, the flood risk equation changes. Actuarial Review, 2018. [6] MS Bartlett. Some aspects of the time-correlation problem in regard to tests of signifi- cance. Journal of the Royal Statistical Society, 98(3):536–543, 1935. [7] CB Bell and EP Smith. Infrence for non-negative autoregressive schemes. Communica- tions in Statistics-Theory and Methods, 15(8):2267–2293, 1986. [8] Ocean Studies Board, National Research Council, et al. Sea-level rise for the coasts of California, Oregon, and Washington: past, present, and future. National Academies Press, 2012. [9] David Bolin, Peter Guttorp, Alex Januzzi, Daniel Jones, Marie Novak, Harry Pod- schwit, Lee Richardson, Aila S¨arkk¨a, Colin Sowder, and Aaron Zimmerman. Statistical prediction of global sea level from global temperature. Statistica Sinica, pages 351–367, 2015. [10] John D Boon. Evidence of sea level acceleration at us and canadian tide stations, atlantic coast, north america. Journal of Coastal Research, 28(6):1437–1445, 2012. [11] John D Boon, John M Brubaker, and David R Forrest. Chesapeake bay land subsidence and sea level change: An evaluation of past and present trends and future outlook. 2010. [12] George Edward Pelham Box. Time series analysis forecasting and control. 1970. 110 [13] George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung. Time series analysis: forecasting and control. John Wiley & Sons, 2015. [14] George EP Box and David A Pierce. Distribution of residual autocorrelations in autoregressive-integrated moving average time series models. Journal of the American statistical Association, 65(332):1509–1526, 1970. [15] Peter J Brockwell, Richard A Davis, and Stephen E Fienberg. Time series: theory and methods: theory and methods. Springer Science & Business Media, 1991. [16] J Peter Byrne. The cathedral engulfed: Sea-level rise, property rights, and time. La. L. Rev., 73:69, 2012. [17] Niamh Cahill, Andrew C Kemp, Benjamin P Horton, Andrew C Parnell, et al. Modeling sea-level change using errors-in-variables integrated gaussian processes. The Annals of Applied Statistics, 9(2):547–571, 2015. [18] Richard Chandler and Marian Scott. Statistical methods for trend detection and analysis in the environmental sciences. John Wiley & Sons, 2011. [19] Junesang Choi and HM Srivastava. Evaluation of higher-order derivatives of the gamma function. Publikacije Elektrotehniˇckog fakulteta. Serija Matematika, pages 9–18, 2000. [20] JA Church, PU Clark, A Cazenave, JM Gregory, S Jevrejeva, A Levermann, et al. Sea level change in: Stocker tf, qin d, plattner gk, tignor m, allen sk, boschung j, et al., editors. climate change 2013: The physical science basis. contribution of working group i to the fifth assessment report of the intergovernmental panel on climate change. cambridge, united kingdom and new york, ny, 2013. [21] John A Church and Neil J White. A 20th century acceleration in global sea-level rise. Geophysical research letters, 33(1), 2006. [22] John A Church and Neil J White. Sea-level rise from the late 19th to the early 21st century. Surveys in geophysics, 32(4-5):585–602, 2011. [23] Denis Conniffe and John E Spencer. When moments of ratios are ratios of moments. Journal of the Royal Statistical Society: Series D (The Statistician), 50(2):161–168, 2001. [24] Harald Cramir. Mathematical methods of statistics. Princeton U. Press, Princeton, page 500, 1946. 111 [25] Nelville Davies, CM Triggs, and Paul Newbold. Significance levels of the box-pierce portmanteau statistic in finite samples. Biometrika, 64(3):517–522, 1977. [26] Neville Davies, Trevor Spedding, and William Watson. Autoregressive moving average processes with non-normal residuals. Journal of Time Series Analysis, 1(2):103–109, 1980. [27] L Dewald and P Lewis. A new laplace second-order autoregressive time-series model– nlar (2). IEEE Transactions on Information theory, 31(5):645–651, 1985. [28] Bruce C Douglas. Global sea level acceleration. Journal of Geophysical Research: Oceans, 97(C8):12699–12706, 1992. [29] Gilles R Ducharme and Pierre Lafaye de Micheaux. Goodness-of-fit tests of normality for the innovations in arma models. Journal of Time Series Analysis, 25(3):373–395, 2004. [30] Philip A Ernst, Lawrence D Brown, Larry Shepp, and Robert L Wolpert. Stationary gaussian markov processes as limits of stationary autoregressive time series. Journal of Multivariate Analysis, 155:180–186, 2017. [31] Tal Ezer, Larry P Atkinson, William B Corlett, and Jose L Blanco. Gulf stream’s induced sea level rise and variability along the us mid-atlantic coast. Journal of Geo- physical Research: Oceans, 118(2):685–697, 2013. [32] Tal Ezer and William Bryce Corlett. Analysis of relative sea level variations and trends In 2012 in the chesapeake bay: Is there evidence for acceleration in sea level rise? Oceans, pages 1–5. IEEE, 2012. [33] Tal Ezer and William Bryce Corlett. Is sea level rise accelerating in the chesapeake bay? a demonstration of a novel new approach for analyzing sea level data. Geophysical Research Letters, 39(19), 2012. [34] Jianqing Fan and Wenyang Zhang. Generalised likelihood ratio tests for spectral density. Biometrika, 91(1):195–209, 2004. [35] Donald P Gaver and PAW Lewis. First-order autoregressive gamma sequences and point processes. Advances in Applied Probability, 12(3):727–745, 1980. [36] GK Grunwald, RJ Hyndman, and LM Tedesco. A unified view of linear ar (1) models. 1995. 112 [37] John A Hall, Christopher P Weaver, Jayantha Obeysekera, Mark Crowell, Radley M Horton, Robert E Kopp, John Marburger, Douglas C Marcy, Adam Parris, William V Sweet, et al. Rising sea levels: Helping decision-makers confront the inevitable. Coastal management, 47(2):127–150, 2019. [38] James Hamilton. D.(1994), time series analysis, 1994. [39] Keith W Hipel and A Ian McLeod. Time series modelling of water resources and environmental systems. Elsevier, 1994. [40] Wassily Hoeffding, Herbert Robbins, et al. The central limit theorem for dependent random variables. Duke Mathematical Journal, 15(3):773–780, 1948. [41] SJ Holgate and PL Woodworth. Evidence for enhanced coastal sea level rise during the 1990s. Geophysical research letters, 31(7), 2004. [42] JRM Hosking. A unified derivation of the asymptotic distributions of goodness-of-fit statistics for autoregressive time-series models. Journal of the Royal Statistical Society: Series B (Methodological), 40(3):341–349, 1978. [43] M Imani, R-J You, and C-Y Kuo. Caspian sea level prediction using satellite altime- try by artificial neural networks. International Journal of Environmental Science and Technology, 11(4):1035–1042, 2014. [44] H Bˆaki Iz, L Berry, and M Koch. Modeling regional sea level rise using local tide gauge data. Journal of Geodetic Science, 2(3):188–199, 2012. [45] Patricia A Jacobs and Peter AW Lewis. Discrete time series generated by mixtures. i: Correlational and runs properties. Journal of the Royal Statistical Society: Series B (Methodological), 40(1):94–105, 1978. [46] Patricia A Jacobs and Peter AW Lewis. Discrete time series generated by mixtures ii: asymptotic properties. Journal of the Royal Statistical Society: Series B (Methodologi- cal), 40(2):222–228, 1978. [47] S Jevrejeva, T Frederikse, RE Kopp, Goneri Le Cozannet, LP Jackson, and RSW van de Wal. Probabilistic sea level projections at the coast by 2100. Surveys in Geophysics, 40(6):1673–1696, 2019. [48] S Jevrejeva, A Grinsted, JC Moore, and S Holgate. Nonlinear trends and multiyear cycles in sea level records. Journal of Geophysical Research: Oceans, 111(C9), 2006. 113 [49] DN Joanes and CA Gill. Comparing measures of sample skewness and kurtosis. Journal of the Royal Statistical Society: Series D (The Statistician), 47(1):183–189, 1998. [50] Karl G J¨oreskog, Ulf H Olsson, and Fan Y Wallentin. Multivariate analysis with LIS- REL. Springer, 2016. [51] Jana Jureˇckov´a, Hira L Koul, and Jan Picek. Testing the tail index in autoregressive models. Annals of the Institute of Statistical Mathematics, 61(3):579–598, 2009. [52] Nobuhisa Kashiwagi and Takemi Yanagimoto. Smoothing serial count data through a state-space model. Biometrics, pages 1187–1194, 1992. [53] Robert E Kopp, Jerry X Mitrovica, Stephen M Griffies, Jianjun Yin, Carling C Hay, and Ronald J Stouffer. The impact of greenland melt on local sea levels: a partially coupled analysis of dynamic and static equilibrium effects in idealized water-hosing experiments. Climatic Change, 103(3-4):619–625, 2010. [54] Hira L Koul. Minimum distance estimation and goodness-of-fit tests in first-order au- toregression. The Annals of Statistics, pages 1194–1213, 1986. [55] Hira L Koul, Shiqing Ling, et al. Fitting an error distribution in some heteroscedastic time series models. The Annals of Statistics, 34(2):994–1012, 2006. [56] Hira L Koul, Winfried Stute, et al. Nonparametric model checks for time series. The Annals of Statistics, 27(1):204–236, 1999. [57] Honghai Li, Lihwa Lin, and Kelly A Burks-Copes. Modeling of coastal inundation, storm surge, and relative sea-level rise at naval station norfolk, norfolk, virginia, usa. Journal of coastal research, 29(1):18–30, 2013. [58] Wai Keung Li. Diagnostic checks in time series. CRC Press, 2003. [59] Wai Keung Li and A Ian McLeod. Arma modelling with non-gaussian innovations. Journal of Time Series Analysis, 9(2):155–168, 1988. [60] WK Li and AI McLeod. Distribution of the residual autocorrelations in multivariate arma time series models. Journal of the Royal Statistical Society: Series B (Method- ological), 43(2):231–239, 1981. [61] Greta M Ljung and George EP Box. On a measure of lack of fit in time series models. Biometrika, 65(2):297–303, 1978. 114 [62] Jenny N Lye and Vance L Martin. Non-linear time series modelling and distributional flexibility. Journal of Time Series Analysis, 15(1):65–84, 1994. [63] Gerald A Meehl, Curt Covey, Thomas Delworth, Mojib Latif, Bryant McAvaney, John FB Mitchell, Ronald J Stouffer, and Karl E Taylor. The wcrp cmip3 multimodel dataset: A new era in climate change research. Bulletin of the American meteorological society, 88(9):1383–1394, 2007. [64] Gerard A Meehl, Thomas F Stocker, William D Collins, Pierre Friedlingstein, T Gaye, Jonathan M Gregory, Akio Kitoh, Reto Knutti, James M Murphy, Akira Noda, et al. Global climate projections. 2007. [65] Mark F Meier, Mark B Dyurgerov, Ursula K Rick, Shad O’neel, W Tad Pfeffer, Robert S Anderson, Suzanne P Anderson, and Andrey F Glazovsky. Glaciers dominate eustatic sea-level rise in the 21st century. Science, 317(5841):1064–1067, 2007. [66] Terence C Mills. A Statistical Biography of George Udny Yule: A Loafer of the World. Cambridge Scholars Publishing, 2017. [67] Jerry X Mitrovica, Mark E Tamisiea, James L Davis, and Glenn A Milne. Recent mass balance of polar ice sheets inferred from patterns of global sea-level change. Nature, 409(6823):1026–1029, 2001. [68] Richard H Moss, Jae A Edmonds, Kathy A Hibbard, Martin R Manning, Steven K Rose, Detlef P Van Vuuren, Timothy R Carter, Seita Emori, Mikiko Kainuma, Tom Kram, et al. The next generation of scenarios for climate change research and assessment. Nature, 463(7282):747–756, 2010. [69] Manfred Mudelsee. Climate time series analysis. Springer, 2013. [70] Robert J Nicholls and Anny Cazenave. Sea-level rise and its impact on coastal zones. science, 328(5985):1517–1520, 2010. [71] PE O’Connell and DA Jones. Some experience with the development of models for the stochastic simulation of daily flows. Inputs for Risk Analysis in Water Resources, pages 287–314, 1979. [72] M Perrette, F Landerer, R Riva, K Frieler, and M Meinshausen. A scaling approach to project regional sea level rise and its uncertainties. 2013. [73] David L Prothero and Kenneth F Wallis. Modelling macroeconomic time series. Journal of the Royal Statistical Society: Series A (General), 139(4):468–486, 1976. 115 [74] Asbury H Sallenger, Kara S Doran, and Peter A Howd. Hotspot of accelerated sea-level rise on the atlantic coast of north america. Nature Climate Change, 2(12):884–888, 2012. [75] Robert H Shumway and David S Stoffer. Time series analysis and its applications. Studies In Informatics And Control, 9(4):375–376, 2000. [76] Chiaw Hock Sim. First-order autoregressive logistic processes. Journal of applied prob- ability, 30(2):467–470, 1993. [77] ABA Slangen, Mark Carson, CA Katsman, RSW Van de Wal, Armin K¨ohl, LLA Ver- meersen, and Detlef Stammer. Projecting twenty-first century regional sea-level changes. Climatic Change, 124(1-2):317–332, 2014. [78] ABA Slangen, CA Katsman, RSW Van de Wal, LLA Vermeersen, and REM Riva. Towards regional projections of twenty-first century sea-level change based on ipcc sres scenarios. Climate dynamics, 38(5-6):1191–1209, 2012. [79] Giorgio Spada, JL Bamber, and RTWL Hurkmans. The gravitationally consistent sea- level fingerprint of future terrestrial ice loss. Geophysical Research Letters, 40(3):482– 486, 2013. [80] William V Sweet, Robert E Kopp, Christopher P Weaver, Jayantha Obeysekera, Radley M Horton, E Robert Thieler, and Chris Zervas. Global and regional sea level rise scenarios for the united states. 2017. [81] William V Sweet and Joseph Park. From the extreme to the mean: Acceleration and tipping points of coastal inundation from sea level rise. Earth’s Future, 2(12):579–600, 2014. [82] A Louise Swift and Gareth J Janacek. Forecasting non-normal time series. Journal of Forecasting, 10(5):501–520, 1991. [83] Karl E Taylor, Ronald J Stouffer, and Gerald A Meehl. An overview of cmip5 and the experiment design. Bulletin of the American Meteorological Society, 93(4):485–498, 2012. [84] Ruey S Tsay. Time series and forecasting: Brief history and future research. Journal of the American Statistical Association, 95(450):638–643, 2000. [85] Aad W Van der Vaart. Time series. VU University Amsterdam, lecture notes, 2010. 116 [86] Manouchehr Vaziri. Predicting caspian sea surface water level by ann and arima models. Journal of waterway, port, coastal, and ocean engineering, 123(4):158–162, 1997. [87] Hans Visser, Soenke Dangendorf, and Arthur C Petersen. A review of trend models ap- plied to sea level data with reference to the “acceleration-deceleration debate”. Journal of Geophysical Research: Oceans, 120(6):3873–3895, 2015. [88] AM Walker. Some properties of the asymptotic power functions of goodness-of-fit tests for linear autoregressive schemes. Journal of the Royal Statistical Society: Series B (Methodological), 14(1):117–134, 1952. [89] W Watson, TG King, Trevor A Spedding, and KJ Stout. The machined surface—time series modelling. Wear, 57(1):195–205, 1979. [90] Mike West, P Jeff Harrison, and Helio S Migon. Dynamic generalized linear models and bayesian forecasting. Journal of the American Statistical Association, 80(389):73–83, 1985. [91] PL Woodworth. A search for accelerations in records of european mean sea level. International Journal of Climatology, 10(2):129–143, 1990. [92] PL Woodworth, Neil J White, S Jevrejeva, SJ Holgate, JA Church, and WR Gehrels. Evidence for the accelerations of sea level on multi-decade and century timescales. International Journal of Climatology: A Journal of the Royal Meteorological Society, 29(6):777–789, 2009. [93] Jianjun Yin, Michael E Schlesinger, and Ronald J Stouffer. Model projections of rapid sea-level rise on the northeast coast of the united states. Nature Geoscience, 2(4):262– 266, 2009. [94] G. Udny Yule. On the time-correlation problem, with especial reference to the variate- difference correlation method. Journal of the Royal Statistical Society, 84(4):497–537, 1921. 117