The law of large numbers, part II: Estimating expected values with Monte Carlo analysis

Photo by Daniel Frank from Pexels

Little rephrasing from the first part about this important theorem:

The average of a sequence of random variables from the same distribution converges to the expected value of that distribution. — “Bayesian Methods for Hackers: Probabilistic Programming and Bayesian Inference” by Cameron Davidson-Pilon

In human language, as I had said in the first part, it means, that in short term we can experience wild swings of lucky or unlucky events, but on the longer term, we most probably will get only the expected value, which, generally, is the average of distribution of all events.

Understanding of psychology in trading and investing is exceptionally important, because our nature fools us, we want to believe what we saw or experienced recently (recency, availability, clustering biases), without calculating the average expectancy of all possibilities.

So, we need as much samples as possible of the event we want to test in order to estimate its probable outcomes, but here is the problem, often we don’t have enough samples to estimate that average or can’t do so much of the required (trading) actions.

It is said, that, for example, for daily returns, it is enough to have minimum of 25–30 samples, a lot of scientists are also doing their research with such a low number of population, but coin toss experiment (probably the simplest form of Monte Carlo method) will show us that it is clearly not enough.

We know that fair coin heads or tails probability is exactly 50%, but we may need 100,000 draws (i.e., large number) to achieve that fairness:

from matplotlib import pyplot as plt
from numpy import random
def main():
    total_draws = 100000
    draws = random.randint(2, size=total_draws)
    cnt = 1
    drawed = 0
    probs = []
    for toss in range(total_draws):
        drawed += draws[toss]
        head_probability = drawed / cnt
        probs.append(head_probability)
        cnt += 1
plt.plot(range(total_draws), probs)
plt.axhline(0.5)
plt.show()

Here we can see how probabilities from the first draws jump around like crazy and true expected value is achieved only after large number of draws. If we just had only few initial samples, our judgement would have been extremely skewed to one or another side without us even realizing the truth until it’s too late. (In trading no one can withstand thousands of steps without knowing anything.)

Discovering the true value through trading is expensive, especially if we need thousands of samples. We need a solution.

The solution

One of solutions for this problem can be Monte Carlo experiment.

Monte Carlo experiments are “a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results.” In our case we will use Monte Carlo in one problem class — drawing many samples from the probability distribution.

Long version steps:

1. Determine the distribution of in-sample daily returns.
2. Get parameters for the selected distribution.
3. Draw many samples for sample paths.
4. Plot paths.
5. Get the expected average.
6. Compare average with the test values.

Here are the functions for the first three steps:

