new-repeated-simulation.py

import json
import os
import sys
from copy import deepcopy

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.backends.backend_pdf import PdfPages

"""
This file runs the same simulation multiple times with the same parameters and averages the results.
The final product is a PDF file with multiple plots.
"""


def usage():
	print(
		"usage: python3 new-repeated-simulation.py (infl | not-infl) [<output-pdf>] [<output-csv>]"
	)
	exit()


if len(sys.argv) < 2:
	usage()

# A very small constant
EPS = np.finfo(float).eps

# Sometimes the price might get too high or too low and Python complains, in
# those chases, either tweak the parameters or reduce the time horizon.
TIME_HORIZON = 3000

# Number of independent runs to perform to average the results
RUNS = 10

# Save the results in the Overleaf folder for automatic upload
home_directory = os.path.expanduser("~")
CSV_DEST_TEMPLATE = (
	home_directory + "/Dropbox/Applicazioni/Overleaf/LOB & RL/plots/learning-simulation.csv"
)

# Write in the CSV one round every {CSV_ROUND_SKIP} to avoid crowded and slow plots
CSV_ROUND_SKIP = 10

# ---- Price noise

NOISE_STD = 0.5

# ---- Define feasibility constraints

MIN_INV = 0
MIN_CASH = 0

# --- A more readable dict


class dotdict(dict):
	"""dot.notation access to dictionary attributes"""

	__getattr__ = dict.get
	__setattr__ = dict.__setitem__
	__delattr__ = dict.__delitem__

	def __add__(self, other):
		return dotdict({k: v + other[k] for (k, v) in self.items()})

	def __sub__(self, other):
		return dotdict({k: v - other[k] for (k, v) in self.items()})

	def __truediv__(self, div):
		return dotdict({k: v / div for (k, v) in self.items()})

	def __pow__(self, p):
		return dotdict({k: v**p for (k, v) in self.items()})


# ---- Choose a strategy

if sys.argv[1] == "infl":
	PHI = 0.7
	KALPHA = 0.5
	KBETA = 0.5
	VALPHA = 0.5
	VBETA = 0.5
elif sys.argv[1] == "not-infl":
	PHI = 0.3
	KALPHA = 0.5
	KBETA = 0.5
	VALPHA = 0.5
	VBETA = 0.5
else:
	usage()


KAPPA = PHI * KALPHA * VALPHA - (1 - PHI) * KBETA * VBETA
MU = PHI * np.log(1 + KALPHA * VALPHA) + (1 - PHI) * np.log(1 - KBETA * VBETA)


def inflationary():
	return MU > 0


# ---- Set initial state

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

# --- Quantity selection


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)


def pick_quantity(state, side):
	return KALPHA**2 * A(state) if side else -(KBETA**2) * B(state)


# ---- Misc


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


def reject_outliers(mean, std, err=6):
	return np.clip(mean, a_max=std * err, a_min=-std * err)


# ---- Price evolution


def noise():
	# return np.random.uniform(0.6, 1.4)
	return np.random.lognormal(0, NOISE_STD)


def delta(state, quantity):
	if quantity > 0:
		return state.price * KALPHA * VALPHA
	return -state.price * KBETA * VBETA


def update(state, quantity):
	dt = delta(state, quantity)
	epsilon = noise()
	next_price = state.price + dt
	price_difference = epsilon * dt - (1 - epsilon) * state.price
	return dotdict(
		{
			"price": next_price * epsilon,
			"tradedcash": next_price * quantity,
			"kappa": dt / state.price,
			"taker": dotdict(
				{
					"inv": state.taker.inv + quantity,
					"cash": state.taker.cash - next_price * quantity,
					"wealth": state.taker.wealth + price_difference * state.taker.inv,
				}
			),
			"maker": dotdict(
				{
					"inv": state.maker.inv - quantity,
					"cash": state.maker.cash + next_price * quantity,
					"wealth": state.maker.wealth + price_difference * state.maker.inv,
				}
			),
		}
	)


# ---- Data visualization

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


