# Parametric Methodologies for Detecting Changes in Maximum Temperature of Tlaxco, Tlaxcala, México

Authors: Silvia Herrera Cortés; Bulmaro Juárez Hernández; Victor Hugo Vázquez Guevara; Hugo Adán Cruz Suárez
Journal: Journal of Probability and Statistics
 (2019)

---

## Abstract

In this paper, comparison results of parametric methodologies of change points, applied to maximum temperature records from the municipality of Tlaxco, Tlaxcala, México, are presented. Methodologies considered are likelihood ratio test, score test, andbinary segmentation (BS), pruned exact linear time (PELT), and segment neighborhood (SN). In order to compare such methodologies, a quality analysis of the data was performed; in addition, lost data were estimated with linear regression, and finally, SARIMA models were adjusted.


---

## Body

## 1. Introduction

Global warming is the most evident manifestation of climate change; it refers to the average increase in global, terrestrial, and marine temperatures [1] and is attributed to the increase in greenhouse gases emissions (mainly CO2, CH4, and N2O). Although greenhouse gases are naturally present in the atmosphere, consumption of fossil fuels for the production of energy and supplies, product of human development, causes emissions of such gases faster than in its natural cycle [2, 3]. The Intergovernmental Panel on Climate Change (IPCC) estimates radical impacts that imply increases in temperature and the rise in sea level (derived from the melting of the poles), as well as the destruction of ecosystemic services [4]; therefore, it is important to monitor and analyze the behaviour of climatological data in order to establish policies, in case of being necessary, that contribute to diminish impacts on the environment and society.Based on the previous remarks, it was considered pertinent to perform an analysis of the behaviour of the maximum temperature in a region of the Mexican state of Tlaxcala, which has a climatological station that registers, among other variables, the maximum temperature; such station is located at municipality Tlaxco-Tlaxcala, and it is the responsibility of the Civil Protection Commission of Tlaxcala State and the Climate Change Research Center of the Autonomous University of Tlaxcala. With aforementioned records, change point methodologies were applied and compared for identifying changes in maximum temperature of the early mentioned municipality; for this, data quality analysis was performed, lost data were completed by linear regression, and SARIMA models were adjusted. In order to contrast the goodness of fit of considered models, the Akaike information criterion (AIC) as well as the Bayesian information criterion (BIC) was considered. In Section2, applied parametric methodologies are described. In Section 3, characterization of the area under study is presented as well as the results obtained from the applied methodologies. In Section 4, the results obtained are discussed.


## 2. Parametric Methodologies for Detecting Change Points

Change point problem started in control theory, and it dates back to 1950s decade [5, 6] where the existence of a single change point was tested. In order to detect a single change point, the first case analyzed dealt with independent random variables considering parametric methodologies [7–13]. To improve criteria for identifying and estimating change points, Bayesian and nonparametric methodologies have been used [7, 8, 14–21].LetXii=1n be a sequence of n independent random vectors each one with probability distribution functions Fi,i=1,2,…,n, respectively. The goal of change point problem is to test(1)H0:F1=F2=⋯=Fn,versus(2)H1:F1=F2=⋯=Fk1≠Fk1+1=⋯=Fk2≠Fk2+1=⋯=Fkq≠Fkq+1=⋯=Fn,where 1<k1<k2<⋯<kq<n, q is the unknown number of change points, and k1<k2<⋯<kq are the respective unknown positions that have to be estimated. When the distributions belong to a common parametric family Fθ, with θ∈ℝp, then the hypothesis tests are about the population parameters:(3)H0:θ1=θ2=⋯=θn,versus(4)H1:θ1=θ2=⋯=θk1≠θk1+1=⋯=θk2≠θk2+1=⋯=θkq≠θkq+1=⋯=θn,where q and k1,k2,…,kq have to be estimated.Since methodologies to detect change points in the context of independent random variables cannot be applied directly in time series [22, 23], parametric, nonparametric, and Bayesian methodologies have also been developed [12, 16, 24–31].Due to the natural disasters recorded in recent decades attributed mainly to anthropogenic activities, there have been an increasing number of studies using change point methods to detect shifts in climate. According to Beaulieu et al. [32], change point detection techniques have been used to detect changes in temperature and in precipitation to detect regime shifts, to detect shifts in aerosol and cloud data, and also to study past changes in the land uptake of carbon. There exist an extensive literature about detecting shifts into environmental data, which apply parametric and nonparametric methodologies [32–38]. The main objective of this work is to compare some parametric methodologies to detect abrupt and multiple changes in the maximum temperature in the municipality of Tlaxco, Tlaxcala, Mexico.


### 2.1. Parametric Methodologies for a Single Change Point in Time Series

For detecting change points in environmental time series, most studies have to apply Bayesian and nonparametric methodologies; various nonparametric tests, including Mann–Kendall test and Pettit’s test, are widely used to detect trend and change point in the historical series of climatic and hydrological variables [35]. In this study, we consider identifying a single change point under the assumption that the data fit an autoregressive model.LetX1, X2,…,Xn be consecutive observations from the model(5)Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt,if−∞<t<τ,α0+α1Xt−1+⋯+αpXt−p+εt,ift≥τ+1,where p∈ℕ, τ∈p,n, εt is a white-noise sequence that has a finite fourth moment, and ϕz=1−ϕ1−⋯−ϕpzp≠0,∀ z≤1.To test(6)H0:τ=nvsHa:τ=k,p<k≤n,under H0 and normality assumption, via likelihood ratio technique conditioned on the first p observations and applying −2ln to the likelihood ratio, is possible to establish(7)Λnk≔−2lnLnLk=minϕ∑p+1≤t≤nXt−ϕ0−∑1≤j≤pϕjXt−j2−minϕ∑p+1≤t≤kXt−ϕ0−∑1≤j≤pϕjXt−j2−minα∑k+1≤t≤nXt−α0−∑1≤j≤pαjXt−j2.For studying the asymptotic behaviour of the distribution of the likelihood ratio statistic underH0, Davis et al. [39] considered a quadratic form and defined the process Qnt=Λnnt over 0,1 and assumed that suptEεt4+δ<∞, for some 0≤δ≤1. The details are in the following theorem.Theorem 2.1.
LetXt be the process defined by(8)Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt,−∞<t<∞,and assume that under H0, εt satisfies Eεiεj=σ2fori=j, with suptEεt4+δ<∞, for some 0<δ≤1, and Xt is strongly related with ρn≪n−1+ε1+4/δ for some ε>0. Then, under H0, we have that ∀ x∈ℝ(9)ℙσ−2maxp<k≤nΛnk−bnp+1anp+1≤x⟶exp−2e−x/2,where bnd=2lnlnn+d/2lnlnlnn−lnΓd/22/2lnlnn and and=bnd/2lnlnn are normalization constants and Γ⋅ is the gamma function.Proof 2.1.
Details of the proof can be consulted in [39].
The methodology proposed by Davis et al. [39] uses a quadratic form and does not allow identifying which parameter of the autoregressive model has changed. Gombay [40] has proposed a method that establishes one- or two-sided hypothesis tests to identify which parameter of an autoregressive model has changed; for this, the proposed methodology uses the components of the efficient score vector as test statistic.
Consider the followingARp process:(10)Xi−μ=ϕ1Xi−1−μ+⋯+ϕpXi−p−μ+εi,i≥p+1,where εi is a white noise sequence with Eε12=σ2.
To test(11)H0:there does not exist a change inξ=μ,σ2,ϕ1,…,ϕpT,versus(12)H1:there exists a change inξ=μ,σ2,ϕ1,…,ϕpT,components of the efficient score vector under Gaussian assumption are needed.
Let▽ξlkX1,…,Xk;ξ=▽ξlkξ, the components are(13)∂∂μlkξ=1−∑j=1pϕjσ2∑i=1kXi−μ−∑j=1pϕjXi−j−μ,∂∂σ2lkξ=−k2σ2+12σ4∑i=1kXi−μ−∑j=1pϕjXi−j−μ2,∂∂ϕslkξ=1σ2∑i=1kXi−μ−∑j=1pϕjXi−j−μXi−s−μ,s=1,…,p,where lk denotes the logarithm function of likelihood based on observations X1,…,Xk. Information matrix of size p+2×p+2 is given by(14)Iξ=Iμ,σ2,ϕ1,…,ϕp=1−∑i=1pϕj2σ20001σ40001σ2V,with V being the covariance matrix of X1,X2,…,Xp. In addition, let μ^n,σ^n2, and ϕ^n be the maximum likelihood estimators of parameters, and consider under H0(15)B^u=n−1/2Iξ^n∂∂μlnuξ^n∂∂σ2lnuξ^n∇ϕlnuξ^n,as test statistic. For testing the change in d parameters with an exact level of significance α, Gombay [40] suggested the use of α∗=1−1−α1/d on each component, for both unilateral and bilateral tests, as follows:Unilateral test: ifsup0≤u≤1B^ju≥C1α∗, then there exists a change in parameter ξjj=1,…,p+2 over the sequence X1,…,Xn. Critic value C1α∗ is obtained from the relation(16)ℙsup0≤u≤1B1u≥x=e−2x.For unilateral test, the estimator time isτ^n=mink∂/∂ξjlkξ^n=max1<m≤n∂/∂ξjlmξ^n.Bilateral test: ifsup0≤u≤1B^ju≥C2α∗, then there exits a change in parameter ξjj=1,…,p+2 over the sequence X1,…,Xn. Critic value C2α∗ is obtained from the relation(17)ℙsup0≤u≤1B1u>x=∑k≠0−1k+1e−2k2x2.
Once a change inξj,j=1,…,p+2 over τ=ρn,0≤ρ≤1, is guaranteed, the estimator for τ is(18)τ^=mink∂∂ξjlkξ^n=max1<m≤n∂∂ξjlmξ^n,orτ^=minn∂∂ξjlkξ^n=min1<m≤n∂∂ξjlmξ^n.Proof 2.2.
Details of the proof can be consulted in [40].


