import numpy as np
from scipy.stats import levy_stable, ks_2samp
from argparse import ArgumentParser
import matplotlib.pyplot as plt


def calculate_ndays_returns(returns_1day_list, nday=10):
    """
    Calculate n-day returns based on 1-day returns
    :param returns_1day_list: 1D array of 1-day returns
    :param nday: default=10;  number of days (==n) for returns calculation
    :return: 1D array n-day returns
    """
    returns_1day_list += 1
    iterations_num = len(returns_1day_list) - (nday - 1)
    returns_ndays = np.zeros(iterations_num)
    for idx in range(iterations_num):
        input_data = returns_1day_list[idx:idx + nday]
        returns_ndays[idx] = np.prod(input_data) - 1
    return returns_ndays


def generate_percentile(distribution_params, size, ndays, percentile_val):
    """
    Generate percentile value
    :param distribution_params:  ordered list with stable distributions parameters (alpha, beta, gamma, delta)
    :param size: size of series of 1-day returns
    :param ndays: number of days (==n) for returns calculation
    :param percentile_val: percentile
    :return: calculated percentile
    """
    returns_1day = levy_stable.rvs(*distribution_params, size=size)
    if abs(returns_1day).max() == np.inf:
        return None
    returns_10day = calculate_ndays_returns(returns_1day, nday=ndays)
    percentile = np.percentile(returns_10day, percentile_val)
    return percentile


def show_distribution_info(data, bins=1000, title='Distribution Info'):
    """
    Show histogram plot and accumulative distribution plot
    :param data: 1D array
    :param bins: number of bins
    :param title: main title
    :return: None
    """
    fig = plt.figure(figsize=(15, 8))
    ax = plt.axes()
    plt.ylabel("Proportion for histogram")
    values, base, _ = plt.hist(data, bins=bins, normed=True, alpha=0.5, color="green", label="Histogram")
    ax.legend()
    ax_bis = ax.twinx()
    values = np.append(values, 0)
    ax_bis.plot(base, np.cumsum(values) / np.cumsum(values)[-1], color='darkorange', marker='o', ms=1,
                label="Cumulative Histogram")
    plt.ylabel("Proportion for accumulation value")
    plt.title(title)
    perc = np.percentile(data, 1)
    plt.xlim(perc, 10)
    ax_bis.legend()
    plt.show()


def make_monte_carlo_simulation(distribution_params,
                                datasize,
                                significance_level=0.50,
                                min_dstat=0.01,
                                percentile_val=1,
                                max_iterations=10 ** 6):
    """
    Make Monte-Carlo simulation for generation of distribution of percentile_val% percentile of 10‐days returns.
    1-day returns are generated by stable distribution
    :param distribution_params: ordered list with stable distributions parameters (alpha, beta, gamma, delta)
    :param datasize: size of series of 1-day returns
    :param significance_level: const for Smirnov-Kolmogorov test
    :param min_dstat: dstat valid value for Smirnov-Kolmogorov test
    :param percentile_val: percentile
    :param max_iterations: maximum number of search iterations
    :return: Flag True if Smirnov-Kolmogorov test is passed or False if - not,
            number of performed iterations ,
            pval- probability of not rejecting the hypothesis for Smirnov-Kolmogorov test,
            dstat,
            array with percentile values
    """
    pval, dstat = 0, np.inf
    cnt = 0
    ndays = 10
    perc_accum1, perc_accum2 = [], []
    success_flag = True

    while (pval < significance_level) or (dstat > min_dstat):
        percentile1 = generate_percentile(distribution_params, datasize, ndays, percentile_val)
        percentile2 = generate_percentile(distribution_params, datasize, ndays, percentile_val)

        if percentile1 is None or percentile2 is None: continue

        perc_accum1.append(percentile1)
        perc_accum2.append(percentile2)

        if (cnt + 1) % 10 == 0:
            dstat, pval = ks_2samp(perc_accum1, perc_accum2)
        cnt += 1
        if cnt >= max_iterations:
            success_flag = False
            break
    result_values_list = perc_accum1 + perc_accum2

    return success_flag, cnt, pval, dstat, result_values_list


if __name__ == '__main__':
    parser = ArgumentParser()
    parser.add_argument('-n', '--observations_num', default=750, required=False, type=int,
                        help='number of 1-day returns observations')
    parser.add_argument('-a', '--alpha', default=1.7, required=False, type=float,
                        help='alpha parameter for stable distribution function. Valid value belongs to range (0; 2]')
    parser.add_argument('-b', '--betta', default=0, required=False, type=float,
                        help='betta (skewness parameter) for stable distribution function is a measure of asymmetry.'
                             'Valid value belongs to range [-1; 1]')
    parser.add_argument('-g', '--gamma', default=1., required=False, type=float,
                        help='gamma - shift  parameter for stable distribution function')
    parser.add_argument('-d', '--delta', default=1., required=False, type=float,
                        help='delta - scale parameter for stable distribution function')
    parser.add_argument('-p', '--percentile', default=1., required=False, type=float,
                        help='percentile value by which the distribution is generated')
    parser.add_argument('-sl', '--significance_level', default=0.5, required=False, type=float,
                        help='probability of not rejecting the hypothesis for Smirnov-Kolmogorov test')
    parser.add_argument('-dst', '--min_dstat_value', default=0.01, required=False, type=float,
                        help='dstat valid value for Smirnov-Kolmogorov test')
    parser.add_argument('-mi', '--max_iter', default=10 ** 6, required=False, type=int,
                        help='maximum number of search iterations')
    args = parser.parse_args()

    distribution_params = args.alpha, args.betta, args.gamma, args.delta

    success_flag, cnt, pval, dstat, result_values_list = make_monte_carlo_simulation(distribution_params,
                                                                                     datasize=args.observations_num,
                                                                                     significance_level=args.significance_level,
                                                                                     min_dstat=args.min_dstat_value,
                                                                                     percentile_val=args.percentile,
                                                                                     max_iterations=args.max_iter)
    print(f'Smirnov-Kolmogorov test was passed: {success_flag}\n'
          f'Number of iterations: {cnt}\n'
          f'dstat = {dstat}, pvalue = {pval}')

    show_distribution_info(result_values_list)


