The overall procedure for NA-MEMD Causal Decomposition is illustrated in Figure 1 by 5 steps. It first adds multichannel auxiliary white noise (shown in step a-2) to original bivariate time series signals (shown in step a-1), followed by NA-MEMD to obtain two corresponding IMFs sets (shown in steps b-1 and b-2) as well as their phase coherence. Next, the main ICC (framed by the purple rectangles in steps c and d) and ICCs sets (framed by the yellow rectangles in steps c and d) are selected according to average frequency and phase difference. Only the IMFs with identical frequency scales and small phase differences can be regarded as ICCs. The main ICC is usually selected as the first ICC in the set. Removing the main ICC from IMFs sets and adding the remaining signals in each IMFs setup, two pairs of signals are formed by one modified signal without the main ICC and the other original one (shown in step d). The phase coherence of the two pairs is calculated. Combined with the value of the phase coherences of two redecomposed signal pairs and the original signal pair, the value of Absolute Causal Strengths (ACSs) and Relative Causal Strengths (RCSs) can be obtained by NA-MEMD Causal Decomposition (shown in step e).

The procedure for function na_memd is illustrated in Figure 2. The function enables to form a matrix, giving Intrinsic Mode Functions (IMFs) decomposed from NA-MEMD with causal-effect time series defined by input matrix input_data. Relevant variables about function na_memd are listed in Table 1. The function na_memd is invoked when variables input_data, ave_noise, level_noise, noise_channel_num, and en_num are set. Then, function size returns the number of row elements of matrix input_data to variable len_inp, and the number of column elements of matrix input_data to variable wid_inp. level_noise and noise_channel_num are used to generate the random Gaussian-White noise time series. After then, function memd is run if vector/matrix inp_noise is appended to matrix input_data, which is named matrix input_cha_noi. It repeats for constant en_num times (i is the loop variable). For each repetition, a three-dimension array imfs is returned which contains all IMFs decomposed from matrix input_cha_noi. In order to facilitate the ensemble process of NA-MEMD, array imfs is reshaped into cell imf_result. Considering repetitions in-between, if the matrix sizes in cell imf_result perhaps are inconsistent, the residual is embedded with a pseudo-monotone variation. In order to guarantee the consistency of the size in cell imf_result, the residuals can be ignored together. Matrix elements of imf_result, sum_imf_result, and imf_result are first initialized to be zero, which stand for reshaped IMF data, the summation of the reshaped IMF, and the average value of the summation respectively. The noise data are generated by random alternating noise components determined by parameter level_noise. This process is repeated by times of en_num. For each loop, the noise is appended to matrix input_data which is matrix input_cha_noi. Then matrix imfs is calculated by performing MEMD on input_cha_noi with function memd. Next, imfs is reshaped into imf_result for further summation which data is stored in sum_imf_result. The average value is also calculated and assigned to ave_imf_result.

The workflow for function PLseires is shown in Figure 3. The function outputs a matrix that shows the average frequency and their phase difference of IMFs decomposed from NA-MEMD by matrices imf1 and imf2 (Table 2), respectively. The row of imf1 and imf2 stand for the index of IMFs and the column of imf1 and imf2 refers to the length of the time series. The output matrix is defined as matrix peakMatrix. The function PLseries is called when matrix imf1 and imf2 are input. The elements in vectors av_fre1, av_fre2, and difference are initialized to be zero and their length of the row is also set to be consistent with that of imf1 and imf2 (their rows are in the same length). After that, for each IMF, Hilbert Transform is used to calculate its instantaneous phase. To alleviate the effect of phase winding, the function unwarp is called to obtain the real phase which is defined as vectors un_ang1 and un_ang2. Then the average value of the unwrapped phase is calculated and stored in vectors av_ang1 and av_ang2, respectively. Their difference is stored in vector difference. According to vectors un_ang1 and un_ang2, the phase information is converted to frequency information which is stored in vectors fre1 and fre2, respectively. Their average values are calculated and are stored in vectors av_fre1 and av_fre2. The final output matrix peakMatrix is the combination of vectors av_fre1, av_fre2, and difference which are in the first, the second, and the third column, respectively.