### 2.2. Parametric Methodologies for Multiple Change Points

Letxt,t=1,…,n be an ordered data sequence, assuming that there are m ordered change points in random times τ=τ1,…,τm, with τ0=1, τm+1=n; hence, such change points induce a partition over the data into m+1 segments, where the ith segment will contain observations xτi−1+1,…,xτi. Something commonly used to identify multiple points is to minimize(19)∑1≤i≤m+1Cxτi−1+1,…,τi+βfm,with respect to m and τ1,τ2,…,τm. Here, Cxτi−1+1,…,τi is a cost function for a segment and βfm is a penalty term in order to avoid over fitting [19]. Some commonly used cost functions are the additive inverse of the log likelihood function, accumulated sums [12], and other methods which include applying entropy information, for the penalization of βfm [24]. The most common choice of penalty is one which is linear in the number of change points, that is, βfm=βm. Examples of such penalties include Akaike information criterion (AIC) and Schwarz information criterion (BIC) [41].For the following methods, the change point hypothesis test is(20)H0:there does not exist a change point in the data,versus(21)H1:there exists a change point in the data.


#### 2.2.1. Pruned Exact Linear Time Method

Pruned exact linear time (PELT) method is based on the method proposed by Jackson et al. [42] and Yao [43]. Let xt,t=1,…,n be a set of observations; the idea of the method is the search for an optimal partition of data given on a time interval I whose elements are sets of data cells, resulting in no significant loss of information.A partitionP of an interval I is a set of τ1,…,τm blocks, PI=Bk,k∈ℳ,ℳ≡1,2,…,m, where the blocks are sets of data cells defined by index sets Nk: Bk=xt,t∈Nk satisfying the usual conditions ∪kBk=I and Bk∩Bk′=∅ if k≠k′, and the number of blocks must satisfy 0≤m≤n. Let P∗ be the (finite) set of all partitions of I into blocks. Take as given an additive fitness function that assigns a value to any partition P∈P∗ in the form(22)VP=∑1≤i≤m+1Cxτi−1+1,…,τi+β,where Cxτi−1+1,…,τi is the fitness of block Bk and has to be minimized in order to detect a change point. According to [19], the PELT method proposes a pruning step, which increases the computational efficiency of the previous method while ensuring a global minimum of the cost function which is linear in n (22). The optimal segmentation is Fn, where(23)Fn=minτ∑1≤i≤m+1Cxτi−1+1,…,xτi+β.The pseudocode of the PELT method can be consulted in [19, 44] and proofs about the optimal partition in [42].


#### 2.2.2. Binary Segmentation Method

In thebinary segmentation (BS) method, consider the model xt=ft+ξt,t=1,…,n, where ft is a deterministic, one-dimensional, piecewise-constant signal with change points whose number N and locations τ1,…,τn are unknown, and the random sequence ξt1n is i.i.d. Gaussian with mean zero and variance one.The method uses the CUSUM statistic defined by the inner product between the vectorxs,…,xe and a particular vector of contrast weights given by(24)x˜s,eb=e−bnb−s+1∑t=sbxt−b−s+1n−e−b∑t=b+1ext,where s≤b<e, with n=e−s+1. In its first step, the BS algorithm computes x˜1,nb and then takes b1,1=argmaxb:1≤b<nx˜1,nb to be the first change point candidate, whose significance is to be judged against a certain criterion. If it is considered significant, the domain 1,n is split into two subintervals to the left and to the right of b1,1, and the recursion continues by computing x˜1,b1,1n and x˜b1,1+1,nb possibly resulting in further splits [44].This method extends any single change point method for multiple change points by iteratively repeating the method on different subsets of the sequence for this. It begins by initially applying the single change point method to the entire data and prove if there existsτ such that(25)Cx1:τ+Cxτ+1:n+β<Cx1:n.If (25) does not hold, it is guaranteed that there is no change point along x=x1,x2,…,xn, and the process ends; otherwise, the data are partitioned into two sets, which consist of the data sequence before and after the change point, τα; the process is applied again in each partition until no more change points are found (see pseudocode in [44, 45]).


#### 2.2.3. Segment Neighborhood Method

Segment neighborhood (SN) method uses dynamic programming to detect change point; it begins by setting an upper limit on the size of the segmentation space (i.e., the maximum number of change points) that is required, and this method continues by computing the cost function for all possible segments.It takes the constrained case which segments the data up tot, for t≥m+1, into m+1 segments (using k change points) and denote the minimum value of the cost by Cm,t. The idea of segment neighborhood search is to derive a relationship between Ct,m and for Cs,m−1 for s<t:(26)Cm,t=minτ∑1≤j≤m+1Cxτj−1:τj,=minτ∈m,…,t−1Cm−1,τ+Cxτ:t.If this is run for all values oft up to n and for j=2,…,m, then the exact segmentations with 1,…,m segments can be acquired. To extract the exact segmentation, we first let τl∗t to denote the optimal position of the last change point if we segment data y1:t using l change points. This can be calculated as(27)τl∗targminτ∈l,…,t−1Cl−1,τ+Cxτ+1:t.The pseudocode of the SN method can be consulted in [14].PELT, BS, and SN methods are implemented inR within the library changepoint; according to [46], the library changepoint creates class objects .cpt. The structure in R to detect change in the mean and variance is as follows:(28)cpt.meandata,method=“method”,test.stat=“Normal or CUSUM”,where data = a vector of objects that contain the data where a change in mean will be tested; method = choice for a single change point (AMOC) or multiple change points (PELT, SN, or BS); test.stat = statistic test, i.e., a distribution is assumed; choice “Normal” or “CUSUM” = see [11] or [6] for the normality assumption or the CUSUM test, respectively.For variance(29)cpt.vardata,know.mean=FALSE,method,test.stat=“Normal or CSS”,where data = vector of objects that contain the data where a change in variance will be tested. know.mean = logical argument required only under the normality assumption; if the logical value is TRUE, then the mean is assumed known, and in the other case, mean is estimated via maximum likelihood. test.stat = choice normal or CUSUM test.


## 2.1. Parametric Methodologies for a Single Change Point in Time Series

