<jupyter_start><jupyter_text>EF3M Implementation Testing AbstractThis notebook tests the implementation of the EF3M algorithm using synthetic data as well as the example used in the source literature. This notebook is intended to provide convincing evidence of the accuracy of this EF3M implementation. IntroductionThe Exact Fit of the first 3 Moments (EF3M) algorithm allows the parameters of a mixture of Gaussian distributions to be estimated given the first 5 moments of the mixture distribution, as well as the assumption that the mixture distribution is composed of a number of Gaussian distributions. The algorithm, its development, and its original implementation are described in <NAME>, M. and <NAME> (2014): "A mixture of Gaussians approach to mathematical portfolio oversight: The EF3M algorithm." _Quantitative Finance_, Vol. 14, No. 5, pp. 913-930. The implementation tested here can be found in the `EF3M\ef3m.py` module located in this directory of the repository. Three test cases are presented. First, the example from the algorithm's source paper is presented to confirm that we can replicate the results of literature. Second, a user-defined mixture of 2 Gaussian distributions is generated and the EF3M implementation is tasked to recover the original mixture parameters from the raw moments of the synthetic distribution. Third, the second case is executed using a series of randomly chosen parameters from which to compose the mixture distribution. The notebook user can alter this at will but should keep in mind that running a large number of trials with hundreds or thousands of rounds each will take hours or days even using the current multiprocessing implementation. ConclusionThe tests run in this notebook conclude that using the EF3M algorithm to generate a series of estimates followed by a kernel density estimation to determine the most likely value for each parameter from said series results in a recovery of the original mixture parameters to within a very close tolerance. There appear to be at least two central solutions when fitting a mixture of two Gaussians, but both solutions appear to reproduce the mixture distribution, thus both can be considered viable and correct solutions in the absence of further evidence.  Next StepsThe implementation tested here supports only distributions assumed to be composed of two Gaussian distributions. A natural next step would be to extend this algorithm to $n$ Gaussian distributions. This implementation would take exponentially more computing power but may prove useful in more complex mixtures. Additionally, extending or encapsulating these methods to make use of a distributed computing resource would greatly speed up simulations that consist of more than 1000 fits. -------- Test CasesTest cases and accompanying code are described below. The next cell sets up some helper functions for the below analysis.<jupyter_code># imports
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde, moment
from IPython.display import display

# import function from our custom module
from mlfinlab.bet_sizing.ef3m import M2N, raw_moment, most_likely_parameters


def mixture_samples(params, n=10_000):
    # returns 'n' samples from a mixture distribution defined by 'params'
    # :param params: (list) mixture parameters [mu1, mu2, sigma1, sigma2, p1]
    # :param n: (int) number of samples to draw
    # :return: (numpy.array) an array of the length 'n' of distribution mixture samples
    # ===================================================
    n1 = int(n*params[4])
    n2 = n - n1
    d1 = np.random.normal(params[0], params[2], n1)
    d2 = np.random.normal(params[1], params[3], n2)
    return np.concatenate((d1, d2))

def plot_fit(true_list, fit_list, n=10_000):
    # plots the true and fit mixture distributions
    # :param true_list: (list) true parameters [mu1, mu2, sigma1, sigma2, p1]
    # :param fit_list: (list) true parameters [mu1, mu2, sigma1, sigma2, p1]
    # :param n: (int) number of samples to draw
    # :return: (None)
    # ===================================================
    d_true = mixture_samples(true_list, n)
    d_fit = mixture_samples(fit_list, n)
    fig, ax = plt.subplots(figsize=(10, 8))
    sns.distplot(d_true, kde=True, bins=200, ax=ax, label=f'True Mixture Distribution {true_list}')
    sns.distplot(d_fit, kde=True, bins=200, ax=ax, label=f'Fitted Mixture Distribution {fit_list}')
    ax.legend(loc='upper left')
    plt.show()
    return None<jupyter_output><empty_output><jupyter_text>---- Case 1: Example data from source literature<NAME>, M. and <NAME> (2014): "A mixture of Gaussians approach to mathematical portfolio oversight: The EF3M algorithm." _Quantitative Finance_, Vol. 14, No. 5, pp. 913-930. The below cell calculates the raw moments from the sample distribution parameters, $(-2, 1, 2, 1, \frac{1}{10})$, as a demonstration of the helper functions.<jupyter_code># moments and parameters from source literature
true_params = [-2, 1, 2, 1, 0.1]

true_samples = mixture_samples(true_params, n=1_000_000)

# calculate the central moments
mmnts = [moment(true_samples, moment=i) for i in range(1, 6)]

