learning-simulation.py

import sys
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt
import tikzplotlib

# A very small constant
# TINY = np.finfo(np.float64).tiny # 2e-308
EPS = np.finfo(np.float64).eps # 2e-16

# Set the seed, if needed
np.random.seed(hash("Ford, you are turning into a penguin, stop!") % 2**31)

DEFAULT_PHI = 0.5

# ---- Length and amount of runs

TIME_HORIZON = 500
RUNS = 10
SAMPLES = 50

# ---- Feasibility and starting position

MIN_INV = 0
MIN_CASH = 0

class dotdict(dict):
    """dot.notation access to dictionary attributes"""
    __getattr__ = dict.get
    __setattr__ = dict.__setitem__
    __delattr__ = dict.__delitem__

INITIAL_STATE = dotdict({
	"price": 1,
	"taker": dotdict({"inv": 1, "cash": 1}),
	"maker": dotdict({"inv": 1, "cash": 1}),
})

# ---- State functions

def A(state):
	return min(state.maker.inv - MIN_INV, (state.taker.cash - MIN_CASH) / state.price)

def B(state):
	return min((state.maker.cash - MIN_CASH) / state.price, state.taker.inv - MIN_INV)

# ---- Misc

def L(strategy, state, quantity, alpha = 0):
	return dotdict({
		"taker": (1 - alpha) * strategy.delta(state, quantity) * state.taker.inv,
		"maker": (1 - alpha) * strategy.delta(state, quantity) * state.maker.inv
	})
		
def sanity_check(state):
	assert state.price > 0, f"negative price ({state.price})"
	assert state.taker.inv > MIN_INV, f"taker inv too low ({state.taker.inv})"
	assert state.maker.inv > MIN_INV, f"maker inv too low ({state.maker.inv})"
	assert state.taker.cash > MIN_CASH, f"taker cash too low ({state.taker.cash})"
	assert state.maker.cash > MIN_CASH, f"maker cash too low ({state.maker.cash})"

def avg(xs):
	return np.cumsum(xs) / np.arange(1, len(xs) + 1)

def t(xs):
	return abs(min(xs[xs != 0], key=abs))

def to_unconstrained(x, a, b):
    # x = np.clip(x, a + EPS, b - EPS)
    return np.log((x - a) / (b - x))

def to_constrained(y, a, b):
    return a + (b - a) / (1 + np.exp(-y))

def dx_dy(y, a, b):
    s = 1 / (1 + np.exp(-y))
    return (b - a) * s * (1 - s)

# w => constrained space
# tilde w => unconstrained space
# w_half => extra gradient intermediate step
# w_one => updated point

def ascent(grad_f, w, par_min, par_max, lr):
	w_one = w + lr * grad_f(w)
	return np.clip(w_one, par_min + EPS, par_max - EPS)

def descent(grad_f, w, par_min, par_max, lr):
	w_one = w - lr * grad_f(w)
	return np.clip(w_one, par_min + EPS, par_max - EPS)

def ascent_in_unconstrained_space(grad_f, w, par_min, par_max, lr):
	tilde_w = to_unconstrained(w, par_min, par_max)
	grad_wrt_tilde_w = grad_f(w) * dx_dy(tilde_w, par_min, par_max)
	tilde_w_one = tilde_w + lr * grad_wrt_tilde_w 
	return to_constrained(tilde_w_one, par_min, par_max)

def descent_in_unconstrained_space(grad_f, w, par_min, par_max, lr):
	tilde_w = to_unconstrained(w, par_min, par_max)
	grad_wrt_tilde_w = grad_f(w) * dx_dy(tilde_w, par_min, par_max)
	tilde_w_one = tilde_w - lr * grad_wrt_tilde_w 
	return to_constrained(tilde_w_one, par_min, par_max)

