Expectation of Maximum and Minimum of Partial Sums of Normal Random Variables

arkon.md

original question

https://stats.stackexchange.com/questions/526581/expectation-of-maximum-and-minimum-of-partial-sums-of-normal-random-variables

distro_normal.py

"""
special cases for normal distribution

https://en.wikipedia.org/wiki/Mills_ratio
https://stats.stackexchange.com/questions/308224/expectation-of-x-given-x-c
https://math.stackexchange.com/questions/1367607/expectation-of-minx-c-for-x-truncated-r-v-and-c-constant
https://stats.stackexchange.com/questions/526581/expectation-of-maximum-and-minimum-of-partial-sums-of-normal-random-variables
https://stats.stackexchange.com/questions/166273/expected-value-of-x-in-a-normal-distribution-given-that-it-is-below-a-certain-v

"""

import math

import numpy as np

from scipy.stats import norm as sci_norm


class DistroNormalALG:
    """
    special cases for normal distribution
    """

    @classmethod
    def expect_vector_X_with_X_GT_C(cls, mean_X:np.array, sigma_X:np.array, level_C:np.array) -> np.array:
        "expectation of X, given X > C"
        """
        E[ X | X > C] = %mu + %sigma pdf(z) / (1 - cdf(z))
        z = (C - %mu) / %sigma
        """
        z = (level_C - mean_X) / sigma_X
        pdf = sci_norm.pdf(z)
        cdf = sci_norm.cdf(z)
        E_with_X_GT = mean_X + sigma_X * (pdf) / (1.0 - cdf)
        return E_with_X_GT

    @classmethod
    def expect_vector_X_with_X_LT_C(cls, mean_X:np.array, sigma_X:np.array, level_C:np.array) -> np.array:
        "expectation of X, given X < C"
        """
        E[ X | X < C] = %mu - %sigma pdf(z) / (cdf(z))
        z = (C - %mu) / %sigma
        """
        z = (level_C - mean_X) / sigma_X
        pdf = sci_norm.pdf(z)
        cdf = sci_norm.cdf(z)
        E_with_X_LT = mean_X - sigma_X * (pdf) / (cdf)
        return E_with_X_LT

    @classmethod
    def expect_vector_min_X_C(cls, mean_X:np.array, sigma_X:np.array, level_C:np.array) -> np.array:
        "expectation of min(X,c)"
        """
        E[ min(X,C) ] = E[ X | X < C ] * Pr( X < C ) + E[ C | X > C] * Pr( X > C )
        """
        z = (level_C - mean_X) / sigma_X
        cdf = sci_norm.cdf(z)
        E_lo = cls.expect_vector_X_with_X_LT_C(mean_X, sigma_X, level_C)
        E_hi = level_C
        E_min_X_C = E_lo * (cdf) + E_hi * (1.0 - cdf)
        return E_min_X_C

    @classmethod
    def expect_vector_max_X_C(cls, mean_X:np.array, sigma_X:np.array, level_C:np.array) -> np.array:
        "expectation of max(X,c)"
        """
        E[ max(X,C) ] = E[ C | X < C ] * Pr( X < C ) + E[ X | X > C] * Pr( X > C )
        """
        z = (level_C - mean_X) / sigma_X
        cdf = sci_norm.cdf(z)
        E_lo = level_C
        E_hi = cls.expect_vector_X_with_X_GT_C(mean_X, sigma_X, level_C)
        E_max_X_C = E_lo * (cdf) + E_hi * (1.0 - cdf)
        return E_max_X_C

    @classmethod
    def expect_provided_min_sum_X(cls, mean_X:float, sigma_X:float, size_sum:int) -> float:
        "expectation of min(S_1,S_2,...,S_N), S_K=sum(X_1,X_2...,X_K), K=1...N, N=size_sum"
        """
        E[ min(S_1, S_2, ..., S_N) ]
        =
        E[ S_N ]
        -
        sum from { K=1 } to { N-1 } { 1 over K } E[ max(0,+S_K) ]
        """
        K = np.arange(1, size_sum)  # 1 ... N-1

        mean_SK = K * mean_X
        sigma_SK = np.sqrt(K) * sigma_X
        E_max_0_pos_SK = cls.expect_vector_max_X_C(+mean_SK, sigma_SK, 0.0)

        E_SN = size_sum * mean_X
        inv_K = 1.0 / K
        E_min_sum = E_SN - np.sum(inv_K * E_max_0_pos_SK)

        return E_min_sum

    @classmethod
    def expect_provided_max_sum_X(cls, mean_X:float, sigma_X:float, size_sum:int) -> float:
        "expectation of max(S_1,S_2,...,S_N), S_K=sum(X_1,X_2...,X_K), K=1...N, N=size_sum"
        """
        E[ max(S_1, S_2, ..., S_N) ]
        =
        E[ S_N ]
        +
        sum from { K=1 } to { N-1 } { 1 over K } E[ max(0,-S_K) ]
        """
        K = np.arange(1, size_sum)  # 1 ... N-1

        mean_SK = K * mean_X
        sigma_SK = np.sqrt(K) * sigma_X
        E_max_0_neg_SK = cls.expect_vector_max_X_C(-mean_SK, sigma_SK, 0.0)

        E_SN = size_sum * mean_X
        inv_K = 1.0 / K
        E_max_sum = E_SN + np.sum(inv_K * E_max_0_neg_SK)

        return E_max_sum

    @classmethod
    def expect_standard_min_sum_X(cls, size_sum:int) -> float:
        "expectation of min(S_1,S_2,...,S_N), S_K=sum(X_1,X_2...,X_K), K=1...N, N=size_sum"
        "when mean=0 and variance=1"
        """
        E[ min(S_1, S_2, ..., S_N) ]
        =
        -
        {1 over sqrt{2 %pi}} sum from { K=1 } to { N-1 } { 1 over sqrt(K) }
        """
        K = np.arange(1, size_sum)  # 1 ... N-1

        inv_sqrt_K = 1.0 / np.sqrt(K)
        factor = 1.0 / math.sqrt(2.0 * math.pi)
        E_min_sum = -factor * np.sum(inv_sqrt_K)

        return E_min_sum

    @classmethod
    def expect_standard_max_sum_X(cls, size_sum:int) -> float:
        "expectation of max(S_1,S_2,...,S_N), S_K=sum(X_1,X_2...,X_K), K=1...N, N=size_sum"
        "when mean=0 and variance=1"
        """
        E[ max(S_1, S_2, ..., S_N) ]
        =
        +
        {1 over sqrt{2 %pi}} sum from { K=1 } to { N-1 } { 1 over sqrt(K) }
        """
        K = np.arange(1, size_sum)  # 1 ... N-1

        inv_sqrt_K = 1.0 / np.sqrt(K)
        factor = 1.0 / math.sqrt(2.0 * math.pi)
        E_max_sum = +factor * np.sum(inv_sqrt_K)

        return E_max_sum

