# Empirical Mode Decomposition Combined with Local Linear Quantile Regression for Automatic Boundary Correction

Authors: Abobaker M. Jaber; Mohd Tahir Ismail; Alssaidi M. Altaher
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Empirical mode decomposition (EMD) is particularly useful in analyzing nonstationary and nonlinear time series. However, only partial data within boundaries are available because of the bounded support of the underlying time series. Consequently, the application of EMD to finite time series data results in large biases at the edges by increasing the bias and creating artificial wiggles. This study introduces a new two-stage method to automatically decrease the boundary effects present in EMD. At the first stage, local polynomial quantile regression (LLQ) is applied to provide an efficient description of the corrupted and noisy data. The remaining series is assumed to be hidden in the residuals. Hence, EMD is applied to the residuals at the second stage. The final estimate is the summation of the fitting estimates from LLQ and EMD. Simulation was conducted to assess the practical performance of the proposed method. Results show that the proposed method is superior to classical EMD.


---

## Body

## 1. Introduction

We consider the following general nonparametric regression model:(1)
y
=
f
(
x
)
+
ε
i
⋯
,
where Y is the response variable, x is a covariate, f
(
x
)
=
E
(
y
∣
x
) is assumed to be a smooth nonparametric function, and ε
i represents independent and identical random errors with mean 0 and variance σ
2.Empirical mode decomposition (EMD) is a form of analysis based on nonparametric methods [1]. This technique is particularly useful for analyzing nonlinear and nonstationary time series. This method has been widely applied over the last few years to analyze data in different disciplines, such as biology, finance, engineering, and climatology. EMD can enhance estimation performance. Applying the capabilities of EMD as a fully adaptive method and its advantages of handling nonlinear and nonstationary signal behaviors leads to better results. However, EMD suffers from boundary extension, curve fitting, and stopping criteria [2]. Such problems may corrupt the entire data and result in a misleading conclusion [3]. Given that finite data are involved, the algorithms must be adjusted to use certain boundary conditions. In EMD, the end points are also considered problems. The influence of the end points propagates into the data range during sifting. Data extension (or data prediction) is a risky procedure for linear and stationary processes and is more difficult for nonlinear and nonstationary processes. The work in [1] indicated that only the values and locations of the next several extrema, and not all extended data, need to be predicted for EMD. Widely used approaches, such as the characteristic wave extending method, mirror extending method [4], data extending method [5], data reconstruction method [6], and similarity searching method [7], were proposed to overcome the problem and generate a more reasonable solution. The work in [8] introduced quantile regression, a significant extension of traditional parametric and nonparametric regression methods. Quantile regression has been largely used in statistics since its introduction because of its ease of interpretation, robustness, and numerous applications in important areas, such as medicine, economics, environment modeling, toxicology, and engineering [9, 10]. A robust version of classical local linear regression (LLR) known as local linear quantile regression (LLQ) by [11, 12] respectively, have increasingly drawn interest. With its robust behavior, LLQ exhibits excellent boundary adjustment. This characteristic can more efficiently distinguish systematic differences in dispersion, tail behavior, and other features with respect to covariates [12, 13].The current study aims to use the advantages of LLQ to automatically reduce the boundary effects of EMD instead of using classical boundary solutions mentioned previously. The proposed method consists of two stages that automatically decrease the boundary effects of EMD. At the first stage, LLQ is applied to the corrupted and noisy data. The remaining series is then expected to be hidden in the residuals. At the second stage, EMD is applied to the residuals. The final estimate is the summation of the fitting estimates from LLQ and EMD. Compared with EMD, this combination obtains more accurate estimates.The remainder of this study is organized as follows. In Section2, we present a brief background of EMD and LLQ. Section 3 introduces the proposed method. Section 4 compares the results of the original EMD algorithm and the proposed new boundary adjustment by simulation experiments. Conclusions are drawn in Section 5.


## 2. Background

### 2.1. History of Boundary Treatment in Nonparametric Estimators

Most nonparametric techniques such as kernel regression, wavelet thresholding, and empirical mode decomposition show a sharp increase in variance and bias at points near the boundary. Lots of works have been reported in the literature in order to reduce the effects of boundary problem. For kernel regression solution, see [14, 15]. For wavelet thresholding, in addition to use of periodic or symmetric assumption, the authors in [16, 17] used polynomial regression to improve the boundary problem. For empirical mode decomposition the authors in [18] provided a new idea about the boundary extension instead of using the traditional mirror extension on the boundary, and they proposed a ratio extension on boundary. The authors in [19] applied neural network to each IMF to restrain the end effect. The work in [2] provided an algorithm based on the sigma-pi neural network which is used to extend signals before applying EMD. The authors in [20] proposed a new approach that couples the mirror expansion with the extrapolation prediction of regression function to solve boundary problem. The algorithm includes two steps: the extrapolation of the signal through support vector (SV) regression at both endpoints to form the primary expansion signal, and then the primary signal is further expanded through extrema mirror expansion and EMD is performed on the resulting signal to obtain reduced end limitations.In this paper we have followed [16] and [17] strategies to handle end effects of boundary problem in EMD. Instead of using classical polynomial nonparametric regression we will replace it by using a more robust nonparametric estimator, called local linear quantile regression. Practical justifications for choosing such estimator will be explained in Section 2.4.


