In [1]:
get_ipython().ast_node_interactivity = 'all'
import os
import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(suppress=True)
import matplotlib
import math
from PIL import Image
import pandas as pd
import random
import tqdm
from collections import defaultdict
from scipy.optimize import *
from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor
from queue import SimpleQueue

import warnings
warnings.filterwarnings('ignore')
In [2]:
targets = (
    #'MATRIKS.AKSA',
    #'MATRIKS.AKSEN',
    #'MATRIKS.AKBNK',
    'MATRIKS.ADEL',
    #'MATRIKS.CCOLA',
    #'MATRIKS.ALKIM',
    #'MATRIKS.ANHYT',
    #'MATRIKS.ALARK',
    #'MATRIKS.KAREL',
    #'MATRIKS.ISMEN',
    #'MATRIKS.ARCLK',
    #'MATRIKS.EREGL',
    #'MATRIKS.SISE',
    #'MATRIKS.PETKM',
    #'MATRIKS.SASA',
    #'MATRIKS.VESTL',
    'MATRIKS.VESBE',
    #'MATRIKS.XAUTRY',
    #'MATRIKS.XAGTRY',
    #'MATRIKS.USDTRY',
    #'MATRIKS.EURTRY',
    #'MATRIKS.GBPTRY',
    #'MATRIKS.BIMAS',
    #'MATRIKS.TAVHL',
    'MATRIKS.MTRKS',
    #'MATRIKS.MGROS',
    'MATRIKS.SOKM',
    'MATRIKS.TCELL',
    'MATRIKS.GARAN',
    'MATRIKS.ENJSA',
    'MATRIKS.ASELS',
    'MATRIKS.SAHOL',
    'MATRIKS.TTKOM',
    'MATRIKS.ENKAI',
    'MATRIKS.AEFES',
    'MATRIKS.GSDHO',
    'MATRIKS.AKSEN',
    'MATRIKS.TKFEN',
    'MATRIKS.ALARK',
    'MATRIKS.CIMSA',
    'MATRIKS.ECILC',
    'MATRIKS.SISE',
    'MATRIKS.TTRAK',
    # test
    #'MATRIKS.DOHOL',
    #'MATRIKS.THYAO',
    #'MATRIKS.KCHOL',
    #'MATRIKS.EKGYO',
    # end test
    #'TEFAS.FYD',
    #'TEFAS.TCD',
    #'TEFAS.TCB',
    #'MATRIKS.SAHOL',
    #'MATRIKS.MAVI',
    #'MATRIKS.BUCIM',
    #'MATRIKS.GLDGR',
    #'MATRIKS.GWIND',
    #'MATRIKS.PGSUS',
)
In [3]:
# Fetch data
import miniclickhouse

ch = miniclickhouse.connect(database="algotrade")

qry = """
select * from daily_price_pivot
where `date` > today() - 365
and `MATRIKS.CCOLA` is not null
order by date asc
"""

symbol_data = defaultdict(lambda: [])

for row in ch.execute(qry):
    for target in targets:
        symbol_data[target].append(row[target])
In [4]:
def start_from_zero(vals):
    return [x - vals[0] for x in vals]
In [5]:
def returns(values):
    values = list(values)
    prev = values[0]

    for val in values[1:]:
        yield (val - prev) / prev
        prev = val
In [6]:
def accumulate(values):
    t = 0
    yield t
    for val in values:
        t += val
        yield t
In [7]:
_ = plt.figure(dpi=200)

for target in targets:
    _ = plt.plot(list(accumulate(returns(symbol_data[target]))), label=target)

_ = plt.legend()
Out:
<Figure size 1280x960 with 1 Axes>
In [8]:
def normalize(x):
    x = np.array(x)
    x = x + np.abs(np.min(x))
    return x / np.sum(x)
In [9]:
def portfolio_to_prices(portfolio):
    assert len(portfolio) == len(targets)
    prices = []
    for row in zip(*[symbol_data[target] for target in targets]):
        p = 0
        for v, w in zip(portfolio, row):
            p += v * w
        prices.append(p)
    return prices