The workflow for function causal_decomposition is demonstrated in Figure 4. The function generates a matrix causal_matrix to provide values of RCSs and ACSs according to the times series inputted by matrix input referred to in Table 3. The raw data of bivariate time series are row-by-row loaded to input which is subsequently defined as vectors s1 and s2. Those two decomposed IMF sets are then assigned as matrices imfs1 and imfs2, respectively. Intrinsic Causal Components (ICCs) sets, vector ICC, are manually reviewed and selected by function plot_and_pick. Based on ICC, the number ICC_main is confirmed to indicate RCSs and ACSs of raw bivariate time series of input. Main ICC is selected manually by function plot_and_pick, and is defined as vector imfss1 and imfss2. Their phase coherence and variance are calculated by function phasefcimf and function nvar, respectively. After obtaining the main ICC, it is removed from the corresponding inputs s1 and s2, and the rest are defined as vectors s1r and s2r. They are used to replace the corresponding signal in input to form two new input matrices input1 and input2. Then they are re-decomposed by calling function na_memd to calculate two IMFs pairs (matrices imfsr11, imfsr12 and matrices imfsr21, imfsr22, respectively). The phase coherence between paired IMFs sets and the distance between the phase coherences of the original IMFs pairs and the redecomposed ones are calculated and stored in (ps12, ps21) and (p12, p21), respectively. The final causal strengths are outputted by the matrix causal_matrix.

To guarantee the experimental efficiency in the real application, variable noise_channel_num and en_num and the data length should be modified. The input data is the Gaussian-White noise with a mean value of 0 and variance of 1, and the experiment is performed on Matlab 2020b on the fully powered laptop.

Since the random noise is involved in NA-MEMD, it is likely that the results of different execution will be slightly different from one another. Therefore, the consistency of the output should be tested. In this test, level_noise, noise_channel_num, en_num, and data length are set to be 0.001, 3, 5, and 61, respectively, and the open-access predator-prey data is provided by Vucetich and Peterson (2012). The result is shown in Table 7.

In Table 7, since the X-to-Y causality value and Y-to-X causality value sum to one, only the X-to-Y causal value is considered in the stability discussion.

From the statistical point of view, the mean value of C_xy is 0.6384 and the variance is 0.004712 which is less volatile. Combined with the actual causal process demonstrated by the causality values by Yang et al. (2018), all of them indicate strong and valid causality since our results are close to that of Yang et al. (2018), which demonstrates the outcome of the procedure is consistent with the algorithm.

Generally, input parameters level_noise, en_num, and noise_channel_num are used to determine the total amount of noise appended to the original signal. The description of the relevant variables can be checked in Table 3. Before setting the parameters, the noise level in the original input data should be evaluated. If the input has already been overwhelmed by noise, both the two parameters are supposed to be set higher to better remove the noise in the original data. However, with the Gaussian-White noise increasing, the effective data in the input is also likely to be attenuated by the appended noise. As a result, it is of significance to choose a proper value of level_noise and noise_channel_num to avoid excessive Gaussian-White noise. One method of removing the original noise without doing damage to the effective data is to reduce the product of level_noise and noise_channel_num and to increase en_num.

For ICCs, the ICCs set and the main ICC need to be selected for further operation. Among all the IMFs, ICCs are the ones whose average frequencies and phases are on the same scale. IMFs in different frequency scales or with large phase differences should be excluded from the ICCs set. In ICCs set, the main ICC is the one with the smallest average phase difference. Besides, the phase diagram of IMFs generated by Matlab can also be used to discriminate ICCs from all IMFs. Assist from both data in peakMatrix and visual aid of signal diagram is necessary for ICCs discrimination.

In this code, it is recommended that only RCSs are used to detect the causal-effect relationship between bivariate time series signals.

Absolute Causal Strengths is defined as the variance weighted Euclidean distance (De Leeuw and Pruzansky, 1978) between the re-decomposed IMFs set and the original IMFs set. Its value can reflect the relative cause-effect relationship between multiple signals. The signals with higher ACS can be regarded as the cause and the ones with smaller ACS can be considered as the effect, respectively. Specifically, if only two signals are considered, RCS can also be calculated to illustrate their cause-effect relationship. It normalizes the ACS into an interval between zero and one, and we use the value of RCS to interfere with the causal-effect relationship. When the value of RCS is larger than 0.5, it represents the relationship of cause, and when the value of RCS is smaller than 0.5, it represents the relationship of effect. When the value is equal to 0.5, the relationship can be either reciprocal causation or an irrelevant relationship. When the ACS of two candidate signals is negligible, RCS will converge to 0.5. The threshold value is normally set as 0.05 (Yang et al., 2018).