### 2.2. Empirical Mode Decomposition (EMD)

EMD [1] has proven to be a natural extension and an alternative technique to traditional methods for analyzing nonlinear and nonstationary signals, such as wavelet methods, Fourier methods, and empirical orthogonal functions [21]. In this section, we briefly describe the EMD algorithm. The main objective of EMD is to decompose the data y
t into small signals called intrinsic mode functions (IMF). An IMF is a function in which the upper and the lower envelopes are symmetric; in addition, the number of zero-crossings and the number of extremes are equal or differ by at most one [22]. The algorithm for extracting IMFs for a given time series y
t is called shifting and consists of the following steps.(I)
Setting initial estimates for the residue asr
0
(
t
)
=
y
t, g
0
(
t
)
=
r
k
-
1
(
t
), i
=
1, and the index of IMF as k
=
1.
(II)
Constructing the lower minimaI
min
⁡
i
-
1 and the upper I
max
⁡
i
-
1 envelopes of the signal by the cubic spline method.
(III)
Computing the mean values,m
i, by averaging the upper envelope and the lower envelope as m
i
-
1
=
[
I
max
⁡
i
-
1
+
I
min
⁡
i
-
1
]
/
2.
(IV)
Subtracting the mean from the original signal, that is,g
i
=
g
i
-
1
-
m
i
-
1 and i
=
i
+
1. Steps II to IV are repeated until g
i becomes an IMF. If so, the kth IMF is given by IMF
K
=
g
i.
(V)
Updating the residue asr
k
=
(
t
)
=
r
k
-
1
(
n
)
-
IMF
K. This residual component is treated as new data and subjected to the process described above to calculate the next IMF
K
+
1.
(VI)
Repeating the steps above until the final residual componentr
(
x
) becomes a monotonic function and then considering the final estimation of residue r
^
(
x
).Many methods have been presented to extract trends from a time series. Freehand and least squares methods are the commonly used techniques; the former depends on the experience of users, and the latter is difficult to use when the original series are very irregular [23]. EMD is another effective method for extracting trends [19].


### 2.3. Local Linear Quantile Regression (LLQ)

The seminal study by [8] introduced the parametric quantile regression, which can be considered an alternative to classical regression in both parametric and nonparametric fields. Many models for the nonparametric approach, including locally polynomial quantile regression by [11] and kernel methods by [24], have since been introduced into the statistical literature. In this paper we adopt local linear regression (LLQ) introduced by [12].Let{
(
x
i
,
y
i
), i
=
1
,
…
,
n
} be bivariate observations. To estimate the τth conditional quantile function of response y, the equation below is defined given X
=
x:
(2)
g
(
x
)
=
Q
y
(
τ
∣
x
)
.
Let K be a positive symmetric unimodal kernel function and consider the following weighted quantile regression problem:
(3)
min
⁡
β
∈
ℝ
2
⁡
∑
i
=
1
n
w
i
(
x
)
ρ
τ
(
y
i
-
β
0
-
β
1
(
x
i
-
x
)
)
,
where w
i
(
x
)
=
k
(
(
x
i
-
x
)
/
h
)
/
h. Once the covariate observations are centered at point x, the estimate of g
(
x
) is simply β
^
0, which is the first component of the minimizer of (2). β
^
1 determines an estimate of the slope of the function g at point x.The higher-order LLQ estimate is the minimizer of the following:(4)
min
⁡
β
∈
ℝ
2
⁡
∑
i
=
1
n
w
i
(
x
)
ρ
τ
(
y
i
-
β
0
-
β
1
(
x
i
-
x
)
-
⋯
-
β
p
(
x
i
-
x
)
p
)
.
The choice of the bandwidth parameter hsignificantly influences all nonparametric estimations. An excessively large h obscures too much local structure by excessive smoothing. Conversely, an excessively small h introduces too much variability by relying on very few observations in the local polynomial fitting [13].


### 2.4. Bandwidth Selection

