# Asymptotically Sufficient Statistics in Nonparametric Regression Experiments with Correlated Noise

Authors: Andrew V. Carter
Journal: Journal of Probability and Statistics
 (2009)

---

## Abstract

We find asymptotically sufficient statistics that could help
simplify inference in nonparametric regression problems with correlated errors. These statistics are derived from a wavelet decomposition that is used to whiten the noise process and to effectively separate high-resolution and low-resolution components. The lower-resolution components contain nearly all the available information about the mean function, and the higher-resolution components can be used to estimate the error
covariances. The strength of the correlation among the errors is related to the speed at which the variance of the higher-resolution components shrinks, and this is considered an additional nuisance parameter in the model. We show that the NPR experiment with correlated noise is asymptotically equivalent to an experiment that observes the mean function in the presence of a continuous Gaussian process that is similar to a fractional
Brownian motion. These results provide a theoretical motivation for some commonly proposed wavelet estimation techniques.


---

## Body

## 1. Introduction

A nonparametric regression (NPR) problem consists of estimating an unknown mean function that smoothly changes between observations at different design points. There aren observations Yi of the form(1.1)Yi=μ(in)+ξifori=1,…,n,
where μ is the unknown smooth mean function on [0,1] and the errors ξi are observations from a zero-mean Gaussian process. For NPR problems that have a particular long memory structure to the covariance of the error terms, we will find a continuous Gaussian experiment approximation to the problem of estimating the mean.Brown and Low [1] showed that the NPR experiment is asymptotically equivalent to the white-noise model where the mean function is observed in the presence of a Brownian motion process. This result paralleled the work by Nussbaum [2] in showing that asymptotic results in nonparametric function estimation problems can be simplified using approximations by the continuous white-noise experiments that Pinsker [3] studied. The original asymptotic equivalence results for NPR experiments were extended by Brown et al. [4] and Carter [5, 6] along with refinements in the approximations from Rohde [7] and Reiss [8].All of these results assume that the errorsξi in (1.1) are all independent, and this assumption is critical in establishing the appropriateness of a white-noise model that also has independent increments. We want to consider the effect of correlation between the observations on these approximations. Presumably, if the correlation is weak, then the effect washes out asymptotically. However, we wish to consider cases where there is sufficient long-range correlation to affect the form of the approximation. In particular, we will show that the appropriate approximation is by a continuous Gaussian process experiment that is no longer white noise but is closer to a fractional Brownian motion.Our approach is motivated by the work by Johnstone and Silverman [9] and Johnstone [10]. They investigated the wavelet decomposition of data of this type and used a fractional Brownian motion approximation in the limit:(1.2)dY(t)=μ(t)dt+n-(β+1)/2dBK(t)t∈[0,1].
They argued that the wavelet decomposition resulted in nearly independent coefficients which simplified the inference significantly. We will assume that the BK(t) process is decorrelated by a wavelet decomposition and then show that this continuous model is asymptotically equivalent to the NPR experiment with the same covariance structure.Theorem 1.1.
The nonparametric regression experimentℱ observes Yi as in (1.1) for an unknown mean function μ from a parameter set ℳ(M,α) defined in Section 1.2 and a known covariance structure as described in Section 1.3. This experiment is asymptotically equivalent to the experiment ℰ that observes
(1.3)dY(t)=μ(t)dt+σn-(β+1)/2dBK(t),
where BK(t) is a Brownian motion with covariance kernel K.This will be proven in two steps. First, Lemma2.1 proves that the first n wavelet coefficients in a decomposition of dY(t) are asymptotically sufficient in ℰ for estimating μ. For the second step, Lemma 3.1 shows that a discrete wavelet transform of the observations from ℱ produces observations with nearly the same distribution as these asymptotically sufficient statistics.Furthermore, in both experiments the lower-frequency terms in the wavelet decomposition are sufficient for estimating the means, allowing the higher-frequency terms to be used to give information about the variance process. This leads to Theorem1.2, which proposes an experiment that allows some flexibility in the error structure.Theorem 1.2.
The NPR experimentℱ̃ observes the Yi as in (1.1), where the covariance structure depends on the parameters β and γ and is such that the variance of the wavelet coefficients is 2γ+β(j+1).
The experimentℰ̂ observes the pair
(1.4)(γ̂β̂)~𝒩((γβ),log22(x⊤Λx)-1),
(where x and Λ are defined in Section 4) and then observes the continuous Gaussian process conditionally on γ̂ and β̂:
(1.5)dY(t)=μ(t)dt+n-(1+β̃)2γ̃dBK̂(t),
where the covariance K̂ is such that Var(BK̂(ψjk))=2β̃(j+1). The estimators β̃ and γ̃ are the same as (β̂,γ̂) but truncated so that -1≤β̃≤0, and γ̃≥-c.
Forμ∈ℳ(M,α), -1<β≤0, and γ>-c for some constant c, the experiments ℱ̃ and ℰ̂ are asymptotically equivalent.This theorem can be seen as an extension of Carter [6] Theorem 1.1 from a case where there is a single unknown variance for all the wavelet coefficients to a case where the variance changes as a log-linear function of the resolution level (or frequency).Wang [11] addressed the issue of asymptotically sufficient statistics in the fractional Brownian motion process. In Section 3 of that article there is an argument that bounds the difference between minimax errors in an NPR experiment with correlated errors and an experiment that observes the mean in the presence of fractional Brownian motion error. This result extends the sort of approximation by Donoho and Johnstone [12] to correlated errors and is very much in the spirit of our Theorem 1.1 here. Our results differ from Wang [11] in that we have made a stronger assumption on the covariance structure of the errors in order to obtain the full asymptotic equivalence of the experiments as discussed in Section 1.1.Lemma2.1 is presented and proven in Section 2. Section 3 presents Lemma 3.1 and the proof of Theorem 1.1. The proof for Theorem 1.2 is in Section 4 with some relevant bounds in Sections 5 and 6.


### 1.1. Asymptotic Sufficiency