For detecting change points in environmental time series, most studies have to apply Bayesian and nonparametric methodologies; various nonparametric tests, including Mann–Kendall test and Pettit’s test, are widely used to detect trend and change point in the historical series of climatic and hydrological variables [35]. In this study, we consider identifying a single change point under the assumption that the data fit an autoregressive model.LetX1, X2,…,Xn be consecutive observations from the model(5)Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt,if−∞<t<τ,α0+α1Xt−1+⋯+αpXt−p+εt,ift≥τ+1,where p∈ℕ, τ∈p,n, εt is a white-noise sequence that has a finite fourth moment, and ϕz=1−ϕ1−⋯−ϕpzp≠0,∀ z≤1.To test(6)H0:τ=nvsHa:τ=k,p<k≤n,under H0 and normality assumption, via likelihood ratio technique conditioned on the first p observations and applying −2ln to the likelihood ratio, is possible to establish(7)Λnk≔−2lnLnLk=minϕ∑p+1≤t≤nXt−ϕ0−∑1≤j≤pϕjXt−j2−minϕ∑p+1≤t≤kXt−ϕ0−∑1≤j≤pϕjXt−j2−minα∑k+1≤t≤nXt−α0−∑1≤j≤pαjXt−j2.For studying the asymptotic behaviour of the distribution of the likelihood ratio statistic underH0, Davis et al. [39] considered a quadratic form and defined the process Qnt=Λnnt over 0,1 and assumed that suptEεt4+δ<∞, for some 0≤δ≤1. The details are in the following theorem.Theorem 2.1.
LetXt be the process defined by(8)Xt=ϕ0+ϕ1Xt−1+⋯+ϕpXt−p+εt,−∞<t<∞,and assume that under H0, εt satisfies Eεiεj=σ2fori=j, with suptEεt4+δ<∞, for some 0<δ≤1, and Xt is strongly related with ρn≪n−1+ε1+4/δ for some ε>0. Then, under H0, we have that ∀ x∈ℝ(9)ℙσ−2maxp<k≤nΛnk−bnp+1anp+1≤x⟶exp−2e−x/2,where bnd=2lnlnn+d/2lnlnlnn−lnΓd/22/2lnlnn and and=bnd/2lnlnn are normalization constants and Γ⋅ is the gamma function.Proof 2.1.
Details of the proof can be consulted in [39].
The methodology proposed by Davis et al. [39] uses a quadratic form and does not allow identifying which parameter of the autoregressive model has changed. Gombay [40] has proposed a method that establishes one- or two-sided hypothesis tests to identify which parameter of an autoregressive model has changed; for this, the proposed methodology uses the components of the efficient score vector as test statistic.
Consider the followingARp process:(10)Xi−μ=ϕ1Xi−1−μ+⋯+ϕpXi−p−μ+εi,i≥p+1,where εi is a white noise sequence with Eε12=σ2.
To test(11)H0:there does not exist a change inξ=μ,σ2,ϕ1,…,ϕpT,versus(12)H1:there exists a change inξ=μ,σ2,ϕ1,…,ϕpT,components of the efficient score vector under Gaussian assumption are needed.
Let▽ξlkX1,…,Xk;ξ=▽ξlkξ, the components are(13)∂∂μlkξ=1−∑j=1pϕjσ2∑i=1kXi−μ−∑j=1pϕjXi−j−μ,∂∂σ2lkξ=−k2σ2+12σ4∑i=1kXi−μ−∑j=1pϕjXi−j−μ2,∂∂ϕslkξ=1σ2∑i=1kXi−μ−∑j=1pϕjXi−j−μXi−s−μ,s=1,…,p,where lk denotes the logarithm function of likelihood based on observations X1,…,Xk. Information matrix of size p+2×p+2 is given by(14)Iξ=Iμ,σ2,ϕ1,…,ϕp=1−∑i=1pϕj2σ20001σ40001σ2V,with V being the covariance matrix of X1,X2,…,Xp. In addition, let μ^n,σ^n2, and ϕ^n be the maximum likelihood estimators of parameters, and consider under H0(15)B^u=n−1/2Iξ^n∂∂μlnuξ^n∂∂σ2lnuξ^n∇ϕlnuξ^n,as test statistic. For testing the change in d parameters with an exact level of significance α, Gombay [40] suggested the use of α∗=1−1−α1/d on each component, for both unilateral and bilateral tests, as follows:Unilateral test: ifsup0≤u≤1B^ju≥C1α∗, then there exists a change in parameter ξjj=1,…,p+2 over the sequence X1,…,Xn. Critic value C1α∗ is obtained from the relation(16)ℙsup0≤u≤1B1u≥x=e−2x.For unilateral test, the estimator time isτ^n=mink∂/∂ξjlkξ^n=max1<m≤n∂/∂ξjlmξ^n.Bilateral test: ifsup0≤u≤1B^ju≥C2α∗, then there exits a change in parameter ξjj=1,…,p+2 over the sequence X1,…,Xn. Critic value C2α∗ is obtained from the relation(17)ℙsup0≤u≤1B1u>x=∑k≠0−1k+1e−2k2x2.
Once a change inξj,j=1,…,p+2 over τ=ρn,0≤ρ≤1, is guaranteed, the estimator for τ is(18)τ^=mink∂∂ξjlkξ^n=max1<m≤n∂∂ξjlmξ^n,orτ^=minn∂∂ξjlkξ^n=min1<m≤n∂∂ξjlmξ^n.Proof 2.2.
Details of the proof can be consulted in [40].


## 2.2. Parametric Methodologies for Multiple Change Points

Letxt,t=1,…,n be an ordered data sequence, assuming that there are m ordered change points in random times τ=τ1,…,τm, with τ0=1, τm+1=n; hence, such change points induce a partition over the data into m+1 segments, where the ith segment will contain observations xτi−1+1,…,xτi. Something commonly used to identify multiple points is to minimize(19)∑1≤i≤m+1Cxτi−1+1,…,τi+βfm,with respect to m and τ1,τ2,…,τm. Here, Cxτi−1+1,…,τi is a cost function for a segment and βfm is a penalty term in order to avoid over fitting [19]. Some commonly used cost functions are the additive inverse of the log likelihood function, accumulated sums [12], and other methods which include applying entropy information, for the penalization of βfm [24]. The most common choice of penalty is one which is linear in the number of change points, that is, βfm=βm. Examples of such penalties include Akaike information criterion (AIC) and Schwarz information criterion (BIC) [41].For the following methods, the change point hypothesis test is(20)H0:there does not exist a change point in the data,versus(21)H1:there exists a change point in the data.


### 2.2.1. Pruned Exact Linear Time Method

Pruned exact linear time (PELT) method is based on the method proposed by Jackson et al. [42] and Yao [43]. Let xt,t=1,…,n be a set of observations; the idea of the method is the search for an optimal partition of data given on a time interval I whose elements are sets of data cells, resulting in no significant loss of information.A partitionP of an interval I is a set of τ1,…,τm blocks, PI=Bk,k∈ℳ,ℳ≡1,2,…,m, where the blocks are sets of data cells defined by index sets Nk: Bk=xt,t∈Nk satisfying the usual conditions ∪kBk=I and Bk∩Bk′=∅ if k≠k′, and the number of blocks must satisfy 0≤m≤n. Let P∗ be the (finite) set of all partitions of I into blocks. Take as given an additive fitness function that assigns a value to any partition P∈P∗ in the form(22)VP=∑1≤i≤m+1Cxτi−1+1,…,τi+β,where Cxτi−1+1,…,τi is the fitness of block Bk and has to be minimized in order to detect a change point. According to [19], the PELT method proposes a pruning step, which increases the computational efficiency of the previous method while ensuring a global minimum of the cost function which is linear in n (22). The optimal segmentation is Fn, where(23)Fn=minτ∑1≤i≤m+1Cxτi−1+1,…,xτi+β.The pseudocode of the PELT method can be consulted in [19, 44] and proofs about the optimal partition in [42].


### 2.2.2. Binary Segmentation Method