def plot_history(history, std):
	time = np.arange(TIME_HORIZON)

	Ps = np.array([s.price for (s, _) in history])
	std_Ps = np.array([s.price for (s, _) in std])
	Ps = reject_outliers(Ps, std_Ps)

	W_Ts = np.array([s.taker.wealth for (s, _) in history])
	std_W_Ts = np.array([s.taker.wealth for (s, _) in std])
	W_Ts = reject_outliers(W_Ts, std_W_Ts)

	W_Ms = np.array([s.maker.wealth for (s, _) in history])
	std_W_Ms = np.array([s.maker.wealth for (s, _) in std])
	W_Ms = reject_outliers(W_Ms, std_W_Ms)

	Qs = np.array([q for (_, q) in history])
	std_Qs = np.array([q for (_, q) in std])
	Qs = reject_outliers(Qs, std_Qs)

	I_Ts = np.array([s.taker.inv for (s, _) in history])
	std_I_Ts = np.array([s.taker.inv for (s, _) in std])
	I_Ts = reject_outliers(I_Ts, std_I_Ts)

	I_Ms = np.array([s.maker.inv for (s, _) in history])
	std_I_Ms = np.array([s.maker.inv for (s, _) in std])
	I_Ms = reject_outliers(I_Ms, std_I_Ms)

	C_Ts = np.array([s.taker.cash for (s, _) in history])
	std_C_Ts = np.array([s.taker.cash for (s, _) in std])
	C_Ts = reject_outliers(C_Ts, std_C_Ts)

	C_Ms = np.array([s.maker.cash for (s, _) in history])
	std_C_Ms = np.array([s.maker.cash for (s, _) in std])
	C_Ms = reject_outliers(C_Ms, std_C_Ms)

	QPs = np.array([s.tradedcash for (s, _) in history])
	std_QPs = np.array([s.tradedcash for (s, _) in std])
	QPs = reject_outliers(QPs, std_QPs)

	As = np.array([A(s) for (s, _) in history])
	# std_As = np.array([A(s) for (s, _) in std])

	Bs = np.array([B(s) for (s, _) in history])
	# std_Bs = np.array([B(s) for (s, _) in std])

	kappas = np.array([s.kappa for (s, q) in history])

	mpl.rcParams.update(
		{
			"font.size": 16,
			"text.usetex": True,
			"text.latex.preamble": r"""
			    \usepackage{libertine}
			    \usepackage[libertine]{newtxmath}
			    \newcommand{\taker}{{\mathbb{T}}}
				\newcommand{\maker}{{\mathbb{M}}}
		    """,
			"axes.linewidth": 0.4,
			"lines.linewidth": 1.5,
			"lines.markersize": 3.5,
			"xtick.direction": "in",
			"ytick.direction": "in",
			"xtick.minor.visible": False,
			"ytick.minor.visible": False,
			"xtick.major.width": 0.4,
			"ytick.major.width": 0.4,
			"grid.linewidth": 0.2,
			"axes.grid": True,
			"grid.linestyle": "--",
			"savefig.pad_inches": 0.02,
			"axes.prop_cycle": plt.cycler(color=plt.cm.tab10.colors),
		}
	)

	taker_color = "C1"
	maker_color = "C2"
	other_color = "C4"
	other_other_color = "C6"
	alt_color = "C9"
	zero_color = "grey"
	# initial_color = "C6"

	pp = PdfPages(sys.argv[2])

	FIG_SIZE = (10, 4)
	fig, axs = plt.subplots(1, 2, figsize=FIG_SIZE, constrained_layout=True, sharex=True)

	axs[0].fill_between(time, Ps + std_Ps, Ps - std_Ps, color=other_color, alpha=0.3)
	axs[0].plot(Ps, c=other_color)
	# axs[0].set_title("$P_t$")
	axs[0].set_yscale("log")

	axs[1].fill_between(time, W_Ms + std_W_Ms, W_Ms - std_W_Ms, color=maker_color, alpha=0.3)
	axs[1].plot(W_Ms, label="$W^\\maker_t$", c=maker_color)
	axs[1].fill_between(time, W_Ts + std_W_Ts, W_Ts - std_W_Ts, color=taker_color, alpha=0.3)
	axs[1].plot(W_Ts, label="$W^\\taker_t$", c=taker_color)
	# axs[1].set_title("Wealth")
	axs[1].set_yscale("log")
	axs[1].legend()

	pp.savefig(fig)
	fig, axs = plt.subplots(1, 2, figsize=FIG_SIZE, constrained_layout=True, sharex=True)

	axs[0].fill_between(time, Qs + std_Qs, Qs - std_Qs, color=other_color, alpha=0.3)
	axs[0].scatter(time, Qs, s=0.5, c=other_color)
	if not inflationary():
		axs[0].plot(avg(Qs), label="avg", c=alt_color, linestyle="--", lw=2)
	axs[0].set_yscale("symlog", linthresh=1e-8)
	axs[0].yaxis.get_major_locator().numticks = 6
	# axs[0].set_title("$Q_t$")

	axs[1].fill_between(time, QPs + std_QPs, QPs - std_QPs, color=other_color, alpha=0.3)
	axs[1].scatter(time, QPs, s=0.5, c=other_color)
	if inflationary():
		axs[1].plot(avg(QPs), label="avg", c=alt_color, linestyle="--", lw=2)
	axs[1].axhline(0, c=zero_color, alpha=0.4)
	# axs[1].set_title("$Q_t P_{t \\text{+} 1}$")
	# axs[1].legend()

	pp.savefig(fig)
	fig, axs = plt.subplots(1, 2, figsize=FIG_SIZE, constrained_layout=True, sharex=True)

	axs[0].fill_between(time, I_Ts + std_I_Ts, I_Ts - std_I_Ts, color=taker_color, alpha=0.3)
	axs[0].plot(I_Ts, c=taker_color, label="$I^\\taker_t$")
	if not inflationary():
		axs[0].plot(avg(I_Ts), c=alt_color, linestyle="--", lw=2)
	axs[0].fill_between(time, I_Ms + std_I_Ms, I_Ms - std_I_Ms, color=maker_color, alpha=0.3)
	axs[0].plot(I_Ms, c=maker_color, label="$I^\\maker_t$")
	if not inflationary():
		axs[0].plot(avg(I_Ms), c=alt_color, linestyle="--", lw=2)
	# axs[0].set_title("Inventory")
	axs[0].legend()

	axs[1].fill_between(time, C_Ts + std_C_Ts, C_Ts - std_C_Ts, color=taker_color, alpha=0.3)
	axs[1].plot(C_Ts, c=taker_color, label="$C^\\taker_t$")
	if inflationary():
		axs[1].plot(avg(C_Ts), c=alt_color, linestyle="--", lw=2)
	axs[1].fill_between(time, C_Ms + std_C_Ms, C_Ms - std_C_Ms, color=maker_color, alpha=0.3)
	axs[1].plot(C_Ms, c=maker_color, label="$C^\\maker_t$")
	if inflationary():
		axs[1].plot(avg(C_Ms), c=alt_color, linestyle="--", lw=2)
	# axs[1].set_title("Cash")
	axs[1].legend()

	pp.savefig(fig)
	fig, axs = plt.subplots(1, 2, figsize=FIG_SIZE, constrained_layout=True, sharex=True)

	# axs[0].fill_between(time, kappas + std_kappas, kappas - std_kappas, color=other_color, alpha=0.3)
	axs[0].scatter(time, kappas, s=0.5, c=other_color)
	axs[0].plot(avg(kappas), c=alt_color, linestyle="--", lw=2)
	axs[0].axhline(KAPPA, label="$\\kappa$", c=other_other_color)
	# axs[0].axhline(0, c=zero_color, alpha=0.4)
	# axs[0].set_title("$\\delta_t / P_t$")
	axs[0].legend()

	axs[1].plot(As, label="$A_t$", c=other_color)
	axs[1].plot(Bs, label="$B_t$", c=other_other_color)
	# axs[1].set_title("Tradeable amounts")
	axs[1].set_yscale("log")
	axs[1].legend()

	pp.savefig(fig)
	pp.close()