# calculate the raw moments
raw_mmnts = raw_moment(central_moments=mmnts, dist_mean=np.mean(true_samples))

# moments caluclated from the paper
moments = [0.7, 2.6, 0.4, 25, -59.8]

# let's look at the moments
display(pd.DataFrame.from_dict({'Moment': [f'E[r^{i+1}]' for i, m in enumerate(mmnts)],
                                'From paper': moments,
                                'Central (about mean)': mmnts,
                                'Raw (about origin)': raw_mmnts,
                                'Significance': ['Mean', 'Variance', 'Skewness', 'Kurtosis', 'Hyperskewness']}))<jupyter_output><empty_output><jupyter_text>The cell below executes a number of fits and collects the results into a table of the most likely parameters (using a kernel density estimation), and the average and standard deviation of the difference between the true value and fit values from each fitting run. The cell after this plots the distributions generated from the true and most likely parameters as an ultimate comparison of the fit. It is especially useful in instances where the most likely set of fitted parameters is one of the multiple solutions that does not exactly match the true parameters.<jupyter_code># 20_000 runs takes approx. 4h30min (Intel i7-7700K)

n_runs = 100
epsilon = 10**-6
factor = 5

m2n = M2N(moments, epsilon=epsilon, factor=factor, n_runs=n_runs, variant=2, max_iter=10_000_000)
df2 = m2n.mp_fit()
df2 = df2.sort_values('error')
results = most_likely_parameters(df2)

# setup comparison DataFrame
para_cols = ['mu_1', 'mu_2', 'sigma_1', 'sigma_2', 'p_1']

for ei, pcol in enumerate(para_cols):
    df2[pcol+'_err'] = true_params[ei] - df2[pcol]

df_results_case1 = pd.DataFrame.from_dict({'Parameter': ['mu1', 'mu2', 'sigma1', 'sigma2', 'p1'],
                                           'Original': true_params,
                                           'Most Likely': list(results.values()),
                                           'Deviation Average': [df2[col].mean() for col in [p+'_err' for p in para_cols]],
                                           'Deviation StDev': [df2[col].std() for col in [p+'_err' for p in para_cols]]})

display(df_results_case1)
# plot the distributions generated from the true and most likely parameters
plot_fit(true_params, list(results.values()), n=100_000)<jupyter_output><empty_output><jupyter_text>---- Case 2: User-defined single exampleDefine the parameters of the two Gaussian distributions to mix in the cell below.<jupyter_code># setup user sample
mu1 = -1
mu2 = 4
sigma1 = 2
sigma2 = 1.5
p1 = 0.45<jupyter_output><empty_output><jupyter_text>Below we generate the mixture, and then calculate the central moments (using `scipy.stats.moment`) and finally the raw momments (about the origin) as required by the EF3M algorithm.<jupyter_code># generate 10_000 samples for our custom mixture distribution
true_params = [mu1, mu2, sigma1, sigma2, p1]
true_samples = mixture_samples(true_params, n=100_000)

# calculate the central moments
mmnts = [moment(true_samples, moment=i) for i in range(1, 6)]

# calculate the raw moments
raw_mmnts = raw_moment(central_moments=mmnts, dist_mean=np.mean(true_samples))

# let's look at the moments
display(pd.DataFrame.from_dict({'Moment':[f'E[r^{i+1}]' for i, m in enumerate(mmnts)],
                                'Central (about mean)': mmnts,
                                'Raw (about origin)': raw_mmnts,
                                'Significance': ['Mean', 'Variance', 'Skewness', 'Kurtosis', 'Hyperskewness']}))<jupyter_output><empty_output><jupyter_text>Now create the M2N object and run a fit. Take a look at the top values, find the most likely values, and plot to compare the original distribution against the distribution made using the recovered parameters. Try adjusting the `n_runs`, or `epsilon` in `m2n.mpFit` to get better results. Variant 2 usually converges faster than variant 1 and typically leads to a lower error, but both are available for use.<jupyter_code>n_runs = 50
epsiolon = 10**-5
factor = 5

m2n = M2N(raw_mmnts, epsilon=epsilon, factor=factor, n_runs=n_runs, variant=2, max_iter=10_000_000)
df_results = m2n.mp_fit()
display("Top values from the EF3M fit:")
display(df_results.head())
results = most_likely_parameters(df_results)
display(pd.DataFrame.from_dict({'Original': true_params,
                                'Estimated': list(results.values())}))