In thebinary segmentation (BS) method, consider the model xt=ft+ξt,t=1,…,n, where ft is a deterministic, one-dimensional, piecewise-constant signal with change points whose number N and locations τ1,…,τn are unknown, and the random sequence ξt1n is i.i.d. Gaussian with mean zero and variance one.The method uses the CUSUM statistic defined by the inner product between the vectorxs,…,xe and a particular vector of contrast weights given by(24)x˜s,eb=e−bnb−s+1∑t=sbxt−b−s+1n−e−b∑t=b+1ext,where s≤b<e, with n=e−s+1. In its first step, the BS algorithm computes x˜1,nb and then takes b1,1=argmaxb:1≤b<nx˜1,nb to be the first change point candidate, whose significance is to be judged against a certain criterion. If it is considered significant, the domain 1,n is split into two subintervals to the left and to the right of b1,1, and the recursion continues by computing x˜1,b1,1n and x˜b1,1+1,nb possibly resulting in further splits [44].This method extends any single change point method for multiple change points by iteratively repeating the method on different subsets of the sequence for this. It begins by initially applying the single change point method to the entire data and prove if there existsτ such that(25)Cx1:τ+Cxτ+1:n+β<Cx1:n.If (25) does not hold, it is guaranteed that there is no change point along x=x1,x2,…,xn, and the process ends; otherwise, the data are partitioned into two sets, which consist of the data sequence before and after the change point, τα; the process is applied again in each partition until no more change points are found (see pseudocode in [44, 45]).


### 2.2.3. Segment Neighborhood Method

Segment neighborhood (SN) method uses dynamic programming to detect change point; it begins by setting an upper limit on the size of the segmentation space (i.e., the maximum number of change points) that is required, and this method continues by computing the cost function for all possible segments.It takes the constrained case which segments the data up tot, for t≥m+1, into m+1 segments (using k change points) and denote the minimum value of the cost by Cm,t. The idea of segment neighborhood search is to derive a relationship between Ct,m and for Cs,m−1 for s<t:(26)Cm,t=minτ∑1≤j≤m+1Cxτj−1:τj,=minτ∈m,…,t−1Cm−1,τ+Cxτ:t.If this is run for all values oft up to n and for j=2,…,m, then the exact segmentations with 1,…,m segments can be acquired. To extract the exact segmentation, we first let τl∗t to denote the optimal position of the last change point if we segment data y1:t using l change points. This can be calculated as(27)τl∗targminτ∈l,…,t−1Cl−1,τ+Cxτ+1:t.The pseudocode of the SN method can be consulted in [14].PELT, BS, and SN methods are implemented inR within the library changepoint; according to [46], the library changepoint creates class objects .cpt. The structure in R to detect change in the mean and variance is as follows:(28)cpt.meandata,method=“method”,test.stat=“Normal or CUSUM”,where data = a vector of objects that contain the data where a change in mean will be tested; method = choice for a single change point (AMOC) or multiple change points (PELT, SN, or BS); test.stat = statistic test, i.e., a distribution is assumed; choice “Normal” or “CUSUM” = see [11] or [6] for the normality assumption or the CUSUM test, respectively.For variance(29)cpt.vardata,know.mean=FALSE,method,test.stat=“Normal or CSS”,where data = vector of objects that contain the data where a change in variance will be tested. know.mean = logical argument required only under the normality assumption; if the logical value is TRUE, then the mean is assumed known, and in the other case, mean is estimated via maximum likelihood. test.stat = choice normal or CUSUM test.


## 2.2.1. Pruned Exact Linear Time Method

Pruned exact linear time (PELT) method is based on the method proposed by Jackson et al. [42] and Yao [43]. Let xt,t=1,…,n be a set of observations; the idea of the method is the search for an optimal partition of data given on a time interval I whose elements are sets of data cells, resulting in no significant loss of information.A partitionP of an interval I is a set of τ1,…,τm blocks, PI=Bk,k∈ℳ,ℳ≡1,2,…,m, where the blocks are sets of data cells defined by index sets Nk: Bk=xt,t∈Nk satisfying the usual conditions ∪kBk=I and Bk∩Bk′=∅ if k≠k′, and the number of blocks must satisfy 0≤m≤n. Let P∗ be the (finite) set of all partitions of I into blocks. Take as given an additive fitness function that assigns a value to any partition P∈P∗ in the form(22)VP=∑1≤i≤m+1Cxτi−1+1,…,τi+β,where Cxτi−1+1,…,τi is the fitness of block Bk and has to be minimized in order to detect a change point. According to [19], the PELT method proposes a pruning step, which increases the computational efficiency of the previous method while ensuring a global minimum of the cost function which is linear in n (22). The optimal segmentation is Fn, where(23)Fn=minτ∑1≤i≤m+1Cxτi−1+1,…,xτi+β.The pseudocode of the PELT method can be consulted in [19, 44] and proofs about the optimal partition in [42].


## 2.2.2. Binary Segmentation Method

In thebinary segmentation (BS) method, consider the model xt=ft+ξt,t=1,…,n, where ft is a deterministic, one-dimensional, piecewise-constant signal with change points whose number N and locations τ1,…,τn are unknown, and the random sequence ξt1n is i.i.d. Gaussian with mean zero and variance one.The method uses the CUSUM statistic defined by the inner product between the vectorxs,…,xe and a particular vector of contrast weights given by(24)x˜s,eb=e−bnb−s+1∑t=sbxt−b−s+1n−e−b∑t=b+1ext,where s≤b<e, with n=e−s+1. In its first step, the BS algorithm computes x˜1,nb and then takes b1,1=argmaxb:1≤b<nx˜1,nb to be the first change point candidate, whose significance is to be judged against a certain criterion. If it is considered significant, the domain 1,n is split into two subintervals to the left and to the right of b1,1, and the recursion continues by computing x˜1,b1,1n and x˜b1,1+1,nb possibly resulting in further splits [44].This method extends any single change point method for multiple change points by iteratively repeating the method on different subsets of the sequence for this. It begins by initially applying the single change point method to the entire data and prove if there existsτ such that(25)Cx1:τ+Cxτ+1:n+β<Cx1:n.If (25) does not hold, it is guaranteed that there is no change point along x=x1,x2,…,xn, and the process ends; otherwise, the data are partitioned into two sets, which consist of the data sequence before and after the change point, τα; the process is applied again in each partition until no more change points are found (see pseudocode in [44, 45]).


## 2.2.3. Segment Neighborhood Method

Segment neighborhood (SN) method uses dynamic programming to detect change point; it begins by setting an upper limit on the size of the segmentation space (i.e., the maximum number of change points) that is required, and this method continues by computing the cost function for all possible segments.It takes the constrained case which segments the data up tot, for t≥m+1, into m+1 segments (using k change points) and denote the minimum value of the cost by Cm,t. The idea of segment neighborhood search is to derive a relationship between Ct,m and for Cs,m−1 for s<t:(26)Cm,t=minτ∑1≤j≤m+1Cxτj−1:τj,=minτ∈m,…,t−1Cm−1,τ+Cxτ:t.If this is run for all values oft up to n and for j=2,…,m, then the exact segmentations with 1,…,m segments can be acquired. To extract the exact segmentation, we first let τl∗t to denote the optimal position of the last change point if we segment data y1:t using l change points. This can be calculated as(27)τl∗targminτ∈l,…,t−1Cl−1,τ+Cxτ+1:t.The pseudocode of the SN method can be consulted in [14].PELT, BS, and SN methods are implemented inR within the library changepoint; according to [46], the library changepoint creates class objects .cpt. The structure in R to detect change in the mean and variance is as follows:(28)cpt.meandata,method=“method”,test.stat=“Normal or CUSUM”,where data = a vector of objects that contain the data where a change in mean will be tested; method = choice for a single change point (AMOC) or multiple change points (PELT, SN, or BS); test.stat = statistic test, i.e., a distribution is assumed; choice “Normal” or “CUSUM” = see [11] or [6] for the normality assumption or the CUSUM test, respectively.For variance(29)cpt.vardata,know.mean=FALSE,method,test.stat=“Normal or CSS”,where data = vector of objects that contain the data where a change in variance will be tested. know.mean = logical argument required only under the normality assumption; if the logical value is TRUE, then the mean is assumed known, and in the other case, mean is estimated via maximum likelihood. test.stat = choice normal or CUSUM test.


## 3. Study Area and Application of Methodologies

### 3.1. Study Area and Data Quality Analysis

