Devenir rentier / effet de levier et stress-test

import numpy as np
import random
import matplotlib.pyplot as plt
from progress.bar import Bar
import time as tm

N = 9*int(365*5/7) ##Market open monday-friday
Nbins = 200
def historical_return(markets):
#Generate a list of daily returns
    import yfinance as yf
    r=[]
    for market in markets :
        df = yf.download(market,period='max',interval='1d')
        close = df['Close']
        for i in range(1,len(close)) : r.append(100*((close[i]-close[i-1])/close[i-1]))
    return r

def gaus(x,a,x0,sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))

def generate_random(on,N):
    H = np.histogram(on,bins=Nbins,density=True)
    y = np.array(H[0])
    x = np.linspace(H[1][0],H[1][-1],len(y))
    R = []
    for i in range(N): ##Méthode de rejet
        search=True
        while search:
            j = random.randint(0,len(x)-1)
            Y=x[j]
            U=random.random()
            if U<=y[j] :
                search=False
                R.append(Y)
    return R

def cumulative(array):
    S = np.sum(array)
    C = [array[0]/S]
    for i in range(1,len(array)):
        C.append(C[-1]+array[i]/S)
    return C

def gaussian_fit(R,array):
    import pylab as plb
    from scipy.optimize import curve_fit
#    N = 200
    H = plt.hist(array,bins=Nbins,label='real data',density=True)

#    Compute gaussian fit
#    x = np.linspace(min(array),max(array),N)
    y = np.array(H[0])
    x = np.linspace(H[1][0],H[1][-1],len(y))
    popt,pcov = curve_fit(gaus,x,y,p0=[1,np.mean(y),np.std(y)])
#    Shell = gaussian_shell(x,popt)

#    Plot setup
    xtick = np.arange(-30,30,step=5)
    xticklabel = []
    for xt in xtick:
        xticklabel.append(str(xt)+' %')
    plt.xticks(xtick)
    plt.xticks(xtick,xticklabel)
    plt.plot(x,gaus(x,*popt),'--',label='gaussian fit')
    plt.yscale('log')
    plt.ylim([0.0001,1])
    plt.xlim([min(x),max(x)])
    plt.xlabel('Daily return')
    plt.ylabel('Normalized')
    plt.hist(R,bins=100,alpha=0.5,label='simulated data',density=True,color='r')
    plt.grid()
    plt.legend()
    plt.title('Comparing gaussian fit to reality')
    plt.show()

def generate_single_scenario(multipliers,R,N):
    daily_returns={}
    daily_cumulated={}
    for m in multipliers :
        daily_returns[m]=[]
        daily_cumulated[m]=[1]
    for i in range(N):
        rtrn=random.choice(R)
        for m in multipliers :
            daily_returns[m].append(1+m*rtrn/100)
            daily_cumulated[m].append(daily_cumulated[m][-1]*daily_returns[m][-1])
    return daily_returns,daily_cumulated

def computeTRI(initialInvest,monthlyInvest,finalValue,nYears) :
    search = True
    order = 1
    new_rate = 1
    old_rate = 1
    old_VAN = np.nan
    new_VAN = np.nan

    while search :
        CF = [-initialInvest]
        for i in range(nYears*12) :
            CF.append(-monthlyInvest*(1+new_rate/100)**(-(i+1)/12))
        CF.append(finalValue*(1+new_rate/100)**(-(nYears+1/12)))
        old_VAN = new_VAN
        new_VAN = np.sum(CF)
        old_rate = new_rate
        if(new_VAN*old_VAN<0) :
            if order == 0.001 : search = False
            order/=10
            if abs(old_VAN) < abs(new_VAN) : new_rate = old_rate
        elif new_VAN<0 : new_rate-=order
        elif new_VAN>0 : new_rate+=order
    return new_rate

def payTaxes(finalValue,initialValue,m):
    tax = 0
    PEA_tax = .172
    CTO_tax = .3
    if m==1 : #PEA
        if finalValue>initialValue : tax = (finalValue-initialValue)*PEA_tax
    else : ##CTO
        if finalValue>initialValue : tax = (finalValue-initialValue)*CTO_tax

    finalValue -= tax
    return finalValue