Instead of focusing on single estimation techniques, we will consider approximations of the entire statistical experiment. For large sample sizes, there is often a simpler statistical experiment that can approximate the problem at hand. One benefit of finding an approximating experiment is that it may have convenient sufficient statistics even when they are not available in the original experiment.Our approximations will therefore be ofexperiments rather than particular distributions. A statistical experiment 𝒫 that observes data X consists of a set of probability distributions {ℙθ} indexed by the parameter set θ∈Θ. We wish to compare the information about θ in 𝒫 to another experiment 𝒬 that observes data Y from among the set of distributions {ℚθ} that are indexed by the same parameter θ. Implicitly, we are concerned with two sequences of experiments 𝒫n and 𝒬n where n roughly denotes the increasing sample size, but generally, we will leave off the subscript n. It will always be understood that the distributions depend on the “sample size.”The NPR experiment will be approximated using Le Cam's notion of asymptotically equivalent experiments [13, 14] and asymptotically sufficient statistics [15]. Asymptotically equivalent experiments have corresponding inference procedures (such as estimators or tests) in each experiment that perform nearly as well. Specifically, if there is an estimator τ(X) in 𝒫 with risk ℙθL(τ(X)), then, for bounded loss functions, there is an estimator σ(Y) such that(1.6)supθℙθL(τ(X))-ℚθL(σ(Y))→0
as n→∞. These asymptotic equivalence results are stronger than the equivalence of minimax rates that is derived under a similar model by, for example, Wang [11]. Our results imply a correspondence over a range of bounded loss functions. Thus, the equivalence holds for a global L2 error as well as local error measurements or other distances.Asymptotic sufficiency is a stronger notion, where ifT(X) is a sufficient statistic for inference about θ in 𝒫, then T(Y) is asymptotically sufficient for 𝒬 when the total-variation distance between ℙθ and ℚθ is negligible. These asymptotically sufficient statistics generate experiments that are all asymptotically equivalent. In particular, 𝒫 and 𝒬 are asymptotically equivalent, and they are also asymptotically equivalent to the experiments generated by the distributions of T(X) and T(Y). As a result, an estimator in 𝒫 should generally be of the form τ(T(X)) and there is a corresponding estimator τ(T(Y)) that performs nearly as well in the 𝒬 experiment. There is a basic transitive property to the asymptotic equivalence that implies if 𝒫 is asymptotically equivalent to 𝒬, and 𝒬 is asymptotically equivalent to ℛ, then 𝒫 is asymptotically equivalent to ℛ.Le Cam's asymptotic equivalence is characterized using the total-variation distanceδ(ℙθ,ℚθ) between the distributions. We will abuse this notation a bit by writing δ(𝒫,𝒬)=supθδ(ℙθ,ℚθ). It will often be more convenient to use the Kullback–Leibler divergence (D(ℙ,ℚ)=ℙlog[dℙ/dℚ]) to bound the total variation distance:(1.7)δ(ℙθ,ℚθ)≤2D(ℙθ,ℚθ)1/2
See [16]. The divergence is convenient for product experiments because D(∏ℙi,∏iℚi)=∑iD(ℙi,ℚi).


### 1.2. Wavelet Basis

We will use orthonormal wavelet bases to characterize the function space and to simplify the covariance structure of the errors.Assuming we are considering periodic functions on the interval[0,1], we can construct periodic wavelet bases as by Daubechies [17, Chapter 9.3]. We start with a space 𝒱j which consists of functions of the form(1.8)fj(t)=∑k=12jajkϕjk(t),t∈[0,1],
where ϕjk is an orthonormal set of periodic functions generated via(1.9)ϕjk(t)=2j/2ϕ(2jt-k).
We will work with a ϕ function having finite support [0,N], and at the boundaries of the interval the ϕjk(t) are given the proper periodic extensions. This space generates wavelet functions ψjk that span the difference between the 𝒱j and 𝒱j-1 and can be written ψjk(t)=2j/2ψ(2jt-k) with the proper periodic adjustment at the boundary (e.g., ψj,2j(ϵ)=2j/2ψ(2jϵ) for a small ϵ). This periodic adjustment has a small effect at the high resolution levels but is a larger factor for small values of j. In particular, the scaling function at level 0 is ϕ0(t)=∑k=0Nϕ(k+t)=1.The mean functionsμ(t) will be assumed to be constructed from this wavelet basis:(1.10)μ(t)=θ0ϕ0(t)+∑j≥0∑k=12jθjkψjk(t).
We will restrict the mean functions to those that belong to a Hölder(α) class of functions. Specifically, the class of periodic mean functions μ(t) is ℳ(M,α):(1.11)supℳ(M,α)supt,s|μ(t)-μ(s)||t-s|α≤M
for some 1/2<α<1 and M>0. This smoothness condition on the functions bounds the rate of growth of the higher-frequency terms in the orthonormal expansion. Originally from Meyer [18], in Daubechies [17, page 299 equation, equation (9.2.1)] gives the bound(1.12)maxk|θjk|≤M22-j(α+1/2)
for a constant M2 related to M. This bounds implies the useful bound on the squared error in the tail of the basis expansion:(1.13)∑j>j*2j∑k=12jθjk2≤M222-εj*(1-2-ε)-1,
where ε=2α-1.


### 1.3. Error Structure

These results rely on a specific structure to the covariance matrix of the errors in the NPR experiment. As by Johnstone [10], the fractional Brownian motion is the motivating example for our continuous Gaussian model. However, this model does not necessarily provide the independent coefficients that would simplify the inference. Instead, an error structure that has roughly some of the properties of the fractional Brownian motion will be considered.Traditionally, the asymptotics of the NPR experiment have assumed independent noise. This white-noise model is especially convenient because all of the eigenvalues of the covariance operator are equal. Thus, any orthonormal basis generates a set of independent standard normal coefficients. With a more general covariance function, the eigenvalues are different and only particular decompositions lead to independent coefficients. Thus there is much less flexibility in the choice of basis, and this basis determines some of the structure of the covariance.Following Johnstone [10], Johnstone and Silverman [9], Zhang and Waiter [19], Wang [11], and Cavalier [20] among others, we will assume a covariance structure that is whitened by a wavelet decomposition. When there is a long-range positive correlation between the observations, the wavelet decomposition tends to decorrelate the error process because the wavelet functions act like bandpass filters.We will assume that there exists an orthonormal basisϕ0 and ψjk for j≥0 and k=0,…,2j-1 such that the decomposition of the error process generates independent normal coefficients. In other words, the error process is a zero-mean Gaussian process that is roughly(1.14)dBK(t)=ξ0ϕ0(t)+∑j≥0‍∑kξjkψjk(t)
in the distributional sense where the ξkj are independent normals. The Var(ξjk) will be assumed to depend on j and not k as a sort of stationarity condition. In particular, we will assume that Var(ξ0)=σ2 and then Var(ξjk)=σ22β(j+1) for some β in the interval (-1,0]. If β=0, then this is the white-noise process.This is a convenient form for the error, but not completely unrealistic. Wavelet decompositions nearly whiten the fractional Brownian motion process. Wornell [21] argued that long-memory processes can be constructed via a wavelet basis with variances at resolution level j shrinking like 2-γj for 0<γ<2. McCoy and Walden [22] showed that the discrete wavelet transform nearly decorrelates the noise in fractionally differenced white-noise processes. Alternatively, Wang [11] used a wavelet—vaguelette decomposition [23] to find a decomposition of the fractional Brownian motion that results in independent coefficients for a nearly orthonormal basis.Section7 demonstrates some properties of the specific Gaussian process generated by using the Haar basis as the wavelet basis. These properties are consistent with the sort of behavior that we want in the covariances of our observations. The correlation between observations decreases like d-(β+1) for β<0 where d measures the distance between the locations of the coefficients.A well-established method for estimating the parameterβ in these long-range dependent models is to fit a linear function to the log of an estimate of the variances of the coefficients at each resolution level. This idea goes back to at least Abry and Veitch [24] and is now a standard approach that has been improved upon in subsequent work; see Veitch and Abry [25], Stoev et al. [26], among others. This motivates the asymptotic sufficient statistics in Theorem 1.2 which are least squares estimates from the fitted line.The assumptions in Theorem1.1 on the covariance structure of the errors are strong and could limit the applicability of the result. However, if we allow the variances at different scales to have a range of linear relationship, we could then have a sufficiently rich class of error models. Theorem 1.2 allows for this somewhat larger class of models, and it seems likely that the changing magnitude of the variances over different resolutions level will have a greater effect on the distribution of the observed errors than the underlying basis.