distro_normal_test.py

"""
"""

import numpy as np

from trader.model.distro_normal import DistroNormalALG


def test_expect_vector_more_less():
    print()

    mean_X = 0.0
    sigma_X = 1.0
    level_C = np.array([-1.0, 0.0, +1.0])

    expect_GT = DistroNormalALG.expect_vector_X_with_X_GT_C(mean_X, sigma_X, level_C)
    print(f"expect_GT={expect_GT}")
    assert np.all((np.isclose(expect_GT, np.array([+0.28759997, +0.79788456, +1.52513528]))))

    expect_LT = DistroNormalALG.expect_vector_X_with_X_LT_C(mean_X, sigma_X, level_C)
    print(f"expect_LT={expect_LT}")
    assert np.all((np.isclose(expect_LT, np.array([-1.52513528, -0.79788456, -0.28759997]))))


def test_expect_vector_min_max():
    print()

    mean_X = 0.0
    sigma_X = 1.0
    level_C = np.array([-1.0, 0.0, +1.0])

    expect_MIN = DistroNormalALG.expect_vector_min_X_C(mean_X, sigma_X, level_C)
    print(f"expect_MIN={expect_MIN}")
    assert np.all((np.isclose(expect_MIN, np.array([-1.08331547, -0.39894228, -0.08331547]))))

    expect_MAX = DistroNormalALG.expect_vector_max_X_C(mean_X, sigma_X, level_C)
    print(f"expect_MAX={expect_MAX}")
    assert np.all((np.isclose(expect_MAX, np.array([+0.08331547, +0.39894228, +1.08331547]))))


def test_expect_sum_min_max():
    print()

    mean_X = 0.0
    sigma_X = 1.0
    size_list = [1, 3, 10, 30, 100]

    for size_sum in size_list:
        print(f"size_sum={size_sum}")

        pro_max_sum = DistroNormalALG.expect_provided_max_sum_X(mean_X, sigma_X, size_sum)
        std_max_sum = DistroNormalALG.expect_standard_max_sum_X(size_sum)
        print(f"pro_max_sum={pro_max_sum:+} std_max_sum={std_max_sum:+}")
        assert np.isclose(pro_max_sum, std_max_sum)

        pro_min_sum = DistroNormalALG.expect_provided_min_sum_X(mean_X, sigma_X, size_sum)
        std_min_sum = DistroNormalALG.expect_standard_min_sum_X(size_sum)
        print(f"pro_min_sum={pro_min_sum:+} std_min_sum={std_min_sum:+}")
        assert np.isclose(pro_min_sum, std_min_sum)