def compareToRegular(R) :
    oneShot = []
    mulShot = []
    n = 2000
    bar = Bar('Single shot vs. Multiple shot', max=n)
    for i in range(n):
        bar.next()
        df,dc = generate_single_scenario([1],R,N)
        df = df[1]
        initialInvest = 10000
        regularInvest = 500
        totalInvest = initialInvest
        old_portfolio = [df[0]]
        new_portfolio = [df[0]*initialInvest]
        for i in range(1,len(df)):
            old_portfolio.append(old_portfolio[-1]*df[i])
            new_portfolio.append(new_portfolio[-1]*df[i])
            if i>0 and i%22==0 :
                new_portfolio[-1]+=regularInvest
                totalInvest+=regularInvest

        old_portfolio[-1] = payTaxes(old_portfolio[-1],1,1)
        new_portfolio[-1] = payTaxes(new_portfolio[-1],totalInvest,1)
        oneShot_tri = 100*(old_portfolio[-1]**(1/9)-1)
        regular_tri = computeTRI(initialInvest,regularInvest,new_portfolio[-1],9)
        oneShot.append(oneShot_tri)
        mulShot.append(regular_tri)
    bar.finish()

    K = 30
    H,xH,yH = np.histogram2d(oneShot,mulShot,bins=K)
    plt.contourf(np.linspace(xH[0],xH[-1],K),np.linspace(yH[0],yH[-1],K),H,levels=100)
    plt.xlabel('One Shot')
    plt.ylabel('Multiple Shot')
    x=np.linspace(min(oneShot),max(oneShot),100)
    print(min(oneShot),max(oneShot))
    print(min(mulShot),max(mulShot))
    plt.plot(x,x,'r--')
    plt.colorbar()
    plt.show()

def simulate_portfolio(R) :
    multipliers = [1,2,3,4]
    brokerMarginRequest = 0.3
    def generate_multiple_scenario(multipliers,R,N):
        V={}
        for m in multipliers : V[m]=[]
        nHisto=1500
        bar = Bar('Generating histograms', max=nHisto)
        for j in range(nHisto):
            if j>0 and j%1000==0 : tm.sleep(10) ##Short break for CPU
            bar.next()
            df,dc=generate_single_scenario(multipliers,R,N)
            lowestReferenceValue = min(dc[1])

            for m,returns in df.items():
                maxLoss = m*(1-lowestReferenceValue)
                if maxLoss > 1-brokerMarginRequest and m>1 : finalValue = 1e-7
                else : finalValue = np.product(returns)

                finalValue = payTaxes(finalValue,1,m)
                V[m].append(finalValue)
        bar.finish()
        return V

    V = generate_multiple_scenario(multipliers,R,N)

    #Plot histograms and cumulative function
    fig,ax = plt.subplots(len(multipliers),2,sharex=True)
    pos=0
    for m in multipliers :
        #Histogram
        Values = V[m]
        hist, bins = np.histogram(Values, bins=Nbins)
        logbins = np.logspace(np.log10(bins[0]),np.log10(bins[-1]),len(bins))
        ax[pos][0].hist(Values,bins=logbins,label=str(m)+"-multiplier")
        ax[pos][0].grid(True)
        ax[pos][0].legend()
        ax[pos][0].set_xscale('log')

        #cumulative
        C=cumulative(hist)
        logx = np.logspace(np.log10(bins[0]),np.log10(bins[-1]),len(C))
        H,B = np.histogram(Values, bins=logbins)
        C = cumulative(H)
        ax[pos][1].plot(logx,C)
        ax[pos][1].grid(True)
        ax[pos][1].set_ylim([0,1])
        ax[pos][1].set_xscale('log')
        pos+=1

        print(np.mean(Values),np.std(Values))

    plt.show()

def conditional_probability(on):
    def bin2label(bin,new_size):
        return np.linspace(bin[0],bin[-1],new_size)

    matrix,xbins,ybins = np.histogram2d(on[:-1],on[1:],bins=50)
    xbins = bin2label(xbins,matrix.shape[0])
    ybins = bin2label(ybins,matrix.shape[1])
    for i in range(len(ybins)):
        S = np.sum(matrix[:,i])
        if S>0 : matrix[:,i]/=S
    plt.contourf(ybins,xbins,matrix,levels=100)
    plt.xlabel('Yesterday daily-return')
    plt.ylabel('Today daily-return')
    plt.colorbar()
    plt.show()

on = historical_return(['^GSPC'])
R = generate_random(on,30000)
#compareToRegular(R)
#gaussian_fit(R,on)
#simulate_portfolio(R)
conditional_probability(on) #on is "real" data, replace by R for simulated data
#conditional_probability(R)