## 1.1. Asymptotic Sufficiency

Instead of focusing on single estimation techniques, we will consider approximations of the entire statistical experiment. For large sample sizes, there is often a simpler statistical experiment that can approximate the problem at hand. One benefit of finding an approximating experiment is that it may have convenient sufficient statistics even when they are not available in the original experiment.Our approximations will therefore be ofexperiments rather than particular distributions. A statistical experiment 𝒫 that observes data X consists of a set of probability distributions {ℙθ} indexed by the parameter set θ∈Θ. We wish to compare the information about θ in 𝒫 to another experiment 𝒬 that observes data Y from among the set of distributions {ℚθ} that are indexed by the same parameter θ. Implicitly, we are concerned with two sequences of experiments 𝒫n and 𝒬n where n roughly denotes the increasing sample size, but generally, we will leave off the subscript n. It will always be understood that the distributions depend on the “sample size.”The NPR experiment will be approximated using Le Cam's notion of asymptotically equivalent experiments [13, 14] and asymptotically sufficient statistics [15]. Asymptotically equivalent experiments have corresponding inference procedures (such as estimators or tests) in each experiment that perform nearly as well. Specifically, if there is an estimator τ(X) in 𝒫 with risk ℙθL(τ(X)), then, for bounded loss functions, there is an estimator σ(Y) such that(1.6)supθℙθL(τ(X))-ℚθL(σ(Y))→0
as n→∞. These asymptotic equivalence results are stronger than the equivalence of minimax rates that is derived under a similar model by, for example, Wang [11]. Our results imply a correspondence over a range of bounded loss functions. Thus, the equivalence holds for a global L2 error as well as local error measurements or other distances.Asymptotic sufficiency is a stronger notion, where ifT(X) is a sufficient statistic for inference about θ in 𝒫, then T(Y) is asymptotically sufficient for 𝒬 when the total-variation distance between ℙθ and ℚθ is negligible. These asymptotically sufficient statistics generate experiments that are all asymptotically equivalent. In particular, 𝒫 and 𝒬 are asymptotically equivalent, and they are also asymptotically equivalent to the experiments generated by the distributions of T(X) and T(Y). As a result, an estimator in 𝒫 should generally be of the form τ(T(X)) and there is a corresponding estimator τ(T(Y)) that performs nearly as well in the 𝒬 experiment. There is a basic transitive property to the asymptotic equivalence that implies if 𝒫 is asymptotically equivalent to 𝒬, and 𝒬 is asymptotically equivalent to ℛ, then 𝒫 is asymptotically equivalent to ℛ.Le Cam's asymptotic equivalence is characterized using the total-variation distanceδ(ℙθ,ℚθ) between the distributions. We will abuse this notation a bit by writing δ(𝒫,𝒬)=supθδ(ℙθ,ℚθ). It will often be more convenient to use the Kullback–Leibler divergence (D(ℙ,ℚ)=ℙlog[dℙ/dℚ]) to bound the total variation distance:(1.7)δ(ℙθ,ℚθ)≤2D(ℙθ,ℚθ)1/2
See [16]. The divergence is convenient for product experiments because D(∏ℙi,∏iℚi)=∑iD(ℙi,ℚi).


## 1.2. Wavelet Basis

We will use orthonormal wavelet bases to characterize the function space and to simplify the covariance structure of the errors.Assuming we are considering periodic functions on the interval[0,1], we can construct periodic wavelet bases as by Daubechies [17, Chapter 9.3]. We start with a space 𝒱j which consists of functions of the form(1.8)fj(t)=∑k=12jajkϕjk(t),t∈[0,1],
where ϕjk is an orthonormal set of periodic functions generated via(1.9)ϕjk(t)=2j/2ϕ(2jt-k).
We will work with a ϕ function having finite support [0,N], and at the boundaries of the interval the ϕjk(t) are given the proper periodic extensions. This space generates wavelet functions ψjk that span the difference between the 𝒱j and 𝒱j-1 and can be written ψjk(t)=2j/2ψ(2jt-k) with the proper periodic adjustment at the boundary (e.g., ψj,2j(ϵ)=2j/2ψ(2jϵ) for a small ϵ). This periodic adjustment has a small effect at the high resolution levels but is a larger factor for small values of j. In particular, the scaling function at level 0 is ϕ0(t)=∑k=0Nϕ(k+t)=1.The mean functionsμ(t) will be assumed to be constructed from this wavelet basis:(1.10)μ(t)=θ0ϕ0(t)+∑j≥0∑k=12jθjkψjk(t).
We will restrict the mean functions to those that belong to a Hölder(α) class of functions. Specifically, the class of periodic mean functions μ(t) is ℳ(M,α):(1.11)supℳ(M,α)supt,s|μ(t)-μ(s)||t-s|α≤M
for some 1/2<α<1 and M>0. This smoothness condition on the functions bounds the rate of growth of the higher-frequency terms in the orthonormal expansion. Originally from Meyer [18], in Daubechies [17, page 299 equation, equation (9.2.1)] gives the bound(1.12)maxk|θjk|≤M22-j(α+1/2)
for a constant M2 related to M. This bounds implies the useful bound on the squared error in the tail of the basis expansion:(1.13)∑j>j*2j∑k=12jθjk2≤M222-εj*(1-2-ε)-1,
where ε=2α-1.


## 1.3. Error Structure

These results rely on a specific structure to the covariance matrix of the errors in the NPR experiment. As by Johnstone [10], the fractional Brownian motion is the motivating example for our continuous Gaussian model. However, this model does not necessarily provide the independent coefficients that would simplify the inference. Instead, an error structure that has roughly some of the properties of the fractional Brownian motion will be considered.Traditionally, the asymptotics of the NPR experiment have assumed independent noise. This white-noise model is especially convenient because all of the eigenvalues of the covariance operator are equal. Thus, any orthonormal basis generates a set of independent standard normal coefficients. With a more general covariance function, the eigenvalues are different and only particular decompositions lead to independent coefficients. Thus there is much less flexibility in the choice of basis, and this basis determines some of the structure of the covariance.Following Johnstone [10], Johnstone and Silverman [9], Zhang and Waiter [19], Wang [11], and Cavalier [20] among others, we will assume a covariance structure that is whitened by a wavelet decomposition. When there is a long-range positive correlation between the observations, the wavelet decomposition tends to decorrelate the error process because the wavelet functions act like bandpass filters.We will assume that there exists an orthonormal basisϕ0 and ψjk for j≥0 and k=0,…,2j-1 such that the decomposition of the error process generates independent normal coefficients. In other words, the error process is a zero-mean Gaussian process that is roughly(1.14)dBK(t)=ξ0ϕ0(t)+∑j≥0‍∑kξjkψjk(t)
in the distributional sense where the ξkj are independent normals. The Var(ξjk) will be assumed to depend on j and not k as a sort of stationarity condition. In particular, we will assume that Var(ξ0)=σ2 and then Var(ξjk)=σ22β(j+1) for some β in the interval (-1,0]. If β=0, then this is the white-noise process.This is a convenient form for the error, but not completely unrealistic. Wavelet decompositions nearly whiten the fractional Brownian motion process. Wornell [21] argued that long-memory processes can be constructed via a wavelet basis with variances at resolution level j shrinking like 2-γj for 0<γ<2. McCoy and Walden [22] showed that the discrete wavelet transform nearly decorrelates the noise in fractionally differenced white-noise processes. Alternatively, Wang [11] used a wavelet—vaguelette decomposition [23] to find a decomposition of the fractional Brownian motion that results in independent coefficients for a nearly orthonormal basis.Section7 demonstrates some properties of the specific Gaussian process generated by using the Haar basis as the wavelet basis. These properties are consistent with the sort of behavior that we want in the covariances of our observations. The correlation between observations decreases like d-(β+1) for β<0 where d measures the distance between the locations of the coefficients.A well-established method for estimating the parameterβ in these long-range dependent models is to fit a linear function to the log of an estimate of the variances of the coefficients at each resolution level. This idea goes back to at least Abry and Veitch [24] and is now a standard approach that has been improved upon in subsequent work; see Veitch and Abry [25], Stoev et al. [26], among others. This motivates the asymptotic sufficient statistics in Theorem 1.2 which are least squares estimates from the fitted line.The assumptions in Theorem1.1 on the covariance structure of the errors are strong and could limit the applicability of the result. However, if we allow the variances at different scales to have a range of linear relationship, we could then have a sufficiently rich class of error models. Theorem 1.2 allows for this somewhat larger class of models, and it seems likely that the changing magnitude of the variances over different resolutions level will have a greater effect on the distribution of the observed errors than the underlying basis.