To detect change points in temperature data applying the methodologies described in the last section, the Civil Protection Commission of the State of Tlaxcala and the Climate Change Research Center of the Autonomous University of Tlaxcala provided information of 6 years (September 2, 2011, to September 2, 2017,) from a climatological station, located at the center of the municipality of Tlaxco, Tlaxcala, México.For the quality analysis, the missing data were estimated with regression models; also, to understanding the behaviour of the maximum temperature of the region studied, a box-plot of the data was made (see Figure1), where the maximum temperature per week was considered; for this, it was taken as the first week of the year, the week 1 of January (the week that begins with the first Sunday of the month) and week 52 as the last week of December. We verify that the assumptions of the methodologies studied were fulfilled to apply the studied theory and finally analyzed the results.Figure 1
Box-plot for weekly maximum temperature registers, period: 2/09/2011–2/09/2017.In Figure1, the increase in maximum temperature is observed until it reaches its maximum values in weeks 11 to 23 (the spring season); with the arrival of rainy season (month of June, week 24), it can be seen a decrease in temperature in the box-plot. It also highlights the effect of the heat or drought intraestival (climatic event which consists of a decrease in the amount of precipitation in the middle of the rainy season which occurs mainly in the warm months of July and August); this event is also characterized by the decrease in rainfall due to an increase in temperature (weeks 28 to 33).


### 3.2. Application of Methodologies

#### 3.2.1. Detection of a Single Change Point

In order to apply the methodologies of Davis et al., as well as Gombay, ARIMA models were fitted applying the Box–Jenkins methodology [47]. To verify stationarity and seasonality, autocorrelation function (acf) and partial autocorrelation function (apf) were analyzed. Figure 2 shows the time series for maximum temperature; according to the behaviour of the autocorrelation function (Figure 2(b)), it is observed that the time series is not stationary, with a difference the acf decreases in a fast way guaranteeing the stationarity of the series.Figure 2
(a) Time series for weekly registers. (b) Autocorrelation function of the time series.
(a)
(b)Via the analysis of theacf, ARIMA models were fitted; we consider as the best model the one with the lowest AIC value. Table 1 shows all fitted models and its corresponding AIC values. According to the AIC criterion, the best model to fit the data was the ARIMA3,1,0 model, and coefficients of this model are shown in Table 2. After fitting the model, a residual analysis was made. Noncorrelation of the residuals was verified with the Ljung-Box test. In addition, the Shapiro test was used for verifying normality (Figure 3).Table 1
Models fitted to the data.
Model
AIC
BIC
ARIMA(3,1,0)
1279.85
1.344602
ARIMA(2,1,1)
1281.01
1.348356
ARIMA(3,1,2)
1281.79
1.38138
ARIMA(3,1,1)
1281.79
1.362887
ARIMA(2,1,2)
1282.04
1.363619Table 2
Coefficients of the ARIMA(3,1,0) model.
ar1
ar2
ar3
Constant
s.e.
−0.4872
−0.2981
−0.1795
0.0036
0.0559
0.0600
0.0559
0.0548
σ
2 estimated: 3.563, AIC = 1279.85Figure 3
Residual analysis of theARIMA3,1,0 model.With presented coefficients in Table2, we have that the ARIMA3,1,0 model for fitting the maximum temperature is(30)Tmaxt=0.0036+Tmaxt−1−0.4872Tmaxt−1−Tmaxt−2−0.2981Tmaxt−2−Tmaxt−3−0.1795Tmaxt−3−Tmaxt−4+εt,fort∈ℝ+,where Tmaxt denotes the maximum temperature at time t and εt is a white noise process.After verifying assumptions of Theorem 2.1, values of statisticΛnk were calculated for k=1,…,313 (Figure 4), and the maximum value Λnk was 29.17452 with k=36. Considering α=0.05, critic value in order to accept or reject a change in parameters of model expressed by (3.2.1) is 15.76829; therefore, there exists sufficient evidence to guarantee a change in some parameters.Figure 4
Values of statisticΛnk, for k=1,…,313.For applying Gombay’s methodology, each component of vectorB^u was estimated and plotted (Figure 5). In addition, in Table 3, maximum values of B^u with 0<u<1 as well as the value of α were presented; hence, it may be concluded that there is not considered enough evidence to guarantee a change in any parameter of the model given in Section 3.2.1.Figure 5
Values ofB^k,0<k≤n−1, for each parameter of model given in Section 3.2.1: (a) μ. (b) σ2. (c) ϕ1. (d) ϕ2. (e) ϕ3.
(a)
(b)
(c)
(d)
(e)Table 3
Maximum values to test a change in parametersμ, σ2, ϕ1, ϕ2, and ϕ3.
To test a single change in parameters
k
sup
0
≤
u
≤
1
B
^
u
C
1
0.01021614
μ
245
0.3325948
1.5139
σ
2
300
0.1475097
1.5139
ϕ
1
48
0.5952468
1.5139
ϕ
2
52
0.258343
1.5139
ϕ
3
263
0.4006972
1.5139Applying the librarychangepoint of R, was detected a single change point (AMOC) in the mean and variance of database (Table 4 Figure 6).Table 4
Single change point location for the mean and the variance for maximum temperature using the librarychangepoint of R.
Estimated change point in the mean
Week
Estimated mean before change point
Estimated mean after change point
250
22.64560
21.50152
Estimated change point in the variance
Week
Estimated variance before change point
Estimated variance after change point
27
4.078351
9.970822Figure 6
Location of change points for the (a) mean and (b) variance.
(a)
(b)


#### 3.2.2. Detection of Multiple Change Points

With the librarychangepoint of R, change point analysis of our data was performed with the PELT, BS, and SN methods. Tables 5 and 6 present the number of change points as well as their estimate for mean and variance of the maximum temperature, respectively.Table 5
Number of change points for the mean of maximum temperature.
Method
Number of change points
Estimated change point
Estimated mean
PELT
48
See Table6
See Table6
SN
4
237, 249, 270, 286
22.48312, 25.79167, 21.23333, 17.75000, 23.75333
BS
5
237, 250, 270, 286, 306,
22.48312, 25.60769, 21.12500, 17.75000, 24.48500, 22.29000Table 6
Location of change points in the mean of maximum temperature, determined with the PELT method.
Week
8
16
23
24
25
34
36
42
52
Average
21.38750
23.20000
20.94286
16.90000
20.30000
23.56667
26.90000
24.70000
21.14000
Week
61
62
74
83
93
122
125
129
133
Average
22.32222
18.60000
22.05000
23.94444
25.73000
21.84138
18.96667
22.65000
24.47500
Week
135
141
148
152
164
165
171
172
175
Average
27.45000
25.46667
22.50000
19.35000
22.00000
18.00000
22.06667
17.20000
22.06667
Week
176
181
184
188
194
227
229
234
235
Average
16.90000
20.28000
24.20000
22.07500
25.60000
22.54848
18.00000
21.38000
25.40000
Week
237
249
259
260
270
286
288
289
291
Average
20.60000
25.79167
21.78000
12.20000
21.59000
17.75000
22.70000
25.70000
21.10000
Week
299
302
306
Average
24.80000
26.90000
24.32500
22.29000Tables7 and 8 present all mean and variance change points detected with the PELT method. To have a better conception of the location of change identified points by PELT, BS, and SN methods, and in order to interpret the location of change points, Figures 7 and 8 are presented.Table 7
Number of change points in the variance of maximum temperature.
Method
Number of change points
Estimated change points
Estimated variances
PELT
17
See Table8
See Table8
SN
4
18 194 227 306
1.847222, 4.804260 1.454451, 11.119685, 1.043222
BS
4
18, 194, 227, 306
1.847222, 4.804260, 1.454451, 11.119685, 1.043222Table 8
Location of change points determined with the PELT method.
Week
18
42
82
93
122
152
Variance
1.847222222
6.236503623
2.305583333
1.078545455
1.564655172
7.142298851
Week
158
194
218
222
254
257
Variance
0.092000000
5.700571429
1.654184783
0.056666667
7.507006048
0.003333333
Week
260
270
286
288
306
Variance
25.270000000
1.105444444
0.548000000
0.020000000
3.610882353
1.043222222Figure 7
Location of change points for the mean of maximum temperature. (a) PELT method. (b) SN method. (c) BS method.
(a)
(b)
(c)Figure 8
Location of change points for the variance of maximum temperature. (a) PELT method. (b) SN method. (c) BS method.
(a)
(b)
(c)


### 3.3. Fitted Models after Detecting Change Points in the Mean of Maximum Temperature