def extra_ascent_in_unconstrained_space(grad_f, w, par_min, par_max, lr):
	tilde_w = to_unconstrained(w, par_min, par_max)
	grad_wrt_tilde_w = grad_f(w) * dx_dy(tilde_w, par_min, par_max)

	tilde_w_half = tilde_w + lr * grad_wrt_tilde_w
	w_half = to_constrained(tilde_w_half, par_min, par_max)
	grad_wrt_tilde_w_half = grad_f(w_half) * dx_dy(tilde_w_half, par_min, par_max)

	tilde_w_one = tilde_w + lr * grad_wrt_tilde_w_half
	return to_constrained(tilde_w_one, par_min, par_max)

def extra_descent_in_unconstrained_space(grad_f, w, par_min, par_max, lr):
	tilde_w = to_unconstrained(w, par_min, par_max)
	grad_wrt_tilde_w = grad_f(w) * dx_dy(tilde_w, par_min, par_max)

	tilde_w_half = tilde_w - lr * grad_wrt_tilde_w
	w_half = to_constrained(tilde_w_half, par_min, par_max)
	grad_wrt_tilde_w_half = grad_f(w_half) * dx_dy(tilde_w_half, par_min, par_max)

	tilde_w_one = tilde_w - lr * grad_wrt_tilde_w_half
	return to_constrained(tilde_w_one, par_min, par_max)

 # ---- Strategy