The practical performance ofQ
^
α
(
x
) depends strongly on selected of bandwidth parameter. In this study we adopt the strategy of [12]. In summary, we have the automatic bandwidth selection strategy for smoothing conditional quantiles as follows.(1)
Use ready-made and sophisticated methods to selecth
mean; we use the technique of [25].
(2)
Useh
τ
=
h
mean
{
τ
(
1
-
τ
)
/
ϕ
(
Φ
-
1
(
τ
)
)
2
}
1
/
5  to obtain all other h
τ
s from h
mean.Here, ϕ and Φ are standard normal density and distribution function and h
mean is a bandwidth parameter for regression mean estimation with various existing methods. As it can be seen, this procedure leads to identical bandwidth for τ and (
1
-
τ
) quantiles.


### 2.5. The Behavior of Local Linear Quantile Estimator at Boundary Region

To examine the asymptotic the asymptotic behavior of the local linear quantile estimators at the boundaries, we offer this theorem which has been discussed in detail; see [26]. Here we omitted the proofs and summarized the key points. Without loss of generality, one can consider only the left boundary point u
0
=
c
h, 0
<
c
<
1, if U
t takes value only from [
0,1
]. However, a similar result holds for the right boundary point u
0
=
1
-
c
h.Define(5)
u
j
,
c
=
∫
-
c
1
u
j
K
(
u
)
d
u
,
v
j
,
c
=
∫
-
c
1
u
j
K
2
(
u
)
d
u
.Theorem 1 (see [26]).
Consider the following assumptions.(1)
a
(
u
) is twice continuously differentiable in a neighborhood of u
0 for any u
0.
(2)
f
u
(
u
) is continuous and f
u
(
u
0
)
>
0.
(3)
f
y
∣
u
,
x
(
y
) is bounded and satisfies the Lipschitz condition.
(4)
The kernel functionK
(
·
) is symmetric and has a compact support, say [
-
1,1
].
(5)
{
(
X
t
,
Y
t
,
U
t
)
} is a strictly α-mixing stationary process with mixing coefficient which satisfies ∑
t
≥
1
∞
t
l
α
(
δ
-
2
)
/
δ
(
t
)
<
∞ for some positive real number δ
>
2 and l
>
(
δ
-
2
)
/
δ.
(6)
E
∥
X
t
∥
2
δ
<
∞ with δ
*
>
δ.
(7)
Ω
(
u
0
) is positive-definite and continuous in a neighborhood of u
0.
(8)
Ω
*
(
u
0
) is continuous and positive-definite in a neighborhood of u
0.
(9)
The bandwidthh satisfies h
→
0 and h
→
∞.
(10)
f
(
u
,
v
∣
x
0
,
x
s
;
s
)
≤
M
<
∞ for s
≥
1 where f
(
u
,
v
∣
x
0
,
x
s
;
s
) is the conditional density of (
U
0
,
U
s
) given (
X
0
=
x
0
,
X
s
=
x
s
).
(11)
n
1
/
2
-
δ
/
4
h
δ
/
δ
*
-
1
/
2
-
δ
/
4
=
O
(
1
).The asymptotic normality of the local linear quantile estimator at the left boundary point is given by
(6)
n
h
[
a
^
(
c
h
)
-
a
(
c
h
)
-
h
2
b
c
2
a
′
′
(
0
+
)
+
o
p
(
h
2
)
]
⟶
N
{
o
,
τ
(
1
-
τ
)
v
c
∑
α
0
+
}
,
where
(7)
b
c
=
u
2
,
c
2
-
u
1
,
c
u
3
,
c
u
2
,
c
u
0
,
c
-
u
1
,
c
2
,
v
c
=
u
2
2
,
c
v
0
,
c
-
2
u
1
,
c
u
2
,
c
v
1
,
c
+
u
1
,
c
2
v
2
,
c
[
u
2
,
c
u
0
,
c
-
u
1
,
c
2
]
2
.
Further, the asymptotic normality of the local constant quantile estimator at the left boundary point u
0
=
c
h for 0
<
c
<
1 is
(8)
n
h
[
a
~
(
c
h
)
-
a
(
c
h
)
-
b
~
c
+
o
p
(
h
2
)
]
⟶
N
{
0
,
τ
(
1
-
τ
)
v
0
,
c
∑
α
(
0
+
)
u
0
,
c
2
}
,
where
(9)
b
~
c
=
(
[
h
u
1
,
c
a
′
(
0
+
)
+
h
2
u
2
,
c
2
×
{
a
′
′
(
0
+
)
+
2
a
′
(
0
+
)
f
u
′
(
0
+
)
f
u
(
0
+
)
+
2
Ω
*
-
1
(
0
+
)
Ω
*
′
(
0
+
)
a
′
(
0
+
)
2
a
′
(
0
+
)
f
u
′
(
0
+
)
f
u
(
0
+
)
}
]
)
×
(
u
0
,
c
)
-
1
.From the above theorem, one can deduce that, at the boundaries, the asymptotic bias term for the local constant quantile estimate is of the orderh, compared to the order h
2 for the local linear quantile estimate. Hence, the local linear estimation possesses good behavior at boundaries and there is no need for any boundary correction. In other words, the local linear quantile estimate does not suffer from boundary effects but the local constant quantile estimate does. Therefore, local linear quantile is preferable in practice.