In [10]:
_ = plt.figure(dpi=200)
_ = plt.plot(list(accumulate(returns(portfolio_to_prices([1 / len(targets)] * len(targets))))))
Out:
<Figure size 1280x960 with 1 Axes>
In [11]:
def random_weights(N):
    w = np.random.uniform(0, 1, N)
    return normalize(w)

def random_portfolio():
    return random_weights(len(targets))

_ = plt.figure(dpi=200)
for _ in range(512):
    _ = plt.plot(list(accumulate(returns(portfolio_to_prices(random_portfolio())))), linewidth=0.2)
Out:
<Figure size 1280x960 with 1 Axes>

As you can see, there is a quite a bit of variance on the final results.

Let’s try to plot random portfolios to see if we can see the expected shape.

In [12]:
def volatility_flat(values):
    n = len(values)
    # values = np.interp(values, (np.min(values), np.max(values)), (0, 1))
    ideal = np.linspace(values[0], values[-1], num=n)
    # return np.sum(np.abs(ideal - np.array(values))) / n
    return np.sum((ideal - np.array(values)) ** 2) / n
In [13]:
def hillclimb(f, start, iters, directions, lr):
    d = len(start)
    best = start
    
    for _ in range(iters):
        news = [np.random.normal(best, lr) for _ in range(directions)]
        best = min(news, key=f)
    return best
In [14]:
def coord_descent(fn, nparams, iters, step):
    x0 = [random.uniform(0, 1) for _ in range(nparams)]
    for iteration in range(iters):
        for param in range(nparams):
            candidates = []
            for val in np.arange(0.01, 1, step):
                n = list(x0)
                n[param] = val
                candidates.append(np.array(n))
            x0 = min(candidates, key=fn)
    return x0
In [15]:
def uniform_sphere(dimensions, radius):
    pn = np.random.normal(0, 1, dimensions)
    r = np.linalg.norm(pn)
    return radius / r * pn
In [16]:
def pattern_search(fn, nparams, ncandidates=16, decay=0.75, startradius=1):
    x0 = normalize([0.5 for _ in range(nparams)])
    f = fn(x0)
    radius = startradius
    
    while radius > 0.001:
        candidates = []
        for _ in range(ncandidates):
            candidate = x0 + uniform_sphere(nparams, radius)
            candidate = normalize(candidate)
            candidates.append(candidate)
        scores = list(map(fn, candidates))
        ind = min(range(len(candidates)), key=lambda x: scores[x])
        best = candidates[ind]
        bf = scores[ind]
        if bf < f:
            f = bf
            x0 = best
        else:
            radius *= decay
    return x0
In [17]:
target_returns = ["75.0", "85.0", "95.0", "100.0", "110.0"]
target_returns = []
In [18]:
Xs = []
Ys = []
Cs = []

minvol_at_target = {}

Q = SimpleQueue()

def minvol_from_queue():
    while True:
        try:
            rpct, std, params = Q.get(False)
            
            if "min" not in minvol_at_target:
                minvol_at_target["min"] = (std, params)
            else:
                mvs, mvp = minvol_at_target["min"]
                if std < mvs:
                    minvol_at_target["min"] = (std, params)

            if rpct not in minvol_at_target:
                minvol_at_target[rpct] = (std, params)
            else:
                mvs, mvp = minvol_at_target[rpct]
                if std < mvs:
                    minvol_at_target[rpct] = (std, params)
        except Exception:
            break
            

def record_portfolio(params, save=True):
    params = normalize(params)
    values = portfolio_to_prices(params)
    rets = list(returns(values))
    std = np.std(rets)
    # std = volatility_flat(list(accumulate(rets)))
    mean = np.mean(rets)
    pct = list(accumulate(rets))[-1] * 100
    rpct = str(round(pct * 10) / 10)
    
    Q.put((rpct, std, params))
    
    if save and pct > 0 and pct < 200 and std < 1:
        Xs.append(std)
        Ys.append(pct)
        Cs.append(mean / std * np.sqrt(249))
    return values, std, mean, pct