After applying methods AMOC, BS, PELT, and SN for mean, we observe that PELT method realizes various segmentations on data; hence, an ARIMA fit is not possible. Method BS realizes six partitions to the database; however, last segmentation has only seven data (from weeks 307 to 313) and then is not possible a SARIMA fit in this segment.In order to model the time series partitions and compare their results, this paper considered partitions from methods AMOC and SN. In the SN method, we considered, like one segmentation the partitions of weeks 276 to 286 (see Figure7), three ARIMAs were fitted for the maximum temperature in lapses: weeks 1 to 237, 238 to 286, and 287 to 313 applying the methodology of Box-Jenkins.Table9 shows the ARIMA models fitted; Figure 9 presents the partitioned time series and its autocorrelation function for each segment.Table 9
Models fitted after applying the AMOC method.
Fitted modelFirst data block
AIC
Fitted modelSecond data block
AIC
AMOC method
ARIMA(2,1,1)
987.06
ARIMA(1,1,1)
273.61
ARIMA(3,1,0)
987.74
ARIMA(2,1,0)
274.79
ARIMA(2,1,2)
988.47
ARIMA(2,1,1)
275.61
ARIMA(3,1,1)
988.85
ARIMA(2,1,2)
277.52
ARIMA(3,1,2)
990.47
ARIMA(2,2,2)
991.34
ARIMA(2,1,0)
994.72Figure 9
Partitioned time series after applying the AMOC method. (a) First data block. (b) Second data block.
(a)
(b)Table10 presents fitted ARIMA models for three partitions considered with the SN method; Figure 10 shows each partition for the time series.Table 10
Models fitted after applying the SN method.
Fitted modelFirst data block
AIC
Fitted modelSecond data block
AIC
Fitted modelThird data block
AIC
ARIMA(2,1,1)
914.29
ARIMA(1,1,1)
191.99
MA(1)
108.29
Método
ARIMA(3,1,2)
917.13
ARIMA(2,1,2)
194.44
ARMA(1,1)
110.21
SN
ARIMA(2,1,2)
924.6
ARIMA(1,1,0)
194.84
MA(2)
110.23
ARIMA(3,1,1)
924.64
ARIMA(1,1,2)
197.58
ARMA(1,2)
111.71
ARMA(2,2)
113.5Figure 10
Partitioned time series after applying the SN method. (a) First data block. (b) Second data block. (c) Third data block.
(a)
(b)
(c)To choose the best model for the partitioned time series, the Akaike information criteria were considered, and after applying the AMOC method, the fit to the time series is presented in (31); equation (32) represents the fit after applying the SN method,(31)Tmaxt=Tmaxt−1+0.0213−0.7750Tmaxt−1−Tmaxt−2+0.0732Tmaxt−2−Tmaxt−3−0.2464εt−1+0.5789εt−2+εt,if1≤t≤250,0.0102+Tmaxt−1+0.0419Tmaxt−1−Tmaxt−2+0.5162εt−1+εt,ift≥250,(32)Tmaxt=−0.0008+Tmaxt−1+0.4547Tmaxt−1−Tmaxt−2+0.1698Tmaxt−2−Tmaxt−3−εt−1+εt,if1≤t≤237,−0.2044+Tmaxt−1+0.3740Tmaxt−1−Tmaxt−2+εt−1+εt,if238≤t≤286,23.8714+0.3968εt−1,ift≥287.


## 3.1. Study Area and Data Quality Analysis

To detect change points in temperature data applying the methodologies described in the last section, the Civil Protection Commission of the State of Tlaxcala and the Climate Change Research Center of the Autonomous University of Tlaxcala provided information of 6 years (September 2, 2011, to September 2, 2017,) from a climatological station, located at the center of the municipality of Tlaxco, Tlaxcala, México.For the quality analysis, the missing data were estimated with regression models; also, to understanding the behaviour of the maximum temperature of the region studied, a box-plot of the data was made (see Figure1), where the maximum temperature per week was considered; for this, it was taken as the first week of the year, the week 1 of January (the week that begins with the first Sunday of the month) and week 52 as the last week of December. We verify that the assumptions of the methodologies studied were fulfilled to apply the studied theory and finally analyzed the results.Figure 1
Box-plot for weekly maximum temperature registers, period: 2/09/2011–2/09/2017.In Figure1, the increase in maximum temperature is observed until it reaches its maximum values in weeks 11 to 23 (the spring season); with the arrival of rainy season (month of June, week 24), it can be seen a decrease in temperature in the box-plot. It also highlights the effect of the heat or drought intraestival (climatic event which consists of a decrease in the amount of precipitation in the middle of the rainy season which occurs mainly in the warm months of July and August); this event is also characterized by the decrease in rainfall due to an increase in temperature (weeks 28 to 33).


## 3.2. Application of Methodologies

### 3.2.1. Detection of a Single Change Point

In order to apply the methodologies of Davis et al., as well as Gombay, ARIMA models were fitted applying the Box–Jenkins methodology [47]. To verify stationarity and seasonality, autocorrelation function (acf) and partial autocorrelation function (apf) were analyzed. Figure 2 shows the time series for maximum temperature; according to the behaviour of the autocorrelation function (Figure 2(b)), it is observed that the time series is not stationary, with a difference the acf decreases in a fast way guaranteeing the stationarity of the series.Figure 2
(a) Time series for weekly registers. (b) Autocorrelation function of the time series.
(a)
(b)Via the analysis of theacf, ARIMA models were fitted; we consider as the best model the one with the lowest AIC value. Table 1 shows all fitted models and its corresponding AIC values. According to the AIC criterion, the best model to fit the data was the ARIMA3,1,0 model, and coefficients of this model are shown in Table 2. After fitting the model, a residual analysis was made. Noncorrelation of the residuals was verified with the Ljung-Box test. In addition, the Shapiro test was used for verifying normality (Figure 3).Table 1
Models fitted to the data.
Model
AIC
BIC
ARIMA(3,1,0)
1279.85
1.344602
ARIMA(2,1,1)
1281.01
1.348356
ARIMA(3,1,2)
1281.79
1.38138
ARIMA(3,1,1)
1281.79
1.362887
ARIMA(2,1,2)
1282.04
1.363619Table 2
Coefficients of the ARIMA(3,1,0) model.
ar1
ar2
ar3
Constant
s.e.
−0.4872
−0.2981
−0.1795
0.0036
0.0559
0.0600
0.0559
0.0548
σ
2 estimated: 3.563, AIC = 1279.85Figure 3
Residual analysis of theARIMA3,1,0 model.With presented coefficients in Table2, we have that the ARIMA3,1,0 model for fitting the maximum temperature is(30)Tmaxt=0.0036+Tmaxt−1−0.4872Tmaxt−1−Tmaxt−2−0.2981Tmaxt−2−Tmaxt−3−0.1795Tmaxt−3−Tmaxt−4+εt,fort∈ℝ+,where Tmaxt denotes the maximum temperature at time t and εt is a white noise process.After verifying assumptions of Theorem 2.1, values of statisticΛnk were calculated for k=1,…,313 (Figure 4), and the maximum value Λnk was 29.17452 with k=36. Considering α=0.05, critic value in order to accept or reject a change in parameters of model expressed by (3.2.1) is 15.76829; therefore, there exists sufficient evidence to guarantee a change in some parameters.Figure 4
Values of statisticΛnk, for k=1,…,313.For applying Gombay’s methodology, each component of vectorB^u was estimated and plotted (Figure 5). In addition, in Table 3, maximum values of B^u with 0<u<1 as well as the value of α were presented; hence, it may be concluded that there is not considered enough evidence to guarantee a change in any parameter of the model given in Section 3.2.1.Figure 5
Values ofB^k,0<k≤n−1, for each parameter of model given in Section 3.2.1: (a) μ. (b) σ2. (c) ϕ1. (d) ϕ2. (e) ϕ3.
(a)
(b)
(c)
(d)
(e)Table 3
Maximum values to test a change in parametersμ, σ2, ϕ1, ϕ2, and ϕ3.
To test a single change in parameters
k
sup
0
≤
u
≤
1
B
^
u
C
1
0.01021614
μ
245
0.3325948
1.5139
σ
2
300
0.1475097
1.5139
ϕ
1
48
0.5952468
1.5139
ϕ
2
52
0.258343
1.5139
ϕ
3
263
0.4006972
1.5139Applying the librarychangepoint of R, was detected a single change point (AMOC) in the mean and variance of database (Table 4 Figure 6).Table 4
Single change point location for the mean and the variance for maximum temperature using the librarychangepoint of R.
Estimated change point in the mean
Week
Estimated mean before change point
Estimated mean after change point
250
22.64560
21.50152
Estimated change point in the variance
Week
Estimated variance before change point
Estimated variance after change point
27
4.078351
9.970822Figure 6
Location of change points for the (a) mean and (b) variance.
(a)
(b)