## 2. Approximate Sufficiency in the Gaussian Sequence Experiment

The first step in the proof of Theorem1.1 is to establish that a truncated wavelet decomposition is asymptotically sufficient for the continuous Gaussian experiment.Lemma 2.1.
The experimentℰ is a Gaussian process experiment that observes
(2.1)dY(t)=μ(t)dt+σn-(β+1)/2dBK(t),
where BK is a zero-mean continuous Gaussian process with covariance K. There are asymptotically sufficient statistics for estimating μ∈ℳ(M,α):
(2.2)y0,yjk~𝒩(θjk,σ2n-β-12βj)
for 0≤j≤j* as long as
(2.3)j*>log2n2α.In the Gaussian sequence experimentℰ where only the mean μ(t) is to be estimated, the likelihood is(2.4)dℙμdℙ0(y)=exp[2y0θ0-θ022σ2+σ-2n(1+β)∑j≥0σ-22-β(j+1)∑k(θjkyjk-12θkj2)],
where ℙμ is the distribution of Y(t) and ℙ0 is the distribution of the version with mean 0 which would just be σn-(β+1)/2BK(t).We want to approximate this experimentℰ by a similar experiment ℱ where the mean is projected onto the first j* resolution levels, that is; μ is replaced by μ̅:(2.5)μ̅(t)=θ0ϕ0(t)+∑j=0j*∑kθjkψjk(t).
The likelihood becomes(2.6)dℚμdℚ0(y)=exp[2y0θ0-θ022σ2+σ-2n(1+β)∑j=0j*2-βj∑k(θjkyjk-12θkj2)].
Therefore, this experiment ℱ has sufficient statistics y0 and yjk for 0≤j≤j*. These observations are approximately sufficient in the ℰ experiment if the distance between the distributions in the two experiments is small.By (5.3), the distance between these two sets of experiments is(2.7)δ(ℙ,ℚ)≤σ-1n(β+1)/2(∑j>J2-βj∑kθjk2)1/2.
For the parameter space ℳ(M,α), (5.4) bounds the distance between the two experiments as(2.8)δ(ℙ,ℚ)≤σ-1n(β+1)/2M2-εJ≤Mσ-12(J-j*)(1+β)/2-εj*,
which is negligible as n→∞ when the dimension of the sufficient statistics increases like(2.9)j*>(1+ββ+2α)log2n
for -1<β<0. The worst case is when β=0, and thus we have proved Lemma 2.1.


## 3. Approximating the NPR Experiment

Theorem1.1 can now be proven by approximating the sufficient statistics from Lemma 2.1 using the observations from the NPR experiment ℱ.We suppose that we haven observations from the NPR experiment as in (1.1) where the ξi are Gaussian random variables with a specified covariance function. Specifically, let W be the n×n matrix that performs the discrete wavelet transform, and let W⊤ be its inverse. The vector of random wavelet coefficients from Lemma 2.1, y=(y0,y00,y01,…,yJ-1,2J-1-1)⊤ where J=log2n, can be transformed via the discrete wavelet transform to create ỹJ=W⊤y.The expected value of this transformed vector is(3.1)𝔼ỹJi=μ̃Ji=∫2J/2ϕ(2Jt-i)μ(t)dt.
For a μ(t) function that is nearly constant around i/n, μ̃Ji≈2-J/2μ(i/n); so we can approximate ỹJ by (1/n)(Y1,Y2,…,YN)⊤.In the original NPR experiment, the variances areVar(Yi)=Var(ξi)=Cσ2 for a constant C that depends on β and the basis we will be using. The covariance matrix for these Yi will be assumed to be Σ=WDW⊤ where D is the diagonal matrix of Var(yjk)=σn22β(j+1). The variance of the ỹJk should be the same as that of Yin-1/2, and in the model described, Var(ỹJk)∝σn2nβ. Therefore, σn2 should be set to σ2n-1-β.Lemma 3.1.
If the mean functionμ(t) is in ℳ(M,α), then the total variation distance between the distribution ℙμ of the vector (1/n)(Y1,Y2,…,YN)⊤ and the distribution ℙ̃μ of the vector ỹJ is
(3.2)supμ∈ℳ(M,α)δ(ℙμ,ℙ̃μ)≤(MNα+1α+1)n1/2-α+β/2.This lemma essentially goes back to the original work of Mallat [27] and parallels some of what was done by Brown and Low [1], Johnstone [10], Rohde [7], and Carter [5].The NPR observations are such that the covariance matrix ofξi is also Σ, and therefore the total variation distance between the distributions is bounded in (5.2) by(3.3)δ(ℙ,ℙ̃)≤12πΔ1/2
with(3.4)Δ=(μ(in)-μ̃Ji)⊤Σ-1(μ(in)-μ̃Ji).A standard calculation bounds the difference between the means whenϕ(t)<M with support on [0,N] and the μ(t) are Hölder (α) for α<1:(3.5)|μ(in)-μ̃Ji|≤n-1/2-αMα+1Nα+1.The covariance matrix is a positive definite matrix such thatΣ-1=MD-1M⊤. The first column of the wavelet transform matrix is n1 where 1 is the vector of 1's. Therefore,(3.6)Δ≤(MNα+1α+1)n-1-2αnn1+β=(MNα+1α+1)n1-2α+β,
which is negligible for large n and α>1/2.


### 3.1. Proof of Theorem1.1