# Lowest volatility
import cma

for _ in range(4096):
    _ = record_portfolio(random_portfolio())
minvol_from_queue()

def fitness(params):
    values, std, mean, pct = record_portfolio(params)
    return std

#print("Initial pattern search")
#_ = pattern_search(fitness, len(targets), 128, 0.7)
print("Initial least-squares")
_ = least_squares(fitness, minvol_at_target["min"][1])
print("Initial CMA")
_ = cma.fmin(fitness, minvol_at_target["min"][1], 0.001, options={"verb_disp": 0})
print("Initial direct")
_ = direct(fitness, [(0, 1)] * len(targets), locally_biased=False)

minvol_from_queue()

values, std, mean, pct = record_portfolio(minvol_at_target["min"][1])
print(f"Min volatility of initial search: {std} @ {pct:.2f}%")

# Min vol at target

def fitness(params):
    values, std, mean, pct = record_portfolio(params)
    pct = round(pct)
    return std + (pct - target) ** 2

for target in tqdm.tqdm(list(np.arange(60, 100, 1))):
    globals()["target"] = target
    # _ = minimize(fitness, [0.3] * len(targets))
    _ = least_squares(fitness, [0.3] * len(targets), jac="3-point")
    # _ = coord_descent(fitness, len(targets), 8, 8)
    # _ = pattern_search(fitness, len(targets), 15, 0.8)
    # _ = direct(fitness, [(0, 1)] * len(targets), locally_biased=False)
    # _ = cma.fmin(fitness, [0.3] * len(targets), 0.001, options={"verb_disp": 0, "bounds": [0, 1]})
    minvol_from_queue()

for tgt in target_returns:
    globals()["target"] = float(tgt)
    print("Target:", tgt)
    _ = cma.fmin(fitness, [0.3] * len(targets), 0.001, options={"verb_disp": 0, "bounds": [0, 1]})
    minvol_from_queue()

#_ = plt.scatter(Xs, Ys, s=3, c=Cs, cmap="gray")
Xs.clear()
Ys.clear()
Cs.clear()

values, std, mean, pct = record_portfolio(minvol_at_target["min"][1])
print(f"Min volatility after search: {std} @ {pct:.2f}%")
Out:
Initial least-squares
Initial CMA
Initial direct
Min volatility of initial search: 0.015215287173937646 @ 102.25%
Out:
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 40/40 [00:41<00:00,  1.04s/it]
Out:
Min volatility after search: 0.015215287173937646 @ 102.25%
Out:
 
In [19]:
Xs = []
Ys = []
Cs = []

for target in dict(minvol_at_target):
    mvs, mvp = minvol_at_target[target]
    values, std, mean, pct = record_portfolio(mvp, save=False)
    Xs.append(std)
    Ys.append(pct)
    Cs.append(mean / std * np.sqrt(249))

_ = plt.figure(dpi=200)
_ = plt.scatter(Xs, Ys, s=3, c=Cs, cmap="copper")
_ = plt.colorbar()


values, std, mean, pct = record_portfolio(minvol_at_target["min"][1])
_ = plt.scatter(std, pct, marker="*", label=f"Min volatility ({pct:.2f}%)")

#values, std, mean, pct = record_portfolio(minvol_at_target["min"][1])
#pct = int(round(pct))
#for target in target_returns:
#    values, std, mean, pct = record_portfolio(minvol_at_target[target][1])
#    _ = plt.scatter(std, pct, marker="*", label=f"Min volatility ({target}%)")
_ = plt.legend()
Out:
<Figure size 1280x960 with 2 Axes>
In [20]:
def print_weights(params):
    weights = []
    for x, y in zip(targets, params):
        weights.append((x, y))
    for x, y in sorted(weights, key= lambda x: x[1], reverse=True):
        print(f"{x} - {y * 100:.2f}")