## 2.1. History of Boundary Treatment in Nonparametric Estimators

Most nonparametric techniques such as kernel regression, wavelet thresholding, and empirical mode decomposition show a sharp increase in variance and bias at points near the boundary. Lots of works have been reported in the literature in order to reduce the effects of boundary problem. For kernel regression solution, see [14, 15]. For wavelet thresholding, in addition to use of periodic or symmetric assumption, the authors in [16, 17] used polynomial regression to improve the boundary problem. For empirical mode decomposition the authors in [18] provided a new idea about the boundary extension instead of using the traditional mirror extension on the boundary, and they proposed a ratio extension on boundary. The authors in [19] applied neural network to each IMF to restrain the end effect. The work in [2] provided an algorithm based on the sigma-pi neural network which is used to extend signals before applying EMD. The authors in [20] proposed a new approach that couples the mirror expansion with the extrapolation prediction of regression function to solve boundary problem. The algorithm includes two steps: the extrapolation of the signal through support vector (SV) regression at both endpoints to form the primary expansion signal, and then the primary signal is further expanded through extrema mirror expansion and EMD is performed on the resulting signal to obtain reduced end limitations.In this paper we have followed [16] and [17] strategies to handle end effects of boundary problem in EMD. Instead of using classical polynomial nonparametric regression we will replace it by using a more robust nonparametric estimator, called local linear quantile regression. Practical justifications for choosing such estimator will be explained in Section 2.4.


## 2.2. Empirical Mode Decomposition (EMD)

EMD [1] has proven to be a natural extension and an alternative technique to traditional methods for analyzing nonlinear and nonstationary signals, such as wavelet methods, Fourier methods, and empirical orthogonal functions [21]. In this section, we briefly describe the EMD algorithm. The main objective of EMD is to decompose the data y
t into small signals called intrinsic mode functions (IMF). An IMF is a function in which the upper and the lower envelopes are symmetric; in addition, the number of zero-crossings and the number of extremes are equal or differ by at most one [22]. The algorithm for extracting IMFs for a given time series y
t is called shifting and consists of the following steps.(I)
Setting initial estimates for the residue asr
0
(
t
)
=
y
t, g
0
(
t
)
=
r
k
-
1
(
t
), i
=
1, and the index of IMF as k
=
1.
(II)
Constructing the lower minimaI
min
⁡
i
-
1 and the upper I
max
⁡
i
-
1 envelopes of the signal by the cubic spline method.
(III)
Computing the mean values,m
i, by averaging the upper envelope and the lower envelope as m
i
-
1
=
[
I
max
⁡
i
-
1
+
I
min
⁡
i
-
1
]
/
2.
(IV)
Subtracting the mean from the original signal, that is,g
i
=
g
i
-
1
-
m
i
-
1 and i
=
i
+
1. Steps II to IV are repeated until g
i becomes an IMF. If so, the kth IMF is given by IMF
K
=
g
i.
(V)
Updating the residue asr
k
=
(
t
)
=
r
k
-
1
(
n
)
-
IMF
K. This residual component is treated as new data and subjected to the process described above to calculate the next IMF
K
+
1.
(VI)
Repeating the steps above until the final residual componentr
(
x
) becomes a monotonic function and then considering the final estimation of residue r
^
(
x
).Many methods have been presented to extract trends from a time series. Freehand and least squares methods are the commonly used techniques; the former depends on the experience of users, and the latter is difficult to use when the original series are very irregular [23]. EMD is another effective method for extracting trends [19].


## 2.3. Local Linear Quantile Regression (LLQ)

