The Local Partial Autocorrelation Function and Some Applications

The classical regular and partial autocorrelation functions are powerful tools for stationary time series modelling and analysis. However, it is increasingly recognized that many time series are not stationary and the use of classical global autocorrelations can give misleading answers. This article introduces two estimators of the local partial autocorrelation function and establishes their asymptotic properties. The article then illustrates the use of these new estimators on both simulated and real time series. The examples clearly demonstrate the strong practical benefits of local estimators for time series that exhibit nonstationarities.


Introduction
Much work has been undertaken to develop both theory and methods for the use of the autocorrelation and partial autocorrelation for mean zero second-order stationary time series. See, for example, Priestley (1983), Brockwell and Davis (1991) or Chatfield (2003). For stationary time series, both autocorrelations are fundamental for eliciting second-order structure and are particularly useful for subsequent modelling and prediction. Unfortunately, in many applied situations, for example neurophysiology (Fiecas and Ombao, 2016) or biology (Hargreaves et al., 2018), the stationarity assumption is not tenable and, hence, use of the classical stationary-based autocorrelations is highly questionable. Indeed, it is not possible for a time-varying parameter to be adequately summarised by a single coefficient. Before practical analysis, one should therefore attempt to assess whether the series is stationary or not. Many techniques and software packages exist that enable such assessment, see reviews in Dahlhaus (2012) or Cardinali and Nason (2018) or newer techniques that measure, rather than test, the degree of nonstationarity, e.g. Das and Nason (2016).
A large literature on nonstationary time series modelling has developed since the 1950s. See, for example, Page (1952), Silverman (1957), Whittle (1963), Priestley (1965), Tong (1974) and Dahlhaus (1997). Alternative model forms including the piecewise stationary time series of Adak (1998); the wavelet models of Nason et al. (2000); and the SLEX models of Ombao et al. (2002) have been proposed. A comprehensive review of locally stationary series can be found in Dahlhaus (2012). As part of these developments, the local autocovariance, for non-or locally stationary processes, has been studied in the literature and details on specific estimators can be found in Hyndman and Wand (1997), Nason (2013c), Cardinali (2014) and Zhao (2015), for example. However, to date, little attention seems to have been paid to local partial autocorrelation and the benefits it could bring. An exception is Degerine and Lambert (1996) and Degerine and Lambert-Lacroix (2003), who extended the classical partial autocorrelation to encompass nonstationary processes. Their seminal work mentions estimation, including the windowing idea that we use in Section 3, but provides no theory for their estimator nor evaluation via simulation or on real time series. More recently, Yang et al. (2016) use a hierarchical Bayesian modelling approach to estimate process time-frequency structure, linking the time-dependent partial autocorrelations to the coefficients of a time-varying autoregressive process.
Autocorrelation and partial autocorrelation are intimately related, presenting complementary views on the underlying structure within a time series. For example, arguably, partial autocorrelation provides direct information on the order and underlying structure of autoregressive-type processes (see Appendix A for additional background on its interpretation). As in the stationary case, for real-life statistical analysis one needs both local autocorrelation and partial autocorrelation. This article fills the gap for the latter. We introduce two new estimators of the local partial autocorrelation function, supplying new results on their theoretical properties. We further exhibit our estimators on a simulated series and three real time series that demonstrate the importance of using a local approach. In addition, our work also provides a freeware R software package, lpacf, for local partial autocorrelation that complements existing software for local autocorrelations, such as lacf in the locits package.
2 The Local Partial Autocorrelation Function 2.1 The (process) local partial autocorrelation function, q T , for a locally stationary process Let {X t,T } be a zero-mean locally stationary process such as the locally stationary Fourier process, Dahlhaus (1997, Definition 2.1), or the locally stationary wavelet process, Nason et al. (2000, Definition 1) (for ease of reference, these definitions can also be found in Appendix B). Locally stationary process theory supports short-memory processes and often has quantities of interest such as the time-varying spectrum, f (z, ω) at (rescaled) time z ∈ (0, 1) and frequency ω, or local autocovariance c(z, τ ) at location z and lag τ , which are estimated via a process quantity (f T or c T ), which depends on the sample size T and asymptotically approaches to the quantity of interest as T → ∞. Consider, for example, f T (z, ω) from Neumann and von Sachs (1997) or c T (z, τ ) from Nason et al. (2000). We follow this paradigm by first introducing the process local partial autocorrelation, q T . The (process) partial autocorrelation function, q T (z, τ ), of a zero-mean locally stationary process can be understood informally as q T (z, τ ) = corr X [zT ],T , X [zT ]+τ,T |"in-between" data , where [x] denotes the integer part of the real number x. A formal definition follows.
The next proposition shows an alternative useful representation of q T .
Proposition 2.2. Let {X t,T } be a zero-mean locally stationary process. Then the process local partial autocorrelation, q T , can be expressed as where ϕ [zT ],τ,τ ;T is from projecting X Proof. See Section I.1.
Formulae (1) and (2) are natural generalisations of their stationary equivalents, compare for example with Definitions 3.4.1 and 3.4.2 from Brockwell and Davis (1991). 2.2 Equivalent expressions for the process local partial autocorrelation function, q T As a step to estimation, we will express q T by exploiting a well-known connection between partial autocorrelation and linear prediction. We introduce the following notation P [zT ],τ (X [zT ],T ) =X (b) [zT ],T and P [zT ],τ (X [zT ]+τ ,T ) =X (f ) [zT ]+τ,T . These are simply the respective linear predictors of X [zT ],T (back-casted), and X [zT ]+τ ,T (forecasted), using the predictor set X [zT ]+1,T , . . . , X [zT ]+τ −1,T . The numerator and denominator in (2) can be re-expressed as a Mean Squared Prediction Error (MSPE). Consequently, we can rewrite q T (z, τ ) as (4) For details see Section I.2. For stationary processes the square root term in (4) equals one and q T (z, τ ) coincides with the classical q(τ ). In general, given t observations of a zero-mean locally stationary process, X 0,T , . . . , X t−1,T , the mean squared prediction error of a linear predictor of X t,T ,X t,T = t−1 s=0 b t−1−s,T X s,T , can be written as where b t = (b t−1,T , . . . , b 0,T , −1) T and Σ t,T is the covariance of X 0,T , . . . , X t,T , see, e.g., Fryzlewicz et al. (2003, Section 3.3). In our case, the back-casted and forecasted values of X [zT ],T and X [zT ]+τ ,T are also linear predictors using the window of observations X [zT ]+1,T , . . . , X [zT ]+τ −1,T , and can be expressed aŝ , b (f ) coefficient vectors are obtained through minimisation of the corresponding mean squared prediction error using the same principle as in the stationary case. We next give a proposition that paves the way towards a natural definition of the local partial autocorrelation function q in Section 2.3. Proposition 2.3. Let {X t,T } be a zero-mean locally stationary process. Then q T can also be expressed as where ϕ [zT ],τ,τ ;T is as in (3), and b [zT ]+τ ;T on τ , even though it is still present. The τ × τ covariance matrices Σ (b) [zT ];T and Σ (f ) [zT ]+τ ;T are given in Appendix C.