# ---- Misc


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})"


# ---- Play the game


def main():
	if inflationary():
		print(f"The proposed strategy IS inflationary (mu: {MU:.5f})")
	else:
		print(f"The proposed strategy IS NOT inflationary (mu: {MU:.5f})")

	states = [[deepcopy(INITIAL_STATE) for _ in range(RUNS)] for _ in range(TIME_HORIZON)]
	quantities = [[0 for _ in range(RUNS)] for _ in range(TIME_HORIZON)]

	# Set the seed
	rng = np.random.default_rng()  # 11235813

	for run in range(RUNS):
		state = deepcopy(INITIAL_STATE)
		sides = rng.random(TIME_HORIZON) < PHI

		for round, side in enumerate(sides):
			quantity = pick_quantity(state, side)
			state = update(state, quantity)
			sanity_check(state)

			states[round][run] = deepcopy(state)
			quantities[round][run] = quantity

	avg_state = [sum(ss, start=INITIAL_STATE) / RUNS for ss in states]

	# Hand-compute std
	avg_2_state = [(sum(ss, start=INITIAL_STATE) / RUNS) ** 2 for ss in states]
	avg_state_2 = [sum([s**2 for s in ss], start=INITIAL_STATE) / RUNS for ss in states]
	std_state = [(a - b) ** 1 / 2 for (a, b) in zip(avg_state_2, avg_2_state)]

	avg_quantity = [np.sum(qs) / RUNS for qs in quantities]
	std_quantity = [np.std(qs) for qs in quantities]

	avg_history = list(zip(avg_state, avg_quantity))
	std_history = list(zip(std_state, std_quantity))

	print("Final average state")
	print(json.dumps(avg_history[-1][0], indent=2))

	plot_history(avg_history, std_history)


if __name__ == "__main__":
	main()