class Strategy:
	def __init__(self, kalpha, kbeta, valpha, vbeta, lr=5):
		self.phi = DEFAULT_PHI
		self.kalpha = kalpha
		self.kbeta = kbeta
		self.valpha = valpha
		self.vbeta = vbeta
		self.check()

		self.lr = lr # learning rate


	def kalpha_grad(self, kalpha):
		c = 2 * np.sqrt(2) / 3
		return self.phi * c * self.valpha
    
	def kbeta_grad(self, kbeta):
		c = 2 * np.sqrt(2) / 3
		return - (1 - self.phi) * c * self.vbeta

	def valpha_grad(self, valpha):
		c = 2 * np.sqrt(2) / 3
		return self.phi * c * self.kalpha

	def vbeta_grad(self, vbeta):
		c = 2 * np.sqrt(2) / 3
		return - (1 - self.phi) * c * self.kbeta
	

	def kalpha_ascent(self):
		return ascent_in_unconstrained_space(self.kalpha_grad, self.kalpha, 0, 1/2, self.lr)
		
	def kbeta_ascent(self):
		return ascent_in_unconstrained_space(self.kbeta_grad, self.kbeta, 0, 1/2, self.lr)
	
	def valpha_ascent(self):
		return ascent_in_unconstrained_space(self.valpha_grad, self.valpha, 0, 1, self.lr)

	def vbeta_ascent(self):
		return ascent_in_unconstrained_space(self.vbeta_grad, self.vbeta, 0, 1, self.lr)


	def kalpha_descent(self):
		return descent_in_unconstrained_space(self.kalpha_grad, self.kalpha, 0, 1/2, self.lr)
		
	def kbeta_descent(self):
		return descent_in_unconstrained_space(self.kbeta_grad, self.kbeta, 0, 1/2, self.lr)
	
	def valpha_descent(self):
		return descent_in_unconstrained_space(self.valpha_grad, self.valpha, 0, 1, self.lr)

	def vbeta_descent(self):
		return descent_in_unconstrained_space(self.vbeta_grad, self.vbeta, 0, 1, self.lr)

	def kalpha_extra_ascent(self):
		return extra_ascent_in_unconstrained_space(self.kalpha_grad, self.kalpha, 0, 1/2, self.lr)
		
	def kbeta_extra_ascent(self):
		return extra_ascent_in_unconstrained_space(self.kbeta_grad, self.kbeta, 0, 1/2, self.lr)
	
	def valpha_extra_ascent(self):
		return extra_ascent_in_unconstrained_space(self.valpha_grad, self.valpha, 0, 1, self.lr)

	def vbeta_extra_ascent(self):
		return extra_ascent_in_unconstrained_space(self.vbeta_grad, self.vbeta, 0, 1, self.lr)

	def kalpha_extra_descent(self):
		return extra_descent_in_unconstrained_space(self.kalpha_grad, self.kalpha, 0, 1/2, self.lr)
		
	def kbeta_extra_descent(self):
		return extra_descent_in_unconstrained_space(self.kbeta_grad, self.kbeta, 0, 1/2, self.lr)
	
	def valpha_extra_descent(self):
		return extra_descent_in_unconstrained_space(self.valpha_grad, self.valpha, 0, 1, self.lr)

	def vbeta_extra_descent(self):
		return extra_descent_in_unconstrained_space(self.vbeta_grad, self.vbeta, 0, 1, self.lr)

	def mu(self):
		return self.phi * np.log(1 + 2/3 * np.sqrt(2) * self.kalpha * self.valpha) \
			+ (1 - self.phi) * np.log(1 - 2/3 * np.sqrt(2) * self.kbeta * self.vbeta)

	def eta(self):
		return 1 + self.phi * 2/3 * np.sqrt(2) * self.kalpha * self.valpha \
			- (1 - self.phi) * 2/3 * np.sqrt(2) * self.kbeta * self.vbeta

	def inflationary(self):
		return self.mu() > 0

	def quantity(self, state, sign = None):
		if sign is not None:
			return self.kalpha ** 2 * A(state) if sign else -self.kbeta ** 2 * B(state)
			
		if np.random.random() < self.phi:
			return self.kalpha ** 2 * A(state) 
		return -self.kbeta ** 2 * B(state)

	def delta(self, state, quantity):
		if quantity > 0:
			return state.price * 2/3 * np.sqrt(2) * self.kalpha * self.valpha
		return -state.price * 2/3 * np.sqrt(2) * self.kbeta * self.vbeta

	def generate(inflationary = None):
		# phi = np.random.uniform(0, 1)
		kalpha = np.random.uniform(0, 1/2)
		kbeta = np.random.uniform(0, 1/2)
		valpha = np.random.uniform(0, 1)
		vbeta = np.random.uniform(0, 1)

		s = Strategy(kalpha, kbeta, valpha, vbeta)
		if inflationary is None or s.inflationary() == inflationary: return s
		return Strategy.generate(inflationary)

	def check(self):
		assert 0 < self.phi < 1, f"invalid phi: {self.phi}"
		assert 0 < self.kalpha < 1/2, f"invalid kalpha: {self.kalpha}"
		assert 0 < self.kbeta < 1/2, f"invalid kbeta: {self.kbeta}"
		assert 0 < self.valpha < 1, f"invalid valpha: {self.valpha}"
		assert 0 < self.vbeta < 1, f"invalid vbeta: {self.vbeta}"

	def dist(self, strategy):
		return (self.kalpha - strategy.kalpha) ** 2 + \
			(self.kbeta - strategy.kbeta) ** 2 + \
			(self.valpha - strategy.valpha) ** 2 + \
			(self.vbeta - strategy.vbeta) ** 2

	def __str__(self):
		return f"phi={self.phi:.2f} kalpha={self.kalpha:.2f} kbeta={self.kbeta:.2f} valpha={self.valpha:.2f} vbeta={self.vbeta:.2f} (μ={self.mu():.2f})"

# ---- Price evolution

def update(strategy, state, quantity):
	next_price = state.price + strategy.delta(state, quantity)
	return dotdict({
		"price": next_price,
		"taker": dotdict({
			"inv": state.taker.inv + quantity,
			"cash": state.taker.cash - next_price * quantity,
		}),
		"maker": dotdict({
			"inv": state.maker.inv - quantity,
			"cash": state.maker.cash + next_price * quantity,
		}),
	})

# ---- Data visualization