the coefficient vectors
[zT ] are obtained as the solution to the forecasting and back-casting prediction equations, or equivalently through minimisation of the MSPE. See Section 3.1 and Proposition 3.1 from Fryzlewicz et al. (2003) for details.
Next, Proposition 2.5 shows that the (process) local partial autocorrelation, q T , converges to the local partial autocorrelation, q, defined by (6).

Wavelet local partial autocorrelation estimation
We now consider the important problem of local partial autocorrelation estimation. We begin by first noting that all the quantities on the right-hand side of (6) for q(z, τ ) are based on the local autocovariance c(z, τ ). A natural estimator of q can thus be obtained by replacing all occurrences of c(z, τ ) by the wavelet-based estimatorĉ(z, τ ) from Nason (2013c, Section 3.3) as follows.
Definition 2.6. The wavelet-based local partial autocorrelation estimator is defined asq where the matrix estimates,B [zT ] ,B [zT ]+τ , and vector estimatesb [zT ]+τ , are obtained from their population quantities in Sections 2.2 and 2.3 by plugging in the wavelet-based local autocovariance estimatorĉ from Nason (2013c). Similarly, the vectorφ [zT ],τ is obtained as the solution to the local Yule-Walker equations in Definition 2.4 again replacing c byĉ. We next establish the consistency ofq for q.
Proposition 2.7. Let {X t,T } be a zero-mean locally stationary wavelet process under the assumptions given in Definition 2.4. The local partial autocorrelation estimatorq (z, τ ) from (7) is consistent for the true local partial autocorrelation Proof. See Section I.4.
Our wavelet-based estimator,q (z, τ ), develops earlier work on forecasting by Fryzlewicz (2003) in a new direction. However, the estimator is not simple to implement and, as we will see later, does not perform as well as the following alternative approach, which applies a window to the classical partial autocorrelation.
3 Windowed Estimation of Local Partial Autocorrelation