The seminal study by [8] introduced the parametric quantile regression, which can be considered an alternative to classical regression in both parametric and nonparametric fields. Many models for the nonparametric approach, including locally polynomial quantile regression by [11] and kernel methods by [24], have since been introduced into the statistical literature. In this paper we adopt local linear regression (LLQ) introduced by [12].Let{
(
x
i
,
y
i
), i
=
1
,
…
,
n
} be bivariate observations. To estimate the τth conditional quantile function of response y, the equation below is defined given X
=
x:
(2)
g
(
x
)
=
Q
y
(
τ
∣
x
)
.
Let K be a positive symmetric unimodal kernel function and consider the following weighted quantile regression problem:
(3)
min
⁡
β
∈
ℝ
2
⁡
∑
i
=
1
n
w
i
(
x
)
ρ
τ
(
y
i
-
β
0
-
β
1
(
x
i
-
x
)
)
,
where w
i
(
x
)
=
k
(
(
x
i
-
x
)
/
h
)
/
h. Once the covariate observations are centered at point x, the estimate of g
(
x
) is simply β
^
0, which is the first component of the minimizer of (2). β
^
1 determines an estimate of the slope of the function g at point x.The higher-order LLQ estimate is the minimizer of the following:(4)
min
⁡
β
∈
ℝ
2
⁡
∑
i
=
1
n
w
i
(
x
)
ρ
τ
(
y
i
-
β
0
-
β
1
(
x
i
-
x
)
-
⋯
-
β
p
(
x
i
-
x
)
p
)
.
The choice of the bandwidth parameter hsignificantly influences all nonparametric estimations. An excessively large h obscures too much local structure by excessive smoothing. Conversely, an excessively small h introduces too much variability by relying on very few observations in the local polynomial fitting [13].


## 2.4. Bandwidth Selection

The practical performance ofQ
^
α
(
x
) depends strongly on selected of bandwidth parameter. In this study we adopt the strategy of [12]. In summary, we have the automatic bandwidth selection strategy for smoothing conditional quantiles as follows.(1)
Use ready-made and sophisticated methods to selecth
mean; we use the technique of [25].
(2)
Useh
τ
=
h
mean
{
τ
(
1
-
τ
)
/
ϕ
(
Φ
-
1
(
τ
)
)
2
}
1
/
5  to obtain all other h
τ
s from h
mean.Here, ϕ and Φ are standard normal density and distribution function and h
mean is a bandwidth parameter for regression mean estimation with various existing methods. As it can be seen, this procedure leads to identical bandwidth for τ and (
1
-
τ
) quantiles.


## 2.5. The Behavior of Local Linear Quantile Estimator at Boundary Region

To examine the asymptotic the asymptotic behavior of the local linear quantile estimators at the boundaries, we offer this theorem which has been discussed in detail; see [26]. Here we omitted the proofs and summarized the key points. Without loss of generality, one can consider only the left boundary point u
0
=
c
h, 0
<
c
<
1, if U
t takes value only from [
0,1
]. However, a similar result holds for the right boundary point u
0
=
1
-
c
h.Define(5)
u
j
,
c
=
∫
-
c
1
u
j
K
(
u
)
d
u
,
v
j
,
c
=
∫
-
c
1
u
j
K
2
(
u
)
d
u
.Theorem 1 (see [26]).
Consider the following assumptions.(1)
a
(
u
) is twice continuously differentiable in a neighborhood of u
0 for any u
0.
(2)
f
u
(
u
) is continuous and f
u
(
u
0
)
>
0.
(3)
f
y
∣
u
,
x
(
y
) is bounded and satisfies the Lipschitz condition.
(4)
The kernel functionK
(
·
) is symmetric and has a compact support, say [
-
1,1
].
(5)
{
(
X
t
,
Y
t
,
U
t
)
} is a strictly α-mixing stationary process with mixing coefficient which satisfies ∑
t
≥
1
∞
t
l
α
(
δ
-
2
)
/
δ
(
t
)
<
∞ for some positive real number δ
>
2 and l
>
(
δ
-
2
)
/
δ.
(6)
E
∥
X
t
∥
2
δ
<
∞ with δ
*
>
δ.
(7)
Ω
(
u
0
) is positive-definite and continuous in a neighborhood of u
0.
(8)
Ω
*
(
u
0
) is continuous and positive-definite in a neighborhood of u
0.
(9)
The bandwidthh satisfies h
→
0 and h
→
∞.
(10)
f
(
u
,
v
∣
x
0
,
x
s
;
s
)
≤
M
<
∞ for s
≥
1 where f
(
u
,
v
∣
x
0
,
x
s
;
s
) is the conditional density of (
U
0
,
U
s
) given (
X
0
=
x
0
,
X
s
=
x
s
).
(11)
n
1
/
2
-
δ
/
4
h
δ
/
δ
*
-
1
/
2
-
δ
/
4
=
O
(
1
).The asymptotic normality of the local linear quantile estimator at the left boundary point is given by
(6)
n
h
[
a
^
(
c
h
)
-
a
(
c
h
)
-
h
2
b
c
2
a
′
′
(
0
+
)
+
o
p
(
h
2
)
]
⟶
N
{
o
,
τ
(
1
-
τ
)
v
c
∑
α
0
+
}
,
where
(7)
b
c
=
u
2
,
c
2
-
u
1
,
c
u
3
,
c
u
2
,
c
u
0
,
c
-
u
1
,
c
2
,
v
c
=
u
2
2
,
c
v
0
,
c
-
2
u
1
,
c
u
2
,
c
v
1
,
c
+
u
1
,
c
2
v
2
,
c
[
u
2
,
c
u
0
,
c
-
u
1
,
c
2
]
2
.
Further, the asymptotic normality of the local constant quantile estimator at the left boundary point u
0
=
c
h for 0
<
c
<
1 is
(8)
n
h
[
a
~
(
c
h
)
-
a
(
c
h
)
-
b
~
c
+
o
p
(
h
2
)
]
⟶
N
{
0
,
τ
(
1
-
τ
)
v
0
,
c
∑
α
(
0
+
)
u
0
,
c
2
}
,
where
(9)
b
~
c
=
(
[
h
u
1
,
c
a
′
(
0
+
)
+
h
2
u
2
,
c
2
×
{
a
′
′
(
0
+
)
+
2
a
′
(
0
+
)
f
u
′
(
0
+
)
f
u
(
0
+
)
+
2
Ω
*
-
1
(
0
+
)
Ω
*
′
(
0
+
)
a
′
(
0
+
)
2
a
′
(
0
+
)
f
u
′
(
0
+
)
f
u
(
0
+
)
}
]
)
×
(
u
0
,
c
)
-
1
.From the above theorem, one can deduce that, at the boundaries, the asymptotic bias term for the local constant quantile estimate is of the orderh, compared to the order h
2 for the local linear quantile estimate. Hence, the local linear estimation possesses good behavior at boundaries and there is no need for any boundary correction. In other words, the local linear quantile estimate does not suffer from boundary effects but the local constant quantile estimate does. Therefore, local linear quantile is preferable in practice.


