Monte_Carlo.py

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


def compute_ndays_returns(returns_1day_list, nday=10):
    """
    Calculate n-day returns based on 1-day returns.

    Args:
        returns_1day_list: 1D array of 1-day returns
        nday: number of days for returns calculation (default 10)

    Returns:
        returns_ndays: 1D array n-day returns
    """
    returns_1day_list += 1
    iterations_num = len(returns_1day_list) - (nday - 1)
    returns_ndays = np.zeros(iterations_num)
    for i in range(iterations_num):
        input_data = returns_1day_list[i:i + nday]
        returns_ndays[i] = np.prod(input_data) - 1
    return returns_ndays


def create_percentile(distribution_params, size, ndays, percentile_val):
    """
    Generate percentile value.

    Args:
        distribution_params: list of stable distributions parameters (alpha, beta, gamma, delta)
        size: size of series of 1-day returns
        ndays: number of days for returns calculation
        percentile_val: percentile

    Returns:
        percentile: calculated percentile
    """
    returns_1day = levy_stable.rvs(*distribution_params, size=size)
    returns_10day = compute_ndays_returns(returns_1day, nday=ndays)
    percentile = np.percentile(returns_10day, percentile_val)
    return percentile


def show_distribution(data, bins=1000):
    """
    Show histogram plot and Emperical distribution function plot.

    Args:
        data: 1D array of percentile values
        bins: number of bins

    Returns:
        None
    """
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.ylabel("Count")
    plt.xlabel("Percentile")
    values, base, _ = plt.hist(data, bins=bins, color="blue", label="Histogram of procentiles")
    ax2 = ax.twinx()
    values = np.append(values, 0)
    ax2.plot(base, np.cumsum(values) / np.cumsum(values)[-1], color="red", label="Emperical distribution function")
    perc = np.percentile(data, 1)
    plt.xlim(perc, 10)
    fig.legend(loc="upper left", bbox_to_anchor=(0, 1), bbox_transform=ax.transAxes)
    plt.show()


def do_monte_carlo_simulation(distribution_params,
                              datasize,
                              significance_level,
                              min_err,
                              percentile_val,
                              max_iterations):
    """
    Do Monte-Carlo simulation for generation of distribution of percentiles of 10‐days returns.
    1-day returns are generated by stable distribution.

    Args:
        distribution_params: list of stable distributions parameters (alpha, beta, gamma, delta)
        datasize: size of series of 1-day returns
        significance_level: const for Smirnov-Kolmogorov test
        min_err: err valid value for Smirnov-Kolmogorov test
        percentile_val: percentile
        max_iterations: maximum number of search iterations

    Returns:
        success_indicator: indicator True if Smirnov-Kolmogorov test is passed or False if - not,
        cnt: number of performed iterations
        pval: probability of not rejecting the hypothesis for Smirnov-Kolmogorov test
        err: err valid value for Smirnov-Kolmogorov test
        result: 1D array of percentile values
    """
    pval, err = 0, np.inf
    cnt = 0
    ndays = 10
    perc_list1, perc_list2 = [], []
    success_indicator = True

    while (pval < significance_level) or (err > min_err):
        perc1 = create_percentile(distribution_params, datasize, ndays, percentile_val)
        perc2 = create_percentile(distribution_params, datasize, ndays, percentile_val)
        perc_list1.append(perc1)
        perc_list2.append(perc2)

        if (cnt + 1) % 10 == 0:
            err, pval = ks_2samp(perc_list1, perc_list2)
        cnt += 1
        if cnt >= max_iterations:
            success_indicator = False
            break
    result = perc_list1 + perc_list2

    return success_indicator, cnt, pval, err, result


distribution_params = 1.7, 0.0, 1.0, 1.0

success_indicator, cnt, pval, err, result = do_monte_carlo_simulation(distribution_params,
                                                                      datasize=750,
                                                                      significance_level=0.5,
                                                                      min_err=0.02,
                                                                      percentile_val=1.0,
                                                                      max_iterations=10 ** 6)
print(f'Smirnov-Kolmogorov test: {success_indicator}\n'
      f'Number of iterations: {cnt}\n'
      f'error = {err}, pvalue = {pval}')

show_distribution(result)