The integrated local wavelet periodogram
We introduce an alternative estimator,q W (z, τ ), that is simpler to implement thanq (z, τ ), and turns out to perform better. This new estimator is constructed by windowing the classical partial autocorrelation (designed for stationary processes) over an interval centred at time [zT ] with length L(T ), where L(T ) → ∞ and L(T )/T → 0, as T → ∞. Proposition 3.4, in Section 3.2, establishes the asymptotic behaviour ofq W (z, τ ) by approximating the integrated local wavelet periodogram of a (zero-mean) locally stationary wavelet process by its equivalent stationary version at a fixed rescaled time (see Theorem 1). The proof of the theorem introduces new bounds for quantities involving cross-correlation wavelets, as well as a new exact formula for cross-correlation Haar wavelets. Key definitions and results are presented below, while full proofs are provided in Section J.
Definition 3.1. Let {X t,T } be a locally stationary wavelet process as in Definition 1 from Nason et al. (2000) with evolutionary wavelet spectrum {S j (z)} ∞ j=1 for z ∈ (0, 1), Lipschitz constants {L j } ∞ j=1 , process constants {C j } ∞ j=1 and underlying discrete nondecimated wavelets {ψ j,k }. The integrated local periodogram on the interval [[zT ] − L(T )/2 + 1, [zT ] + L(T )/2] is given by Here {φ j } ∞ j=1 ∈ Φ and Φ is a set of complex-valued bounded sequences equipped with uniform norm ||φ|| ∞ := sup j |φ j |, z ∈ (0, 1) and, for j ∈ N, I * L(T ) (z, j) is the uncorrected, tapered local wavelet periodogram given by We next approximate the integrated local (wavelet) periodogram, J L(T ) (z, φ), by the corresponding statistics of a stationary process {Y t } with the same local corresponding statistics at t = zT , for fixed z. Conceptually, this is a common approach useful in establishing asymptotic properties for functions of locally stationary processes (Dahlhaus and Giraitis, 1998), which in this work we advance to include wavelet-based expansions. Specifically, define is the wavelet periodogram on the segment [zT ] − L(T )/2 + 1, . . . , [zT ] + L(T )/2 of the stationary process Here ψ j,k is the same wavelet sequence as previously, {ξ j,k } a set of independent identically distributed random variables with mean zero and unit variance and W j (z) is such that W 2 j (z) = S j (z) for all z ∈ (0, 1) and j ∈ N. The next theorem is the key result establishing the asymptotic properties of the integrated local wavelet periodogram.
Theorem 1. Let {X t,T } be a zero-mean Gaussian locally stationary wavelet process as defined by Definition 3.1. Suppose j C 2 j 2 2j < ∞, {W j } j is Lipschitz continuous with Lipschitz constants L j such that j L 2 j 2 2j < ∞ and j W 2 j (z)2 2j < ∞ at any rescaled time z, L(T )/T → 0 as T → ∞, and φ ∈ Φ is a sequence of bounded variation. Further, assume h is a rectangular kernel. Then, using the family of discrete Haar wavelets, we have Proof. Section J.5 contains the full proof.
An important difference between earlier literature in this area and our work is the introduction of windowing. We provide new results on windowed versions of the cross-correlation wavelets, which we denote i N,z , where, to simplify notation, we replace L(T ) by N and sometimes omit z. To prove Theorem 1 we need bounds on quantities involving i N,z which we can obtain via their connection with cross-correlation wavelets and, in particular, our new closed form expression for the cross-correlation Haar wavelet. For completeness, we define the truncated cross-correlation wavelet here and some of the key bounds.
Definition 3.2. For N ∈ N, scales j, ∈ N and rescaled time z ∈ (0, 1), the windowed cross-scale autocorrelation wavelets i N,z (j, , · ) over the interval [[zT ] where {ψ j,m } j,m is a family of discrete wavelets and k ∈ Z.
The similarity between the cross-scale autocorrelation wavelets Ψ j, (· ), defined in Fryzlewicz (2003, Definition 5.4.2) as Ψ j, (τ ) = k∈Z ψ j,k ψ ,k+τ for j, ∈ N and τ ∈ Z, and their windowed version, i N,z (j, , · ) defined above, is key to how we subsequently bound quantities involving i N,z . The exact new formulae for Haar cross-scale autocorrelation wavelets are established in Appendix D, along with a pictorial description in Figure 6 in Appendix F.
As bounds for i N (j, , · ) are a key component of the proof of Theorem 1, these are provided by the next three results. The first bound for i N is valid for all discrete wavelets based on Daubechies (1992) compactly supported wavelets, although we later only use it for Haar wavelets.
Lemma 1. Using previous notation and assumptions, let b 1 = [zT ] + N/2 + 1 and b 2 = [zT ] + N/2 + N − 1. Then holds when [zT ] > N j − 2 and k < b 1 or when for integers k, z ∈ (0, 1), j, ∈ N and N j is the length of the discrete wavelet ψ j,· for all Daubechies compactly supported wavelets. When b 1 ≤ k ≤ b 2 we have (i) for Daubechies' wavelets with two or more vanishing moments: where γ is the Euler-Mascheroni constant and (ii) for Haar wavelets we have: Proof. See Section J.3.
We use Lemma 1 to prove the next two useful results about i N .
Lemma 2. Using previous notation and assumptions, and assuming {ψ j,k } are discrete Haar wavelets Proof. See Section J.4.
These properties of the integrated local wavelet periodogram allow us to establish the asymptotic behaviour ofq W (z, τ ) in the following section.

Windowed local partial autocorrelation estimation
We now define a local partial autocorrelation estimator by using the classical (stationary) partial autocorrelation computed on a window of length L(T ) centred at time [zT ]. The theoretical properties of this windowed estimator are derived and we investigate its empirical behaviour.
Definition 3.3. Letq be the usual partial autocorrelation estimator as defined by Brockwell and Davis (1991, Definition 3.4.3) for example. Define the window I(z, L) := [z − L(T )/2T, z + L(T )/2T ] for some interval length function L(T ) and location z ∈ (0, 1). Define the windowed estimator,q W (z, τ ), of the local partial autocorrelation function at rescaled time z and lag τ , to be the classical partial autocorrelation function evaluated on observations contained in I(z, L) and denoted byq W (z, τ ) =q I(z,L) (τ ).
Our definition uses a rectangular window, but some of our applications later use an Epanechnikov window. Other variants could also be substituted.
The integrated wavelet periodogram approximation derived in Theorem 1 ensures that our windowed estimator can benefit from the established asymptotic distributional properties of the partial autocovariance estimator in the stationary setting, including its standard deviation, relevant for practical tasks.
Proposition 3.4. Let {X t,T } be a zero-mean Gaussian locally stationary wavelet process under the conditions set out by Theorem 1. Then, for the windowed local partial autocorrelation estimatorq W (z, τ ) from Definition 3.3, assuming L(T ) → ∞ and L(T )/T → 0, as T → ∞, we have thatq W (z, τ ) converges in distribution toq Y (τ ), where Y is a stationary process with the same characteristics at rescaled time z as the process {X t,T } (constructed as in equation (9)) andq Y (τ ) =φ Y τ,τ is the classical Yule-Walker partial autocorrelation function estimator.
Proof. Section J.6 contains the proof, which relies on the integrated wavelet periodogram approximation from Theorem 1.
When dealing with processes that can be locally well modelled by an autoregressive structure, the result above amounts to establishing the asymptotic normality of our windowed local partial autocovariance estimator for large lags (see next corollary).
Corollary 1. Under the assumptions from Proposition 3.4 and assuming that {X t,T } can be locally well modelled by an autoregressive structure of order say p, then for lags τ larger than p we have that L(T ) 1/2q W (z, τ ) converges in distribution to a standard normal random variable.
Proof. The proof follows directly from Proposition 3.4 and classical theory on the asymptotic behaviour of Yule-Walker estimates for stationary autoregressive processes (see for instance Theorem 8.1.2 from Brockwell and Davis (1991)).

Choice of Control Parameters
As with many nonparametric estimation methods in the literature, we have to make various choices in an attempt to obtain good estimatorsq W (z, τ ). Unfortunately there is no universal automatic best choice, at least in the real world. For the wavelet estimator,q, we have to specify an underlying wavelet, a method for handling boundaries and also a smoothing parameter, e.g. s in Section 3.3 of Nason (2013c). However, a further advantage of the windowed estimator is that we really only have to choose the window width L(T ) and the window kernel. Dahlhaus and Giraitis (1998) show that the Epanechnikov window is a good choice, which we also advocate here.
Unfortunately, rates of convergence of the estimator, although providing theoretical insight, do not really help with the practical selection of the window width. A promising direction for practical bandwidth selection might be via methods such as the locally stationary process bootstrap for pre-periodogramlike quantities, as proposed by Kreiss and Paparoditis (2014), but development of this is beyond the scope of the current paper.
Below, we use a manually-selected window width, by observing choices that achieve a good balance between estimates that are too rough, and those that appear too smooth (and change little on further smoothing). Section 4.3 and Appendix H provide some empirical evidence that the window width choice is not too hard, and the results are not particularly sensitive to it. Such manuallyselected procedures are well-acknowledged in the literature,e.g. Chaudhuri and Marron (1999), although a cross-validation method for bandwidth selection is available in our associated software at increased computational cost. This crossvalidation combines a series of dyadic cross-validations, each a simple extension of the even/odd dyadic cross-validation for wavelet shrinkage found in Nason (1996). 4 Local partial autocorrelation estimates in practice

Simulated nonstationary autoregressive examples
We illustrate our local partial autocorrelation function estimators on two simple, well-understood examples: (a) simulated time-varying autoregressive process TVAR(1) and(b) piecewise AR(p). Consider a single T = 512 realization from a time-varying autoregressive process with lag one coefficient linearly changing from 0.9 to −0.9 over the series. Figure 1 shows the partial autocorrelation function estimators, under the classical assumption of process stationarity (top left plot) and our two (timedependent) estimators (top right and bottom plots). The 95% confidence bands are constructed under the null hypothesis of white noise and are the standard ones as displayed by, e.g., established R software. The red dotted lines show the true partial autocorrelation, a linear function of time at lag 1, and constant (0) through time from lag 2 onwards.
Unsurprisingly, the classical partial autocorrelation is misleading, indicating a significant incorrect strong lag two structure, and entirely failing to detect the existing (true) lag 1 dependence. By contrast, our two developed local partial autocorrelation estimators correctly track the true time-dependent autoregressive parameters, thereby showing the importance of not using techniques designed for stationary series on nonstationary ones. Amongst our two proposals, the wavelet-based estimate seems a bit worse, particularly for the lag two partial autocorrelation after about time 350. This was confirmed by a small simulation study, based on 100 realizations drawn from the TVAR process. The average root-mean-square error for the wavelet estimator (times 10 2 , standard errors in parantheses) at lags one and two was 2.4(0.76) and 18.0(3.9), respectively, whereas for the windowed estimator it was 1.5(0.70) and 27.9(4.7) respectively. Both estimators are less accurate near the ends of the series, which is a com-  mon problem with such estimators, see Cheng and Hall (2003), for example. However, the windowed estimator usually appears less affected, and thus is the estimator we propose to use in practice. The TVAR process used in Figure 1 exhibits a large range of time-varying parameter values from −0.9 to +0.9. However, we repeated the example for less extreme parameter changes. Unless the parameter change is very small, and the process is close to stationary, the classical partial autocorrelation still misleads. For smaller parameter changes, the classical partial autocorrelation often gets the process order correct, but gives a partial autocorrelation value that is often close to the average of the local partial autocorrelations.
Our second example considers a piecewise stationary AR(p) process of length T = 256. The first and last segments (each of length 85) are realizations of an AR(1) process with φ = −0.2, and the middle segment (of length 86) follows an AR(2) process with φ = (0.5, 0.2). Note the middle segment has a significantly different structure to the first and last. Our estimators correctly identify the process structure, otherwise invisible to classical approaches. This is verified by performing a small simulation study and drawing from this process 100 times. The average root-mean-square error for the wavelet estimator at lags one and two (times 10 2 , standard errors in parentheses) is 11 (2) and 3 (2), respectively, whereas for the windowed estimator it is 7 (1) and 3 (1) respectively. The process lag 2 structure is closer to stationarity (with corresponding true pacf 0, 0.2, 0 in the three segments) and this is reflected in the similar results for the two estimators.

U.K. National Accounts data
The ABML time series obtained from the Office for National Statistics contains values of the U.K. gross-valued added (GVA), which is a component of gross domestic product (GDP). Our ABML series is recorded quarterly from quarter one 1955 to quarter three 2010 and consists of T = 223 observations. As with many economic time series, ABML exhibits a clear trend, which we removed using second-order differences; these are shown in Figure 2. Naturally, other methods for removing the trend could be tried. The second-order differences strongly suggest that the series is not second-order stationary, with the series variance increasing markedly over time. Use of methods from Nason (2013c) show that the autocorrelation also changes over time. In particular, the lag one autocorrelation undergoes a major and rapid shift around 1991.
Much of the increase in variance observed in Figure 2 is probably due to inflation. However, we have also analysed two different inflation-corrected versions of ABML, one provided by the U.K. Office of National Statistics, and both of these are also not second-order stationary, as determined by tests of stationarity in Priestley and Subba Rao (1969) and Nason (2013c).
these estimates suggest significant dependencies up to lag τ = 4. There are times, such as the 1970s, when the higher-order partial autocorrelations are not outside of the approximate significant bands, indicating that a lower lag, 2, might be appropriate. These results are (i) economically interesting as the local variance, autocorrelation and partial autocorrelation all change over time, (ii) highlight the concerns with having no access to second-order conditional information (as was the case until now) and (iii) further pose the challenge of accurately forecasting such data. Although the topic of time series forecasting is outside the scope of this paper, many authors acknowledge the superiority of wavelet-based forecasting (Aminghafari and Poggi, 2007;Schlüter and Deuschle, 2010) and we envisage the proposed local partial autocorrelation estimator could further improve results.

Precipitation in Eastport, U.S.
Understanding precipitation patterns is important for detecting climate change indications and for policy decisions. The left panel in Figure 4 shows monthly precipitation in millimetres from January 1887 until December 1950 (768 observations) at a location in Eastport. The data can be found in Hipel and McLeod (1994) and have been analysed in many publications including Rao et al. (2012);Dhakal et al. (2015). Our windowed local partial autocorrelation estimate of the Eastport data is given in the right panel of Figure 4 and shows clear nonstationarity at lags one through three. Some authors analyse this series as if it were stationary and our analysis suggests that this is inappropriate. Indeed, if one applies a formal hypothesis test of nonstationarity on appropriate lengths of the series, such as that proposed by Cardinali and Nason (2018), there is strong evidence for nonstationarity. From a modelling point of view, the estimated local partial autocorrelation behaviour might support fitting a time-varying AR(3) model. To provide some empirical support to the notion that window width is not visually critical to the interpretation of the local partial autocorrelation, Appendix H shows the smoothed local partial autocorrelation plots similar to that in the right-hand plot of Figure 4, but at three smaller window widths of 160, 80 and 40. The plot at window width of 160 is not that different to the one above at L = 250 and, indeed, the L = 80 plot is not that dissimilar. However, the L = 40 plot almost certainly contains too much 'noise' and should be disregarded.

Euro-Dollar exchange rate
Following the introduction of the Euro currency in 1999 several authors, including Ahamada and Boutahar (2002) and Garcin (2017) properties of this series, which have an influence on setting monetary policy in various jurisdictions. We analyze log returns of the monthly Euro-Dollar exchange rate as provided by EuroStat at http://ec.europa.eu/eurostat/web/products-datasets/-/ei_mfrt_m from January 1999 until October 2017. The log returns and corresponding local partial autocorrelation function estimates are given in Figure 5. This demonstrates that the log returns do not appear to be time varying (outside of the boundary locations) and exhibit only lag one partial autocorrelation. This apparent stationarity is confirmed with formal tests using the locits (Nason (2013b), Nason (2013c)) and fractal (Constantine and Percival (2016), Priestley and Subba Rao (1969)) packages in R. Interestingly this relationship holds throughout the financial crisis, from 2008 to 2011.
These examples highlight the versatility of our method and its potential use to identify stationary behaviour, manifest through local partial autocorrelation estimates that are constant through time. In addition, it highlights how the approach can identify departures from stationarity, evident through explicit time-dependent profiles at particular lags.

Discussion
This article develops two new estimators of the local partial autocorrelation function and studied their theoretical properties when applied to a locally stationary wavelet process. We established consistency for the wavelet-based estimator and asymptotic distribution for the windowed estimator. The latter result relied on new results on the integrated local wavelet periodogram, the (windowed) Haar cross-correlation wavelets and related quantities. For practical reasons, we promote the use of the windowed estimator. We demonstrated the utility of these estimators for eliciting local second-order structure on simulated data, the U.S. Eastport precipitation time series and the U.K. ABML time series. We also demonstrated the versatility of our method in the (desirable) presence of stationarity for the Euro-Dollar exchange rates. On a practical note, should the practitioner believe that higher process powers also display a locally stationary behaviour, the proposed local partial autocorrelation function could then be used to additionally uncover higher-order dependency structures. Most of the theoretical results relating to the generic local partial autocorrelation function estimator presented here are based on Haar wavelets, but many results and definitions also apply to other Daubechies' compactly supported wavelets. The associated software package lpacf contains functionality to compute the estimators for all such wavelets, up to ten vanishing moments as contained within the wavethresh package (Nason, 2013a), as well as a cross-validation method for automatic bandwidth selection. The lpacf package will be released on to the Comprehensive R Archive Network (CRAN) in due course. A Résumé: partial autocorrelation for stationary series Let {X t } t∈Z be a zero-mean second-order stationary process with autocovariance function γ(τ ). Loosely speaking, the partial autocorrelation function at lag τ is the correlation between X 1 and X τ +1 whilst adjusting for the "in-between" observations, X 2 , . . . , X τ . Brockwell and Davis (1991, p. 54) define the closed span sp(X t , t ∈ H) of any subset {X t , t ∈ H} of a Hilbert space H to be the smallest closed subspace of H which contains each X t , t ∈ H. Then, following Brockwell and Davis (1991, p. 98), the lag τ partial autocorrelation function where P 1,τ (· ) denotes the projection operator onto sp(X 2 , . . . , X τ ). See also Fan and Yao (2003, p. 43) Alternatively, if γ(0) > 0 and γ(h) → 0 as h → ∞, then the partial autocorrelation function, q(τ ), can be obtained as the final entry of the vector ϕ τ which is the solution to the well-known Yule-Walker equations Γ τ ϕ τ = γ τ . Here is a vector of covariances. Equivalently, q(τ ) = ϕ τ,τ where ϕ τ,τ is the coefficient of X 1 when projecting X τ +1 on the space spanned by X 1 , . . . , X τ , i.e. the projection For a sampled series {X t } T t=1 , the sample partial autocorrelation at lag τ is often estimated by solvingΓ τφτ =γ τ , whereγ are the usual sample autocovariances, and takingq(τ ) :=q [1,T ] (τ ) =φ τ,τ . Here we use the index notation [1, T ] in order to indicate the range of observations on which the estimation ofΓ τ andγ τ is based. The properties ofq(· ) are well-known, see Brockwell and Davis (1991, Section 8.10). In particular, T 1/2 {q(τ ) − q(τ )} has a limiting Gaussian distribution, as T → ∞, with mean zero and variance proportional to the last term on the diagonal of Γ −1 τ .

B Definitions of locally stationary processes
Definition 2.1 of Dahlhaus (1997) is as follows.
"A sequence of stochastic processes X t,T (t = 0, . . . , T − 1) is called locally stationary with transfer function A 0 and trend µ if there exists a representation where the following holds.
where cum{· · · } denotes the cumulant of kth order, g 1 = 0, is the periodic 2π extension of the Dirac delta function.
(ii) There exists a constant K and a 2π-periodic function A : for all T ; A(u, λ) and µ(u) are assumed to be continuous in u." Definition 1 of Nason et al. (2000), including improvements from Fryzlewicz (2003), is as follows.
"The locally stationary wavelet processes are a sequence of doubly indexed stochastic processes {X t,T } t=0,...,T −1 , T = 2 J ≥ 1, having the representation in the mean-square sense where ξ j,k is a random orthonormal increment sequence and where {ψ j,k (t)} j,k is a discrete non-decimated family of wavelets based on a mother wavelet ψ(t) of compact support. The quantities in (23) have the following properties: (c) There exists, for each j ≥ 1 a Lipschitz continuous function W j (z) for z ∈ (0, 1) which fulfils the following properties: the Lipschitz constants L j are uniformly bounded in j and there exists a sequence of constant C j such that for each T where

C Miscellaneous covariance matrices
The τ × τ covariance matrices are given by

D Cross-scale autocorrelation Haar wavelets
First note that by substituting s = [zT ] − t in (13) and by denoting the rectangular kernel by h, we obtain: The above expression is by no means restricted to h being a rectangular kernel, and other kernels may be used, as explained in the article main text. This new formulation of i N,z is very similar to that of the cross-correlation wavelet Ψ j, , except that the summation limits are [zT −N +1] and [zT ] instead of −∞ and ∞. This similarity is true for all Daubechies' compactly supported wavelets. We shall use this similarity to bound i N,z using Lemma 1.
Corollary 2. The regions on the right-hand side of (28) correspond exactly to the regions III a to III j in (77). For completeness (and usefulness in working out derived quantities) we can write down Ψ j, (τ ) for > j explicitly as E Subsidiary result used in the proof of Lemma 1 Lemma 3. For a, b such that 2a, 2b ∈ N and a, b > 0 it is the case that Proof. See Section J.2.1.

F Additional Results Required For the Proofs from Section 3
Lemma 4. The core function, Ω i (u) for Haar wavelets is given by for u ∈ R and i ∈ N ∪ {0}. Figure 6 shows a depiction of Ω i (u).
Lemma 5. Under the conditions and notations set out so far, for the nondecimated family of discrete Haar wavelets we have Proof. See Section J.7.2.
Lemma 6. Under the conditions and notations set out so far, for the nondecimated family of discrete Haar wavelets we have that the order of the cross terms is: Proof. See Section J.7.3.  Figure 6: Depiction of Ω i (u). The function is symmetric about 2 (i−1) − 1 2 and the extent of the function is from −1 on the left to 2 i on the right. The width of all of the triangles is always 1 for all i. As i increases the function gets stretched to the right (but also anchored on the left at u = −1), the peaks decrease in size like 2 −i/2 . Lemma 7. Under the conditions and notations set out so far, for the nondecimated family of discrete Haar wavelets we have that Proof. See Section J.7.4.

G Some Exact Formulae for Haar wavelets
Fourth-order absolute value wavelet cross-correlations for Haar wavelets. In what follows we demonstrate new results on the the fourth-order absolute value wavelet cross-correlations, for Haar wavelets which were used in showing the previous results in Appendix J.5.
The B products are symmetric in their arguments, Note that for r = 0, for ease of notation, these B (0) quantities appeared as B in the previous proofs.
Proposition G.1. For Haar wavelets. (Part A) For i, j > : (Part B) For i, j < : (Part C) For i < < j: For i < we have the following bound: (Part E) Finally, when all indices are equal we can use (34) from Nason et al. (2000) to show for > 0 and A is the matrix from Nason et al. (2000). The symmetry of B permits evaluation of B (0) (j, i) for other orderings of (i, j). An overall bound for all i, j, is B (0) (j, i) ≤ K2 −(j+i)/2 2 2 for some positive constant K.
Proof. See Section K.

H LPACF of Eastport Precipitation Data at Different Window Widths
The plots below were produced by the following functions executed using the lpacf package with binwidths of 160, 80 and 40.

I.2 Proof of Proposition 2.3
Proof. AsX Using these expressions we can rewrite q T (z, τ ) from (2) as . Now use the fact that the MSPE of a linear predictor of X t,T can be written as where b t = (b t−1,T , . . . , b 0,T , −1) T and Σ t,T is the covariance of X 0,T , . . . , X t,T ( where, as above, the τ × 1 coefficient vectors are [zT ]+τ ;T appear in Appendix C. Therefore, on combining equation (4) with (36) and (37) we obtain as desired

I.3 Proof of Proposition 2.5
Proof. The proof treats the convergence of ϕ [zT ],τ,τ ;T and the quotient that forms the square-root in equation (4) ,τ,j /T 1, where K j,i is the (j, i)th constant, and in what follows we shall seek to bound this quantity.
From the Cauchy-Schwarz inequality ϕ [zT ],τ 1 ≤ |τ | 1/2 ϕ [zT ],τ 2 = C τ as τ is fixed, and by standard properties of the spectral norm Putting parts A and B together: For the last equality the first term is asymptotically zero since or, more concisely, q (z, τ ) < ∞. The second term is O(T −1 ) as each expectation is finite. This concludes the proof.

I.4 Proof of Proposition 2.7
Proof. First recall that we defined the local partial autocorrelation as where the coefficient ϕ [zT ],τ,τ ;T is obtained in a manner akin to the (stationary) partial autocorrelation coefficient by expressing X [zT ]+τ,T as an AR(τ ) process and solving the associated Yule-Walker equations. The fraction under the square root quantifies the ratio between the backward and forward variances associated to the AR(τ ) process. The Yule-Walker equations here are localized at the rescaled time z, in the sense that they involve observations over the interval Recall that, in estimating the local partial autocorrelation, we use theĉ(z, τ ) estimator of Nason et al. (2000), which was shown there to be consistent for the (true) local autocovariance c(z, τ ). By the classical stationary theory, it follows that the estimated Yule-Walker coefficients of the AR(τ ) process (solution vector to the local Yule-Walker equations) are consistent estimators of the true coefficients, henceφ [zT ],τ,τ ;T P −→ ϕ [zT ],τ,τ , and the forward and backward variances are also estimated consistently.
Using the continuous mapping theorem (Billingsley, 1999) and assuming that the variance is non-zero, it follows that the square-root of the ratio of estimated backward and forward variances is a consistent estimator of the true ratio of variances This together with the consistency ofφ [zT ],τ,τ ;T , yields Proof. For completeness, the definition of (continuous-time) Haar wavelets is otherwise. Nason et al. (2000) show that Ψ j (τ ) = Ψ H (2 −j |τ |) where Ψ j (τ ) is the regular discrete autocorrelation wavelet and Ψ H (u) is the continuous Haar autocorrelation wavelet given by Hence, we can derive the following integral equation for Ψ j (τ ) for τ ≥ 0: from the definition of Ψ H (u). Then make the substitution x = 2 −j y to obtain: where ψ j,0 (y) = 2 −j/2 ψ H (2 −j y). Hence, by a similar argument it is the case that For Haar wavelets, since we know the precise form of ψ H we should be able to obtain an analytical formula for Ψ j, . To do this we consider < j and see that Now let x = 2 − y and we obtain Hence, it makes sense to introduce the following core function: for integers i = 0, 1, 2, . . .. Clearly, for < j. Also Ω 0 (u) = Ψ H (u). Using Lemma 4 and (39), we can now specify an exact formula for Ψ j, (τ ). For < j the result is shown in (28). Corollary 2 shows the formula for > j.
Proposition D.1 shows that the support of the cross-correlation wavelet is {k ∈ Z : −2 <= k < 2 j } for < j.
J.2 Subsidiary result used in the proof of Lemma 1

J.2.1 Proof of Lemma 3
Proof. The result is obtained by combining known results on the Fejér and Dirichlet kernels as follows. The Fejér kernel can be defined by: for ω ∈ [−π, π], see Walter and Shen (2000) Section 4.2, for example. The Fejér kernel can also be written in the following alternative form where D k (ω) = π −1 1 2 + k p=1 cos pω is the Dirichlet kernel, see Section 1.2.1 of Walter and Shen (2000).

J.3 Proof of Lemma 1
Proof. It is obvious that inequality (14) holds when i N,z (j, , k) = 0. This occurs when the lower limit in the sum (13), plus the extra k − 2[zT ] + N/2 − 1 exceeds the support of ψ ,· . In other words, i N,z (j, , k) = 0 when: It can also be shown that i N,z (j, , k) = 0 when k < [zT ] − N/2 + 1 but this inequality is not of interest in this proof . For the inequalities in (15) we decompose Ψ into three terms: Clearly, the inequality (14) is satisfied when i N,z (j, , k) = Ψ j, (k − 2[zT ] + N/2 − 1) which occurs when L = U = 0. We now investigate the conditions when L = U = 0.
(A) When is U = 0? When the lower limit of the sum defining U in (50) exceeds the support of ψ j,· , i.e.
[zT ] or when the lower limit exceeds the support of ψ ,· , i.e.
[zT ] (52) (B) When is L = 0? When the upper limit of the sum defining L in (49) is less than the lower support bound of ψ j,· , which is zero, i.e.
or when the upper limit is less than the support of ψ ,· , i.e.
Hence, U = L = 0 when inequalities (51) and (54) are satisfied. Note: we are not particularly interested in inequalities (52) and (53). For the former, the inequality (52) would have to be allied with (53) (as (54) would be contradictory to (52)) and, asymptotically (53) will not hold (as we expect the rate of increase of T to be much bigger than N ).
So far we have demonstrated the Lemma up to inequalities (15) and (16) and now we look to establish the second part of the Lemma.
To establish (17) it can be shown that, for Daubechies' wavelets with two or more vanishing moments, where H n is the nth Harmonic number, and K is a constant (maximum absolute value of the wavelet). Now using the following approximation for H n where γ is the Euler-Mascheroni constant, we can obtain the result in (17).
For Haar wavelets N j = 2 j for j ∈ N. Next we will require the discrete Fourier transform of the Haar wavelet given by: for ω ∈ (−π, π). The inverse of this transform is: for s ∈ Z. Now let us work out the precise form of the Fourier transform of the discrete Haar wavelet:ψ We now directly examine formula (27) with a rectangular kernel, as discussed in the main body of the paper. To simplify notation, we let B = [zT ] − N and r = k − 2[zT ] + N/2 − 1. In (27) replace the discrete wavelets ψ j,s and ψ ,s+r by their Fourier inverse representations given by (56) to obtain: where We now examine what happens to G j, Case-II: Suppose 0 < B + 1 but B + 1 < N j − 1, that is right-hand end of the wavelet support overlaps [B + 1, B + N ] but the left-hand end does not. Then: where Case-IV: Suppose that 0 < B + N and B + N < N j − 1, that is the left-hand end of the wavelet support overlaps [B + 1, B + N ] but the right-hand end does not. And B + 1 < 0 which is not of interest as it means that [zT ] − N + 1 < 0 which should not happen, for large T , as T increases faster than N . Now let us derive G j,B (ν) from (65) for B > 0: (66) Computing the sums in (66) gives: Now returning to the main formula (57). Case-I: suppose [0, N j − 1] ⊆ [B + 1, B + N ] then G j,B,N (ν) =ψ j (ν) as given by (62). Hence, substituting into (57) gives: We now bound i N,z (j, , k) by the integral of the absolute value of its integrand, i.e.
Case-IIb. Consider the case when N j /2 − 1 < B ≤ N j − 1. Again using (64) we have G j,B,N (ν) =ψ j (ν) − G j,B (ν) and from the corresponding value of (67), we obtain (based on the same logic as above in Cases I and IIa): We now bound i N,z (j, , k) by the integral of the absolute value of its integrand, i.e.

J.4 Proof of Lemma 2
This lemma has two parts, hence we next prove the first part.
Proof. In what follows we use the two bounds (14) for k < b 1 and k > b 2 and (18) For the first sum: For the second sum: where B (j, p) is the fourth-order cross-correlation wavelet absolute value product of order r = 0, defined as B (r) (j, i) = ∞ p=−∞ |p| r |Ψ j, (p)Ψ i, (p)| for r = 0, 1 and scales , j, i ∈ N.
Splitting the sum of j and using Proposition G.1 leads to the following: Now we prove the second part of Lemma 2.

Proof. Consider
The term T B B = k ∈B n ∈B is the case where i N,z can be bounded by Ψ and is addressed in detail in Lemma 7. The term T BB = k∈B n∈B is dealt with in Lemma 5 using the bound (18) for i N,z from Lemma 1 and the cross term is dealt with in Lemma 6. Each of these lemmas (below) show that each of the product terms is of order no worse than O{2 2( +m) }.

J.5 Proof of Theorem 1
Proof. First recall that we are under the zero-mean locally stationary wavelet process framework as described in Appendix B, with {X t,T } T −1 t=0 a doubly-index stochastic process with representation given by The integrated local periodogram was defined as where {φ j } ∞ j=1 ∈ Φ, with Φ a set of complex-valued bounded sequences equipped with uniform norm ||φ|| ∞ := sup j |φ j | and in order to avoid notational clutter N replaces the interval length notation L(T ) present in the main body of the paper. The quantity I * N (z, j) denotes the uncorrected tapered local wavelet periodogram with h : [0, 1] → R + a data taper, H N := N −1 j=0 h 2 (j/N ) ∼ N 1 0 h 2 (x) dx the normalizing factor and h(· ) is assumed symmetric and with a bounded second derivative.
As in Dahlhaus and Giraitis (1998) we approximate J N (z, φ) by the corresponding statistics of a stationary process with the same local corresponding statistics at t = zT , z fixed. Let is the wavelet periodogram on the segment [zT ] − N/2 + 1, . . . , [zT ] + N/2 of the stationary process Note: The next section uses sequences of bounded variation The total variation of a sequence {φ j } ∞ j=1 is defined by TV({φ j }) = ∞ j=1 |φ j+1 − φ j | and the space of all sequences of finite total variation is denoted by bv, see Dunford and Schwartz (1958) for example.
From equations (11) and (12), we obtain and using equation (10) it follows that which reveals the approximation we make and should be compared to equation (4.4) in Dahlhaus and Giraitis (1998), where a term O N T appears instead of O N −1 .
Using the uncorrected tapered local periodogram expression in equation (8) and the LSW definition in equation (23), by rearranging formulae we can write the integrated wavelet periodogram: and from (27) i N,z (j, , k) = Using the properties of the {ξ ,k } ,k field, we obtain where the assumption j W j (z) 2 2 2j < ∞ (at any rescaled time z) ensures that E(J Y N (φ)) is finite. Therefore, using the LSW property that sup k w 0 ,k;T − W (k/T ) ≤ C /T , leads to In order to further bound these quantities, we use Lemma 2 that proves hence we obtain using the Lipschitz continuity of {W j } j , T −1 |k − [zT ]| ∈ (0, 1) and the Hölder inequality coupled with the assumptions in the theorem. We now need to bound can be bounded using Hölder's inequality based on the assumptions in the theorem and recalling that we assumed {Y t } to be stationary.
Using Isserlis, we can decompose cov(ξ ,k ξ m,n , ξ ,k ξ m ,n ) = E(ξ ,k ξ m,n ξ ,k ξ m ,n ) Let us now expand the above Using an inequality of the type (a + b + c) 2 ≤ 3(a 2 + b 2 + c 2 ), the above quantity is upper bounded by a linear combination of a finite number of terms, of the following types In order to bound the above quantities, we use Lemma 2 which proves that Using this result we can bound each term in turn, as follows.
The first term where we have used the Lipschitz continuity of W m , | n T −z| 2 ∈ (0, 1), C 2 2 2 < ∞ and m L 2 m 2 2m < ∞. Using the same set of arguments, we bound The term as | k T − z| 2 ∈ (0, 1), L 2 2 2 < ∞ and m W 2 m (z)2 2m < ∞ at a set z (recall the process {Y t } was assumed stationary). Similarly, We therefore obtain that The result for the process {Y t } follows similarly, which concludes the proof of Theorem 1.
The elements of the covariance matrix (Γ) and vector (γ) areĉ(z, τ ) = jŜ j (z)Ψ j (τ ) and thus can be written as integrated periodograms J L(T ) , sincê Using the result in equation (75), in the manner of Dahlhaus and Giraitis (1998), it follows thatq W (z, τ ) has the same asymptotic distribution as in the stationary case.
J.7 Additional Results Required For the Proofs from Section 3 J.7.1 Proof of Lemma 4 Proof. Using the substitution y = x − u we can break the integral for Ω i into two pieces as follows: where I, II are the integrals in the final line of (76). Let us consider integral I first, making the substitution x = 2 −i (y + u) So, integral I is the result of integrating the product of ψ(x) with the moving window [2 −i u, 2 −i (u + 1 2 )]. To derive integral I we first note that I = 0 if 2 −i u ≥ 1 or 2 −i (u + 1 2 ) ≤ 0 which translates into I = 0 if u ≥ 2 i or u ≤ − 1 2 . We break down the remainder of the case − 1 2 < u < 2 i into five subregions.
In other words, we have already done the work to evaluate II(u). Hence, Now we need to put the two results together into (76) and work out the regions of overlap in the two sets of intervals. Hence, define the following regions: Putting together the two integrals with these new domains gives the result in (31).

J.7.3 Proof of Lemma 6
Proof. We bound T BB ( , m) and include the indices , m explicitly. For < m we have where Now, for n ∈ B we can apply inequality (18)  × 2 min(j,m)+min(p,m) + 2 m+min(j,m) + 2 m+min(p,m) + 2 2m . Now let us examine the sum over k in (80) using the Ψ bound for i N,z from (14).

J.7.4 Proof of Lemma 7
Proof. Using the definition of V from (81),

K Proof of Proposition G.1
Proof. Part A: for i, j > . We can work out the exact formula for B (j, j) by direct application of the formula (28) for Ψ j, (p) for j > .
Hence, after some algebra, Next, we examine: for i = j + 1. Examining Table 1 shows that Range I for both the j and j + 1 cross-correlation wavelets always overlap and Range II for i = j + 1 overlaps with Range III for j. Hence: Finally for |i − j| > 1 we examine: Part B: for the case i, j < . First, for j = i we have where L = j and J = and we used Ψ j, (p) = Ψ ,j (−p) since we will use the formula for Ψ ,j where > j in (28). Since L = j < = J this puts us into the situation of (83) which gives: B (0) (j, j) = 2 −J (2 2L−1 + 1) = 2 − (2 2j−1 + 1), as required. For i < j < we use the formula for Ψ for j < given by (29). Due to the form of Ψ j, (p) for i < j < we can split the sum into three parts corresponding to the the three non-zero parts of the autocorrelation wavelet given in (29). The condition i < j < is helpful as the interval associated with i nests within that of j, and the j interval nests within that associated with . First, we will deal with the last two terms of (29) which do not depend on ("front part").

front =
Note how the first term in each of these sums is | − 2 −j p| because this is the formula for Ψ j, (p) over the interval [1, 2 i ] because this interval is always contained within [1, 2 j−1 ] since i < j. We also ensure that the last absolute value term in the sum is positive (if it was negative then we'd switch signs of the contents of the absolute value as |x| = −x if x is negative). Continuing The middle term is constructed in a similar way except that enters into the equation. However, the concept that the j interval can only ever overlap the i interval on its first half still remains. Hence, The back part uses precisely the same part as the middle back = Adding the front, middle and back components together and multiplying through by the constant that appears at the front of (29) gives: B (0) (j, i) = 2 −( −j)/2 2 −( −i)/2 (2 −j 2 −2 2 2i−1 + 2 −j 2 2i−1 + 2 −2 2 −j 2 2i−1 ) = 2 − 2 i/2 2 j/2 × 2 −j 2 2i−1 (2 −2 + 1 + 2 −2 ) = 3 2 2 − 2 −j/2 2 5i/2−1 , as required. Part C: for the case i < < j. When considering the formula for B (0) (j, i) we use formula (28) for Ψ j, (p) because < j and formula (29) for Ψ i, (p) because i < .
Let us first consider positive p first. In this case, the only values of Ψ i, (p) which are nonzero are for 0 < p ≤ 2 i and the only values of Ψ j, (p) are for 2 j−1 − 2 ≤ p. So, if 2 i < 2 j−1 − 2 then there is no overlap between these two cross-correlation wavelets, under what conditions can this occur (subject to i < < j)? Suppose j = + 1 then 2 j−1 − 2 = 0 and 2 i is never less than or equal to zero. Suppose j = + 2, then 2 j−1 − 2 = 2 which is always greater than 2 i because we know that > i. Hence, we have to consider two cases (i) j = + 1 where the positive parts of the cross-correlation functions overlap and j > + 1 where they do not.
We will now work out this 'overlap' contribution where j = + 1: The limits of the sum correspond to where the term Ψ i, (p) is non-zero for p ≥ 0. Now we consider which parts of Ψ +1, (p) are relevant to the sum. With j = + 1 the fifth and sixth ranges in (28) respectively. Since i < we have 2 i < 2 and so the two ranges in (90) and (91) are the only ones we need to consider for Ψ +1, (p) for positive p. Now the highest that i can be is i = − 1 so that maximum p is 2 i = 2 −1 so, in actuality, it is only range (90) that is active for any p > 0. Hence, without the constant 2 −(j− )/2 2 −( −i)/2 = 2 −j/2 2 i/2 at the front of (28) and (29), we can write Note that 1 − 2 −i p is positive on the range of p in the second sum. Hence, for i < and j = + 1.
For p < 0, i < < j the cross-correlation wavelet Ψ j, (p) splits into two parts random each of those two parts the cross-correlation wavelet Ψ i, (p) splits into two further parts. Hence, the sum on the p < 0 terms (the 'back' bit) is given by: We can now split each of these two sums into the two sets of two ranges dictated by the domain of Ψ i, (p)| as follows: For all of the eight components in the four sums in (95) all except the fourth, fifth and seventh are always positive over their respective sum's range of p values. Hence, we replace those terms in the absolute values by their negative (e.g. |x| = −x for x < 0) and obtain: (1 + 2 − p)(p + 2 ) (1 + 2 − p)(2 i − p − 2 ) Now multiplying by the 2 −j/2 2 i/2 which we omitted earlier gives: as required for the second equation in (32) for i < and + 1 < j.
For the first equation in (32) we have to add 2 −j/2 2 i/2 times equation (93) to obtain:
Part D: First consider j = < i. Then: where Ψ (p) is the ordinary autocorrelation wavelet from Nason et al. (2000). The domain of Ψ (p) is from −2 < p2 . The range of Ψ (p) and Ψ i, (p) agrees for p < 0 but the two wavelets only overlap for p > 0 if the lower end of the nonzero range of Ψ i, (p), namely 2 i−1 − 2 is smaller than the upper end of the autocorrelation wavelet 2 , recalling that i > . This only occurs when i = + 1 and for i ≥ + 2 there is no overlap for p > 0.
For the case of i < we use the following bound: as |Ψ (p)| takes its maximum value of 1 at p = 0. We can work out this sum directly taking care to discover when the argument inside the absolute value sign is negative or positive giving: (2 i − 2 − p) Hence as required.