## 3. Proposed Method

This section elaborates on the proposed method. This technique combines EMD and LLQ (EMD-LLQ). Since local linear quantile regression produces excellent boundary treatment [27], it is expected that the addition of this component to empirical mode decomposition will result in equally well-boundary properties. Results from numerical experiments extremely support this claim.The basic idea behind the proposed method is to estimate the underlining functionf with the sum of a set of EMD functions, f
EMD, and an LLQ function, f
LLQ. That is,
(10)
f
^
LLQ
·
EMD
=
f
^
EMD
+
f
^
LLQ
.
We need to estimate the two components f
EMD and f
LLQ to obtain our proposed estimate, f
^
EMD
·
LLQ, by the following steps.(1)
Applying LLQ to the corrupted and noisy data,y
i and obtaining the trend estimatef
^
LLQ.
(2)
Determining the residualse
i from LLQ; that is, e
i
=
y
i
-
f
^
LLQ.
(3)
Applying EMD toe
i, given that the remaining series is expected to be hidden in the residuals. This step is accomplished by performing the following substeps.
(I)
Setting initial estimates for the residue asr
0
(
t
)
=
e, g
o
(
t
)
=
r
k
-
1
(
t
), i
=
1, and the index of IMF as k
=
1.
(II)
Constructing the lower minimaI
min
⁡
i
-
1  and I
max
⁡
i
-
1  envelopes of the signal by the cubic spline method.
(III)
Calculating the mean values by averaging the upper envelope and the lower envelope. Settingm
i
-
1
=
[
I
max
⁡
i
-
1
+
I
min
⁡
i
-
1
]
/
2.
(IV)
Subtracting the mean from the original signal asg
i
=
g
i
-
1
-
m
i
-
1 and i
=
i
+
1. Steps I to IV are repeated until g
i becomes an IMF. The kth IMF is then given as IMF
K
=
g
i.
(V)
Updating the residuer
k
(
x
)
=
r
k
-
1
(
n
)
-
IM
F
K. This residual component is regarded as a new datum and is subjected to the process described above to calculate the next IMF
K
+
1.
(VI)
The steps above are repeated until the final residual componentr
(
x
) becomes a monotonic function. The final estimation of the residue r
^
(
x
) is then considered.
(4)
The final estimate is the summation of the fitting estimates from LLQ and EMD, as follows:(11)
f
^
LLQ
·
EMD
=
f
^
EMD
+
r
^
(
t
)
.


## 4. Simulation Study