The theorem follows from the fact that the observationsy0,{yjk} for j=0,…,J-1 are asymptotically sufficient for the continuous process in (1.2). Then a linear function of these sufficient statistics ỹ=W⊤y is still approximately sufficient. Thus, the experiment that seeks to draw inference about μ from the observations ỹJi is asymptotically equivalent to the experiment that observes (1.2) by Lemma 2.1.Furthermore, by Lemma3.1, the original NPR experiment that has the same covariance structure as ỹJi is asymptotically equivalent to that experiment and thus, by transitivity, to the experiment that observes the process (1.2) as well. This proves Theorem 1.1.


### 3.2. Remarks on the Covariance Structure

This result is restrictive in that it requires a specific known covariance structure. We are working under the assumption that the covariance matrix has eigenfunctions that correspond to a wavelet basis. This does not generally lead to a typical covariance structure. It does not even necessarily lead to a stationary Gaussian process; see the Haar basis example below.The difficulty is that the requirement for having asymptotically equivalent experiments is quite strict, and the total variation distance between the processes with even small differences in the structure of the covariance is not negligible. For two multivariate Guassian distributions with the same means but where one covariance matrix isΣ and the other is D, a diagonal matrix with the same diagonal elements as Σ, the Kullback–Leibler divergence between the distributions is log|Σ|-log|D|.If the correlation between the highest level coefficients isCorr(ξj*,k,ξj*,k+1)=ρ, then the contribution to the difference of the log determinants is on the order of ρ2j*. The dimension of the problem is growing while the correlations are generally not going to 0 significantly quickly. For instance, in a typical wavelet basis decomposition of the true fractional Brownian motion Corr(ξj*,k,ξj*,k+1)=cβ where cβ is a constant that depends on β but not j* or n.Thus, the differencelog|Σ|-log|D| will not go to 0 as the sample size increases. Therefore, for the sort of long-range correlation structures that we are considering here, the eigenfunctions of the kernel K need to be known or else the experiments will not be asymptotically equivalent.


## 3.1. Proof of Theorem1.1

The theorem follows from the fact that the observationsy0,{yjk} for j=0,…,J-1 are asymptotically sufficient for the continuous process in (1.2). Then a linear function of these sufficient statistics ỹ=W⊤y is still approximately sufficient. Thus, the experiment that seeks to draw inference about μ from the observations ỹJi is asymptotically equivalent to the experiment that observes (1.2) by Lemma 2.1.Furthermore, by Lemma3.1, the original NPR experiment that has the same covariance structure as ỹJi is asymptotically equivalent to that experiment and thus, by transitivity, to the experiment that observes the process (1.2) as well. This proves Theorem 1.1.


## 3.2. Remarks on the Covariance Structure

This result is restrictive in that it requires a specific known covariance structure. We are working under the assumption that the covariance matrix has eigenfunctions that correspond to a wavelet basis. This does not generally lead to a typical covariance structure. It does not even necessarily lead to a stationary Gaussian process; see the Haar basis example below.The difficulty is that the requirement for having asymptotically equivalent experiments is quite strict, and the total variation distance between the processes with even small differences in the structure of the covariance is not negligible. For two multivariate Guassian distributions with the same means but where one covariance matrix isΣ and the other is D, a diagonal matrix with the same diagonal elements as Σ, the Kullback–Leibler divergence between the distributions is log|Σ|-log|D|.If the correlation between the highest level coefficients isCorr(ξj*,k,ξj*,k+1)=ρ, then the contribution to the difference of the log determinants is on the order of ρ2j*. The dimension of the problem is growing while the correlations are generally not going to 0 significantly quickly. For instance, in a typical wavelet basis decomposition of the true fractional Brownian motion Corr(ξj*,k,ξj*,k+1)=cβ where cβ is a constant that depends on β but not j* or n.Thus, the differencelog|Σ|-log|D| will not go to 0 as the sample size increases. Therefore, for the sort of long-range correlation structures that we are considering here, the eigenfunctions of the kernel K need to be known or else the experiments will not be asymptotically equivalent.


## 4. Estimating the Covariance of the Increments

The key limitation of Theorem1.1 is that it supposes that the covariance structure of the errors is known to the experimenter. To make the approximation more useful, it would help if the covariance structure was more flexible. A strategy similar to that used by Carter [6] can be used to estimate the variances of the coefficients.By Carter [6], I showed that a model with a variance that changes slowly over time can still be approximated by the Gaussian process as long as all of the observations are independent. Our result here is that for correlated observations, if the variance is a linear function of the frequency, then a similar technique can be used to establish a set of asymptotically sufficient statistics.Flexibility with regard to the covariance structure is added by allowing the magnitude of theVar(yjk) to depend on the resolution level j. The variances will be described by two parameters γ and β, which characterize the size of the error and the speed that it shrinks at higher resolution levels. These nuisance parameters can be estimated using part of the data, and then the inference can be carried out conditionally on the estimates.Specifically, the experimentℰn observes independent components y0~𝒩(θ0,n-(1+β)2γ) and(4.1)yjk~𝒩(θjk,n-(1+β)2γ+(j+1)β)for0≤j≤J,
where the n-(1+β) factor is included to match up with the scaling functions at the Jth resolution level. These observations form a new experiment with a parameter set that includes (μ,γ,β) where μ(t)∈ℳ(M,ε), -1<β<0, and γ is bounded below by a constant -c.This experimentℰn with the parametric variances is no longer an approximately sufficient statistic for the experiment that observes all of the θjk. That experiment has too much information about the variance. If we observed the entire sequence at all resolution levels j, then γ and β could be estimated exactly. We need to adopt another approximating experiment as in Carter [6] Many of the bounds in this section follow arguments from that paper.


### 4.1. Proof of Theorem1.2