### 3.2.2. Detection of Multiple Change Points

With the librarychangepoint of R, change point analysis of our data was performed with the PELT, BS, and SN methods. Tables 5 and 6 present the number of change points as well as their estimate for mean and variance of the maximum temperature, respectively.Table 5
Number of change points for the mean of maximum temperature.
Method
Number of change points
Estimated change point
Estimated mean
PELT
48
See Table6
See Table6
SN
4
237, 249, 270, 286
22.48312, 25.79167, 21.23333, 17.75000, 23.75333
BS
5
237, 250, 270, 286, 306,
22.48312, 25.60769, 21.12500, 17.75000, 24.48500, 22.29000Table 6
Location of change points in the mean of maximum temperature, determined with the PELT method.
Week
8
16
23
24
25
34
36
42
52
Average
21.38750
23.20000
20.94286
16.90000
20.30000
23.56667
26.90000
24.70000
21.14000
Week
61
62
74
83
93
122
125
129
133
Average
22.32222
18.60000
22.05000
23.94444
25.73000
21.84138
18.96667
22.65000
24.47500
Week
135
141
148
152
164
165
171
172
175
Average
27.45000
25.46667
22.50000
19.35000
22.00000
18.00000
22.06667
17.20000
22.06667
Week
176
181
184
188
194
227
229
234
235
Average
16.90000
20.28000
24.20000
22.07500
25.60000
22.54848
18.00000
21.38000
25.40000
Week
237
249
259
260
270
286
288
289
291
Average
20.60000
25.79167
21.78000
12.20000
21.59000
17.75000
22.70000
25.70000
21.10000
Week
299
302
306
Average
24.80000
26.90000
24.32500
22.29000Tables7 and 8 present all mean and variance change points detected with the PELT method. To have a better conception of the location of change identified points by PELT, BS, and SN methods, and in order to interpret the location of change points, Figures 7 and 8 are presented.Table 7
Number of change points in the variance of maximum temperature.
Method
Number of change points
Estimated change points
Estimated variances
PELT
17
See Table8
See Table8
SN
4
18 194 227 306
1.847222, 4.804260 1.454451, 11.119685, 1.043222
BS
4
18, 194, 227, 306
1.847222, 4.804260, 1.454451, 11.119685, 1.043222Table 8
Location of change points determined with the PELT method.
Week
18
42
82
93
122
152
Variance
1.847222222
6.236503623
2.305583333
1.078545455
1.564655172
7.142298851
Week
158
194
218
222
254
257
Variance
0.092000000
5.700571429
1.654184783
0.056666667
7.507006048
0.003333333
Week
260
270
286
288
306
Variance
25.270000000
1.105444444
0.548000000
0.020000000
3.610882353
1.043222222Figure 7
Location of change points for the mean of maximum temperature. (a) PELT method. (b) SN method. (c) BS method.
(a)
(b)
(c)Figure 8
Location of change points for the variance of maximum temperature. (a) PELT method. (b) SN method. (c) BS method.
(a)
(b)
(c)


## 3.2.1. Detection of a Single Change Point

In order to apply the methodologies of Davis et al., as well as Gombay, ARIMA models were fitted applying the Box–Jenkins methodology [47]. To verify stationarity and seasonality, autocorrelation function (acf) and partial autocorrelation function (apf) were analyzed. Figure 2 shows the time series for maximum temperature; according to the behaviour of the autocorrelation function (Figure 2(b)), it is observed that the time series is not stationary, with a difference the acf decreases in a fast way guaranteeing the stationarity of the series.Figure 2
(a) Time series for weekly registers. (b) Autocorrelation function of the time series.
(a)
(b)Via the analysis of theacf, ARIMA models were fitted; we consider as the best model the one with the lowest AIC value. Table 1 shows all fitted models and its corresponding AIC values. According to the AIC criterion, the best model to fit the data was the ARIMA3,1,0 model, and coefficients of this model are shown in Table 2. After fitting the model, a residual analysis was made. Noncorrelation of the residuals was verified with the Ljung-Box test. In addition, the Shapiro test was used for verifying normality (Figure 3).Table 1
Models fitted to the data.
Model
AIC
BIC
ARIMA(3,1,0)
1279.85
1.344602
ARIMA(2,1,1)
1281.01
1.348356
ARIMA(3,1,2)
1281.79
1.38138
ARIMA(3,1,1)
1281.79
1.362887
ARIMA(2,1,2)
1282.04
1.363619Table 2
Coefficients of the ARIMA(3,1,0) model.
ar1
ar2
ar3
Constant
s.e.
−0.4872
−0.2981
−0.1795
0.0036
0.0559
0.0600
0.0559
0.0548
σ
2 estimated: 3.563, AIC = 1279.85Figure 3
Residual analysis of theARIMA3,1,0 model.With presented coefficients in Table2, we have that the ARIMA3,1,0 model for fitting the maximum temperature is(30)Tmaxt=0.0036+Tmaxt−1−0.4872Tmaxt−1−Tmaxt−2−0.2981Tmaxt−2−Tmaxt−3−0.1795Tmaxt−3−Tmaxt−4+εt,fort∈ℝ+,where Tmaxt denotes the maximum temperature at time t and εt is a white noise process.After verifying assumptions of Theorem 2.1, values of statisticΛnk were calculated for k=1,…,313 (Figure 4), and the maximum value Λnk was 29.17452 with k=36. Considering α=0.05, critic value in order to accept or reject a change in parameters of model expressed by (3.2.1) is 15.76829; therefore, there exists sufficient evidence to guarantee a change in some parameters.Figure 4
Values of statisticΛnk, for k=1,…,313.For applying Gombay’s methodology, each component of vectorB^u was estimated and plotted (Figure 5). In addition, in Table 3, maximum values of B^u with 0<u<1 as well as the value of α were presented; hence, it may be concluded that there is not considered enough evidence to guarantee a change in any parameter of the model given in Section 3.2.1.Figure 5
Values ofB^k,0<k≤n−1, for each parameter of model given in Section 3.2.1: (a) μ. (b) σ2. (c) ϕ1. (d) ϕ2. (e) ϕ3.
(a)
(b)
(c)
(d)
(e)Table 3
Maximum values to test a change in parametersμ, σ2, ϕ1, ϕ2, and ϕ3.
To test a single change in parameters
k
sup
0
≤
u
≤
1
B
^
u
C
1
0.01021614
μ
245
0.3325948
1.5139
σ
2
300
0.1475097
1.5139
ϕ
1
48
0.5952468
1.5139
ϕ
2
52
0.258343
1.5139
ϕ
3
263
0.4006972
1.5139Applying the librarychangepoint of R, was detected a single change point (AMOC) in the mean and variance of database (Table 4 Figure 6).Table 4
Single change point location for the mean and the variance for maximum temperature using the librarychangepoint of R.
Estimated change point in the mean
Week
Estimated mean before change point
Estimated mean after change point
250
22.64560
21.50152
Estimated change point in the variance
Week
Estimated variance before change point
Estimated variance after change point
27
4.078351
9.970822Figure 6
Location of change points for the (a) mean and (b) variance.
(a)
(b)


## 3.2.2. Detection of Multiple Change Points