In this simulation, the software packageR was employed to evaluate classical EMD by [1] and the proposed combined method, EMD-LLQ. The following conditions were set.(1)
Three different test functions (Table1).
(2)
Three different values of quantileτ (0.25, 0.50, and 0.75).
(3)
Three different kinds of noise structure errors, namely:
(a)
normal distribution with zero mean and unity variance,
(b)
correlated noise from the first-order autoregressive model AR (1) with parameter (0.5),
(c)
heavy-tailed noise fromt distribution with three degrees of freedom.Table 1
Formula of the test functions used in the simulation.
Test function
Formula
1
f
(
x
)
=
sin
(
π
x
)
-
sin
(
2
π
x
)
+
0.5
x
2
f
(
x
)
=
{
10
e
-
10
x
+
2
i
f
x
≤
0.5
3cos
(
10
π
x
)
0.5
<
x
<
1
3
f
(
x
)
=
5
e
(
-
10
(
x
-
0.5
)
2
)
+
e
(
-
10
x
)Datasets were simulated from each of the three test functions with a sample size ofn
=
100 (Figure 1). For each simulated dataset, the above two methods were applied to estimate the test function. In each case, 1,000 replications of the sample size n
=
100 were made. The mean squared error (MSE) was used as the numerical measure to assess the quality of the estimate. The MSE was calculated for those observations that were at most 10 sample points away from the boundaries of the test functions:(12)
MS
E
Δ
(
f
^
)
=
1
2
Δ
∑
i
∈
N
(
Δ
)
{
f
(
x
i
)
-
f
^
(
x
i
)
}
2
a
a
a
a
a
a
a
a
(
Δ
=
1,2
,
…
,
[
n
2
]
;
x
i
=
i
n
)
,
where N
(
Δ
)
=
{
1
,
…
,
Δ
,
n
-
Δ
+
1
,
…
,
n
}.Three test functions used in the simulation.
(a)
(b)
(c)To compare the methods, Tables2, 3, and 4 present the numerical results of the classical EMD with respect to the proposed method.Table 2
The MSE of the classical EMD and proposed method under variety of boundary solution noise structure, different values of quantileτ (0.25, 0.50, and 0.75), and sample size 100.
Test function 1
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (none)
0.269175
0.310785
0.32917
0.261690
0.3027736
0.3190445
0.2552796
0.306369
0.3223
EMD-LLQ
0.11712
0.122119
0.145831
0.06343
0.058571
0.010130
0.070048
0.067575
0.10510
Wilcoxon test
V
=
451927
459306
449112
49589
497928
480711
493926
494605
479475
Test function 1
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (periodic)
1.38123
1.40508
1.45777
1.42545
1.43629
1.47529
1.40223
1.41718
1.44426
EMD-LLQ
0.11287
0.12451
0.14150
0.06375
0.05967
0.09826
0.06509
0.06284
0.10475
Wilcoxon test
V
=
500496
50076
500427
500500
500076
500477
500493
500490
500476
P
-value
<
2.2
e
-
16
Test function 1
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (symmetric)
0.86598
0.91897
1.005607
0.87603
0.916301
1.005788
0.8820793
0.92558
0.96602
EMD-LLQ
0.11656
0.12304
0.14698
0.8760
0.05816
0.09760
0.070185
0.06459
0.10498
Wilcoxon test
V
=
500500
499701
500456
500498
499540
500497
500500
500497
500492
P
-value
<
2.2
e
-
16
Test function 1
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
Errors
N~(0,1
)
EMD (wave)
1.20811
1.18084
1.250195
1.19274
1.20855
1.226365
1.199136
1.20274
1.22829
EMD-LLQ
0.11318
0.11967
0.14845
0.06456
0.06173
0.098486
0.066864
0.06339
0.10583
Wilcoxon test
V
=
500500
500500
500500
500500
500500
500500
500500
500500
500500Table 3
The MSE of the classical EMD and proposed method under variety of boundary solution noise structure, different values of quantileτ (0.25, 0.50, and 0.75), and sample size = 100.
Test function 2
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (none)
13.71622
8.517261
7.290762
7.818736
8.094024
6.974746
7.168348
8.445187
7.364002
EMD-LLQ
2.200657
2.148026
2.063819
2.147824
0.8982332
0.847309
1.68474
1.69877
1.57035
Wilcoxon test
V
=
467459
463117
500260
499774
498875
489150
482567
483765
489150
P
-value
<
2.2
e
-
16
Test function 2
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (periodic)
7.430924
6.976339
6.77826
7.135594
6.982454
6.827251
7.226472
6.865929
6.725184
EMD-LLQ
2.118629
2.145495
2.052055
0.904618
0.9031342
0.8702984
1.672161
1.691775
1.598264
Wilcoxon test
V
=
498562
497328
495251
500441
500350
500203
498068
496900
496766
P
-value
<
2.2
e
-
16
Test function 2
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (symmetric)
8.693953
8.800117
8.678484
8.679718
8.902223
8.76185
8.623087
8.756977
8.718655
EMD-LLQ
2.121942
2.142339
2.059669
0.918506
0.8823274
0.8637991
1.655741
1.670038
1.571081
Wilcoxon test
V
=
500500
500500
500500
500500
500500
500500
500500
500500
500500
P
-value
<
2.2
e
-
16
Test function 2
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (wave)
7.243865
7.493569
7.341059
7.138806
7.430882
7.449
7.124459
7.374468
7.302789
EMD-LLQ
2.12519
2.149727
2.053073
0.908747
0.9056965
0.8560205
1.653102
1.704716
1.552656
Wilcoxon test
V
=
500500
500500
500500
500500
500500
500500
500500
500497
500500
P
-value
<
2.2
e
-
16Table 4
The MSE of the classical EMD and proposed method under variety of boundary solution noise structure, different values of quantileτ (0.25, 0.50, and 0.75), and sample size = 100.
Test function 3
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (none)
0.034744
0.037490
0.06867
0.035888
0.035473
0.0708678
0.03540
0.0375
0.07260
EMD-LLQ
0.01494
0.01827
0.04637
0.0142052
0.0168427
0.04288602
0.014913
0.0197
0.048020
Wilcoxon test
V
=
453610
426890
371062
740664
435835
415746
455217
420619
393951
P
-value
<
2.2
e
-
16
Test function 3
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (periodic)
0.01495
0.01601
1.46184
0.014931
0.017672
1.475389
0.01513
0.0173
1.44361
EMD-LLQ
0.01201
0.01425
0.14555
0.010031
0.016859
0.0967748
0.01176
0.0157
0.10338
Wilcoxon test
V
=
142798
167295
500426
159630
193475
500491
133946
170435
500474
P
-value
<
2.2
e
-
16
Test function 3
Quantiles
τ
=
0.25
τ
=
0.50
τ=0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (symmetric)
0.856264
0.938382
1.01019
0.853808
0.928291
0.9873542
0.86854
0.9322
0.96783
EMD-LLQ
0.11334
0.12452
0.14649
0.065446
0.056100
0.09569597
0.06553
0.06180
0.104910
Wilcoxon-test
V
=
500489
499304
500420
500500
500498
500474
500500
500485
500463
P
-value
<
2.2
e
-
16
Test function 3
Quantiles
τ
=
0.25
τ
=
0.50
τ
=
0.75
Errors
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
N~(0,1
)
T
(
100,3
)
AR(0.5)
EMD (wave)
1.19549
1.19853
1.22884
1.180275
1.200136
1.217799
1.193706
1.20841
1.236941
EMD-LLQ
0.11508
0.12732
0.14624
0.063552
0.055484
0.094623
0.06860
0.20362
0.110907
Wilcoxon test
V
=
500500
500500
500500
500500
500500
500500
500500
499500
500500
P
-value
<
2.2
e
-
16


### 4.1. Results

From the simulation results, reported in Tables2, 3, and 4, we have observed the following. Regardless of the boundary assumptions, test functions, noise structures, and different values of quantile, the proposed method is constantly superior to the classical EMD under periodic, symmetric (Mirror) and wave conditions. Tables 2, 3, and 4 summarize the results.To ensure that the improvement in mean squared error is due to our proposed treatment, not to something else, we evaluated the classical method and our proposed one when no boundary treatment has been set up at all. From simulation result, we observed that even though the classical solutions help improve the mean squared error, our improvement is much better. Then, at the end, to get rid of some suspicions that the differences might not be significant, we used rank Wilcoxon test. This provided us evidence that our proposed method still achieves a better performance near the boundaries than EMD. AllP value for Wilcoxon test are less than 0.05.


## 4.1. Results

From the simulation results, reported in Tables2, 3, and 4, we have observed the following. Regardless of the boundary assumptions, test functions, noise structures, and different values of quantile, the proposed method is constantly superior to the classical EMD under periodic, symmetric (Mirror) and wave conditions. Tables 2, 3, and 4 summarize the results.To ensure that the improvement in mean squared error is due to our proposed treatment, not to something else, we evaluated the classical method and our proposed one when no boundary treatment has been set up at all. From simulation result, we observed that even though the classical solutions help improve the mean squared error, our improvement is much better. Then, at the end, to get rid of some suspicions that the differences might not be significant, we used rank Wilcoxon test. This provided us evidence that our proposed method still achieves a better performance near the boundaries than EMD. AllP value for Wilcoxon test are less than 0.05.


## 5. Conclusions

In this study, a new two-stage method is introduced to decrease the effects of the boundary problem in EMD. The proposed method is based on a coupling of LLQ at the first stage and classical EMD at the second stage. The empirical performance of the proposed method was tested on different numerical experiments by simulation and real data application. The results of these experiments illustrate the improvement of the EMD estimation in terms of MSE.



---

*Source: 731827-2014-03-25.xml*