values, std, mean, pct = record_portfolio(minvol_at_target["min"][1])
pct = int(round(pct))
for target in ("min", *target_returns):
    print(f"{target}%")
    portfolio = minvol_at_target[target][1]
    print_weights(portfolio)
    print()
Out:
min%
MATRIKS.VESBE - 29.51
MATRIKS.SOKM - 13.21
MATRIKS.TCELL - 9.43
MATRIKS.AKSEN - 7.33
MATRIKS.GARAN - 6.96
MATRIKS.ENKAI - 6.68
MATRIKS.ECILC - 6.54
MATRIKS.ADEL - 5.55
MATRIKS.AEFES - 4.76
MATRIKS.MTRKS - 3.73
MATRIKS.ALARK - 3.28
MATRIKS.TKFEN - 1.65
MATRIKS.CIMSA - 1.06
MATRIKS.TTRAK - 0.30
MATRIKS.SISE - 0.00
MATRIKS.ENJSA - 0.00
MATRIKS.GSDHO - 0.00
MATRIKS.SAHOL - 0.00
MATRIKS.TTKOM - 0.00
MATRIKS.ASELS - 0.00

In [21]:
_ = plt.figure(dpi=200)

for target in ("min", *target_returns):
    values, std, mean, pct = record_portfolio(minvol_at_target[target][1])
    _ = plt.plot(list(accumulate(returns(values))), label=f"{target}% return")

_ = plt.legend()
Out:
<Figure size 1280x960 with 1 Axes>
In [22]:
_ = plt.figure(dpi=200)

for sym, weight in zip(targets, minvol_at_target['min'][1]):
    _ = plt.plot([x * weight for x in accumulate(returns(symbol_data[sym]))], label=f"{sym} ({round(weight * 100)}%)")

_ = plt.legend()
Out:
<Figure size 1280x960 with 1 Axes>

Allocate a given budget to the weights.

In [23]:
BUDGET = 45_000
target = 'min'

weights = minvol_at_target[target][1]

total = 0

for symbol, weight in zip(targets, weights):
    total += BUDGET * weight

for symbol, weight in sorted(zip(targets, weights), key=lambda x: x[1], reverse=True):
    print(f"{symbol} - {weight * 100:.02f}% - {total * weight} TRY - {total * weight / symbol_data[symbol][-1]}")