import warnings
warnings.simplefilter('ignore')
from scipy import stats
from pandas import DataFrame
cdfs = [
'norm',            #Normal (Gaussian)
'alpha',           #Alpha
'anglit',          #Anglit
'arcsine',         #Arcsine
'beta',            #Beta
'betaprime',       #Beta Prime
'bradford',        #Bradford
'burr',            #Burr
'cauchy',          #Cauchy
'chi',             #Chi
'chi2',            #Chi-squared
'cosine',          #Cosine
'dgamma',          #Double Gamma
'dweibull',        #Double Weibull
'erlang',          #Erlang
'expon',           #Exponential
'exponweib',       #Exponentiated Weibull
'exponpow',        #Exponential Power
'fatiguelife',     #Fatigue Life (Birnbaum-Sanders)
'foldcauchy',      #Folded Cauchy
'f',               #F (Snecdor F)
'fisk',            #Fisk
'foldnorm',        #Folded Normal
'frechet_r',       #Frechet Right Sided, Extreme Value Type II
'frechet_l',       #Frechet Left Sided, Weibull_max
'gamma',           #Gamma
'gausshyper',      #Gauss Hypergeometric
'genexpon',        #Generalized Exponential
'genextreme',      #Generalized Extreme Value
'gengamma',        #Generalized gamma
'genlogistic',     #Generalized Logistic
'genpareto',       #Generalized Pareto
'genhalflogistic', #Generalized Half Logistic
'gilbrat',         #Gilbrat
'gompertz',        #Gompertz (Truncated Gumbel)
'gumbel_l',        #Left Sided Gumbel, etc.
'gumbel_r',        #Right Sided Gumbel
'halfcauchy',      #Half Cauchy
'halflogistic',    #Half Logistic
'halfnorm',        #Half Normal
'hypsecant',       #Hyperbolic Secant
'invgamma',        #Inverse Gamma
'invnorm',         #Inverse Normal
'invweibull',      #Inverse Weibull
'johnsonsb',       #Johnson SB
'johnsonsu',       #Johnson SU
'laplace',         #Laplace
'logistic',        #Logistic
'loggamma',        #Log-Gamma
'loglaplace',      #Log-Laplace (Log Double Exponential)
'lognorm',         #Log-Normal
'lomax',           #Lomax (Pareto of the second kind)
'maxwell',         #Maxwell
'mielke',          #Mielke's Beta-Kappa
'nakagami',        #Nakagami
'ncx2',            #Non-central chi-squared
'ncf',             #Non-central F
'nct',             #Non-central Student's T
'norm'           # Normal
'pareto',          #Pareto
'powerlaw',        #Power-function
'powerlognorm',    #Power log normal
'powernorm',       #Power normal
'rdist',           #R distribution
'reciprocal',      #Reciprocal
'rayleigh',        #Rayleigh
'rice',            #Rice
'recipinvgauss',   #Reciprocal Inverse Gaussian
'semicircular',    #Semicircular
't',               #Student's T
'triang',          #Triangular
'truncexpon',      #Truncated Exponential
'truncnorm',       #Truncated Normal
'tukeylambda',     #Tukey-Lambda
'uniform',         #Uniform
'vonmises',        #Von-Mises (Circular)
'wald',            #Wald
'weibull_min',     #Minimum Weibull (see Frechet)
'weibull_max',     #Maximum Weibull (see Frechet)
'wrapcauchy',      #Wrapped Cauchy
'ksone',           #Kolmogorov-Smirnov one-sided (no stats)
'kstwobign']       #Kolmogorov-Smirnov two-sided test for Large N
def distr_finder(data):
    res = []
    for cdf in cdfs:
        try:
            params = eval('stats.'+cdf+'.fit(data)')
            D, p = stats.kstest(data, cdf, args=params)
            D = round(D, 5)
            p = round(p, 5)
            res.append([cdf.ljust(16), D, p])
        except Exception as err:
            continue
    return DataFrame(res, columns=['Distribution', 'D', 'p']).sort_values('D')

def get_params(cdf, data):
    return eval('stats.'+cdf+'.fit(data)')

def gen_data(cdf, params, size):
    roll = eval('stats.'+cdf+'.rvs(params[0], size=size)')
    return roll

Paths generation:

from matplotlib import pyplot as plt
from pandas import DataFrame
from app.models.distribution import distr_finder, gen_data, get_params
from app.data import get_pickle
from app.utils import train_test_split
def main():
    data = get_pickle('tiingo', 'SPY')['SPY_AdjClose'].pct_change().dropna()
    train, test = train_test_split(data, 0.3)
    distr = distr_finder(data=train.values)
    paths = []
    params = get_params(cdf=distr.iloc[0]['Distribution'], data=train.values)
    plor_paths = True
    for path in range(10000):
        samples = gen_data(cdf='dgamma', params=params, size=len(test))
        paths.append(samples)
        if plor_paths:
            df = DataFrame(samples)
            plt.plot(df.cumsum())
    total = DataFrame(paths).T.sum(axis=1)
    total = total / len(paths)
    plt.plot(range(len(test)), test.cumsum().values, lw=3, color='r')
    plt.plot(range(len(test)), total.cumsum().values, lw=3, color='g')
    plt.show()

Let’s see what we got. We got 10,000 paths from dgamma distribution, strange fit for my liking:

Then we averaged those paths and got the following expected mean (green) for test data (red):

As we can see, for the first few out of sample years Monte Carlo average follows pretty closely the true data. For the far future, of course, estimated information fades due to fundamental regime change.

So, let’s review everything. We had 30% of SPY adjusted close daily return events. We found the distribution for those events. Then we got parameters for this distribution, generated 10,000 paths with size equal to test data size (remaining 70% of SPY adjusted close data), averaged paths and (visually) compared average to true out of sample data. All this, of course, involves many other details, for example, in quantitative finance, we don’t estimate correctness of predictions visually, but we still can derive good mots probable estimates of the future without any knowledge of that future.