MAX_INV = INITIAL_STATE.taker.inv + INITIAL_STATE.maker.inv
MAX_CASH = INITIAL_STATE.taker.cash + INITIAL_STATE.maker.cash

def plot_trajectory(ax, params):
	# Keep only a fixed amount of samples
	coarse_trajectory = params[::len(params) // SAMPLES]
	x, y = zip(*coarse_trajectory)
	
	# Add invisible points to size the plot correctly
	# ax.scatter(x, y, s=0)
	
	# Plot arrows inbetween points
	for i in range(len(x) - 1):
		ax.annotate("",
            xy=(x[i+1], y[i+1]),
            xytext=(x[i], y[i]),
            arrowprops=dict(arrowstyle="->", ls="--", color="tab:blue", lw=1)
		)

def plot_all_runs(history):
	starting_strategy = history[0][0][0]
	
	fig, axs = plt.subplots(3, 2, figsize=(10, 15), constrained_layout=True)

	# You get a grid, and you get a grid, and you get a grid
	for ax in axs.flatten(): ax.grid(ls="--")
	
	# Set scatter marker size
	# mpl.rcParams["lines.markersize"] = 0.5

	# history := [[(strategy, state, quantity) for t in time] for i in run]	
	time = np.arange(TIME_HORIZON)
	
	# The evolution of the strategy is deterministic, we can just pick the first run
	mu = [strategy.mu() for (strategy, _, _) in history[0]]
	taker_params = [(strategy.kalpha, strategy.kbeta) for (strategy, _, _) in history[0]]
	maker_params = [(strategy.valpha, strategy.vbeta) for (strategy, _, _) in history[0]]

	# The distance between the optimum and the trajectory of strategies
	optimum = Strategy(1/2 - EPS, 0 + EPS, 1 - EPS, 0 + EPS)
	error = [strategy.dist(optimum) for (strategy, _, _) in history[0]]

	price_avg, price_std = np.zeros((TIME_HORIZON,)), np.zeros((TIME_HORIZON,))
	
	w_taker_avg, w_taker_std = np.zeros((TIME_HORIZON,)), np.zeros((TIME_HORIZON,))
	w_maker_avg, w_maker_std = np.zeros((TIME_HORIZON,)), np.zeros((TIME_HORIZON,))

	for t in range(TIME_HORIZON):
		across_runs = [history[i][t] for i in range(RUNS)]

		price = np.log([state.price for (_, state, _) in across_runs])
		price_avg[t] = np.average(price)
		price_std[t] = np.std(price)
		
		w_taker = np.log([state.price * state.taker.inv + state.taker.cash
		    for (_, state, _) in across_runs])
		w_taker_avg[t] = np.average(w_taker)
		w_taker_std[t] = np.std(w_taker)

		w_maker = np.log([state.price * state.maker.inv + state.maker.cash
		    for (_, state, _) in across_runs])
		w_maker_avg[t] = np.average(w_maker)
		w_maker_std[t] = np.std(w_maker)

	axs[0, 0].plot(mu)
	# axs[0, 0].axhline(optimum.mu(), ls="--", label="optimum")
	axs[0, 0].set_title("μ")
	# axs[0, 0].legend()

	axs[0, 1].plot(price_avg)
	axs[0, 1].fill_between(time, price_avg - price_std, price_avg + price_std, alpha=.1)
	axs[0, 1].set_title("log P_t")
	axs[1, 0].plot(w_taker_avg, label="taker")
	axs[1, 0].fill_between(time, w_taker_avg - w_taker_std, w_taker_avg + w_taker_std, alpha=.1)
	axs[1, 0].plot(w_maker_avg, label="maker")
	axs[1, 0].fill_between(time, w_maker_avg - w_maker_std, w_maker_avg + w_maker_std, alpha=.1)
	# axs[1, 0].set_yscale("log")
	axs[1, 0].set_title("log W_t")
	axs[1, 0].legend()

	axs[1, 1].plot(error, "--")
	# axs[1, 1].fill_between(time, 0, [1/t for t in range(1, TIME_HORIZON + 1)], alpha=.2, label="1/t")
	# axs[1, 1].legend()
	axs[1, 1].set_title("Distance from the optimum")

	plot_trajectory(axs[2, 0], taker_params)
	axs[2, 0].plot(starting_strategy.kalpha, starting_strategy.kbeta, "o", color="tab:blue", label="start")
	axs[2, 0].plot(1/2, 0, "*", color="orange", label="optimum")
	axs[2, 0].set_title("Taker parameters")
	axs[2, 0].set_xlabel("k alpha")
	axs[2, 0].set_ylabel("k beta")
	# axs[2, 0].set_xlim([-0.1, 1.1])
	# axs[2, 0].set_ylim([-0.1, 1.1])
	axs[2, 0].legend()

	plot_trajectory(axs[2, 1], maker_params)
	axs[2, 1].plot(starting_strategy.valpha, starting_strategy.vbeta, "o", color="tab:blue", label="start")
	axs[2, 1].plot(1, 0, "*", color="orange", label="optimum")
	axs[2, 1].set_title("Maker parameters")
	axs[2, 1].set_xlabel("v alpha")
	axs[2, 1].set_ylabel("v beta")
	# axs[2, 1].set_xlim([-0.1, 1.1])
	# axs[2, 1].set_ylim([-0.1, 1.1])
	axs[2, 1].legend()

	title = f"""
        Starting strategy {starting_strategy}
        Starting positions [taker I: {INITIAL_STATE.taker.inv}, C: {INITIAL_STATE.taker.cash}] [maker inv: {INITIAL_STATE.maker.inv}, cash: {INITIAL_STATE.maker.cash}]
        Time horizon: {TIME_HORIZON}    Runs: {RUNS}     Learning rate: {starting_strategy.lr:.4f}
    """
	fig.suptitle(title)

	if len(sys.argv) > 1: plt.savefig(sys.argv[1], dpi=200)
	else: plt.show()

	tikzplotlib.save("learning-simulation.tex")

# ---- Play the game

def main():
	# Draw a random strategy
	# starting_strategy = Strategy.generate()

	# Draw a random inflationary strategy
	# starting_strategy = Strategy.generate(inflationary=True)

	# Draw a random non-inflationary strategy
	# starting_strategy = Strategy.generate(inflationary=False)
	
	# Play this specific strategy
	starting_strategy = Strategy(
		kalpha=0.25,
		kbeta=0.25,
		valpha=0.5,
		vbeta=0.5,
	)

	print(f"""phi={starting_strategy.phi:.2f},
kalpha={starting_strategy.kalpha},
kbeta={starting_strategy.kbeta},
valpha={starting_strategy.valpha},
vbeta={starting_strategy.vbeta},
η={starting_strategy.eta()},
μ={starting_strategy.mu()}""")

	histories = []
	for run in range(RUNS):
		# Pick the side of the trade for the whole horizon
		signs = np.random.rand(TIME_HORIZON) < DEFAULT_PHI

		# Pick which player is updating their parameters on a given round
		turns = [True, False] * (TIME_HORIZON // 2)
		
		strategy = deepcopy(starting_strategy)
		state = INITIAL_STATE
		history = []

		for sign, turn in zip(signs, turns):
			quantity = strategy.quantity(state, sign)
			history.append((strategy, state, quantity))
			state = update(strategy, state, quantity)
			sanity_check(state)

			if run == 0 and len(history) % 10 == 1:
				print(strategy)
			
			# Update strategy alternately
			strategy = Strategy(
				strategy.kalpha_ascent(),
				strategy.kbeta_ascent(),
				strategy.valpha_ascent(),
				strategy.vbeta_ascent(),
			)

		histories.append(history)
		
	plot_all_runs(histories)

if __name__ == "__main__":
	main()