\documentclass{CUP-JNL-AAS}

\usepackage{placeins}
\pdfoutput=1 
\usepackage[utf8]{inputenc} % allow utf-8 input
\usepackage[T1]{fontenc}    % use 8-bit T1 fonts
\usepackage{hyperref}       % hyperlinks
\usepackage{url}            % simple URL typesetting
\usepackage{booktabs}       % professional-quality tables
\usepackage{amsfonts}       % blackboard math symbols
\usepackage{nicefrac}       % compact symbols for 1/2, etc.
\usepackage{graphicx}
\usepackage{bm}
\usepackage{float}
\usepackage{amsmath}
\usepackage{subcaption}
\usepackage{enumitem}
\usepackage[english]{babel}
\usepackage{placeins}
\usepackage{algorithm2e}
\usepackage{caption}
\captionsetup[table]{skip=10pt}
\usepackage{amsthm}
\usepackage{xcolor}
\usepackage{titlesec} % custom formatting Appendix
\usepackage{environ} % custom formatting Appendix
\usepackage{amsmath} % added for cases environment
\usepackage{bbm} % added for indicator function
\usepackage{amssymb} % added for real number symbol
\usepackage{relsize} % added to increase size of sum symbol
\usepackage{tikz} % drawing diagrams
\usepackage{algorithm} % pseudo-code
\usepackage{algorithmic} % pseudo-code
\usepackage{pifont} % table checks and crosses
\usepackage{array} % text wrapping in table
\usepackage{multirow} % text wrapping in table
\usepackage{tabu} % text wrapping in table
\usepackage{placeins} % for \Floatbarrier and white space management

\newcommand{\cmark}{\color{teal}\ding{51}} % green check
\newcommand{\xmark}{\color{red}\ding{55}} % red check