plot_fit(true_params, list(results.values()), n=100_000)  # adjust 'n' to get a smoother distribution<jupyter_output>|#########################| Completed 50 of 50 fitting rounds.<jupyter_text>---- Case 3: Series of randomly generated examplesIn <NAME> and Foreman (2014), the authors used randomly generated Gaussian mixture distributions to test that the EF3M algorithm performed well in general. The parameter tuple desribing the observed distribution, $(\tilde{\mu}_1, \tilde{\mu}_2, \tilde{\sigma}_1, \tilde{\sigma}_2, \tilde{p}_1)$, was drawn from the uniform distribution with boundaries $-1 \leq \mu_1 \leq 0 \leq \mu_2 \leq 1$, $0 < \sigma_1 \leq 1$, $0 < \sigma_2 \leq 1$, $0 < p < 1$. This is reflected in the code below. Change the number of rounds, number of fits per round, and mixture parameters to test additional scenarios. Authors used 10,000 test rounds and 10,000 fits per round, but note that running this on a single developer workstation will take days or weeks to complete, so smaller numbers are advised.<jupyter_code>test_rounds = 5  # number of random tests to conduct
fit_runs = 10  # M2N fitting rouns per round
epsiolon = 10**-5
factor = 5

# define parameter ranges, (low, high)
mu1_range = (-1, 0)
mu2_range = (0, 1)
sigma1_range = (1e-5, 1)
sigma2_range = (1e-5, 1)
p1_range = (1e-5, 1-1e-5)

param_ranges = [mu1_range, mu2_range, sigma1_range, sigma2_range, p1_range]
true_collection, fit_collection = [], []
mu1_err, mu2_err, sigma1_err, sigma2_err, p1_err = [], [], [], [], []
m1_err, m2_err, m3_err, m4_err, m5_err = [], [], [], [], []

for round_i in range(test_rounds):
    true_params = [np.random.uniform(r[0], r[1]) for r in param_ranges]  # generate random parameters
    true_samples = mixture_samples(true_params, n=10_000)
    mmnts = [moment(true_samples, moment=i) for i in range(1, 6)]  # calculate the central moments
    raw_mmnts = raw_moment(central_moments=mmnts, dist_mean=np.mean(true_samples))  # calculate the raw moments
    m2n = M2N(raw_mmnts, epsilon=epsilon, factor=factor, n_runs=fit_runs, variant=2, max_iter=10_000_000)
    df_results = m2n.mp_fit()
    results = list(most_likely_parameters(df_results).values())
    fit_mmnts = m2n.get_moments(results, return_result=True)
    
    # store results
    true_collection.append(true_params)
    fit_collection.append(results)
    
    # store parameter deviations
    for ai, a_err in enumerate([mu1_err, mu2_err, sigma1_err, sigma2_err, p1_err]):
        a_err.append(true_params[ai] - results[ai])
    
    # store moment deviation
    for bi, m_err in enumerate([m1_err, m2_err, m3_err, m4_err, m5_err]):
        m_err.append(mmnts[bi] - fit_mmnts[bi])
    

df_collection = pd.DataFrame.from_dict({'Original': true_collection,
                                        'Fit': fit_collection,
                                        'mu1_err': mu1_err, 'mu2_err': mu2_err, 'sigma1_err': sigma1_err, 'sigma2_err': sigma2_err, 'p1_err': p1_err,
                                        'm1_err': m1_err, 'm2_err': m2_err, 'm3_err': m3_err, 'm4_err': m4_err, 'm5_err': m5_err})

err_cols = ['mu1_err', 'mu2_err', 'sigma1_err', 'sigma2_err', 'p1_err']
m_cols = ['m1_err', 'm2_err', 'm3_err', 'm4_err', 'm5_err']

df_results_case2 = pd.DataFrame.from_dict({'Errors': [1, 2, 3, 4, 5],
                                           'Errors Average': [df_collection[col].mean() for col in m_cols],
                                           'Errors StDev': [df_collection[col].std() for col in m_cols],
                                           'Deviation': ['mu1', 'mu2', 'sigma1', 'sigma2', 'p1'],
                                           'Deviation Average': [df_collection[col].mean() for col in err_cols],
                                           'Deviation StDev': [df_collection[col].std() for col in err_cols]})

print("Compare results to Figure 4 in López de Prado and Foreman (2014)")
display(df_results_case2)<jupyter_output>|####################     | Completed 8 of 10 fitting rounds..
{'fraction_non_alphanumeric': 0.10466906105695228, 'fraction_numerical': 0.030198636663490434, 'mean_word_length': 3.8693790149892933, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '149971', 'license': 'MIT', 'n_tokens_mistral': 4254, 'n_tokens_neox': 3972, 'n_words': 1512, 'path': 'Chapter10/ef3m_testing.ipynb', 'repo_name': 'jorses/research', 'size': '13675'}