The theorem can be proven by applying Lemma3.1 and then a version of Lemma 2.1 that uses only a small proportion of the low-frequency wavelet coefficients. The rest of the coefficients can be used to fix the parameters in the covariance of the observations.The first step is to decompose the nonparametric regression into a set of wavelet coefficients. Then NPR observations Yi can be transformed by dividing by n and then performing the discrete wavelet transformation as in Lemma 3.1. The result is that a sequence of n wavelet coefficients y0 and yjk for j=0,…,J-1 is equivalent to the original NPR observations with a total-variation distance between the distributions of(4.2)δ(ℙμ,γ,β,ℙ̃μγβ)≤C2-γn1-2α+β.
The supremum of this bound over all γ>-c and β<0 is(4.3)δ(𝒫,𝒫̃)≤C2cn1-2α,
which will be negative for α>1/2.The key strategy is to break the observations from this wavelet composition into pieces starting at levelj*, where observations on j≤j* are assumed to be informative about the means, and the higher resolution levels are used to estimate the covariance structure.For each resolution level withj>j*, we generate the approximately sufficient statistics Vj=∑kyjk2. Along with the yjk for j≤j*, the collection of Vj is exactly sufficient if the means are θjk=0 for j>j*, because if there is no information about the means in the higher-frequency terms, then we have a piece of the experiment that is like a normal scale family. This new experiment ℰv is asymptotically equivalent to our ℰn.The error in approximatingℰn by ℰv, where the means of the higher-frequencies coefficients are 0, is bounded by (5.3):(4.4)δ(ℰv,ℰn)≤(∑j=j*+1J-1n1+β2-(γ+(j+1)β)∑kθjk2)1/2.
For θjk in ℳ(M,ε) space, (5.4) bounds the distance as bound(4.5)δ(ℰv,ℰn)≤M2-γ/22(J-j*-1)(1+β)/2-β/2-εj*≤M2c/2+1/22J/2-(1/2+ε)j*,
which is negligible when j*>J/(1+2ε).Thisℰv experiment has sufficient statistics y0, yjk for j≤j*, and(4.6)Vj~Γ(2j,n-(1+β)2γ+β(j+1)-j)forj=j*+1,…,J-1.
Furthermore, there are approximately sufficient statistics in this experiment (yjk,γ̂,β̂) where γ̂ and β̂ are the weighted least-squares estimates of γ and β from the data logVj. These are exactly sufficient statistics in the experiment ℰℓ that observes the y0 and yjk for the lower resolution levels j≤j* as before, in addition to the observations 2Wj for j*<j<J where(4.7)Wj~𝒩(-(1+β)J+γ+β(j+1),2-j+1(log2)).The distance betweenℰℓ and ℰv depends on the distance between the distribution of the log of the Gamma variables and the normal approximation to this distribution. The calculation in of [6, Section 10.1] gives a bound on the Kullback–Leibler divergence of D(ℚj,ℚ̆j)≤2-j where ℚj is the distribution of Vj, and ℚ̆j is the distribution of 2Wj. Therefore, the total error between the two experiments is(4.8)δ(ℰℓ,ℰv)≤(∑j=j*+1J-12-j)1/2≤2-j*/2.
Therefore, the observations in ℰℓ are asymptotically sufficient for ℰv and thus also ℰn (as long as j*→∞ with n).In the experimentℰℓ, the sufficient statistics for estimating γ and β are the weighted least-squares estimators γ̂ and β̂:(4.9)(γ̂β̂)=(x⊤Λx)-1x⊤Λ(W+J),
where Λ is the diagonal matrix with 2j for j=j*+1,…,J-1 along its diagonal, x is the design matrix with rows (1,j-J+1), and W+J is the column of observations Wj-J. The vector of estimators is normal with mean (γβ)⊤ and covariance (1/2)(log2)(x⊤Λx)-1.Therefore, we can compare this experimentℰℓ to an experiment ℰ̂ that observes the same (γ̂β̂)⊤, but the y0 and yjk for j≤j* are replaced by Gaussian random variables ŷjk with variances (conditional on (γ̂,β̂)) that are Var(ŷjk)=2γ̂+(j-J+1)β̂-J. The error in this approximation depends on the distance between the two sets of independent normal experiments with different variances. Letting ℙjk be the distribution of yjk and ℙ̂jk the distribution of ŷjk, the bound (6.12) in Section 6 gives(4.10)D(ℚjk,ℚ̂jk)≤(logn)2n-2j*+O((logn)4(n-2j*+1)2),D(ℚ0,ℚ̂0)≤(logn)2n-2j*+O((logn)4(n-2j*+1)2).
There are 2j*+1 independent normals yjk for j≤j* so that the total divergence is(4.11)δ(ℰℓ,ℰ̂)2≤2(logn)22J-j*-1+O(2j*(logn)4(n-2j*+1)2).
Therefore, the experiments ℰℓ and ℰ̂ are asymptotically equivalent for(4.12)j*=J-2log2J-ηn
for some ηn→∞.We can improve this approximation by replacing the estimatorsβ̂ and γ̂ in ℰ̂ by using(4.13)β̃={-1,β̂≤-1,β̂,-1<β̂<0,0,β̂≥0,
and γ̃=γ̂∨c to match up with the bounds on the parameter space. The new version of this experiment therefore observes (γ̂,β̂) and the normal coordinates yjk~𝒩(θjk,n1+β̃2γ̃2β̃(j+1)) for 0≤j≤j*. The error between ℰℓ and this new version of ℰ̂ is smaller because |γ-γ̃+(β-β̃)(j-J)|≤|γ-γ̂+(β-β̂)(j-J)|, which makes the bound in (6.2) uniformly smaller.Finally, we create a continuous Gaussian version of theℰ̂ experiment. This approximation ℰ̃ observes all the ỹkj for k≥0 with means θjk and variances n-(1+β̃)2γ̃+β̃(j+1). The ℰ̂ are actually sufficient statistics for an experiment that observes (γ̂,β̂) and yjk for 0≤j≤j* and for j>j*:(4.14)yjk~𝒩(0,n-(1+β̃)2γ̃+β̃(j+1)).
The difference between the experiments ℰ̂ and ℰ̃ conditional on γ̃ and β̃ is as in Section 2 and (5.4) less than M2-γ̃2(J-j*)(1+β̃)/2-β̃/2-εj*. The expectation of this bound when averaged over the possible values of (γ̃,β̃) is a bound on the unconditional error. Furthermore, this expectation is less than the minimum over possible values of (γ̃,β̃) (this is the real advantage that comes from going from (γ̂,β̂) to (γ̃,β̃)). Thus,(4.15)δ(ℰ̂,ℰ̃)≤2c-12(J-j*)/2-εj*.
As before, this is asymptotically negligible when j*>J/(1+2ε).All that is left to do is to choose the levelj* at which to split the data so that the requirements from (4.5) and (4.15) that j*>J/(1+2ε) and from (4.12) that j*=J-2log2J-ηn are all fulfilled. We could choose ηn=εJ/(1+2ε)-log2J so that(4.16)j*=J1+2ε+εJ1+2ε-log2J,
which is greater than J/(1+2ε) for ε>0. This choice of j* plugged into the bound in (4.5) gives us(4.17)δ(ℰv,ℰn)≤N2(c+1)/2J1/2+ε2-εJ/2→0
as J→∞ for ε>0. At the same time, the bound in (4.11) becomes(4.18)δ(ℰℓ,ℰ̂)2≤2(log2)2J2Jε/(1+2ε)-J-1→0.
Thus, Theorem 1.2 is established.


## 4.1. Proof of Theorem1.2

The theorem can be proven by applying Lemma3.1 and then a version of Lemma 2.1 that uses only a small proportion of the low-frequency wavelet coefficients. The rest of the coefficients can be used to fix the parameters in the covariance of the observations.The first step is to decompose the nonparametric regression into a set of wavelet coefficients. Then NPR observations Yi can be transformed by dividing by n and then performing the discrete wavelet transformation as in Lemma 3.1. The result is that a sequence of n wavelet coefficients y0 and yjk for j=0,…,J-1 is equivalent to the original NPR observations with a total-variation distance between the distributions of(4.2)δ(ℙμ,γ,β,ℙ̃μγβ)≤C2-γn1-2α+β.
The supremum of this bound over all γ>-c and β<0 is(4.3)δ(𝒫,𝒫̃)≤C2cn1-2α,
which will be negative for α>1/2.The key strategy is to break the observations from this wavelet composition into pieces starting at levelj*, where observations on j≤j* are assumed to be informative about the means, and the higher resolution levels are used to estimate the covariance structure.For each resolution level withj>j*, we generate the approximately sufficient statistics Vj=∑kyjk2. Along with the yjk for j≤j*, the collection of Vj is exactly sufficient if the means are θjk=0 for j>j*, because if there is no information about the means in the higher-frequency terms, then we have a piece of the experiment that is like a normal scale family. This new experiment ℰv is asymptotically equivalent to our ℰn.The error in approximatingℰn by ℰv, where the means of the higher-frequencies coefficients are 0, is bounded by (5.3):(4.4)δ(ℰv,ℰn)≤(∑j=j*+1J-1n1+β2-(γ+(j+1)β)∑kθjk2)1/2.
For θjk in ℳ(M,ε) space, (5.4) bounds the distance as bound(4.5)δ(ℰv,ℰn)≤M2-γ/22(J-j*-1)(1+β)/2-β/2-εj*≤M2c/2+1/22J/2-(1/2+ε)j*,
which is negligible when j*>J/(1+2ε).Thisℰv experiment has sufficient statistics y0, yjk for j≤j*, and(4.6)Vj~Γ(2j,n-(1+β)2γ+β(j+1)-j)forj=j*+1,…,J-1.
Furthermore, there are approximately sufficient statistics in this experiment (yjk,γ̂,β̂) where γ̂ and β̂ are the weighted least-squares estimates of γ and β from the data logVj. These are exactly sufficient statistics in the experiment ℰℓ that observes the y0 and yjk for the lower resolution levels j≤j* as before, in addition to the observations 2Wj for j*<j<J where(4.7)Wj~𝒩(-(1+β)J+γ+β(j+1),2-j+1(log2)).The distance betweenℰℓ and ℰv depends on the distance between the distribution of the log of the Gamma variables and the normal approximation to this distribution. The calculation in of [6, Section 10.1] gives a bound on the Kullback–Leibler divergence of D(ℚj,ℚ̆j)≤2-j where ℚj is the distribution of Vj, and ℚ̆j is the distribution of 2Wj. Therefore, the total error between the two experiments is(4.8)δ(ℰℓ,ℰv)≤(∑j=j*+1J-12-j)1/2≤2-j*/2.
Therefore, the observations in ℰℓ are asymptotically sufficient for ℰv and thus also ℰn (as long as j*→∞ with n).In the experimentℰℓ, the sufficient statistics for estimating γ and β are the weighted least-squares estimators γ̂ and β̂:(4.9)(γ̂β̂)=(x⊤Λx)-1x⊤Λ(W+J),
where Λ is the diagonal matrix with 2j for j=j*+1,…,J-1 along its diagonal, x is the design matrix with rows (1,j-J+1), and W+J is the column of observations Wj-J. The vector of estimators is normal with mean (γβ)⊤ and covariance (1/2)(log2)(x⊤Λx)-1.Therefore, we can compare this experimentℰℓ to an experiment ℰ̂ that observes the same (γ̂β̂)⊤, but the y0 and yjk for j≤j* are replaced by Gaussian random variables ŷjk with variances (conditional on (γ̂,β̂)) that are Var(ŷjk)=2γ̂+(j-J+1)β̂-J. The error in this approximation depends on the distance between the two sets of independent normal experiments with different variances. Letting ℙjk be the distribution of yjk and ℙ̂jk the distribution of ŷjk, the bound (6.12) in Section 6 gives(4.10)D(ℚjk,ℚ̂jk)≤(logn)2n-2j*+O((logn)4(n-2j*+1)2),D(ℚ0,ℚ̂0)≤(logn)2n-2j*+O((logn)4(n-2j*+1)2).
There are 2j*+1 independent normals yjk for j≤j* so that the total divergence is(4.11)δ(ℰℓ,ℰ̂)2≤2(logn)22J-j*-1+O(2j*(logn)4(n-2j*+1)2).
Therefore, the experiments ℰℓ and ℰ̂ are asymptotically equivalent for(4.12)j*=J-2log2J-ηn
for some ηn→∞.We can improve this approximation by replacing the estimatorsβ̂ and γ̂ in ℰ̂ by using(4.13)β̃={-1,β̂≤-1,β̂,-1<β̂<0,0,β̂≥0,
and γ̃=γ̂∨c to match up with the bounds on the parameter space. The new version of this experiment therefore observes (γ̂,β̂) and the normal coordinates yjk~𝒩(θjk,n1+β̃2γ̃2β̃(j+1)) for 0≤j≤j*. The error between ℰℓ and this new version of ℰ̂ is smaller because |γ-γ̃+(β-β̃)(j-J)|≤|γ-γ̂+(β-β̂)(j-J)|, which makes the bound in (6.2) uniformly smaller.Finally, we create a continuous Gaussian version of theℰ̂ experiment. This approximation ℰ̃ observes all the ỹkj for k≥0 with means θjk and variances n-(1+β̃)2γ̃+β̃(j+1). The ℰ̂ are actually sufficient statistics for an experiment that observes (γ̂,β̂) and yjk for 0≤j≤j* and for j>j*:(4.14)yjk~𝒩(0,n-(1+β̃)2γ̃+β̃(j+1)).
The difference between the experiments ℰ̂ and ℰ̃ conditional on γ̃ and β̃ is as in Section 2 and (5.4) less than M2-γ̃2(J-j*)(1+β̃)/2-β̃/2-εj*. The expectation of this bound when averaged over the possible values of (γ̃,β̃) is a bound on the unconditional error. Furthermore, this expectation is less than the minimum over possible values of (γ̃,β̃) (this is the real advantage that comes from going from (γ̂,β̂) to (γ̃,β̃)). Thus,(4.15)δ(ℰ̂,ℰ̃)≤2c-12(J-j*)/2-εj*.
As before, this is asymptotically negligible when j*>J/(1+2ε).All that is left to do is to choose the levelj* at which to split the data so that the requirements from (4.5) and (4.15) that j*>J/(1+2ε) and from (4.12) that j*=J-2log2J-ηn are all fulfilled. We could choose ηn=εJ/(1+2ε)-log2J so that(4.16)j*=J1+2ε+εJ1+2ε-log2J,
which is greater than J/(1+2ε) for ε>0. This choice of j* plugged into the bound in (4.5) gives us(4.17)δ(ℰv,ℰn)≤N2(c+1)/2J1/2+ε2-εJ/2→0
as J→∞ for ε>0. At the same time, the bound in (4.11) becomes(4.18)δ(ℰℓ,ℰ̂)2≤2(log2)2J2Jε/(1+2ε)-J-1→0.
Thus, Theorem 1.2 is established.


## 5. Bounding the Total Variation Distance

We need a bound on the distance between two multivariate normal distributions with different means in order to bound the error in many of our approximations.For shifted Gaussian processes, the total-variation distance between the distributions is(5.1)δ(ℙμ1,ℙμ2)=1-2Φ(-Δ1/22),
where Δ=(μ1-μ2)⊤Σ-1(μ1-μ2). The expression in (5.1) for the total variation distance is concave for positive Δ, so a simple expansion gives(5.2)δ(ℙμ1,ℙμ2)≤12πΔ1/2.For the Gaussian process with correlated components, we will assume that the variance of each wavelet coefficient is of the formVar(ψjk)=σ2n-(1+β)2β(j+1) where the variance is calibrated so that Var(Yi)=σ2=nVar(BK(ϕJ,ℓ)). A bound on the error in the projection onto the span of ψjk for j>j* comes from (5.2) which depends on(5.3)Δ=n1+βσ2∑j>j*2-β(j+1)∑kθjk2(5.4)≤n1+βσ22-(1+β)j*-β∑j>j*2j∑jθjk2≤M22σ-22(J-j*)(1+β)-εj*(2β-2ε+β)-1,
where the upper bound in (5.4) follows from n=2J, the definition of ℳ(M,α) and the bound in (1.13), and -1<β<0. This error is negligible as J→∞ whenever(5.5)j*>J2α.


## 6. Bounds from the Estimated Variances

In order to expand our asymptotically sufficient statistics out into a continuous Gaussian experiment, we need a bound on the total-variation distance betweenℰn which, for 0≤j≤j*, observes a sequence of normals with variances n-(1+β)2γ+β(j+1) and ℰg, which observes a similar set of normals with variances n-(1+β)2γ̂+β̂(j+1).For two normal distributions with the same meansμ, and variances σ12 and σ22, respectively, the Kullback–Leibler divergence is(6.1)D(N1,N2)=12[σ12σ22-1-log(σ12σ22)].
Thus, for ℚkj (the distribution of the ykj) and ℚ̂kj (the distribution of the ŷkj), the divergence between the conditional distributions given γ̂ and β̂ is(6.2)D(ℚkj,ℚ̂kj∣γ̂,β̂)=122γ-γ̂+(j+1-J)(β-β̂)-12-12[γ-γ̂+(j+1-J)(β-β̂)]log2,(6.3)D(ℚ0,ℚ̂0∣γ̂,β̂)=122γ-γ̂-12-(γ-γ̂)log22.
This divergence between conditional distributions can be used to bound the joint divergence:(6.4)D(ℚkj,ℚ̂kj)=𝔼D(ℚkj,ℚ̂kj∣γ̂,β̂),
where the expectation is taken over the estimators γ̂ and β̂.To bound the expected value of the divergence in (6.2), we need the distribution of the estimators:(6.5)(γ̂β̂)~𝒩((γβ),log22(x⊤Λx)-1).
This implies that(6.6)𝔼exp[log2(γ-γ̂+(j+1-J)(β-β̂))]=exp[(log2)22Var(γ̂+(j+1-J)β̂)],
and therefore(6.7)𝔼D(ℚkj,ℚ̂kj∣γ̂,β̂)=exp[(log2)22Var(γ̂+(j+1-J)β̂)]-12.Via elementary linear algebra calculations we get that for(6.8)ζ=(∑i=j*+1J-12i)(∑i=j*+1J-1(J-i-1)22i)+(∑i=j*+1J-1(J-i-1)2i)2,Var(γ̂)=log22(∑i=j*+1J-1(J-i-1)22i)ζ-1,Var(β̂)=log22(∑i=j*+1J-12i)ζ-1,Cov(γ̂,β̂)=log22(∑i=j*+1J-1(J-i-1)2i)ζ-1.
As a result,(6.9)Var(γ̂+(j+1-J)β̂)=log22×[∑i(J-i-1)22i+(J-j-1)2∑i2i-2(J-j-1)∑i(J-i-1)2i]ζ-1.
To simplify this expression, let X be a random variable with probability mass function proportional to 2i for i=j*+1,…,J-1. Thus, the ζ is equal to (∑2i)2[𝔼(J-X-1)2-[𝔼(J-X-1)]2]=(∑2i)2Var(X). Similarly, the main factor in (6.9) is equal to(6.10)(∑2i)[𝔼(J-X-1)2+(J-j-1)2-2(J-j-1)𝔼(J-X-1)]=(∑2i)𝔼(X-j)2.A simple bound of0<X-j<J leads to 𝔼(X-j)2<J2. Furthermore, the variance of X is decreasing in j*, and thus it is greater than 2/9 when j*<J-2. Therefore,(6.11)Var(γ̂+(j+1-J)β̂)≤9log24J2∑2i≤J22J-1-2j*
because 9log2/4<2.Thus,(6.12)D(ℚjk,ℚ̂jk)=12exp[(log2)22Var(γ̂+(j+1-J)β̂)]-12≤12exp[(log2)2J22J-2j*+1]-12≤(log2)2J22J-2j*+1+CJ4(2J-2j*)2=(logn)2n-2j*+1+O((logn)4(n-2j*+1)2).
Analogously, the expected divergence between ℚ0 and ℚ̂jk is bounded by(6.13)𝔼D(ℚ0,ℚ̂0∣γ̂,β̂)=12exp[-(log2)22Var(γ̂)]-12=12exp[-(log2)22((log2)𝔼(J-X-1)22(∑2i)Var(X))]-12≤(log2)2J22J-2j*+1+CJ4(2J-2j*)2.If we add up these errors over the2j* observations in the experiment, we get that the error in the approximation is less than C(logn)2/(n2-j*-1), which is negligible for j* sufficiently small.


## 7. Haar Basis Covariance

The Haar basis is a simple enough wavelet basis by which we can make some explicit calculations of the properties of the error distribution. We will show that the resulting errorsξi will have variances of approximately nβ as we expected, and the correlation between ξi and ξj will decrease at about a rate of |i-j|-(1+β).The scaling functions for the Haar basis are constant on2j dyadic intervals at the resolution level j. The assumption is that we have a single scaling function coefficient with Var(y0)=1, and then every wavelet coefficient yjk is independent and has variance 2β(j+1). Then the covariances can be calculated from the synthesis formula for the Haar basis.The formula for synthesizing the scaling function coefficientsỹJk from the wavelet decomposition is(7.1)ỹJk=2-J/2y0+∑j=0J-1ζj,J,k2(j-J)/2yjk*,
where k* is the index such that ψjk* has support that includes the support of ϕJk. The ζj,J,k is either 1 or -1 depending on whether ϕJk sits in the positive or negative half of the ψjk* function.Using the covariance structure described above, the variance ofỹJk is(7.2)Var(ỹJk)=2-J+∑j=0J-12(j-J)+β(j+1)=2-J+2β-J[2(1+β)J-121+β-1]=2βJ[12-2-β]-2-J[2-β-12-2-β]
for -1<β<0. For β=0, the variance of each scaling function coefficient is 1 as in white noise. For β=-1, direct calculation leads to a variance of 2-J(1+J/2).To find the covariance between two variablesỹJk1 and ỹJk2, we need j* which is the highest resolution level such that the support of ψj*k* includes the support of both scaling functions ϕJk1 and ϕJk2. The covariance is thus(7.3)Cov(ỹJk1,ỹJk2)=2-J+∑j=0j*-12(j-J)+β(j+1)-2j*-J+β(j*+1)=2-J-2j*-J+β(j*+1)+2β-J[2(1+β)j*-121+β-1]=2j*-J+β(j*+1)[2-21+β21+β-1]-2-J[2-β-12-2-β]=2βJ(2-(1+β)(J-j*))[1-2β1-2-β-1]-2-J[2-β-12-2-β].
For large J the correlation is on the order of d-(1+β) where d=2J-j* is a proxy for the distance between the observations. For β=0, all of these covariances are 0. For β=-1, the correlation is j*/(J+2).



---

*Source: 275308-2010-01-26.xml*