\newcommand{\pkg}[1]{{\fontseries{m}\fontseries{b}\selectfont #1}}

\renewcommand\abstractname{}

    
\newenvironment{advanced}{\titleformat{name=\section}{\bfseries}{}{0pt}{\large}}{}  % custom formatting Appendix

\usepackage{fancyvrb} % handles verbatim code
\DefineVerbatimEnvironment{Code}{Verbatim}{}

% for custom list items "Condition 1 & Condition 2"
\makeatletter
\def\namedlabel#1#2{\begingroup
	#2%
	\def\@currentlabel{#2}%
	\phantomsection\label{#1}\endgroup
}
\makeatother

% for pseudo code
\usepackage{algorithm}
\usepackage{algorithmic}

% for larger math symbols
\usepackage{relsize}

% customise section number
\renewcommand\thesection{\ifcase\value{section}?\or 10\or S\or D\or E\else ?\fi}

% for custom author output when citing
% concise author in text
% full author in bibliography
\DeclareRobustCommand{\firstsecond}[2]{#1}

\DeclareMathSizes{10}{8}{5}{2}
\DeclareMathOperator*{\tr}{\text{Tr}}
\DeclareMathOperator*{\myvec}{\text{vec}}
\DeclareMathOperator*{\diag}{diag}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator*{\argmax}{argmax}
\allowdisplaybreaks

\theoremstyle{remark}
\newtheorem{remark}{Remark}
\theoremstyle{lemma}
\newtheorem{lemma}{Lemma}[section]

\providecommand{\keywords}[1]{\textbf{\textit{Index terms---}} #1}

\graphicspath{{Images/}}

\begin{document}

\lefttitle{Supplement to: Package \textbf{AdvEMDpy}}
\righttitle{Annals of Actuarial Science}

\begin{authgrp}
\author{Cole van Jaarsveldt$^{1}$, Matthew Ames$^{2}$, Gareth W. Peters$^{3}$, and Mike Chantler$^{4}$}
\affiliation{
$^{1}$School of Mathematical and Computer Sciences, Heriot-Watt University, Edinburgh, United Kingdom of Great Britain and Northern Ireland (UK), \texttt{cv25@hw.ac.uk}
$^{2}$ResilientML, Melbourne, Australia, \texttt{matt.ames@resilientml.com}
$^{3}$Department of Statistics \& Applied Probability, University of California, Santa Barbara, California, United States of America (USA), \texttt{garethpeters@ucsb.edu}
$^{4}$School of Mathematical and Computer Sciences, Heriot-Watt University, Edinburgh, UK, \texttt{m.j.chantler@hw.ac.uk}}
\end{authgrp}

\title{Supplement to: Package \textbf{AdvEMDpy}: Algorithmic Variations of Empirical Mode Decomposition in \textsf{Python}}

\begin{abstract}
    This work serves as a formal supplement to `Package \textbf{AdvEMDpy}: Algorithmic Variations of Empirical Mode Decomposition in \textsf{Python}' with additional synthetic and real-world examples. All of Section 7, from the `Package \textbf{AdvEMDpy}: Algorithmic Variations of Empirical Mode Decomposition in \textsf{Python}', will be repeated here verbatim with other viable algorithmic variations and further details on the extensions to the algorithm already alluded to in Section 7. This is done to make the original work more concise, whilst simultaneously providing those interested in further reading, examples, and algorithmic variations with further compelling literature. \textbf{AdvEMDpy} will be shown to be more accurate than its \textsf{Python} competitors in resolving the underlying driving function of the Duffing Equation, before it is used to isolate different frequency structures present in Carbon ETF data. An annual fluctuation will be extracted and possibly causally linked to the seasonal trend of the Carbon Dioxide concentration in the atmosphere. These examples are by no means exhaustive and merely serve as demonstrations of \textbf{AdvEMDpy}'s usage and superiority.
\end{abstract}

\begin{keywords} {Empirical Mode Decomposition (EMD), Statistical EMD (SEMD), Enhanced EMD (EEMD), Ensemble EMD, Hilbert transform, time series analysis, filtering, graduation, Winsorization, downsampling, splines, knot optimisation, \textsf{Python}, \textsf{R}, \textsf{MATLAB}}
\end{keywords}

\maketitle

\vspace{-2em}
\section*{Software Availability}The software accompanying this paper is available on GitHub at:
\vspace{-1.0em}
\begin{center}
    \url{https://github.com/Cole-vJ/AdvEMDpy}.
\end{center}
\vspace{-1.0em}
Instructions on how to install, extensive worked examples, as well as the package versions required for complete reproducibility can all be found in the repository.

\vspace{-1.5em}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ Core Details of AdvEMDpy Package} \label{sec:Package}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In this section, the details of each of the components of \textbf{AdvEMDpy} package are outlined both in how to interact with the package through specific functionalities and what expectations one should have on outputs created and features that can be customised.

\vspace{-1.0em}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ Base Implementation of AdvEMDpy Package} \label{sec:Base_implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Before discussing the nuances of the algorithmic variations, one would benefit from some choice remarks concerning non-essential, but helpful outputs of the \textbf{AdvEMDpy} package. These play no direct part in the implementation of the algorithm but assist in iterative understanding and assessment of the EMD decompositions performance. The base implementation of the algorithm may be seen below. More detailed descriptions may be found in the associated scripts as well as a base implementation in the associated \path{README.md} file.

\newpage

\begin{Code}
# instantiate EMD class with time (optional 
# keyword argument) and time series
emd = EMD(time=time, time_series=time_series)

# base implementation with optional helpful inputs shown
# knots and knot_time are optional keyword arguments
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots, 
                                                  knot_time=knot_time, 
                                                  debug=False,
                                                  dft=`envelopes',
                                                  verbose=True,
                                                  stopping_criterion=`sd',
                                                  output_coefficients=False,
                                                  output_knots=False)[:3]
\end{Code}

The details of each field in the input to the function are detailed in the following subsections. If time is not specified as a keyword argument, then the following code is implemented where time is simply an index set of the same length as the original time series.

\vspace{-0.5em}
\begin{Code}
self.time = np.arange(len(self.time_series))
\end{Code}
\vspace{-0.5em}

If the knot sequence or the associated knot time are also unspecified keyword arguments, then the following code is implemented that sets knot time equal to time and spaces the knots $ 10 $ times further apart than the time points. This is problematic if the default knot sequence is insufficient to capture the highest frequency structures which is demonstrated in \texttt{new\_user.ipynb}.

\vspace{-0.5em}
\begin{Code}
knots = np.linspace(0, self.time[-1], int(len(self.time) / 10 + 1))
knot_time = self.time  
\end{Code}
\vspace{-0.5em}

The only required argument is the time series. This is detailed further in \texttt{new\_user.ipynb} which is intended for new users to get comfortable with the package and the available options.


\vspace{-1.0em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Debug flag in AdvEMDpy Package}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
If true, each iteration of the local mean estimation through the chosen detrended mean threshold technique will be plotted for analysis. The debugging output displayed depends on the choice of \texttt{dft}. If \texttt{dft=`envelopes'} (base implementation), the extrema, the extrema envelopes, and the calculated local mean will be plotted. If \texttt{dft=`inflection\_points'} the extrema, inflection points and the associated spline fitted through the inflection points will be displayed. If \texttt{dft=`binomial\_average'} the extrema, the points used for the binomial average, the binomial average of these points, and the associated spline fitted will be displayed. Finally, if \texttt{dft=`enhanced'}, the extrema, the optimal 'extrema', the associated splines, and the resulting local mean are displayed.   

\vspace{-1.0em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Verbose flag in AdvEMDpy Package}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
If true, each iteration of the algorithm will output text describing if the stopping criterion chosen or the mean threshold is met. The intermediate numbering of the potential IMFs is output for consistency and the iteration counter is printed if the maximum number of iterations is met. Different text is displayed for each of the optional stopping criteria. If, for example, \texttt{stopping\_criterion=`sd'} then the text displayed for the fifth candidate for the second IMF, should the IMF candidate not meet the stopping criterion, will be \texttt{IMF\_25 Standard deviation STOPPING CRITERION NOT MET with sd = 1000.0} or it will be \texttt{IMF\_25 Standard deviation STOPPING CRITERION MET with sd = 0.05 < sd threshold = 0.1} with the stopping criterion threshold being another possible input such as \texttt{stopping\_criterion\_threshold=0.1}, if the stopping criterion condition is met. For \texttt{stopping\_criterion\_threshold}, a value of $ 0.2 - 0.3 $ is recommended in the paper that originally proposed the method, \cite{Hua98}. This value should be time series dependent as the value is very much dependent on the level of noise in the time series as well as the number of time points available, but more on this (and other stopping criteria) in Section \ref{sec:Stopping_criteria}.

\vspace{-1em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Output Coefficients flag in AdvEMDpy Package}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
If \texttt{output\_coefficients=True}, cubic B-spline coefficients corresponding to each IMF output (including initial smoothed time series and trend) are output. One must remember the six additional coefficients corresponding to the three additional knots on either edge of the time series as demonstrated in Figure \ref{fig:boundary_bases}. These additional bases result in implicit smoothing through a non-natural spline which does not tend to zero at the edges in both position and the higher-order derivatives. It would be unnecessarily restrictive to impose a natural spline on the sifting algorithm and the resulting IMFs.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{boundary_bases_supp.png}
	\caption{Figure demonstrating incomplete bases at the boundary of the time series to create non-natural cubic B-splines that can accommodate non-zero edge values, derivatives, and curvatures.}
	\label{fig:boundary_bases}
\end{figure}

\vspace{-2em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Output Knots flag in AdvEMDpy Package}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
If \texttt{output\_knots=True}, knot points are output. If knot points are not optimised, that is if \texttt{optimise\_knots=0}, then original knot points are repeated at each iteration of the sifting procedure when extracting recursively each IMF. In other words, all IMF's will have the same knot sequence. If the user provides a uniform knot sequence and chooses not to optimise the knot locations, then a uniform knot sequence will be used throughout as demonstrated in Figure \ref{fig:knot_uniform}.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.75\textwidth]{Images/knot_uniform.png}
	\caption{Figure demonstrating predefined uniform knot placement and the resulting IMFs.}
	\label{fig:knot_uniform}
\end{figure}

Alternatively, one can optimize the knots on an original spline (\texttt{optimise\_knots=1}) representation of the input time series signal using optimal knot placement methodology as outlined in Section \ref{sec:Knot_optimisation} where two options are available using simple bisection and serial bisection, but more about this will follow in the dedicated section. In this way the knot points will still be universal for all IMF's however, they will be optimized for the given signal input and will not dynamically adjust as the residuals become increasingly less complex. For this reason, this can be referred to as the statically optimised knots and an example demonstrating this can be seen in Figure \ref{fig:knot_1} where the optimised non-uniform knots can be seen to be used throughout the algorithm despite the decreasing complexity.

\newpage

\begin{figure}[h]
\vspace{-1.5em}
	\centering
	\includegraphics[width=0.75\textwidth]{Images/knot_1.png}
	\caption{Figure demonstrating statically optimised knot placement, which is optimised once at outset and used throughout the sifting, and the resulting IMFs.}
	\label{fig:knot_1}
\end{figure}

\vspace{-1.0em}

In addition, one could also optimize the knot points per iteration (\texttt{optimise\_knots=2}) of the EMD sifting procedure once extracting an IMF, the next iteration of EMD on the residual signal could first have the knot point optimisation performed when fitting the spline to the residual to proceed with the next rounds of sifting to extract the next IMF. This results in different numbers of knot points per IMF and different placements, in general one would expect to require far fewer knots for later IMFs extracted which have lower frequency content compared to those first few IMF bases extracted. An example demonstrating this can be seen in Figure \ref{fig:knot_2} where there are far fewer non-uniform knots for the second IMF than was required for the first IMF.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.75\textwidth]{Images/knot_2.png}
	\caption{Figure demonstrating dynamically optimised knot placement, which is optimised at the beginning of each internal sifting routine, and the resulting IMFs.}
	\label{fig:knot_2}
\end{figure}

Performing the third option may increase the speed of the algorithm if a large number of siftings need to take place and a large number of IMFs need to be extracted, despite the added time required to find the optimised knots at each stage. Progressively fewer knots will be required throughout the algorithm with the added advantage of having a more parsimonious representation with fewer parameters required for higher-order IMFs which have guaranteed reducing oscillation in their representation.

\vspace{-2.0em}

\subsubsection{ Recommendations for Base Implementation of AdvEMDpy Package}

For the initial application of the algorithm to the desired time series, it is recommended that the user display the results of the selected stopping criterion (\texttt{verbose=True}) and run to completion. The user can then inspect the outputs and should anything be irregular or if the user wants to observe each iterative step of the algorithm one can use \texttt{debug=True} which will cause the method to plot each iteration of the algorithm.

\vspace{-0.5em}

\begin{Code}
# recommended initial base implementation
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots, 
                                                  knot_time=knot_time, 
                                                  debug=True,
                                                  verbose=True)[:3]
\end{Code}
\vspace{-2.5em}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ Preprocessing flag in AdvEMDpy Package} \label{sec:Preprocessing}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The raw time series signals that can be passed to the EMD package may contain a wide variety of structures. This could include different non-linear structures, non-stationarity features in trend, volatility etc. as well as corruption by observation noise that may be present in the collection of the data. Therefore, in the \textbf{AdvEMDpy} package there is the option to pre-process the time series data input to reduce the effect of these features if they may be adverse to the user's analysis. This is particularly the case in the context of signals observed in various noise environments.

\newpage

Therefore, before the sifting algorithm is implemented, there are several choices available for the pre-processing of the raw time series. This was a necessary inclusion in the algorithm owing to the wide variety of time series and their potentially vastly different statistical characteristics that would need to be accommodated in the EMD basis decomposition. These preprocessing techniques may be broadly grouped into filtering and smoothing. The filtering methods developed out of a necessity to mitigate the corruption of the IMFs that would otherwise arise as a result of the permeation of error that would occur due to the presence in the input time series signal of noise corruptions such as heavy-tailed noise, mixed noise (Gaussian noise with different variances), and Poisson noise. Smoothing developed out of the broader field of trend extraction amongst cyclical components where within EMD defining an individual cyclical component amongst others and within a non-stationary setting becomes increasingly challenging. The base implementation is as follows:

\vspace{-1.0em}
\begin{Code}
imfs, hts, ifs = \
    emd.empirical_mode_decomposition(knots=knots,
                                     knot_time=knot_time,
                                     initial_smoothing=True,
                                     preprocess=`none',
                                     preprocess_window_length=51,
                                     preprocess_quantile=0.9,
                                     preprocess_penalty=1,
                                     preprocess_order=13,
                                     preprocess_norm_1=2,
                                     preprocess_norm_2=1,
                                     downsample_window=`hamming',
                                     downsample_decimation_factor=20,
                                     downsample_window_factor=20)[:3].
\end{Code}
The following subsections will outline the different functionalities these flags provide in the pre-processing options that users can select in the \textbf{AdvEMDpy} package.

\vspace{-1.0em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Filtering}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Filtering uses a localising window to filter extrema out of the time series before the sifting procedure is performed. The preprocessing window length must be an odd integer as the window centres on a particular time point and calculates the new preprocessed time series accordingly. As an example, if \texttt{preprocess=`median\_filter'} and \texttt{preprocess\_window\_length=51}, the central time series point is replaced with the median of the $ 51 $ time series points. All filtering techniques are displayed in Figure \ref{fig:preprocess_filter} for quick reference with quantile filters used by Winsorization filtering included assisting with understanding.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Mean Filter - }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Each point in the time series is replaced with the mean to mitigate the effects of outliers (\texttt{preprocess=`mean\_filter'}) over the specified window width centred on the specific point. This is extremely susceptible to outliers and is therefore far less robust than the other filters, but is included for classical usage and completeness-sake. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Median Filter - }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Each point in the time series is replaced with the median of the points over the specified window width (\texttt{preprocess=`median\_filter'}). This technique is far less susceptible to outliers with it forming the most basic technique in a toolbox of robust statistics.

\paragraph{Winsorization - }

Both this method and the one below are based on the ground-breaking work in \cite{Has47}. A local quantile window is created for every point such that extreme values are restricted by making them equal to the boundary values (\texttt{preprocess=`winsorize'}). Should the point be above the maximum quantile it is replaced with the maximum quantile value at that point and if it is below the minimum quantile it is replaced with the minimum quantile value at that point. This is a slightly more complex robust statistical technique in that it will leave the majority of the time series untouched and simply restrict the extreme values. The level of the boundary is controlled with \texttt{preprocess\_quantile=0.9} - this would result in the maximum boundary being such that $ 95\% $ of the time series within the window is below the maximum boundary and $ 5\% $ of the time series within the window is below minimum boundary - i.e. $ 90\% $ of the time series is within the boundaries.

\paragraph{Winsorization Interpolation - }

As above, a local quantile window is created for every point, except the extreme values are treated in a slightly different manner than being made equal to the boundary values (\texttt{preprocess=`winsorize\_interpolate'}). All the points above or below the quantile windows are removed and they are linearly interpolated - this further removes extremes from the time series.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/preprocess_filter.png}
	\caption{Figure demonstrating family of filtering methods available to the user with quantile envelope included for clarifying Winsorization calculation.}
	\label{fig:preprocess_filter}
\end{figure}

\vspace{-3.0em}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Smoothing}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Rather than explicitly limiting the effects of extrema by filtering the original time series, one can fit a smoothing spline directly to the time series that implicitly limits the effect of extrema, whilst simultaneously removing less extreme noise present in the time series. The smoothing preprocessing options are displayed in Figure \ref{fig:preprocess_smooth} where both an undecimated and decimated downsampled time series is shown.

\paragraph{Smoothing - }

If an initial smoothing (\texttt{initial\_smoothing=True}) is not done, then the first IMF extracted will not be explainable in terms of cubic B-spline coefficients. It will mostly consist of random noise in the system with little-to-no physical significance as corroborated by numerous papers. This smoothing is implemented separately from other preprocessing techniques to allow smoothing to be performed on a time series once some preprocessing has been done. A natural example promoting this necessity would be a time series with mixed-noise.

\paragraph{Generalised Hodrick-Prescott - }

This is a generalised form of Hodrick-Prescott filtering (\texttt{preprocess=`HP'}), introduced in \cite{Hod97}, but in this context it is classified under smoothing:

\begin{equation}
    min_{z}\Bigg(\sum_{t=0}^{T}(y_t - z_t)^{n_1} + \lambda\sum_{t=0}^{T-d}(\nabla^d{z_t})^{n_2}\Bigg),
\end{equation}

where in the original Hodrick-Prescott filter assumes values $ n_1 = 2 $ (\texttt{preprocess\_norm\_1=2}), $ n_2 = 2 $ (\texttt{preprocess\_norm\_2=1}), the penalty term $ \lambda = 1 $ (\texttt{preprocess\_penalty=1}) was not originally present, and $ d = 2 $ (\texttt{preprocess\_order=2}) is the order of differencing with $ \nabla $ being the differencing function such that $ \nabla{z_t} = z_{t+1} - z_t $, $ \nabla^2{z_t} = \nabla(z_{t+1} - z_t) =  z_{t+2} - 2z_{t+1} + z_t $, and so forth.

\paragraph{Henderson-Whittaker Graduation - }

The well-known classical Whittaker-Henderson graduation method (\texttt{preprocess=`HW'}), introduced in \cite{Hen16}, \cite{Whi22}, and \cite{Hen24}, uses appropriately chosen window widths and weights to graduate the time series to approximate a cubic spline. This may also be generalised to a higher or lower order of smoothing with appropriately chosen constraints (such as in Equation \eqref{HW_objective} and Equation \eqref{HW_constraints} below for a cubic spline). In the interest of being able to generalise, the classical Henderson-Whitaker weighting parameters are derived using the following objective function with the order, $ d $, being such that $ D = d + 6 $:

\vspace{-1.0em}
\begin{equation} \label{HW_objective}
     min_{\omega}\Bigg(\sum_{t=0}^{D-4}(\nabla^3{\omega_t})^2\Bigg),
\end{equation}
subject to the constraints:

\vspace{-1.0em}
\begin{equation} \label{HW_constraints}
\begin{split}
    \sum_{t=0}^{D-1}{\omega_t} = 1,\text{ }\sum_{t=0}^{D-1}{\Bigg(t - \frac{D-1}{2}\Bigg)^2\omega_t} = 0,\text{ }\sum_{t=0}^{D-1}{\Bigg(t - \frac{D-1}{2}\Bigg)\omega_t} = 0,\\
    \omega_0 = 0,\text{ }\omega_1 = 0,\text{ }\omega_2 = 0,\text{ }\omega_{D-3} = 0,\text{ }\omega_{D-2} = 0,\text{ and }\omega_{D-1} = 0,
\end{split}
\end{equation}
with the final cubic spline inducing weighting parameters then being the appropriately chosen subset $ \{\omega_3, \dots, \omega_{D-4}\} $. The derivation of the asymmetric weights at the edges of the times series have been researched extensively by \cite{Mus64a}, \cite{Mus64b}, \cite{Dag06}, \cite{Bia06}, and \cite{Dag08}. This problem is related to the present issue of dealing with the edges discussed in Section \ref{sec:Edge_effects}. This problem goes back much further as evidenced by \cite{DeF77}.

\paragraph{Downsampling - }

Certain time series may benefit from varying degrees of downsampling (\texttt{preprocess=`downsample'} or \texttt{preprocess=`downsample\_decimate'}). Downsampling and the associated process of decimation are thoroughly reviewed in \cite{Cro83}. This process involves two separate, but related transforms. The first transform uses a normalised sinc filter augmented with an appropriate windowing or tapering function - the most common of which are Hamming and Hann windows - to downsample or truncate the frequency content of the time series. A sinc function is an idealised filter that removes all frequency content above a certain frequency, but since this sinc function is tapered using a Hamming window or a Hann window, it is merely a close approximation of this idealised filter. Once the frequency content has been truncated, the time series should undergo a decimation step where an appropriate subset of the time series can be sampled (\texttt{preprocess=`downsample\_decimate'}) to reflect the lower frequency content of the downsampled time series. In the smoothing context, to keep the lengths of the time series consistent with the original unpreprocessed time series, it is advised to perform only the downsampling (\texttt{preprocess=`downsample'}).

Having incorporated downsampling as a preprocessing technique into EMD, Compressive Sampling, introduced in \cite{Can06a}, \cite{Can06b}, and \cite{Can08}, should also be integrated into EMD to better isolate IMFs from sparsely sampled time series. This technique relies on the time series being sparse in some domain, such as the frequency domain. A simple sinusoid is exceptionally sparse in the frequency domain, where it is represented as two Dirac functions at both the appropriate negative and positive frequencies. EMD would most certainly benefit from the integration of Compressive Sampling. The sparsity problem becomes more complicated when the IMFs are amplitude and frequency modulated. This technique is called Compressive Sampling Empirical Mode Decomposition (CSEMD) where amplitude and frequency modulated waves are translated into a sparse domain to fit sparsely sampled time series from which IMFs can be extracted using EMD. This technique is experimental and is included with the package in \texttt{emd\_experimental.py} for completeness sake.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/preprocess_smooth.png}
	\caption{Figure demonstrating family of smoothing methods available to the user with both decimated and undecimated downsampling being displayed.}
	\label{fig:preprocess_smooth}
\end{figure}

\vspace{-2.0em}
\subsubsection{ Recommendations for Preprocessing flag in AdvEMDpy Package}

The preprocessing required depends upon the level of noise present in the time series and is very time series and process dependent. There is no preprocessing to suit all time series, but the most robust version of preprocessing would be median filtering the noise out of the time series (\texttt{preprocess=`median\_filter'} with an appropriate window length. The window length is dependent upon the highest frequency present in the time series and the sampling rate. Following the median filter, the user should then perform an initial smoothing (\texttt{initial\_smoothing=True}) to remove the discontinuities introduced by the median filter.

\begin{Code}
imfs, hts, ifs = \
    emd.empirical_mode_decomposition(knots=knots,
                                     knot_time=knot_time,
                                     initial_smoothing=True,
                                     preprocess=`median_filter',
                                     preprocess_window_length=51)[:3]
\end{Code}

\vspace{-2.0em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ Edge Effects in AdvEMDpy Package} \label{sec:Edge_effects}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The most prevalent challenge in EMD is the propagation of errors throughout the IMF and the permeation of the errors throughout the sifting process. Owing to the iterative nature of the algorithm, incorrectly estimated edge leads to errors being ubiquitous throughout all the IMFs. There are several techniques used to estimate the extrema beyond the edges of the signal - this \textbf{AdvEMDpy} package attempts to give the reader as many reasonable and researched options as possible. The base implementation of the edge effects code is as follows:

\newpage

\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  edge_effect=`symmetric', 
                                                  sym_alpha=0.1,
                                                  nn_m=200, 
                                                  nn_k=100, 
                                                  nn_method=`grad_descent',
                                                  nn_learning_rate=0.01, 
                                                  nn_iter=100)[:3].
\end{Code}

\vspace{-2.0em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Symmetric Methods} \label{symmetric_methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Without loss of generality, the procedure will be explained for the right edge maxima of the signal. Figure \ref{fig:edge_effect_symmetry_anti} demonstrates examples of the three possible symmetric edge effects. In no particular order, the technique demonstrated in \cite{Ril03} and \cite{Wu08} is referred to in this paper as the Symmetric Discard technique owing to the discarding of the end of time series in the new extrema approximation.

\paragraph{Symmetric Discard - } \label{symmetric_discard}

In Figure \ref{fig:edge_effect_symmetry_anti} below, an implementation of the Symmetric Discard technique (\texttt{edge\_effect=`symmetric\_discard'}) maxima is shown with a purple dot. The end of the time series between the last blue dot and the orange dot is disregarded when approximating the next extreme. Taking the last maxima value, $ X(t^{max}_M) $, the associated time point, $ t^{max}_M $, the last minima value, $ X(t^{min}_m) $, and the associated time point, $ t^{min}_m $, the next extreme is calculated as:

\vspace{-1.0em}
\begin{equation}
	X(t^{max}_{M+1}) = X(t^{max}_M),
\end{equation}
with the associated time point being calculated as:

\vspace{-1.0em}
\begin{equation}
	t^{max}_{M+1} = t^{min}_{m} + (t^{min}_{m} - t^{max}_{M}).
\end{equation}

\paragraph{Symmetric Anchor - } \label{symmetric_anchor}

In Figure \ref{fig:edge_effect_symmetry_anti} below, an implementation of the Symmetric Anchor technique (\texttt{edge\_effect=`symmetric\_anchor'}) maxima is shown with an orange dot. In \cite{Zha01} and \cite{Zen04} the technique creates an extreme at the endpoint - this is why this technique is referred to in this paper as the Symmetric Anchor technique. In this paper, this technique has been generalised to conditionally create an extreme depending on the difference in vertical displacement between the last two extrema and the difference in vertical displacement between the last extrema and the end of the signal - this can be referred to as the Conditional Symmetric Anchor technique. The Conditional Symmetric Anchor is calculated as follows - if $ \beta{L} \geq (1 - \alpha)L $ where $ \beta $ is the ratio of the vertical height between the last extrema and the end of the time series to the vertical height between the last two extrema, $ L $ is the vertical height between the last two extrema, and $ \alpha $ (\texttt{sym\_alpha=0.1}) is the significance level input, then:

\vspace{-1.0em}
\begin{equation}
	X(t^{max}_{M+1}) = X(t_N),
\end{equation}
with the associated time point being calculated as:

\vspace{-1.0em}
\begin{equation}
	t^{max}_{M+1} = t_N.
\end{equation}
The above is the method followed (without conditions) in \cite{Zha01} and \cite{Zen04}). If, however, $ \beta{L} < (1 - \alpha)L $, then:

\vspace{-1.0em}
\begin{equation}
X(t^{max}_{M+1}) = X(t^{max}_M),
\end{equation}
with the associated time point being calculated as:

\vspace{-1.0em}
\begin{equation}
t^{max}_{M+1} = t_N + (t_N - t^{max}_{M}).
\end{equation}

\newpage

The Conditional Symmetric Anchor technique collapses to the Symmetric Anchor technique when $ \alpha = 1 $. The other extreme where $ \alpha = -\infty $ leads to the following method.

\paragraph{Symmetric - } \label{symmetric}

The Symmetric technique (\texttt{edge\_effect=`symmetric'}) does not anchor the extrema envelope to the ends of the signal under any condition and is equivalent to the Conditional Symmetric Anchor technique where $ \alpha = -\infty $. The values are calculated as follows:

\vspace{-1.0em}
\begin{equation}
X(t^{max}_{M+1}) = X(t^{max}_M),
\end{equation}
with the associated time point being calculated as:

\vspace{-1.0em}
\begin{equation}
t^{max}_{M+1} = t_N + (t_N - t^{max}_{M}).
\end{equation}
This point is denoted with a gray dot in Figure \ref{fig:edge_effect_symmetry_anti}.

\paragraph{Anti-Symmetric - } \label{anti_symmetric}

The Anti-Symmetric (\texttt{edge\_effect=`anti-symmetric'}) approach reflects the uni-variate signal about both axes - the approximated maximum will be the reflected minimum. It is formally calculated as:

\vspace{-1.0em}
\begin{equation}
X(t^{max}_{M+1}) = X(t_N) + (X(t_N) - X(t^{min}_m)),
\end{equation}
with the associated time point being calculated as:

\vspace{-1.0em}
\begin{equation}
t^{max}_{M+1} = t_N + (t_N - t^{min}_{m}).
\end{equation}
This point is denoted with a green dot in Figure \ref{fig:edge_effect_symmetry_anti}. This technique is a more practical variation of that proposed in \cite{Zen04}, where the points are reflected about the axis rather than about the endpoints.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth, height = 8cm]{edge_effects_symmetry_anti.png}
	\caption{Example time series demonstrating four (five if Conditional Symmetric Anchor is included) different symmetric edge effect techniques with axes of symmetry included.}
	\label{fig:edge_effect_symmetry_anti}
\end{figure}

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Slope-Based Methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The slope-based methods offer a different approach to the edge effect problem. These methods use the calculated slopes and difference in temporal location between the extrema nearest to the edges to approximate the subsequent slopes and temporal locations and therefore the location of the subsequent points.

\paragraph{Slope-Based - }

This technique (\texttt{edge\_effect=`slope\_based\_method'}), which relies on the four nearest extrema, is introduced in \cite{Dat04}. The slopes ($ s_1 $ and $ s_2 $) can be calculated as follows:

\vspace{-1.0em}
\begin{equation}
	s_1 = \frac{X(t^{max}_{M}) - X(t^{min}_{m-1})}{t^{max}_{M} - t^{min}_{m-1}},
\end{equation}
and

\begin{equation}
s_2 = \frac{X(t^{max}_{M}) - X(t^{min}_{m})}{t^{max}_{M} - t^{min}_{m}}.
\end{equation}

With $ \Delta{t^{max}_M} = t^{max}_M - t^{max}_{M-1} $, the slope-based maximum ($ X(t^{max}_{M+1}) $) is calculated as follows:

\vspace{-1.0em}
\begin{equation}
	t^{max}_{M+1} = t^{max}_{M} + \Delta{t^{max}_M},
\end{equation}
and
\begin{equation}
	X(t^{max}_{M+1}) = X(t^{min}_{m}) + s_1\times(t^{max}_{M+1} - t^{min}_{m}).
\end{equation}
With $ \Delta{t^{min}_m} = t^{min}_m - t^{min}_{m-1} $, the slope-based minimum ($ X(t^{min}_{m+1}) $) is calculated as follows:

\vspace{-0.5em}
\begin{equation}
t^{min}_{m+1} = t^{min}_{m} + \Delta{t^{min}_m},
\end{equation}
and
\begin{equation}
X(t^{min}_{m+1}) = X(t^{max}_{M+1}) + s_2\times(t^{min}_{m+1} - t^{max}_{M+1}).
\end{equation}

\paragraph{Improved Slope-Based - }

This (\texttt{edge\_effect=`improved\_slope\_based\_method'}) is introduced in \cite{Wu08} as an improvement upon the slope-based method that takes into account the end of the signal and is therefore a conditional edge effect much like the Symmetric Anchor edge effect. With $ s_1 $, $ s_2 $, $ \Delta{t^{max}_M} $, $ \Delta{t^{min}_m} $, and $ X(t^{max}_{M+1}) $ defined as above the Improved Slope-Based maximum is calculated as follows: If $ X(t_N) < X(t^{max}_{M+1}) $, then proceed as in Slope-Based method, otherwise, set: 

\vspace{-1.0em}
\begin{equation}
	t^{max}_{M+1} = t_N,
\end{equation}
and
\begin{equation}
	X(t^{max}_{M+1}) = X(t_N).
\end{equation}
The Improved Slope-Based minimum is then calculated as:

\vspace{-0.5em}
\begin{equation}
t^{min}_{m+1} = t^{min}_{m} + \Delta{t^{min}_m},
\end{equation}
and
\begin{equation}
X(t^{min}_{m+1}) = X(t^{max}_{M+1}) + s_2\times(t^{min}_{m+1} - t^{max}_{M+1}).
\end{equation}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{edge_effects_slope_based.png}
	\caption{Example time series demonstrating two different slope-based edge effect techniques.}
	\label{fig:edge_effects_slope_based}
\end{figure}

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Characteristic Wave Methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This family of methods uses the extrema at the edge of the time series to estimate amplitudes and periods to approximate the next extrema using a sinusoidal wave. The original Huang Characteristic Wave approach is discussed in \cite{Hua98} and \cite{Wu08} - the Modified Huang Characteristic Wave method takes the trend of the edge implicitly into account when approximating the next extrema. The Coughlin Characteristic Wave approach is discussed in \cite{Cou04} and \cite{Wu08} - this method uses only the last two extrema and approximates the next extrema using a simple sinusoidal extension.

\paragraph{Modified Huang Characteristic Wave - } \label{huang_wave}

The descriptions of the original technique in \cite{Hua98} and \cite{Wu08} are very nebulous and are open to interpretation, but by reviewing \cite{Hua98}, \cite{Cou04}, and \cite{Wu08} a reasonable method (\texttt{edge\_effect=`characteristic\_wave\_mod\_Huang'}) can be inferred and modified with some confidence. By calculating $ p_2 = 2\times(t^{max}_{M-1} - t^{min}_{m-1}) $, $ p_1 = 2\times(t^{max}_{M} - t^{min}_{m}) $,  $ a_2 = \frac{(X(t^{max}_{M-1}) - X(t^{min}_{m-1}))}{2} $, and $ a_1 = \frac{(X(t^{max}_{M}) - X(t^{min}_{m}))}{2} $, as seen in Figure \ref{fig:edge_effects_characteristic_wave}, the next maximum may be estimated as follows:

\vspace{-0.5em}
\begin{equation}
	t^{max}_{M+1} = t^{min}_{m} + \frac{p_1}{p_2}\times(t^{max}_{M} - t^{min}_{m-1}),
\end{equation}
and
\begin{equation}
	X(t^{max}_{M+1}) = X(t^{min}_{m}) + \frac{a_1}{a_2}\times(X(t^{max}_{M}) - X(t^{min}_{m-1})),
\end{equation}
with the next minimum being estimated as:

\vspace{-0.5em}
\begin{equation}
t^{min}_{m+1} = t^{max}_{M+1} + \frac{p_1}{p_2}\times(t^{min}_{m} - t^{max}_{M}),
\end{equation}
and
\begin{equation}
X(t^{min}_{m+1}) = X(t^{max}_{M+1}) + \frac{a_1}{a_2}\times(X(t^{min}_{m}) - X(t^{max}_{M})).
\end{equation}

\paragraph{Coughlin Characteristic Wave - } \label{coughlin_wave}

The Coughlin Characteristic Wave is a less intensive technique (\texttt{edge\_effect=`characteristic\_wave\_Coughlin'}) and only uses $ t^{max}_{M} $, $ t^{min}_{m} $, $ X(t^{max}_{M}) $, and $ X(t^{min}_{m}) $ - the last two extrema. This method disregards any trend in the structure of the edges, but is simpler to calculate. With $ a_1 $ and $ p_1 $ as before, the next maximum is calculated as follows:

\vspace{-0.5em}

\begin{equation}
t^{max}_{M+1} = t^{min}_{m} + \frac{p_1}{2},
\end{equation}
and
\begin{equation}
X(t^{max}_{M+1}) = X(t^{max}_{M}),
\end{equation}
with the next minimum being estimated as:
\begin{equation}
t^{min}_{m+1} = t^{max}_{M+1} + \frac{p_1}{2},
\end{equation}
and
\begin{equation}
X(t^{min}_{m+1}) = X(t^{min}_{m}).
\end{equation}
This is equivalent to a simple sinusoidal wave with amplitude $ a_1 $ and period $ p_1 $. The maximum derived in this manner is equivalent to the maximum derived using the Symmetric Discard technique, but this is where the equivalency ends as the minimums are significantly different.

\paragraph{Average Characteristic Wave - }

The Average Characteristic Wave plot is not included in Figure \ref{fig:edge_effects_characteristic_wave} as there is no amplitude defined for either the maxima characteristic wave or the minima characteristic wave, but the calculation of the maximum and minimum rely on the characteristic periods and extrema (and hence indirectly on the amplitude) and should be included in this class of method. Unlike the other characteristic wave methods (\texttt{edge\_effect=`average'}), the Average Characteristic Wave method implicitly uses two separate characteristic waves for the maximum and minimum calculation respectively - this can be seen in the calculations below as the maximum calculation does not refer to the minima and vice versa. This method is used in \cite{Chi05}. The average maximum is calculated as follows:

\vspace{-0.5em}

\begin{equation}
t^{max}_{M+1} = t^{max}_{M} + (t^{max}_{M} - t^{max}_{M-1}),
\end{equation}
and
\begin{equation}
X(t^{max}_{M+1}) = \frac{(X(t^{max}_{M}) + X(t^{max}_{M-1}))}{2},
\end{equation}
with the average minimum being calculated as follows:

\vspace{-0.5em}

\begin{equation}
t^{min}_{m+1} = t^{min}_{m} + (t^{min}_{m} - t^{min}_{m-1}),
\end{equation}
and
\begin{equation}
X(t^{min}_{m+1}) = \frac{(X(t^{min}_{m}) + X(t^{min}_{m-1}))}{2}.
\end{equation}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{edge_effects_characteristic_wave.png}
	\caption{Example time series demonstrating three different characteristic wave edge effect techniques with average characteristic waves excluded.}
	\label{fig:edge_effects_characteristic_wave}
\end{figure}

\vspace{-2.0em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Explicit Methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

All of the above edge effect techniques may be further sub-classified under implicit edge effect methods as the time series itself is not explicitly extrapolated to find the next extrema from the extrapolated time series. Instead, only the nearest two or four extrema are used to approximate where the next two will occur. Another approach is to explicitly extrapolate the time series and extract the resulting extrema. Only one method in this family will be explored, but there are many more.

\newpage

\vspace{-1.0em}

\paragraph{Single Neuron Neural Network - }

The method (\texttt{edge\_effect=`neural\_network'}) demonstrated here follows the work done in \cite{Den01}. The objective function is:

\vspace{-0.75em}
\begin{equation} \label{neural_network}
    min_{\bar{w}}(\bar{w}P - \bar{y}),
\end{equation}

with

\vspace{-1.5em}
\begin{equation}
\bar{w} = 
\begin{bmatrix}
    w_1 & w_2 & \dots & w_k
\end{bmatrix},
\end{equation}

\vspace{-1.0em}

\begin{equation}
    P = 
\begin{bmatrix}
    y_{-(m+k)} & y_{-(m+k-1)} & \cdots & y_{-(k+1)}\\
    y_{-(m+k-1)} & y_{-(m+k-2)} & \cdots & y_{-(k)}\\
    \vdots & \vdots & \ddots & \vdots\\
    y_{-(m+1)} & y_{-m} & \cdots & y_{-2}
\end{bmatrix},
\end{equation}

and

\vspace{-1.5em}
\begin{equation}
    \bar{y} = 
\begin{bmatrix}
    y_{-m} & y_{-(m-1)} & \dots & y_{-1}
\end{bmatrix}.
\end{equation}

In the absence of an activation function or the use of the identity function as an activation function a single neuron neural network simply reduces to linear regression as can be seen in Equation \eqref{neural_network}. The naming convention is kept for consistency with the method proposed in \cite{Den01}. In this problem $k$ (\texttt{nn\_k=100}) is the number of previous points (number of parameters) needed to estimate the next point in the time series and $m$ (\texttt{nn\_m=200}) is the number of target outputs used to estimate the parameters. Two different techniques may be used to estimate $ \bar{w} $ - either gradient descent (\texttt{nn\_method=`grad\_descent'}) or steepest descent (\texttt{nn\_method=`steep\_descent'}). The optimisation problem is arranged in the same manner for both techniques. The weights are initialised ($\bar{w_0}$), a learning rate is stipulated ($l$) (\texttt{nn\_learning\_rate=0.01}), and a fixed number of iterations is set (\texttt{nn\_iter=100}). The algorithm proceeds as follows:

\begin{Code}
w = w_0

for iter in iterations:

    y_0 = np.matmul(w, P)
    e = (y - y_0)
    grad = - e * P
    av_grad = np.mean(grad, axis=1)
    if method == `grad_descent':
        adj = - l * av_grad
    elif method == `steep_descent':
        max_grad = av_grad * (np.abs(av_grad) == max(np.abs(av_grad)))
        adj = - l * max_grad
    w += adj.
\end{Code}

\vspace{-0.5em}

The gradient descent optimisation method requires a lower learning rate than the steepest descent optimisation method as all the weights are optimised simultaneously in gradient descent in each iteration, whereas only a single weight is optimised in each iteration during steepest descent. Depending on the complexity of the time series, the weighting vector obtained using steepest descent may be very sparse and this could be desirable in certain situations.

\vspace{-1.00em}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{neural_network.png}
	\caption{Example time series demonstrating single neuron neural network edge effect with $m=200$ and $k=100$ with forecasted time series in green.}
	\label{fig:neural_network}
\end{figure}

\vspace{-3.0em}

\subsubsection{ Recommendations for Edge Effects in AdvEMDpy Package}

The edge effect is the most diverse method in this paper with many options available. No method suits every time series. The most robust methods would be the Conditional Symmetric Anchor method or the Coughlin Characteristic Wave. Both of these methods only use the last two extrema to forecast the next extrema, but the Conditional Symmetric Discard method explicitly takes into account the end of the time series. The Conditional Symmetric Anchor method is therefore recommended. The Single Neuron Neural Network has not been as extensively studied and robustified as the other method and it is advised to be used with caution.

\begin{Code}
imfs, hts, ifs = \
    emd.empirical_mode_decomposition(knots=knots,
                                     knot_time=knot_time,
                                     edge_effect=`symmetric_anchor', 
                                     sym_alpha=0.1)[:3]
\end{Code}

\vspace{-1.0em}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ Detrended Fluctuation Analysis in AdvEMDpy Package} \label{sec:DFA}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

At EMD's core, the algorithm iterates by progressively removing the local mean from the IMF candidate. Owing to the numerous ways of defining and estimating a local mean, local mean estimation should be referred to as detrended fluctuation analysis. Detrended fluctuation analysis was originally introduced in \cite{Pen94} and intended for a different purpose where data was partitioned into equal sets and detrended using discontinuous linear approximations of the local trends. The average of the variances of the detrended data in these partitioned sets was calculated before the subsets were increased in size and the process was repeated. This was done to determine the presence of long-memory processes or rather long-term correlation structures. For our purposes, each set of data (data points between two consecutive knots) is detrended using continuous cubic B-splines to approximate the local trend or local mean. The trends are extracted using local windows defined using the knot sequences. The time series can therefore be estimated as a series of bases and thought of as a sequence of locally defined segments.

In Figure \ref{fig:DFA_different_trends} detrended fluctuation analysis is demonstrated using the following time series:

\vspace{-1.5em}

\begin{equation} \label{eq:DFA}
    f(t) = \text{cos}(2t) + \text{cos}(4t) + \text{cos}(8t) + \epsilon(t),
\end{equation}
where $ \epsilon(t) \in \mathcal{N}(0,1) $ for $ t \in \{t_0, \dots, t_{1000}\} $. One can observe in Figure \ref{fig:DFA_different_trends} that as the distance between the uniformly placed knots is increased, the fitted splines cannot detect the higher-frequency structures.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/DFA_different_trends.png}
	\caption{Example time series demonstrating detrended fluctuation analysis of time series with different knot sequences resulting in different trend estimation.}
	\label{fig:DFA_different_trends}
\end{figure}

In Figure \ref{fig:DFA_different_trends} one can observe the relationship between detrended fluctuation analysis and the iterative IMF extraction method known as EMD. In the top image of Figure \ref{fig:DFA_different_trends}, one can observe the trend as the sum of all the IMFs extracted with a sufficient knot sequence. In the middle image, the trend is approximated with a knot sequence that is insufficient to capture the highest frequency structure. One can also see how closely this trend resembles the sum of IMF 2 and IMF 3 for EMD when the knot sequence is sufficient. Finally, in the bottom image, one can see how the trend is estimated with a wholly insufficient knot sequence compared with the lowest order IMF of each previous sifting procedure.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/DFA_different_trends_zoomed.png}
	\caption{Example time series demonstrating detrended fluctuation analysis of time series with different knot sequences resulting in different trend estimation.}
	\label{fig:DFA_different_trends_zoomed}
\end{figure}

In Figure \ref{fig:DFA_different_trends_zoomed} the different frequency structures are more easily visible. It can be shown that:

\begin{equation} \label{eq: DFA_limit}
    \text{lim}_{n\rightarrow{\infty}}\sqrt{\frac{1}{n-1}\sum_{k=1}^{n}\text{cos}^2\bigg(\frac{2\pi{fk}}{n}\bigg)} = \frac{1}{\sqrt{2}} \approx 0.707,
\end{equation}
with $ f \in \mathbbm{N} $ and $ f << n $. With this in mind, and noting the structure of Equation \eqref{eq:DFA}, one can understand the following results:

\begin{equation} \label{eq: DFA_SD}
\begin{split}
    SD\bigg(f(t) - \sum_{i=1}^3{IMF^{51}_i(t)}\bigg) &= 1.005 \approx 1\\
    SD\bigg(f(t) - \sum_{i=1}^2{IMF^{31}_i(t)}\bigg) &= 1.558\\
    SD\bigg(f(t) - IMF^{11}_1(t)\bigg) &= 2.072,
\end{split}
\end{equation}
with $ IMF^{j}_i $ being the $ i^{th} $ IMF with knot sequence $ j $. The first equality in Equation \eqref{eq: DFA_SD} most closely calculates the true underlying noise present in the system as a result of the random fluctuations caused by the Standard Normally distributed Gaussian Noise. In the other equalities, the standard deviation calculations are confounded by undetected underlying high-frequency structures with their individual "standard deviations" being approximated by Equation \eqref{eq: DFA_limit}. The parallels between detrended fluctuation analysis and EMD should be noted by the user. By noting Figure \ref{fig:DFA_different_trends}, Figure \ref{fig:DFA_different_trends_zoomed}, and Figure \ref{fig:DFA_hilbert_spectrum} one can observe that as a result of the defined knot sequence that fluctuation is measured relative to some frequency band. Each IMF exists in some frequency range and the fluctuations calculated in Equation \eqref{eq: DFA_SD} measure the local fluctuation relative to some implicit frequency boundary as a result of the knot sequence defined.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/DFA_hilbert_spectrum.png}
	\caption{Hilbert spectrum of example time series demonstrating the frequencies of the three IMFs present when sufficient knots are used.}
	\label{fig:DFA_hilbert_spectrum}
\end{figure}

In Figure \ref{fig:DFA_hilbert_spectrum} the frequencies of the three constituent structures (apart from the added standard normal Gaussian noise) are visible. The highest frequency structure is perturbed by the noise as expected. EMD can be viewed as a generalisation of detrended fluctuation analysis whereby the trend is decomposed into separate frequency structures in descending order of IF using a defined knot sequence (in general).

The local mean can be estimated using the framework in Section 3.3 without any explicit or implicit smoothing. Explicit smoothing that deals simultaneously with smoothing and the edge effects is referred to in the literature as Statistical EMD (SEMD). The base implementation of the detrended fluctuation analysis is as follows:

\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  smooth=True,
                                                  smoothing_penalty=0.1,
                                                  dft=`envelopes', 
                                                  order=15, 
                                                  increment=10)[:3].
\end{Code}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Statistical EMD (P-Splines)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This method (\texttt{smooth=True}) is introduced in \cite{Kim12}. SEMD, in the cubic B-spline envelope setting, can be implemented by introducing a smoothing parameter into the objective function, Equation (14) of `Package AdvEMDpy: Algorithmic
Variations of Empirical Mode Decomposition in Python'. The specific introduction of a penalty into B-splines is discussed in \cite{Eil96} and is referred to as P-splines. The matricification of the P-spline objective function can be seen below. Second-order smoothing is done on the coefficients, but this can be generalised to higher-order smoothing. The second-order smoothing of the coefficients is incorporated using the $ \boldsymbol{D} $ matrix seen below:

\vspace{-1.0em}
\begin{equation}
\boldsymbol{D} = \begin{bmatrix}
1 & -2 & 1 & 0 & 0 & \cdots & 0\\
0 & 1 & -2 & 1 & 0 & \cdots & 0\\
\vdots & & \ddots & \ddots & \ddots & & \vdots\\
0 & \cdots & 0 & 1 & -2 & 1 & 0\\
0 & \cdots & 0 & 0 & 1 & -2 & 1\\
\end{bmatrix},
\end{equation}
and with $ \boldsymbol{P} $ defined as below,

\vspace{-1.0em}
\begin{equation} \label{P_matrix}
    \boldsymbol{P} = \boldsymbol{D}^T\boldsymbol{D},
\end{equation}
Equation Equation (14) of `Package AdvEMDpy: Algorithmic
Variations of Empirical Mode Decomposition in Python', with $ \boldsymbol{s} $, $ \boldsymbol{B} $, and $ \boldsymbol{c} $ defined as before and with discrete penalty term $ \lambda $, becomes:

\vspace{-1.0em}
\begin{equation} \label{DPMSE}
    DPMSE(\boldsymbol{c}|\boldsymbol{s}) = (\boldsymbol{s} - \boldsymbol{B}\boldsymbol{c})^T(\boldsymbol{s} - \boldsymbol{B}\boldsymbol{c}) + \lambda\boldsymbol{c}^T\boldsymbol{P}\boldsymbol{c}.
\end{equation}
The magnitude of the penalty term, $ \lambda $ (\texttt{smoothing\_penalty=0.1}), determines the amount of smoothing. SEMD is implemented separately to the other detrended fluctuation techniques as smoothing can be done irrespective of the detrended fluctuation technique used.

It is highly recommended that \texttt{smooth=True} and the smoothing penalty is non-zero when \texttt{dft=`envelopes'} as the extrema are not guaranteed to satisfy the Schoenberg–Whitney Conditions (SWC) for an arbitrary knot sequence. This will allow envelopes to be fitted even when there are no extrema between successive knot points, otherwise nonsensical envelopes may result. For an arbitrary set of extrema $ \mathbf{y} = \{y_1, y_2, \dots, y_j\} $, a cubic B-spline knot sequence $ \boldsymbol{\tau} = \{\tau_1, \tau_2, \dots, \tau_{j+4}\} $, the SWC can be stated as:

\vspace{-2.0em}
\begin{equation}
    \tau_i \leq y_i, \leq \tau_{i+1}\text{ }\forall\text{ }i \in \{1, 2, \dots, j-1\}.
\end{equation}

In Figure \ref{fig:Schoenberg_Whitney_Conditions} the following time series is plotted to demonstrate the necessity for either the SWC to be satisfied or for the envelopes to be smoothed:

\vspace{-1.0em}
\begin{equation}
    g(t) = \text{cos}(t) + \text{cos}(5t),
\end{equation}
with $ t \in [0, 5\pi] $. With $ \lambda = 0 $ in Equation \eqref{DPMSE} and with the SWC not being satisfied by either the maxima or the minima, the envelopes are stretched towards zero. This is, unfortunately, unavoidable without either the SWC being satisfied or some form of smoothing.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/Schoenberg_Whitney_Conditions.png}
	\caption{Example time series demonstrating unsmoothed extrema envelopes being fitted when SWC are not satisfied resulting in nonsensical envelopes.}
	\label{fig:Schoenberg_Whitney_Conditions}
\end{figure}

\vspace{-1.0em}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Enhanced EMD}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In \cite{Kop07b} Enhanced EMD (\texttt{dft=`enhanced'}) is introduced. In \cite{Kop07a} and \cite{Kop08} a genetic search algorithm is used to optimise the interpolation point allocation to optimise the local mean estimation when compared against the known local mean but in \cite{Kop07b}, the EMD is performed on the derivative of the time series to isolate the highest frequency component. Using the highest frequency component of the derivative, the optimal interpolation point allocation can be estimated for the estimation of the local mean.

\newpage
\vspace{-1.0em}
The only difference to the above is the preprocessing (in a sense) of the time series to estimate the optimal extrema points. The derivative of the time series can be calculated using the first forward difference:

\vspace{-1.0em}
\begin{equation}
	\boldsymbol{s'} = \begin{bmatrix}
	s'(t_0) = \frac{s(t_1) - s(t_0)}{t_1 - t_0}\\
	s'(t_1) = \frac{s(t_2) - s(t_1)}{t_2 - t_1}\\
	\vdots\\
	s'(t_{N-1}) = \frac{s(t_N) - s(t_{N-1})}{t_N - t_{N-1}}
	\end{bmatrix}.
\end{equation}

Given the B-spline framework developed in Section 3.3 and using a result stated in \cite{Che06} and proved in \cite{Boo78} the derivative of a B-spline basis (and therefore a spline fitted using B-splines) can be calculated as:

\vspace{-1.0em}
\begin{equation}
	B'_{j,k,\tau}(t) = \dfrac{k-1}{\tau_{j+k-1} - \tau_j}\\B_{j,k-1,\tau}(t) - \dfrac{k-1}{\tau_{j+k} - \tau_{j+1}}B_{j+1,k-1,\tau}(t).
\end{equation} 

Therefore, with

\vspace{-2.0em}
\begin{equation}
\boldsymbol{B}' = \begin{bmatrix}
B'_{0,4}(t_0) & \cdots & B'_{(M-4),4}(t_0)\\
B'_{0,4}(t_1) & \cdots & B'_{(M-4),4}(t_1)\\
\vdots & \ddots & \vdots\\
B'_{0,4}(t_{N}) & \cdots & B'_{(M-4),4}(t_{N})\\
\end{bmatrix},
\end{equation}
and $ \boldsymbol{c} $ optimised in either Equation (14) of `Package AdvEMDpy: Algorithmic Variations of Empirical Mode Decomposition in Python' or Equation \eqref{DPMSE}, the derivative of of a B-spline curve can be calculated as:

\vspace{-1.0em}
\begin{equation}
\boldsymbol{s'} = \boldsymbol{B}'\boldsymbol{c}.
\end{equation}

Once the derivative of the spline curve has been obtained (by whichever method), the EMD method is then performed on the derivative until the first IMF (highest frequency structure) is isolated. The derivative is therefore dichotomised such that:

\begin{equation}
	s'(t) = s'_h(t) + s'_l(t).
\end{equation}
Using the same technique as above such that:
\vspace{0.00em}
\begin{equation}
\begin{split}
B''_{j,k,\tau}(t) &= \dfrac{k-1}{\tau_{j+k-1} - \tau_j}B'_{j,k-1,\tau}(t) - \dfrac{k-1}{\tau_{j+k} - \tau_{j+1}}B'_{j+1,k-1,\tau}(t)\\
&= \dfrac{(k-1)(k-2)}{(\tau_{j+k-1} - \tau_j)(\tau_{j+k-2} - \tau_j)}B_{j,k-2,\tau}(t)\\ 
&- \dfrac{(k-1)(k-2)}{(\tau_{j+k-1} - \tau_j)(\tau_{j+k-1} - \tau_{j+1})}B_{j+1,k-2,\tau}(t)\\
&- \dfrac{(k-1)(k-2)}{(\tau_{j+k} - \tau_{j+1})(\tau_{j+k-1} - \tau_{j+1})}B_{j+1,k-2,\tau}(t)\\
&+ \dfrac{(k-1)(k-2)}{(\tau_{j+k} - \tau_{j+1})(\tau_{j+k} - \tau_{j+2})}B_{j+2,k-2,\tau}(t)\\
\end{split},
\end{equation}
the derivative of $ s'_h(t) $, $ s''_h(t) $, can be calculated as:

\vspace{-1.0em}
\begin{equation}
\boldsymbol{s''_h} = \boldsymbol{B}''\boldsymbol{c_h},
\end{equation}
with

\vspace{-1.0em}
\begin{equation}
\boldsymbol{B}'' = \begin{bmatrix}
B''_{0,4}(t_0) & \cdots & B''_{(M-4),4}(t_0)\\
B''_{0,4}(t_1) & \cdots & B''_{(M-4),4}(t_1)\\
\vdots & \ddots & \vdots\\
B''_{0,4}(t_{N}) & \cdots & B''_{(M-4),4}(t_{N})\\
\end{bmatrix},
\end{equation}
and $ \boldsymbol{c_h} $ being the coefficients corresponding to $ s'_h(t) $. The optimised maxima for the extraction of the first IMF can then be calculated as the points in $ s(t) $ such that:

\vspace{-1.0em}
\begin{equation}
	s'_h(t) = 0 \text{ and } s''_h(t) < 0,
\end{equation}
with the optimised minima being calculated as the points in $ s(t) $ such that:

\vspace{-1.0em}
\begin{equation}
s'_h(t) = 0 \text{ and } s''_h(t) > 0.
\end{equation}

The EMD algorithm proceeds as before with the extrema not changing until an IMF is extracted. The new optimised extrema are then calculated and the algorithm continues.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Inflection Point Interpolation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The detrended fluctuation analysis of a time series to estimate the local mean by interpolating through the interpolation points (\texttt{dft=`inflection\_points'}) was first proposed in \cite{Kop07a}. The inflection points are calculated using:

\begin{equation}
\boldsymbol{s''} = \boldsymbol{B}''\boldsymbol{c},
\end{equation}
with the inflection points being such that:

\begin{equation}
	s''(t) = 0.
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{ Binomial Average Interpolation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

All the above methods are well suited to smoothed data. If the above methods were to be directly applied to noisy data there would be a proliferation of extrema and any potentially meaningful information in the first IMF would be obscured by the noise. This is a double-edged sword as over-smoothing would also result in the loss of meaningful high-frequency information.

In \cite{Che06}, estimating the local mean structure is done by taking a binomial average (\texttt{dft=`binomial\_average'}) of the surrounding points. This may be classified as a form of preprocessing as smoothing is then required to create a local mean to be extracted, otherwise high-frequency content may be present in lower-order IMFs.

The binomial average is calculated as:

\vspace{-1.00em}
\begin{equation}
\mu(\tau_j) = \dfrac{1}{2^{k-1}}\mathlarger{\mathlarger{\sum}}_{h = j-\frac{(k-1)}{2}}^{j+\frac{(k-1)}{2}}\binom{k-1}{h-\Big(j-\frac{(k-1)}{2}\Big)}s(\tau_h),
\end{equation}
with $ \mu(\tau_j) $ being the binomial average of the time series at time point $ \tau_j $, k (\texttt{order=15}) being the order of the binomial averaging such that $ k = 2n + 1 $ with $ n \in \mathbb{N} $, and $ s(\tau_h) $ being the value of time series at time point $ \tau_j $ with $ j $ being every $ 10^{th} $ (\texttt{increment=10}) point.

All of the discussed detrended fluctuation techniques can be seen in Figure \ref{fig:detrended_fluctuation_analysis} - when the time series is well-behaved or smooth all the methods converge, but nuances still exist.

\vspace{-0.50em}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{detrended_fluctuation_analysis_resized_ssrn.png}
	\caption{Example time series demonstrating five different local mean estimation techniques through detrended fluctuation analysis.}
	\label{fig:detrended_fluctuation_analysis}
\end{figure}

\vspace{-2.50em}
\subsubsection{ Recommendations for Detrended Fluctuation Analysis in AdvEMDpy Package}

The most studied and utilised local mean estimation technique is the standard envelope technique (\texttt{dft=`envelopes'}) originally put forward in \cite{Hua98}, \cite{Hua99}, and \cite{Hua99pat} and, as already mentioned, should be performed with smoothing such that \texttt{smooth=True} and \texttt{smoothing\_penalty=0.1}. This technique is recommended and is intended to be applied after an initial preprocessing or smoothing, otherwise, the first few IMFs may be nonsensical as the noise in the time series will result in the proliferation of extrema and confound our IMFs.

\begin{Code}
imfs, hts, ifs = \
    emd.empirical_mode_decomposition(knots=knots, knot_time=knot_time,
                                     initial_smoothing=True, smooth=True,
                                     smoothing_penalty=0.1, 
                                     dft=`envelopes')[:3]
\end{Code}
\vspace{-3.00em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ Stopping Criteria in AdvEMDpy Package} \label{sec:Stopping_criteria}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

To prevent over-sifting resulting in physically meaningless IMFs and little discernible information about the process under observation, several stopping criteria are provided to prevent over-sifting. The validity of the various stopping criteria warrants further study as some are more related to algorithmic steps than others. The base implementation of the stopping criteria follows: 
\vspace{-2.00em}
\begin{Code}
imfs, hts, ifs = \ 
    emd.empirical_mode_decomposition(knots=knots, knot_time=knot_time,
                                     stop_crit=`S_stoppage',
                                     stop_crit_threshold=10,
                                     mft_theta_1=0.05, mft_theta_2=0.5,
                                     mft_alpha=0.05, mean_threshold=10, 
                                     max_internal_iter=30, max_imfs=10)[:3].
\end{Code}
\vspace{-0.75em}

\ref{con:Condition_1} and \ref{con:Condition_2} are restated here for ease of reference:
\vspace{-1.00em}
\begin{description}
	\item[\namedlabel{con:Condition_1}{Condition 1}] $ \text{abs}\bigg(\Big|\Big\{\frac{d\gamma_k(t)}{dt} = 0 : t \in (0, T)\Big\}\Big| - \Big|\Big\{\gamma_k(t) = 0 : t \in (0, T)\Big\}\Big|\bigg) \leq 1, $
    \item[\namedlabel{con:Condition_2}{Condition 2}] $ \tilde{\gamma}^{\mu}_k(t) = \Big(\frac{\tilde{\gamma}^M_k(t) + \tilde{\gamma}^m_k(t)}{2}\Big) = 0 \text{ }\forall\text{ }t \in [0, T] $ with,\\
	$ \gamma_k(t) = \tilde{\gamma}^M_k(t) \text{  if  } \frac{d\gamma_k(t)}{dt} = 0 \text{  and  } \frac{d^2\gamma_k(t)}{dt^2} < 0 $,\\
	$ \gamma_k(t) \leq \tilde{\gamma}^M_k(t)\text{ }\forall\text{ }t \in [0, T], $\\
	$ \gamma_k(t) = \tilde{\gamma}^m_k(t) \text{  if  } \frac{d\gamma_k(t)}{dt} = 0 \text{  and  } \frac{d^2\gamma_k(t)}{dt^2} > 0 $, and\\
	$ \gamma_k(t) \geq \tilde{\gamma}^m_k(t)\text{ }\forall\text{ }t \in [0, T]. $\\
\end{description}

\vspace{-3.50em}

\subsubsection{ Modified Mean Threshold} \label{modified_mean}

In order to be classified as an IMF, \ref{con:Condition_1} and \ref{con:Condition_2} need to both be satisfied. In practice, \ref{con:Condition_1} is easily achievable and takes relatively few iterations, but \ref{con:Condition_2} is computationally very expensive and the large number of iterations required for the local mean of an IMF candidate to approach zero would, unfortunately, remove a large portion of meaningful content from the time series. This necessitates the statement of a modified version of \ref{con:Condition_2} as stated in Section 6 and repeated here for ease of reference.
\vspace{-1.00em}
\begin{description}
    \item[\namedlabel{con:remod_Condition_2}{Modified Condition 2}] $ \sum_t\big|\tilde{\gamma}^{\mu}_k(t) \big| \leq \epsilon $, 
\end{description}
\vspace{-1.00em}
for some chosen threshold value $ \epsilon $ (\texttt{mean\_threshold=10}).

\vspace{-1.50em}
\subsubsection{ Fixed Iteration} \label{fixed_iteration}

A simple stopping criterion puts a hard limit on the number of iterations allowed before the internal iteration loop stops and the IMF is extracted (\texttt{max\_internal\_iter=10}). The fixed iteration stopping criterion and the modified mean threshold criterion are implemented in addition to other stopping criteria with them functioning as a hard limit on the number of siftings and a modification of the classical IMF condition, respectively. 

\vspace{-1.50em}
\subsubsection{ S Stoppage} \label{s_stoppage}

Much like the fixed iteration count stopping criterion, the S Stoppage, as in \cite{Hua08} (\texttt{stop\_crit=`S\_stoppage'}), is not mathematically rigorous but relies on a count variable. Unlike the Fixed Iteration count stopping criterion, the S Stoppage criterion has a conditional count variable. The stopping criterion is reached when the difference between the number of zero crossings and the number of extrema remains constant over $ S $ iterations with the number of iterations as another input (\texttt{stop\_crit\_threshold=10}).

\subsubsection{ Cauchy-Type Convergence} \label{Cauchy_type}

This stopping criterion (\texttt{stop\_crit=`sd'}) was first proposed in \cite{Hua98}. Its name results from its similarity to the Cauchy convergence of a series in functional analysis. The standard deviation of IMF candidate $ j $ at iteration $ k $ is calculated as:
\vspace{-0.25em}
\begin{equation} \label{sd_stop}
SD_{(j,k)} = \sum_{t=t_0}^{T_{N}} \Bigg[\dfrac{\big|(h_{(j,k-1)}(t) - h_{(j,k)}(t))\big|^2}{h^2_{(j,k-1)}(t)}\Bigg] < \epsilon,
\end{equation}
for some threshold value $ \epsilon $ (\texttt{stop\_crit\_threshold=10}), with $ h_{(j,k)}(t) $ being IMF candidate $ j $ on iteration $ k $. In \cite{Hua98} a threshold value of $ 0.2-0.3 $ is suggested with other works citing a value of $ 0.1 $. These values are subjective and are very dependent on the length of the time series under observation - a better description would be threshold per time point. 

\vspace{-1.00em}
\subsubsection{ Cauchy-Type Convergence 11a} \label{cauchy_11a}

The first variation (\texttt{stop\_crit=`sd\_11a'}) of the above stopping criterion proposed in \cite{Hua08} as Equation (11a) is calculated as:
\vspace{-0.50em}
\begin{equation}
SD^{11a}_{(j,k)} = \dfrac{\sum_{t=t_0}^{T_{N}}\big|(h_{(j,k-1)}(t) - h_{(j,k)}(t))\big|^2}{\sum_{t=t_0}^{T_{N}}h^2_{(j,k-1)}(t)} < \epsilon.
\end{equation}

This stopping criterion measures the relative energy difference between IMFs, rather than the Cauchy-type convergence.

\vspace{-1.00em}
\subsubsection{ Cauchy-Type Convergence 11b} \label{cauchy_11b}

The second variation (\texttt{stop\_crit=`sd\_11b'}) of the above Cauchy-type convergence stoppping criterion is also proposed in \cite{Hua08} as Equation (11b) and is calculated as:

\vspace{-1.50em}
\begin{equation} \label{11b}
SD^{11b}_{(j,k)} = \dfrac{\sum_{t=t_0}^{T_{N}}\mu^2_{(j,k)}(t)}{\sum_{t=t_0}^{T_{N}}h^2_{(j,k)}(t)} < \epsilon.
\end{equation}

Equation \eqref{11b} can be seen as a combination of the Modified Mean Threshold and the Cauchy-Type Convergence. It has standardised the allowed mean discrepancy.

\vspace{-1.00em}
\subsubsection{ Mean Fluctuation} \label{mean_fluctuation}

This stopping criterion (\texttt{stop\_crit=`mft'}) is proposed in \cite{Ril03} and modified in \cite{Tab14}. It takes into account locally large fluctuations that would otherwise cause over-sifting throughout the majority of the signal. With $ \mu(t) $ being the calculated local mean and $ a(t) $ being the mode amplitude calculated as:

\vspace{-1.00em}
\begin{equation}
	a(t) = \dfrac{|M(t) - m(t)|}{2},
\end{equation}
with $ M(t) $ being the maximum envelope and $ m(t) $ being the minimum envelope the evaluation function can be calculated. The evaluation function is calculated as follows:

\vspace{-1.00em}
\begin{equation}
	\sigma(t) = \bigg|\dfrac{\mu(t)}{a(t)}\bigg|.
\end{equation}

Unlike the other stopping criteria, this stopping criterion requires three parameters, namely $ \theta_1 $ (\texttt{mft\_theta\_1=0.05}), $ \theta_2 $ (\texttt{mft\_theta\_2=0.5}), and $ \alpha $ (\texttt{mft\_alpha=0.05}). The stopping criterion is met (and the internal loop is terminated) when $ \sigma(t) < \theta_1 $ for fraction $ (1 - \alpha) $ of the entire range of the time series, and $ \sigma(t) < \theta_2 $ for the remainder of the time series. In \cite{Ril03} typical values for the parameters are proposed as $ \theta_1 = 0.05 $, $ \theta_2 = 10\theta_1 $, and $ \alpha = 0.05 $.

\subsubsection{ Energy Difference Tracking} \label{energy_difference_tracking}

This method (\texttt{stop\_crit=`edt'}) was proposed in \cite{Jun06}. It takes advantage of the orthogonality of the IMFs and as such the cumulative energy of the IMFs should be equal to the energy of the original time series. The energy of an individual IMF candidate can be calculated as:

\vspace{-1.00em}
\begin{equation}
E_{h_{(j,k)}} = E\big[h_{(j,k)}(t)\big] = \sum_{t=t_0}^{T_{N-1}}\Big|h_{(j,k)}(t)\Big|^2\Delta{t}.
\end{equation}
Should this candidate be selected as an IMF, the energy of the potential remainder of the time series after the removal of the IMF is calculated as:
\begin{equation}
E_{r_j} = E\big[r_j(t)\big] = \sum_{t=t_0}^{T_{N-1}}\Big|r_j(t)\Big|^2\Delta{t}.
\end{equation}
The energy of the residual time series before the removal of the IMF is calculated as:
\begin{equation}
E_{r_{j-1}} = E\big[r_{j-1}(t)\big] = \sum_{t=t_0}^{T_{N-1}}\Big|r_{j-1}(t)\Big|^2\Delta{t}.
\end{equation}
The internal iteration stops when the energy difference falls below a certain threshold value $ \epsilon $:
\begin{equation}
E_{diff} = \Big|E_{r_{j-1}} - \big(E_{h_{(j,k)}} + E_{r_j}\big)\Big| < \epsilon.
\end{equation}

\vspace{-1.50em}
\subsubsection{ Recommendations for Stopping Criteria in AdvEMDpy Package}

There are several recommendations for this particular aspect of the algorithm. The most widely used stopping criterion is the Cauchy-Type Convergence, but this is often unstable and does not converge steadily such as an exponential decaying variable. One should use S Stoppage (\texttt{stop\_crit=`S\_stoppage'}) with a threshold of 10 (\texttt{stop\_crit\_threshold=10}). In addition to this criterion, one should impose a value on the Modified Mean Threshold (\texttt{mean\_threshold=10}), the maximum number of internal iterations (\texttt{max\_internal\_iter=30}), and the maximum allowed number of IMFs (\texttt{max\_imfs=10}). With all these conditions the stopping criteria aspect of the algorithm should be sufficiently managed.

\vspace{-0.50em}
\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  stop_crit=`S_stoppage',
                                                  stop_crit_threshold=10,
                                                  mean_threshold=10, 
                                                  max_internal_iter=30, 
                                                  max_imfs=10)[:3]
\end{Code}
\vspace{-3.00em}
\subsection{ Spline Methods in AdvEMDpy Package} \label{sec:Spline_methods}

All the above techniques are listed with the cubic B-spline implementation of EMD in mind. Other spline techniques are effective when using EMD. The other splines used have a different basis, but because they are both cubic bases, they may be easily mapped onto one another. The base implementation of cubic B-spline EMD is:

\vspace{-0.50em}
\begin{Code}
imfs, hts, ifs = \
    emd.empirical_mode_decomposition(knots=knots,
                                     knot_time=knot_time,
                                     spline_method=`b_spline')[:3].
\end{Code}


EMD using cubic Hermite spline interpolation (CHSI) and Akima spline interpolation (ASI) is performed in \cite{Ega16} to demonstrate their effectiveness in removing eye-blink artifacts from electroencephalograms. The basis for both CHSI and ASI is:

\begin{equation}
\begin{split}
	p(t) &= h_{00}(t^*)y(\tau_k) + h_{10}(t^*)(\tau_{k+1} - \tau_k)m(\tau_k)\\
	&+ h_{01}(t^*)y(\tau_{k+1}) + h_{11}(t^*)(\tau_{k+1} - \tau_k)m(\tau_{k+1}),
\end{split}
\end{equation}
with,

\vspace{-1.00em}
\begin{equation}
\begin{split}
	t^* &= \dfrac{t - \tau_k}{\tau_{k+1} - \tau_k},\\
	h_{00}(t^*) &= 2{t^*}^3 - 3{t^*}^2 + 1,\\
	h_{10}(t^*) &= {t^*}^3 - 2{t^*}^2 + t^*,\\
	h_{01}(t^*) &= -2{t^*}^3 + 3{t^*}^2, \text{ and}\\
	h_{11}(t^*) &= {t^*}^3 - {t^*}^2.
\end{split}
\end{equation}

\vspace{-1.00em}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{comparing_bases.png}
	\caption{Comparison of cubic B-spline bases versus cubic Hermite spline bases between two arbitrary knots of uniform knot sequence.}
	\label{fig:comparing_bases}
\end{figure}

Despite both CHSI and ASI using the same bases, the interpolation techniques differ when calculating the derivative values $ m(\tau_k) $ and $ m(\tau_{k+1}) $. The positional values for both techniques are calculated as $ y(\tau_k) = s(\tau_k) $ and $ y(\tau_{k+1}) = s(\tau_{k+1}) $. One can note the differences between the cubic Hermite bases and the cubic B-spline bases in Figure \ref{fig:comparing_bases}.

\vspace{-1.00em}
\subsubsection{ Cubic Hermite Spline Interpolation} \label{CHSI}

There are a number of methods used to calculate the tangential parameters for CHSI curve-fitting (\texttt{spline\_method=`chsi'}). With $ t_i = \tau_k $ and $ t_j = \tau_{k+1} $ being the time points corresponding to the knot points, a common method, with low computational expense, for calculating the tangential values are:

\begin{equation}
    m_{CHSI}(\tau_k) = \frac{1}{2}\Bigg(\frac{s(t_{i+1}) - s(t_{i})}{t_{i+1} - t_{i}}+\frac{s(t_{i}) - s(t_{i-1})}{t_{i} - t_{i-1}}\Bigg),
\end{equation}
and

\vspace{-1.00em}
\begin{equation}
    m_{CHSI}(\tau_{k+1}) = \frac{1}{2}\Bigg(\frac{s(t_{j+1}) - s(t_{j})}{t_{j+1} - t_{j}}+\frac{s(t_{j}) - s(t_{j-1})}{t_{j} - t_{j-1}}\Bigg).
\end{equation}

\vspace{-1.00em}
\subsubsection{ Akima Spline Interpolation} \label{ASI}

ASI (\texttt{spline\_method=`asi'}) should be viewed as a specific case of CHSI as the bases are the same with the only difference being a defined method for calculating the tangential values. With $ t_i $ and $ t_j $ being defined as above and with:

\vspace{-1.00em}
\begin{equation}
    v(t_i) = \frac{s(t_{i+1}) - s(t_{i})}{t_{i+1} - t_{i}},
\end{equation}
the tangential values for ASI are defined as:
\begin{equation}
    m_{ASI}(\tau_k) = \frac{|v(t_{i+1}) - v(t_{i})| \times v(t_{i-1}) + |v(t_{i-1}) - v(t_{i-2})| \times v(t_{i})}{|v(t_{i+1}) - v(t_{i})| + |v(t_{i-1}) - v(t_{i-2})|},
\end{equation}
and
\vspace{-0.50em}
\begin{equation}
    m_{ASI}(\tau_{k+1}) = \frac{|v(t_{j+1}) - v(t_{j})| \times v(t_{j-1}) + |v(t_{j-1}) - v(t_{j-2})| \times v(t_{j})}{|v(t_{j+1}) - v(t_{j})| + |v(t_{j-1}) - v(t_{j-2})|}.
\end{equation}

\vspace{-1.00em}
\subsubsection{ Recommendations for Spline Methods in AdvEMDpy Package}

Cubic Hermite splines and Akima splines lack second-order continuity and as such are not as smooth as cubic B-splines. This may also lead to overfitting and the higher frequency IMFs being unwanted smoothed noise. In addition to using cubic B-splines interpolation method (\texttt{spline\_method=`b\_spline'}), the user should also smooth the cubic B-splines (\texttt{smooth=True}) for reasons already stated in Section \ref{sec:DFA}.

\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  spline_method=`b_spline',
                                                  smooth=True,
                                                  smoothing_penalty=0.1)[:3]
\end{Code}
\vspace{-2.00em}

\subsection{ Discrete-Time Hilbert Transforms in AdvEMDpy Package} \label{sec:DTHT}

Despite B-splines having closed form solutions to the HT, discrete solutions to the HT may be needed for a number of reasons. Here are two widely known and used methods that are included. The base implementation of the DTHT is as follows:

\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  dtht=False,
                                                  dtht_method=`fft')[:3].
\end{Code}

\newpage

\subsubsection{ Basic DTHT}

This method (\texttt{dtht\_method=`kak'}) is stated and proved in \cite{Kak70} and is often referred to in the literature as the Basis DTHT - it is quick and accurate. The Basic DTHT is calculated as follows: 

\vspace{-1.00em}
\begin{equation}
\bar{\mathcal{H}}(X)(t[k]) \approx \begin{cases}
\sum_{n \text{ is} \text{ odd}}\dfrac{2}{\pi{(k-n)}}X(t[n]),\text{ } k  \text{ is} \text{ even}\\
\sum_{n \text{ is} \text{ even}}\dfrac{2}{\pi{(k-n)}}X(t[n]),\text{ } k  \text{ is} \text{ odd}\\
\end{cases}.
\end{equation}
\vspace{-2.00em}

\subsubsection{ FFT DTHT}

This method (\texttt{dtht\_method=`fft'}), which follows \cite{Python}, relies on the relationship that exists between the Fourier transform and the HT. The fast-Fourier transform is performed on the time series resulting in the below:

\vspace{-0.50em}
\begin{equation}
	\bar{\mathcal{F}}(X)(k) = \sum_{n = 0}^{N - 1}e^{-2\pi{i}(\frac{kn}{N})}X(n).
\end{equation}

The inverse fast-Fourier transform is then performed on twice the positive frequency components, zero times the negative frequency components, and one times the zero frequency component. There is a slight difference in implementation when the time series is of odd or even length. When the length of the time series, $ N $, is even:

\vspace{-0.50em}
\begin{equation}
\begin{split}
	\bar{\mathcal{H}}(X)(t) &\approx \dfrac{1}{N}\sum_{n = 0}^{N - 1}X(n)\\ 
	&+ \dfrac{2}{N}\sum_{k=1}^{\frac{N}{2} - 1}e^{2\pi{i}(\frac{kn}{N})}\bar{\mathcal{F}}(X)(k)\\
	&+ \dfrac{1}{N}(-1)^n\sum_{n = 0}^{N - 1}(-1)^n{X(n)},
\end{split}
\end{equation}
\vspace{-1.50em}

and when the length, $ N $, is odd:

\vspace{-1.00em}
\begin{equation}
\begin{split}
	\bar{\mathcal{H}}(X)(t) &\approx \dfrac{1}{N}\sum_{n = 0}^{N - 1}X(n)\\ 
	&+ \dfrac{2}{N}\sum_{k=1}^{\frac{N - 1}{2}}e^{2\pi{i}(\frac{kn}{N})}\bar{\mathcal{F}}(X)(k).
\end{split}
\end{equation}

\vspace{-1.00em}
\subsubsection{ Recommendations for Discrete-Time Hilbert Transforms in AdvEMDpy Package}

A DTHT is not output by default, but should the user want a DTHT as well as the closed-form cubic B-spline HT for comparison, the FFT DTHT is recommended. The FFT is significantly faster than the Basic DTHT owing to the relationship that exists between the Fourier transform and the HT. Despite this significant increase in computational speed, it is, however, slightly less accurate than the Basic DTHT - especially at the edges of the time series.

\begin{Code}
imfs, hts, ifs, _, _, discrete_time_hts, discrete_time_ifs = \
    emd.empirical_mode_decomposition(knots=knots,
                                     knot_time=knot_time,
                                     dtht=True,
                                     dtht_method=`fft')
\end{Code}

\subsection{ Knot Point Optimisations in AdvEMDpy Package} \label{sec:Knot_optimisation}

Two methods are available for the optimisation of knot point allocation. Both methods are in the bisection family of knot optimisation techniques. The first technique simply bisects the domain iteratively until some error bound is met. The next method is a variation on this that extends the domain or diminishes the domain based on some error bound. The base implementation of the knot point optimisation is as follows:

\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  optimise_knots=0,
                                                  knot_method=`ser_bisect',
                                                  knot_error=10,
                                                  knot_lamda=1,
                                                  knot_epsilon=0.5)[:3].
\end{Code}

\vspace{-2.00em}
\subsubsection{ Bisection}

The bisection method (\texttt{knot\_method=`bisection'}) is the most basic of knot point optimisation techniques in that, given an error term (\texttt{knot\_error=10}), the distance between two knots is iteratively halved until the sum of absolute errors between the spline and the time series over this interval is below the error term. In Figure \ref{fig:bisection}, given knot point $ \tau_k $, the next knot, $ \tau^{0}_{k+1} $, is assumed to be the end of the time series (to start the iterative process). If the error between the spline and the time series is too large the distance between the knots is halved.

The spline is fitted again with the next knot point, $ \tau^{1}_{k+1} $, and if the error is still too large, the distance is halved again. The spline is then fitted between $ \tau_k $ and $ \tau^{2}_{k+1} $. If the calculated error is below the maximum allowed error, then the iterative knot optimisation stops, $ \tau^{2}_{k+1} $ is accepted as an optimised knot point ($ \tau_{k+1} $). The knot optimisation process is repeated for $ \tau_{k+1} = \tau^{2}_{k+1} $ and $ \tau^{0}_{k+2} = \tau^{0}_{k+1} $.

\vspace{+0.50em}
\begin{figure}[h]
	\centering
    \begin{tikzpicture}
    \draw (0,0) -- (12,0);
    \draw (0,0.25) -- (0,-0.25) node[anchor=north] {$ \tau_k $};
    \draw (12,0.25) -- (12,-0.25) node[anchor=north] {$ \tau^{0}_{k+1} $};
    \draw[->] (12,-1) .. controls (9,-2.5) .. (6,-1);
    \draw (6,0.25) -- (6,-0.25) node[anchor=north] {$ \tau^{1}_{k+1} $};
    \draw (9,-0.5) node[anchor=north] {$ \tau^{1}_{k+1} = \frac{\tau_k + \tau^{0}_{k+1}}{2} $};
    \draw[->] (6,-1) .. controls (4.5,-2.5) .. (3,-1);
    \draw (3,0.25) -- (3,-0.25) node[anchor=north] {$ \tau^{2}_{k+1} $};
    \draw (4.5,-0.5) node[anchor=north] {$ \tau^{2}_{k+1} = \frac{\tau_k + \tau^{1}_{k+1}}{2} $};
    \end{tikzpicture}
    \caption{Diagram demonstrating bisection knot point optimisation method where distance between last accepted knot and next potential knot is continuously halved until stopping criteria satisfied.}
	\label{fig:bisection}
\end{figure}

\vspace{-2.00em}
\subsubsection{ Serial Bisection}

The serial bisection knot optimisation method (\texttt{knot\_method=`ser\_bisect'}) is introduced in \cite{Dun17}. As above, the domain is bisected until the calculated error is below the maximum allowed error. Once the error is satisfied, however, the domain between potential knot points is extended again by half the distance between the knot points. If the calculated error is now more than the maximum error, the distance is decreased by half the distance it was increased by previously, if the calculated error is still less than the maximum allowed error the distance is increased by half the distance it was previously increased and so forth until the distance by which the distance is increased or decreased is less than $ \epsilon $ (\texttt{knot\_epsilon=0.5}). 

\begin{figure}[h]
	\centering
    \begin{tikzpicture}
    \draw (0,0) -- (12,0);
    \draw (0,0.25) -- (0,-0.25) node[anchor=north] {$ \tau_k $};
    \draw (12,0.25) -- (12,-0.25) node[anchor=north] {$ \tau^{0}_{k+1} $};
    \draw[->] (12,-1) .. controls (9,-2.5) .. (6,-1);
    \draw (6,0.25) -- (6,-0.25) node[anchor=north] {$ \tau^{1}_{k+1} $};
    \draw (9,-0.5) node[anchor=north] {$ \tau^{1}_{k+1} = \frac{\tau_k + \tau^{0}_{k+1}}{2} $};
    \draw[->] (6,-1) .. controls (4.5,-2.5) .. (3,-1);
    \draw (3,0.25) -- (3,-0.25) node[anchor=north] {$ \tau^{2}_{k+1} $};
    \draw (4.5,-0.5) node[anchor=north] {$ \tau^{2}_{k+1} = \frac{\tau_k + \tau^{1}_{k+1}}{2} $};
    \draw (4.5,-0.25) -- (4.5,0.25) node[anchor=south] {$ \tau^{3}_{k+1} $};
    \draw[->] (3,0.5) .. controls (3.75,2.5) .. (4.5,1);
    \draw (3,1.5) node[anchor=east] {$ \tau^{3}_{k+1} = \frac{\tau^{1}_{k+1} + \tau^{2}_{k+1}}{2} $};
    \end{tikzpicture}
    \caption{Diagram demonstrating serial bisection knot point optimisation method where distance between last accepted knot and next potential knot can also be increased.}
	\label{fig:serial_bisection}
\end{figure}

Since a cubic B-spline is used to optimise the knot allocation, a $ \lambda $ value (\texttt{knot\_lamda=1}) is used to smooth the splines that are fitted to prevent abrupt changes in the coefficients of the B-splines to accommodate the maximum allowed error. Any spline-fitting technique may be used to optimise the knot point allocation, but B-splines are used in this package.

\subsubsection{ Recommendations for Knot Point Optimisations in AdvEMDpy Package}

There are several options and combinations of options available to the user. Depending on the length and complexity of the time series, it is advised that the algorithm is first run with a uniform set of knots (either user-defined or defaults) with \texttt{optimise\_knots=0}. Once the algorithm has run without error, the user may opt for the knots to be optimised once at the outset (\texttt{optimise\_knots=1}). The allowed error (\texttt{knot\_error=10}) and the allowed minimum knot placement distance (\texttt{knot\_epsilon=0.5}) are both very time series specific and should be adjusted with care. The Bisection method (\texttt{knot\_method=`bisection'}) is significantly faster than the Serial Bisection method (\texttt{knot\_method=`ser\_bisect'}), but does result in slightly more knot points. One should first run the Bisection method with one optimisation at the outset. If one wants knots that dynamically adjust to each potential IMF, then one can use \texttt{optimise\_knots=2}.

\vspace{+1.00em}
\begin{Code}
imfs, hts, ifs, _, optimised_knots, _, _ = \
    emd.empirical_mode_decomposition(knots=knots,
                                     knot_time=knot_time,
                                     optimise_knots=1,
                                     knot_method=`bisection',
                                     knot_error=10,
                                     knot_lamda=1,
                                     knot_epsilon=0.5)
\end{Code}

\newpage

\subsection{ Ensemble Empirical Mode Decomposition in AdvEMDpy Package} \label{sec:EEMD}

This is a noise-assisted data analysis (NADA) technique introduced in \cite{Wu09} that utilises white noise. This technique relies on the ability of the algorithm to discern consistent structures in the presence of different sets of randomised white noise. The base implementation of the EEMD is implemented as follows:

\vspace{-1.00em}
\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  ensemble=False,
                                                  ensemble_sd=0.5,
                                                  ensemble_iter=10)[:3].
\end{Code}
\vspace{-1.00em}

The original time series is median filtered to approximate the original level of noise in the system. The median is removed from the time series before the standard deviation is calculated, $ \sigma_s $. The noise added to the system has a mean of zero and a standard deviation of $ \sigma_s \times \sigma_e $, where $ \sigma_e $ (\texttt{ensemble\_sd=0.5}) is the chosen level of standard deviation in the ensemble system as a fraction of the original standard deviation in the system. The IMFs are calculated for this modified time series and stored before a new set of noise is added to the original time series and the sifting process is repeated - the procedure is repeated a predetermined number of iterations (\texttt{ensemble\_iter=10}).

Further applications of this technique will benefit from random additions of other colours of noise such as violet noise, blue noise, pink noise, and red noise (also known as Brownian noise), which all have non-constant power spectral densities, compared with white noise that has a constant power spectral density. This technique, with the random addition of other colours of noise, is called Full-Spectrum Ensemble Empirical Mode Decomposition (FSEEMD). This technique is experimental and is also included with the package in the script entitled \texttt{emd\_experimental.py} for completeness sake.

This method performs a similar task, utilising a different methodology, to ICA-EMD, introduced in \cite{van21}, which is finding the most consistent structures amongst noisy data. EEMD proceeds by introducing noise to the time series and isolating IMFs from the noisy time series. This task is repeated several times and IMFs are simply averaged. A more sophisticated sorting technique such as the Minimum Description Length Principle (introduced in \cite{Fay93}) should be used to ensure structures of similar frequency are grouped as in \cite{van21}. ICA-EMD proceeds by applying Independent Component Analysis (ICA) directly to all the noisy IMFs to isolate the most consistent structure, before EMD is performed on the ICA component to isolate IMFs.

\vspace{-1.00em}
\subsubsection{ Recommendations for Ensemble Empirical Mode Decomposition in AdvEMDpy Package}

This technique shows promise, but it should be implemented with care. The implementation of the EEMD method removes some additional features such as the output of coefficients, knots, as well as the DTHT and IF from the corresponding DTHT. One should apply this technique after applying EMD for comparison.

\vspace{-1.00em}
\begin{Code}
imfs, hts, ifs = emd.empirical_mode_decomposition(knots=knots,
                                                  knot_time=knot_time,
                                                  ensemble=True,
                                                  ensemble_sd=0.5,
                                                  ensemble_iter=10)
\end{Code}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{ Illustrative Examples of EMD Method} \label{main_AAS_Software_Supplement:Worked_example}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In these supplementary examples we provide two additional case studies. The intentions are firstly to demonstrate the application of the toolbox in a real data example and secondly to compare the accuracy of the methods implement in the proposed toolbox versus other existing EMD packages, comparing \textbf{AdvEMDpy} against \textbf{PyEMD 0.2.10} and \textbf{emd 0.3.3}. 

In Section \ref{main_AAS_Software_Supplement:Synthetic}, a synthetic example is created using a specific case of the Duffing equation as provided in \cite{Hua98} to demonstrate using these scripts how to arrive at the desired output. The accuracy of the three methods is quantitatively compared by calculating the absolute deviations of the second IMFs extracted using each method with the underlying driving function.

In Section \ref{main_AAS_Software_Supplement:Real_world}, a well-known real-world data set is used to demonstrate how EMD can be applied to real-world data. The knot sequence would need to be appropriately constructed. Should the data (unlike the Carbon Dioxide data) be noisy, some preprocessing (Package \textbf{AdvEMDpy}: Algorithmic Variations of Empirical Mode Decomposition in \textsf{Python}, Section 7.2) may be appropriate or one would benefit by not performing initial smoothing (\texttt{initial\_smoothing=False}) as all the noise would be captured in a nonsensical first IMF. As there is no true underlying IMF structure as in Section \ref{main_AAS_Software_Supplement:Synthetic}, the three methods are compared by measuring how accurately the IFs of each method replicate the underlying annual structure.

\vspace{-1.00em}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ Supplementary Example One} \label{main_AAS_Software_Supplement:Synthetic}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The motion of a mass attached to a frictionless spring in a vacuum that is perturbed slightly is governed by: 

\vspace{-1.00em}
\begin{equation} \label{spring}
    \frac{d^2x(t)}{dt^2} + {\omega}x(t) = 0,
\end{equation}
with $ \omega^2 $ (squared by convention) being the constant natural frequency of the system and $ 0 $ being the driving force of the system - it oscillates about zero. A specific Duffing equation is governed by:

\vspace{-1.00em}
\begin{equation} \label{Duffing_equation}
    \frac{d^2x(t)}{dt^2} + \Big(1 - x^2(t)\Big)x(t) = \frac{1}{10}\text{cos}\bigg(\frac{1}{25}2{\pi}t\bigg),
\end{equation}
with the natural frequency of the system being a function of the displacement and the driving force of the system no longer being zero. This creates a more complex system than a simple sinusoid (solution to Equation \eqref{spring}) - the displacement and velocity can be viewed in Figure \ref{fig:Duffing_equation}. The simple script below can be used to produce Figure \ref{fig:Duffing_equation}. For further details on how to reproduce this figure, and all others in this section, one can see \texttt{aas\_figures\_replication\_script.py} and \texttt{Worked\_examples.ipynb}.

\begin{Code}
def duffing_equation(XY, t):
    return [XY[1], XY[0] - XY[0] ** 3 + 0.1 * np.cos(0.08 * np.pi * t)]

t = np.linspace(0, 150, 1501)
XY0 = [1, 1]
solution = odeint(duffing_equation, XY0, t)
x = solution[:, 0]
dxdt = solution[:, 1]
\end{Code}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/Duffing_equation.png}
	\caption{Figure demonstrating Duffing equation (Equation \eqref{Duffing_equation}) displacement and velocity for $ t \in [0, 150] $.}
	\label{fig:Duffing_equation}
\end{figure}

In a frictionless environment, the total energy of the system is completely determined by the initial conditions of $ x(0) = 1 $ and $ \frac{dx(0)}{dt} = 1 $. The Hamiltonian frequency of the system (in the absence of a driving function) can be shown to be $ f = 0.124Hz $. This, as well as the frequency of the driving function, $ f = 0.04Hz $, can be seen in Figure \ref{fig:Duffing_equation_ht_pyemd}, Figure \ref{fig:Duffing_equation_ht_emd}, and Figure \ref{fig:Duffing_equation_ht}. One can extract the IMFs from the structure using the code below.

\vspace{-1.00em}
\begin{Code}
emd_duffing = EMD(time=t, time_series=x)
emd_duff, emd_ht_duff, emd_if_duff = \ 
    emd_duffing.empirical_mode_decomposition()[:3]
\end{Code}

In Figure \ref{fig:Duffing_equation_imfs} one can see how accurately the second IMF structure extracted captures the forcing function in the Duffing equation. For comparison, IMF 1 and IMF 2 extracted using \textbf{PyEMD 0.2.10} and \textbf{emd 0.3.3} are also plotted. The various methods are also quantitatively compared by measuring the sum of absolute differences between the second IMFs of each method against the true underlying driving function, the errors are tabulated below:

\vspace{-1.00em}
\begin{equation} \label{eq: imf_error_compare}
\begin{split}
    \sum_{i=0}^{1500}\bigg|IMF^{emd}_{2}(t_i) - 0.1\text{cos}(0.08\pi(t_i))\bigg| &= 21.434\\
    \sum_{i=0}^{1500}\bigg|IMF^{PyEMD}_{2}(t_i) - 0.1\text{cos}(0.08\pi(t_i))\bigg| &= 20.774\\
    \sum_{i=0}^{1500}\bigg|IMF^{AdvEMDpy}_{2}(t_i) - 0.1\text{cos}(0.08\pi(t_i))\bigg| &= 15.978,
\end{split}
\end{equation}
with $ IMF^k_j(t_i) $ being IMF $ j $ isolated using method $ k $ at time point $ t_i $. From Equation \eqref{eq: imf_error_compare} it can be seen that \textbf{AdvEMDpy} most closely captures the underlying driving function.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/Duffing_equation_imfs.png}
	\caption{Figure demonstrating the first three IMFs extracted from Duffing equation displacement along with driving force equation displayed alongside IMF 2.}
	\label{fig:Duffing_equation_imfs}
\end{figure}

The Hilbert spectrum in Figure \ref{fig:Duffing_equation_ht} captures the complex frequency structures present more accurately than \textbf{PyEMD 0.2.10} or \textbf{emd 0.3.3}. The Hilbert spectrum can be calculated using the following code.

\newpage

\begin{Code}
x_hs, y, z = hilbert_spectrum(t, emd_duff, emd_ht_duff, emd_if_duff, 
                              max_frequency=1.3, plot=False)
y = y / (2 * np.pi)  # convert from angular frequency
ax.pcolormesh(x_hs, y, np.abs(z), cmap=`gist_rainbow', 
              vmin=0, vmax=np.abs(z).max())
\end{Code}

For completeness, the Hilbert spectrum of the other two methods can be obtained by using the following code. One can note that in addition to the increased accuracy of \textbf{AdvEMDpy} in resolving the underlying driving function, the resolution of the driving function in the Hilbert spectrum of Figure \ref{fig:Duffing_equation_ht} is not nearly as obscured by the IF of the first IMF as in Figure \ref{fig:Duffing_equation_ht_pyemd} and Figure \ref{fig:Duffing_equation_ht_emd}.

\begin{Code}
pyemd = pyemd0215()
py_emd = pyemd(x)
IP, IF, IA = emd040.spectra.frequency_transform(py_emd.T, 10, `hilbert')
freq_edges, freq_bins = emd040.spectra.define_hist_bins(0, 0.2, 100)
hht = emd040.spectra.hilberthuang(IF, IA, freq_edges)
hht = gaussian_filter(hht, sigma=1)
plt.pcolormesh(t, freq_bins, hht, cmap=`gist_rainbow', 
              vmin=0, vmax=np.max(np.max(np.abs(hht))))
\end{Code}

\begin{Code}
emd_sift = emd040.sift.sift(x)
IP, IF, IA = emd040.spectra.frequency_transform(emd_sift, 10, `hilbert')
freq_edges, freq_bins = emd040.spectra.define_hist_bins(0, 0.2, 100)
hht = emd040.spectra.hilberthuang(IF, IA, freq_edges)
hht = gaussian_filter(hht, sigma=1)
plt.pcolormesh(t, freq_bins, hht, cmap=`gist_rainbow', 
              vmin=0, vmax=np.max(np.max(np.abs(hht))))
\end{Code}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/Duffing_equation_ht_pyemd.png}
	\caption{Hilbert spectrum of Duffing equation IMFs using \textbf{PyEMD 0.2.10}.}
	\label{fig:Duffing_equation_ht_pyemd}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/Duffing_equation_ht_emd.png}
	\caption{Hilbert spectrum of Duffing equation IMFs using \textbf{emd 0.3.3}.}
	\label{fig:Duffing_equation_ht_emd}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/Duffing_equation_ht.png}
	\caption{Hilbert spectrum of Duffing equation IMFs using \textbf{AdvEMDpy} with frequency of driving force equation visible despite being of lower intensity.}
	\label{fig:Duffing_equation_ht}
\end{figure}

\FloatBarrier

\subsection{ Supplementary Example Two} \label{main_AAS_Software_Supplement:Real_world}

One would benefit from an example applied to real-world data. The following well-known data set is from \cite{NOAA20} to demonstrate the effectiveness of \textbf{AdvEMDpy} over \textbf{PyEMD 0.2.10} and \textbf{emd 0.3.3}. One can see the transient nature of the data in Figure \ref{fig:CO2_concentration} which hinders analysis using base-level applications of \textbf{PyEMD 0.2.10} and \textbf{emd 0.3.3}. One can replicate Figure \ref{fig:CO2_concentration} using the below code.

\vspace{-0.50em}
\begin{Code}
CO2_data = pd.read_csv(`../Data/co2_mm_mlo.csv', header=51)
plt.plot(CO2_data[`month'], CO2_data[`decimal date'])
\end{Code}

\vspace{-1.00em}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/CO2_concentration.png}
	\caption{Monthly mean atmospheric concentration of Carbon Dioxide in parts per million from March 1958 until September 2020.}
	\label{fig:CO2_concentration}
\end{figure}

One can apply EMD to the time series with the following code:

\vspace{-0.50em}
\begin{Code}
co2_knots = np.linspace(co2_time[0], co2_time[-1], 200)
emd = EMD(time=co2_time, time_series=co2_time_series)
imfs, hts, ifs = \ 
    emd.empirical_mode_decomposition(knots=co2_knots, 
                                     knot_time=co2_time)[:3].
\end{Code}
\vspace{-0.50em}
and plotting the resulting outputs, one will arrive at Figure \ref{fig:CO2_EMD}. The number of knots is increased above the base-level application to accurately capture the annual cycle. The original time series is plotted as well as the smoothed version of the time series. A single IMF has been extracted as well as a residual after the original time series was smoothed.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/CO2_EMD.png}
	\caption{EMD of monthly mean Carbon Dioxide concentration in parts per million from March 1958 until September 2020.}
	\label{fig:CO2_EMD}
\end{figure}

To plot the Hilbert spectrum, as seen in Figure \ref{fig:CO2_Hilbert}, one can implement the following code:

\begin{Code}
x_hs, y, z = hilbert_spectrum(time, imfs, hts, ifs, 
                              max_frequency=10, which_imfs=[1], 
                              plot=False)
y = y / (2 * np.pi)
fig, ax = plt.subplots()
ax.pcolormesh(x_hs, y, np.abs(z), cmap=`gist_rainbow', 
              vmin=0, vmax=np.abs(z).max()).
\end{Code}

The accuracy of \textbf{AdvEMDpy} versus \textbf{PyEMD 0.2.10} and \textbf{emd 0.3.3} in the resolution of the annual structure can be seen in Equation \eqref{eq: error_compare}.

\begin{equation} \label{eq: error_compare}
\begin{split}
    \sum_{i=0}^{750}\bigg|IF^{PyEMD}_{1}(t_i) - 1\bigg| &= 181.076\\
    \sum_{i=0}^{750}\bigg|IF^{emd}_{1}(t_i) - 1\bigg| &= 174.047\\
    \sum_{i=0}^{750}\bigg|IF^{AdvEMDpy}_{1}(t_i) - 1\bigg| &= 37.616,
\end{split}
\end{equation}
with $ IF^{k}_j(t_i) $ being the IF of IMF $ j $ using method $ k $ at time point $ t_i $. The annual cycle is most accurately resolved in Figure \ref{fig:CO2_Hilbert} compared against Figure \ref{fig:CO2_Hilbert_pyemd} and Figure \ref{fig:CO2_Hilbert_emd}. To reproduce Figure \ref{fig:CO2_Hilbert_pyemd} one can implement the code below. 

\begin{Code}
pyemd = pyemd0215()
py_emd = pyemd(signal)
IP, IF, IA = emd040.spectra.frequency_transform(py_emd[:2, :].T, 
                                                12, `hilbert')
freq_edges, freq_bins = emd040.spectra.define_hist_bins(0, 2, 100)
hht = emd040.spectra.hilberthuang(IF, IA, freq_edges)
hht = gaussian_filter(hht, sigma=1)
fig, ax = plt.subplots()
plt.pcolormesh(time, freq_bins, hht, cmap=`gist_rainbow', 
               vmin=0, vmax=np.max(np.max(np.abs(hht))))
\end{Code}

\vspace{-1.50em}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/CO2_Hilbert_pyemd.png}
	\caption{HT of monthly mean Carbon Dioxide concentration in parts per million from March 1958 until September 2020 using \textbf{PyEMD 0.2.10}.}
	\label{fig:CO2_Hilbert_pyemd}
\end{figure}
\vspace{-1.50em}

To reproduce Figure \ref{fig:CO2_Hilbert_emd} one can implement the code below. By comparing Figure \ref{fig:CO2_Hilbert_pyemd} and Figure \ref{fig:CO2_Hilbert_emd} against Figure \ref{fig:CO2_Hilbert} and by comparing the equalities in Equation \eqref{eq: error_compare}, one can note the increased resolution power of \textbf{AdvEMDpy} over \textbf{PyEMD 0.2.10} and \textbf{emd 0.3.3}.

\vspace{-0.50em}
\begin{Code}
emd_sift = emd040.sift.sift(signal)
IP, IF, IA = emd040.spectra.frequency_transform(emd_sift[:, :1], 
                                                12, `hilbert')
freq_edges, freq_bins = emd040.spectra.define_hist_bins(0, 2, 100)
hht = emd040.spectra.hilberthuang(IF, IA, freq_edges)
hht = gaussian_filter(hht, sigma=1)
fig, ax = plt.subplots()
plt.pcolormesh(time, freq_bins, hht, cmap=`gist_rainbow', 
               vmin=0, vmax=np.max(np.max(np.abs(hht))))
\end{Code}
\vspace{-0.50em}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/CO2_Hilbert_emd.png}
	\caption{HT of monthly mean Carbon Dioxide concentration in parts per million from March 1958 until September 2020 using \textbf{emd 0.3.3}.}
	\label{fig:CO2_Hilbert_emd}
\end{figure}

In \cite{van21}, the resolution of EMD is improved even further by augmenting EMD with X11 and forming the new technique known as EMD-X11. X11 is a simpler trend filtering technique that functions as a post-processing technique in this setting that improves the resolution power of EMD. This will be included in a later version of the software along with many other novel contributions to the field.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/CO2_Hilbert.png}
	\caption{HT of monthly mean Carbon Dioxide concentration in parts per million from March 1958 until September 2020 using \textbf{AdvEMDpy}.}
	\label{fig:CO2_Hilbert}
\end{figure}

In order to minimize carbon emissions in keeping with intergovernmental agreements, nations issue a certain number of allowances which allows the company the right to emit a certain amount of carbon dioxide into the atmosphere in the past business year. One such agreement is the EU's emissions agreement where one European Union Allowance (EUA) allows a company to emit one tonne of carbon dioxide in the previous business year. These EUAs are issued in lots of one thousand EUA which forms a very liquid market. These lots are traded as ECF1 and ECF2 on the EU Emissions Trading Strategy (ETS). The price in Euros of each of these instruments from 22 April 2005 until 17 June 2022 can be seen in Figure \ref{fig:ecf1_and_ecf2}.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/ecf1_and_ecf2.png}
	\caption{ECF1 and ECF2 in Euros from 22 April 2005 until 17 June 2022.}
	\label{fig:ecf1_and_ecf2}
\end{figure}

In Figure \ref{fig:if_imf_1_ecf1} and Figure \ref{fig:if_imf_1_ecf2} the Hilbert transform of the first IMF from ECF1 and ECF2 are presented, respectively. In Figure \ref{fig:if_imf_1_ecf1}, the results are somewhat disturbed from the presence of the discontinuity in ECF1 which is mostly removed in the first IMF. In Figure \ref{fig:if_imf_1_ecf2}, the twelve-to-thirteen day cycle is clearly visible.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/if_imf_1_ecf1.png}
	\caption{IF of IMF 1 for ECF1 in Euros from 22 April 2005 until 17 June 2022 with median and best fits.}
	\label{fig:if_imf_1_ecf1}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/if_imf_1_ecf2.png}
	\caption{IF of IMF 1 for ECF2 in Euros from 22 April 2005 until 17 June 2022 with median and best fits.}
	\label{fig:if_imf_1_ecf2}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/if_ecf1_and_ecf2.png}
	\caption{IF of IMF 2 and IMF 3 for ECF1 and ECF2 in Euros from 22 April 2005 until 17 June 2022 with medians.}
	\label{fig:if_ecf1_and_ecf2}
\end{figure}

In Figure \ref{fig:if_imf_1_ecf2}, the IFs of IMF 2 and IMF 3 from both ECF1 and ECF2 are shown - lower frequency structures are iteratively removed. In Figure \ref{fig:if_ecf1_and_ecf2_45}, the annual structure, first seen in \ref{fig:CO2_Hilbert}, is shown relative to the fifth IMF of both. There is a clear causal link between the two processes, and EMD has assisted in establishing an annual cycle is present in the ECF1 and ECF2 data which can be used in Carbon emission instrument modelling.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Images/if_ecf1_and_ecf2_45.png}
	\caption{IF of IMF 4 and IMF 5 for ECF1 and ECF2 in Euros from 22 April 2005 until 17 June 2022 with medians and annual structure.}
	\label{fig:if_ecf1_and_ecf2_45}
\end{figure}

\FloatBarrier

\bibliographystyle{ACM-Reference-Format}
\bibliography{refs.bib}

\endctxt{C. van Jaarsveldt, M. Ames, G. Peters, and M. Chantler, Supplement to: Package \textbf{AdvEMDpy}: Algorithmic Variations of Empirical Mode Decomposition in \textsf{Python}. {\it Annals of Actuarial Science}, 1--47. \href{10.1017/xxxxx}{10.1017/xxxxx}}

\end{document}