Out:
MATRIKS.VESBE - 29.51% - 13281.345077031805 TRY - 1023.2161076295689
MATRIKS.SOKM - 13.21% - 5946.314963740591 TRY - 212.97689698211286
MATRIKS.TCELL - 9.43% - 4243.0226601318245 TRY - 119.45446678299055
MATRIKS.AKSEN - 7.33% - 3297.696874365242 TRY - 67.52041102303934
MATRIKS.GARAN - 6.96% - 3133.8291980823383 TRY - 104.60044052344253
MATRIKS.ENKAI - 6.68% - 3003.7622855841764 TRY - 106.14001009131366
MATRIKS.ECILC - 6.54% - 2943.619477837192 TRY - 89.20059023749066
MATRIKS.ADEL - 5.55% - 2497.9506723789423 TRY - 26.859684649235938
MATRIKS.AEFES - 4.76% - 2143.6205155356333 TRY - 31.293730153804866
MATRIKS.MTRKS - 3.73% - 1678.2534798657196 TRY - 35.5411579810614
MATRIKS.ALARK - 3.28% - 1473.765197509564 TRY - 17.79909658827976
MATRIKS.TKFEN - 1.65% - 741.2870641534146 TRY - 15.80569433162931
MATRIKS.CIMSA - 1.06% - 478.7728132076663 TRY - 4.797322777631927
MATRIKS.TTRAK - 0.30% - 136.75968483014321 TRY - 0.24539688647073968
MATRIKS.SISE - 0.00% - 1.673258129382668e-05 TRY - 4.085102854938154e-07
MATRIKS.ENJSA - 0.00% - 1.0018461888785301e-05 TRY - 3.0212490617567253e-07
MATRIKS.GSDHO - 0.00% - 8.078640634307484e-06 TRY - 1.8700557023859916e-06
MATRIKS.SAHOL - 0.00% - 4.997802344054734e-07 TRY - 1.161199429380747e-08
MATRIKS.TTKOM - 0.00% - 4.1629242424869473e-07 TRY - 1.9937376640263158e-08
MATRIKS.ASELS - 0.00% - 0.0 TRY - 0.0
In [24]:
sorted(list(minvol_at_target.keys()))
Out [24]:
['101.5',
 '101.9',
 '102.0',
 '102.1',
 '102.2',
 '102.3',
 '102.4',
 '102.5',
 '102.6',
 '102.7',
 '102.8',
 '102.9',
 '103.0',
 '103.1',
 '103.2',
 '103.3',
 '103.4',
 '103.5',
 '103.6',
 '103.7',
 '103.8',
 '103.9',
 '104.0',
 '104.1',
 '104.2',
 '104.3',
 '104.4',
 '104.5',
 '104.6',
 '104.7',
 '104.8',
 '104.9',
 '105.0',
 '105.1',
 '105.2',
 '105.3',
 '105.4',
 '105.5',
 '105.6',
 '105.7',
 '105.8',
 '105.9',
 '106.0',
 '106.1',
 '106.2',
 '106.3',
 '106.4',
 '106.5',
 '106.6',
 '106.7',
 '106.8',
 '106.9',
 '107.0',
 '107.1',
 '107.2',
 '107.3',
 '107.4',
 '107.5',
 '107.6',
 '107.7',
 '107.8',
 '107.9',
 '108.0',
 '108.1',
 '108.2',
 '108.3',
 '108.4',
 '108.5',
 '108.6',
 '108.7',
 '108.8',
 '108.9',
 '109.0',
 '109.1',
 '109.2',
 '109.3',
 '109.4',
 '109.5',
 '109.6',
 '109.7',
 '109.8',
 '109.9',
 '110.0',
 '110.1',
 '110.2',
 '110.3',
 '110.4',
 '110.5',
 '110.6',
 '110.7',
 '110.8',
 '110.9',
 '111.0',
 '111.1',
 '111.2',
 '111.3',
 '111.4',
 '111.5',
 '111.6',
 '111.7',
 '111.8',
 '111.9',
 '112.0',
 '112.1',
 '112.2',
 '112.3',
 '112.4',
 '112.5',
 '112.6',
 '112.7',
 '112.8',
 '112.9',
 '113.0',
 '113.1',
 '113.2',
 '113.3',
 '113.4',
 '113.5',
 '113.6',
 '113.7',
 '113.8',
 '113.9',
 '114.0',
 '114.1',
 '114.2',
 '114.3',
 '114.4',
 '114.5',
 '114.6',
 '114.7',
 '114.8',
 '114.9',
 '115.0',
 '115.1',
 '115.2',
 '115.3',
 '115.4',
 '115.5',
 '115.6',
 '115.7',
 '115.8',
 '115.9',
 '116.0',
 '116.1',
 '116.2',
 '116.3',
 '116.4',
 '116.5',
 '116.6',
 '116.7',
 '116.8',
 '116.9',
 '117.0',
 '117.1',
 '117.2',
 '117.3',
 '117.4',
 '117.5',
 '117.6',
 '117.7',
 '117.8',
 '117.9',
 '118.0',
 '118.1',
 '118.2',
 '118.3',
 '118.4',
 '118.5',
 '118.6',
 '118.7',
 '118.8',
 '118.9',
 '119.0',
 '119.1',
 '119.2',
 '119.3',
 '119.6',
 '119.7',
 '120.1',
 '120.2',
 '120.5',
 '120.6',
 '120.7',
 '120.8',
 '121.0',
 '121.1',
 '121.3',
 '121.4',
 '121.5',
 '121.6',
 '121.9',
 '123.4',
 '123.5',
 'min']