With the librarychangepoint of R, change point analysis of our data was performed with the PELT, BS, and SN methods. Tables 5 and 6 present the number of change points as well as their estimate for mean and variance of the maximum temperature, respectively.Table 5
Number of change points for the mean of maximum temperature.
Method
Number of change points
Estimated change point
Estimated mean
PELT
48
See Table6
See Table6
SN
4
237, 249, 270, 286
22.48312, 25.79167, 21.23333, 17.75000, 23.75333
BS
5
237, 250, 270, 286, 306,
22.48312, 25.60769, 21.12500, 17.75000, 24.48500, 22.29000Table 6
Location of change points in the mean of maximum temperature, determined with the PELT method.
Week
8
16
23
24
25
34
36
42
52
Average
21.38750
23.20000
20.94286
16.90000
20.30000
23.56667
26.90000
24.70000
21.14000
Week
61
62
74
83
93
122
125
129
133
Average
22.32222
18.60000
22.05000
23.94444
25.73000
21.84138
18.96667
22.65000
24.47500
Week
135
141
148
152
164
165
171
172
175
Average
27.45000
25.46667
22.50000
19.35000
22.00000
18.00000
22.06667
17.20000
22.06667
Week
176
181
184
188
194
227
229
234
235
Average
16.90000
20.28000
24.20000
22.07500
25.60000
22.54848
18.00000
21.38000
25.40000
Week
237
249
259
260
270
286
288
289
291
Average
20.60000
25.79167
21.78000
12.20000
21.59000
17.75000
22.70000
25.70000
21.10000
Week
299
302
306
Average
24.80000
26.90000
24.32500
22.29000Tables7 and 8 present all mean and variance change points detected with the PELT method. To have a better conception of the location of change identified points by PELT, BS, and SN methods, and in order to interpret the location of change points, Figures 7 and 8 are presented.Table 7
Number of change points in the variance of maximum temperature.
Method
Number of change points
Estimated change points
Estimated variances
PELT
17
See Table8
See Table8
SN
4
18 194 227 306
1.847222, 4.804260 1.454451, 11.119685, 1.043222
BS
4
18, 194, 227, 306
1.847222, 4.804260, 1.454451, 11.119685, 1.043222Table 8
Location of change points determined with the PELT method.
Week
18
42
82
93
122
152
Variance
1.847222222
6.236503623
2.305583333
1.078545455
1.564655172
7.142298851
Week
158
194
218
222
254
257
Variance
0.092000000
5.700571429
1.654184783
0.056666667
7.507006048
0.003333333
Week
260
270
286
288
306
Variance
25.270000000
1.105444444
0.548000000
0.020000000
3.610882353
1.043222222Figure 7
Location of change points for the mean of maximum temperature. (a) PELT method. (b) SN method. (c) BS method.
(a)
(b)
(c)Figure 8
Location of change points for the variance of maximum temperature. (a) PELT method. (b) SN method. (c) BS method.
(a)
(b)
(c)


## 3.3. Fitted Models after Detecting Change Points in the Mean of Maximum Temperature

After applying methods AMOC, BS, PELT, and SN for mean, we observe that PELT method realizes various segmentations on data; hence, an ARIMA fit is not possible. Method BS realizes six partitions to the database; however, last segmentation has only seven data (from weeks 307 to 313) and then is not possible a SARIMA fit in this segment.In order to model the time series partitions and compare their results, this paper considered partitions from methods AMOC and SN. In the SN method, we considered, like one segmentation the partitions of weeks 276 to 286 (see Figure7), three ARIMAs were fitted for the maximum temperature in lapses: weeks 1 to 237, 238 to 286, and 287 to 313 applying the methodology of Box-Jenkins.Table9 shows the ARIMA models fitted; Figure 9 presents the partitioned time series and its autocorrelation function for each segment.Table 9
Models fitted after applying the AMOC method.
Fitted modelFirst data block
AIC
Fitted modelSecond data block
AIC
AMOC method
ARIMA(2,1,1)
987.06
ARIMA(1,1,1)
273.61
ARIMA(3,1,0)
987.74
ARIMA(2,1,0)
274.79
ARIMA(2,1,2)
988.47
ARIMA(2,1,1)
275.61
ARIMA(3,1,1)
988.85
ARIMA(2,1,2)
277.52
ARIMA(3,1,2)
990.47
ARIMA(2,2,2)
991.34
ARIMA(2,1,0)
994.72Figure 9
Partitioned time series after applying the AMOC method. (a) First data block. (b) Second data block.
(a)
(b)Table10 presents fitted ARIMA models for three partitions considered with the SN method; Figure 10 shows each partition for the time series.Table 10
Models fitted after applying the SN method.
Fitted modelFirst data block
AIC
Fitted modelSecond data block
AIC
Fitted modelThird data block
AIC
ARIMA(2,1,1)
914.29
ARIMA(1,1,1)
191.99
MA(1)
108.29
Método
ARIMA(3,1,2)
917.13
ARIMA(2,1,2)
194.44
ARMA(1,1)
110.21
SN
ARIMA(2,1,2)
924.6
ARIMA(1,1,0)
194.84
MA(2)
110.23
ARIMA(3,1,1)
924.64
ARIMA(1,1,2)
197.58
ARMA(1,2)
111.71
ARMA(2,2)
113.5Figure 10
Partitioned time series after applying the SN method. (a) First data block. (b) Second data block. (c) Third data block.
(a)
(b)
(c)To choose the best model for the partitioned time series, the Akaike information criteria were considered, and after applying the AMOC method, the fit to the time series is presented in (31); equation (32) represents the fit after applying the SN method,(31)Tmaxt=Tmaxt−1+0.0213−0.7750Tmaxt−1−Tmaxt−2+0.0732Tmaxt−2−Tmaxt−3−0.2464εt−1+0.5789εt−2+εt,if1≤t≤250,0.0102+Tmaxt−1+0.0419Tmaxt−1−Tmaxt−2+0.5162εt−1+εt,ift≥250,(32)Tmaxt=−0.0008+Tmaxt−1+0.4547Tmaxt−1−Tmaxt−2+0.1698Tmaxt−2−Tmaxt−3−εt−1+εt,if1≤t≤237,−0.2044+Tmaxt−1+0.3740Tmaxt−1−Tmaxt−2+εt−1+εt,if238≤t≤286,23.8714+0.3968εt−1,ift≥287.


## 4. Comparison of Methodologies and Conclusion

The methodologies of Davis et al. [39] and Gombay [40] allow us to test if there exists a single change point or not in some parameter of the autoregressive model of order p; however, the methodology of Davis et al. does not allow us to estimate the time in which it occurs. When applying these methodologies in the database of maximum temperature, found results show with the methodology of Davis et al., evidence of a change in parameters of the proposed autoregressive model; on the other hand, the methodology of Gombay does not provide sufficient evidence to consider a change; this is attributed to the fact that for fitting the autoregressive model to database, a difference had to be made.The other applied methodologies do not consider the structure of data, apart from being designed to estimate change points in mean and variance. The AMOC method estimates a single change point, while the PELT method estimates 48 change points in the considered database. Segmentation of the latter methodology makes it hard to propose models to fit the database.For change in the mean of maximum temperature, methods AMOC, PELT, SN, and BS coincided with the conclusion of a change in week 250 when comparing this with methodology of Gombay, which proposes a possible change in week 245 is very close to that obtained by the other methods.Gombay’s methodology does not have enough evidence to decide a change in variance; the AMOC method estimates a change in variance in week 27, and the PELT, SN, and BS methods coincide with a change in weeks 18, 194, 227, and 306. With the AMOC method, the change is reported in week 250, which represents the third week of June 2016, reviewing the cold and warm episodes per season reported by the Climate Prediction Center [48]; in the June month of 2016, it is observed that the ocean changed from a heating to a cooling, and this allows us to validate the change reflected with the AMOC method; similarly, for the SN method, estimated change points coincide with the transition from a heating to a cooling of the ocean (week 237, which represents the last week of February month, 2016) and a period without heating of the ocean to a cooling (week 287, which represents June month, 2017).Although PELT and BS methods detect a change in week 306 (Tables5–8) which is almost at the end, this week is in the latest record; therefore, it might be an unreliable change point.The study of the points of change in climate data at the regional level is essential for decision-makers and for the local development of the communities to be analyzed, since the level of impact on their economic growth can be dimensioned at the same time as they will be able to establish measures to adapt to these changes with society and political-economic sectors.For the subsequent studies, we intend to consider seasonal models as well as multivariate models that may be related to describe the climate of some regions.



---

*Source: 3580692-2019-12-16.xml*