Kelly-Optimal Options Positions

As mentioned in the introduction, traders' beliefs about the time-\(T\) state of the market, including the price of some tradable security, can be modeled as a probability measure over \(\mathcal{F}_T\), say \(P_T: \mathcal{F}_T \to [0, 1]\). Derivatives prices also reflect market opinion about a probability measure. According to the famous result of Breeden and Litzenberger, in ideal markets, there is a bijective relationship between the set of possible arbitrage-free² vanilla options chains on an asset \(S\) and the set of equivalent (i.e. mutually absolutely continuous with the "true" underlying measure) local martingale (i.e. under which all tradable assets' \(r\)-discounted payoffs, including premiums, are local martingales) measures over \(\sigma(S)\). In particular, for any particular time \(T\), each options chain with expiry at time \(T\) is compatible with exactly one equivalent martingale measure over \(\sigma (S_T)\), which we'll call \(Q_T\), characterized by the fact that \(S_T\)'s density under \(Q_T\) equals \(-e^{rT} \frac{d^2C^{(0, T)}_K(0)}{dK^2}\), the second derivative of its call premiums with respect to strikes times \(e^{rT}\), almost surely.

As it turns out, the explicit formula for \(Q_T\) in terms of vanilla options prices isn't very useful (at least for our purposes), so we relegate the proof to a footnote³. In markets where the set of tradable strikes are not dense, the second derivative of call premiums with respect to strikes isn't well-defined. In \(\sigma(S_T)\)-complete markets, \(Q_T(A)\) could be directly found as the premium of the European contingent claim associated with \(\chi_A\), the indicator function for any event \(A \in \sigma (S_T)\).

In fact, the latter case is more or less the complement of the former: for open intervals \((x, y)\), each \(S^{-1}((x, y))\)-European contingent claim is the limiting case of constructing iron condors by buying puts at strike \(a - \delta\), selling puts at strike \(a\), selling calls at strike \(b\), and buying calls at \(b + \delta\) as \(\delta \to 0\) (with the number of contracts such that the sum of premiums is 1). Note that "synthetic puts", which mimic the returns of real puts, may be constructed by combining a call \(C_K^{(t, T)}\) and a short position \(-C^{(t, T)}_0\), so all we really need is calls (and vice versa, if we only wanted to use puts). This of course implies such a dense set of calls is enough to replicate any \(\sigma(S_T)\)-measurable European contingent claim arbitrarily-well (in any \(L^p\) sense) due to the definition of the Lebesgue sets and density of such indicators in the space of \(\sigma(S_T)\)-measurable random variables, and the limit of the sum of premiums of such a sequence of approximating strategies for some \(\chi_A\) would likewise equal \(Q_T(A)\). So in markets with a continuum of strikes, by trading "close enough" approximating strategies, all \(\sigma(S_T)\)-measurable European contingent claims are for all practical purposes tradable.

In this sense, an options position is essentially specifying some amount to wager on whether or not events in \(\sigma(S_T)\) occur. The "odds", or the payouts for correct "wagers", are dictated⁴ by the equivalent martingale measure: by definition, every \(\sigma(S_T)\)-measurable European contingent claim's payoff is a martingale under this measure, so the premium for \(\chi_A\)'s European contingent claim is \(\mathbb{E}^{Q_T}[\chi_A] = Q_T(A)\) and its payoff function is \(\chi_A - Q_T(A)\). But unless our trader's measure \(P_T\) is identical to \(Q_T\) on \(\sigma(S_T)\), there is some \(A \in \sigma(S_T)\) whose real probability is believed by the trader to be \(P_T(A) \neq Q_T(A)\), and thus the trader's expected value from going long one unit of the bet is \(\mathbb{E}^{P_T} [ \chi_A - Q_T(A) ] = P_T(A) - Q_T(A)\); the trader should clearly go long or short on this bet according to whether this is positive or negative.

This situation is reminiscent of a card counter betting on blackjack who believes the distribution of undealt cards implies an expected profit, or a professional gambler who believes some horses have better chances of winning a race than their odds as posted by the bookies, with the only difference being that the σ-algebra is not discrete. Indeed, this next section will apply to any game with a set of possible bets of which some are submartingales, and can be directly applied to vanilla options without the above formulation as indicator functions, but we'll see that the indicator formulation provides insight into how \(Q_T\) "sets" the odds.

Moreover, let's for the moment suppose our trader expects to have tradable knowledge fairly often, analogous to a card counter who can play many advantageous hands of blackjack, and that this trade's potential profit is roughly representative of that of future knowledgeable trades – we'll discuss the (un)importance of this assumption later.

One intuitive way for our trader/gambler to proceed is to maximize the expected returns after some number of bets \(n\) on the same game/options chain, i.e. taking \((\alpha, \beta) = \arg \max_{(\tilde{\alpha},\tilde{\beta})} \mathbb{E}^{P_T} \left[ \left( \tilde{\alpha} \cdot C^{(t, T)}(S_T) + \tilde{\beta} \cdot P^{(t, T)}(S_T) \right)^n \mid F_t \right]^{\frac{1}{n}} \). Unfortunately, this is likely to be unsatisfactory for large \(n\) because it views losing one dollar the same whether one has \$1 million or has only that single dollar remaining, and thus is likely to quickly go bankrupt. In fact, readers with an analysis background will recognize that as \(n \to \infty\), this expectation converges to the \(L^\infty\)-norm, or the largest value attained by \(\tilde{\alpha} \cdot C^{(t, T)}(S_T) + \tilde{\beta} \cdot P^{(t, T)}(S_T)\) with a nonzero probability – in other words, for large enough \(n\), the solution to the optimization problem only depends on the best-case scenario, which doesn't seem like a great way to make decisions! The discrete case provides intuition for why this is – in the horse race analogy, it clearly follows from linearity of expectation that if the expected return from betting $1 on one horse is higher than that of all others, the same holds for any dollar amount, and one should therefore bet all one's money on this horse. Of course, if that horse wins with probability \(h < 1\), after \(n\) of the same race, the gambler loses everything with probability \(1 - h^n\), which approaches \(1\) as \(n \to \infty\).

If we're uncomfortable with being overwhelmingly likely to lose our entire investment in order to maximize expected returns, this is a good indication that we have a concave utility function, which is to say that our value of an additional dollar decreases as our wealth increases. One sensible concave utility function is the natural logarithm, which as a utility function is often known as the "Kelly criterion". Logarithmic utility functions have been around for a long time, at least since Daniel Bernoulli's proposed solution to the St. Petersburg Paradox, but a proposal of John Kelly in 1956 to use this utility function and his observation that it results in wager amounts of one's wealth proportional to a bet's expected returns for bets on favorable coin tosses and horse races (we'll later see a continuous analog) birthed the name.

Later authors proved that it is "optimal" in several senses, most notably that among any strategy wagering a fixed fraction of one's wealth on successive iterations of the same event, as the number of iterations approaches infinity, the probability that a gambler following the Kelly-optimal strategy has more wealth than a gambler with the same initial wealth following any other fixed-fraction strategy approaches 1. A survey of other senses in which the Kelly criterion is optimal is given here.

Another reason we choose to use logarithmic utility is that it is isoelastic, which is to say that percentage changes in one's wealth affect utility the same amount regardless of the original magnitude of one's wealth, i.e. that rich and poor people both value percentage changes in their wealth equally. This is a trivial consequence of the "log rule" that \(\log(ab) = \log(a) + \log(b)\): if one's wealth is \(a\), the additive utility of an event causing wealth to be multiplied by \(b\) clearly does not depend on \(a\). Whether this is realistic largely depends on individual investor preferences, but this is an assumption large corporations whose shareholders have diverse wealths are essentially forced to make, and it has the added benefit of making models scale-invariant.

Isoelasticity also implies that the number of times a bet is offered doesn't affect optimal wager sizes: after \(n - 1\) iterations of betting, since our current wealth doesn't affect the proportion of it we bet, we will bet the same fraction as in the 1-iteration case regardless of the outcomes of the earlier bets, and the rest follows by induction. So if we had assumed logarithmic utility from the start, our earlier assumption that we will be offered a large number of roughly similar games would not have mattered, though such an assumption is still necessary for the optimal limiting properties distinguishing the logarithm from other isoelastic utilities to be meaningful. In fact, empirical studies of risk aversion have shown that humans tend to favor more risk-averse isoelastic utility functions. Such functions are not even computationally intractable if used in the same ways as the logarithm below, and these will likely be a topic of future articles on this site, but this article only deals with logarithmic utility for simplicity.

Now that we have a utility function, we can maximize its expected value over all possible options positions. We could of course do this by maximizing \[\mathbb{E}^{P_T} \left[ \log \left( \hat{\alpha} \cdot C^{(t, T)}(S_T) + \hat{\beta} \cdot P^{(t, T)}(S_T) + 1 \right) \: \middle\vert \: F_t \right]\] over all possible positions \(\hat{\alpha}, \hat{\beta}\) directly (the \(+ 1\) is for our current wealth, which we define as one "unit"; our options payoff functions already included premiums and discounted with respect to \(r\), so our positions in those only convey gains/losses, not our starting wealth). However, a position in multiple options is equivalent to a European contingent claim whose payoff is the sum of the payoffs of the options, and we proved above that options at a dense set of strikes allow us to approximate the payoff of any \(\sigma(S_T)\)-measurable European contingent claim arbitrarily well. So this problem is also equivalent to maximizing \[\mathbb{E}^{P_T} \left[ \log \left( X^{(t, T)} + 1 \right) \: \middle\vert \: F_t \right]\] for \(X\) in the set \(\mathcal{X}\) of all \(\sigma(S_T)\)-measurable European contingent claim payoffs such that \(P_T(X > -1 \mid F_t) = 1\).

How can we characterize this set? In complete markets, every \(\sigma(S_T)\)-measurable random variable is the time-\(T\) value of some tradable contingent claim, and all contingent claims are valued under \(Q\) such that they are martingales and their payoffs have 0 expected value at any time. Thus \(\mathcal{X} + 1\) (elementwise) is just the set of positive \(Q\)-martingales whose \(Q_T\)-expected value is 1. This is exactly the set of possible Radon-Nikodym derivatives of probability measures with respect to \(Q\)! Indeed, clearly \(R_X(A) := \int_A X dQ \) is a probability measure for any such \(X\), and for any \(R \ll Q\) and \(t\), we have \(\mathbb{E}^{Q} \left[ \frac{dR}{dQ} \mid \mathcal{F}_t \right] = \mathbb{E}^{R} \left[ 1 \mid \mathcal{F}_t \right] = 1 \ \text{a.s.}\).

This may seem like it comes out of nowhere, but it's actually a very natural corollary of the discrete case and the Lebesgue differentiation theorem. Recall the case of an indicator contingent claim for some event \(A\), and suppose we partition \(\mathbb{R}_{\geq 0}\) into a finite number of half-open intervals \(A_1 \ldots A_n\). Further suppose our trader may allocate our initial wealth 1 among the claims with payoffs \(\frac{\chi_{A_i} - Q_T(A_i)}{Q_T(A_i)} = \frac{\chi_{A_i}}{Q_T(A_i)} - 1\) for \(i = 1 \ldots n\) (indicators normalized so that betting 1 on any claim risks all initial wealth), with "selling" modeled as "buying" \(A_i^C\) rather than its own transaction. Since the allocations sum to one, they can be naturally extended into a probability measure over the set of unions of the \(A_i\)s (where each union's measure is the sum of the measures of the members). Call this measure \(R_T\). Then the trader's final wealth is \[1 + \sum_{i = 1}^{n} \chi_{A_i} \left(\frac{R_T(A_i)}{Q_T(A_i)} - 1 \right) = \sum_{i = 1}^{n} \chi_{A_i} \frac{R_T(A_i)}{Q_T(A_i)}.\] By a trivial consequence of the Lebesgue differentiation theorem, we have that \(y_n, z_n \to 0\) as \(n \to \infty\) implies \(\lim_{n \to \infty} \frac{R([x - y_n, x + z_n))}{Q([x - y_n, x + z_n))} = \frac{dR}{dQ}\). So modeling the final wealth as a Radon-Nikodym derivative is in this sense the limiting case of all-or-nothing bets on a finite number of outcomes.

Thus, the optimization problem can be rewritten as \[ \mathop{\text{arg} \, \text{max}}_{R_T \text{ a probability measure on } \sigma(S_T)} \mathbb{E}^{P_T} \left[ \log \left( \frac{dR_T}{dQ_T} \right) \: \middle\vert \: F_t \right].\] Readers familiar with information theory will recognize that the solution is \(R_T = P_T\), where the value of the objective function equals the Kullback-Leibler divergence \(D_{KL}(P_T \: \Vert \: Q_T)\).⁵ In fact, the solution to the optimization problem (though not the optimal quantity) doesn't even depend on \(Q_T\)! Indeed, \[\log \left( \frac{dR_T}{dQ_T} \right) = \log \left( \frac{\frac{dR_T}{dP_T}}{\frac{dQ_T}{dP_T}} \right) = \log \left( \frac{dR_T}{dP_T} \right) - \log \left( \frac{dQ_T}{dP_T} \right), \] and the second term is constant in \(R_T\), so the optimizer equals \[ \mathop{\text{arg} \, \text{max}}_{R_T \text{ a probability measure on } \sigma(S_T)} \mathbb{E}^{P_T} \left[ \log \left( \frac{dR_T}{dP_t} \right) \: \middle\vert \: F_t \right]. \] The problem is thus equivalent to maximizing \(D_{KL}(R_T \: \Vert \: P_T),\) which is equivalent to \[ \mathop{\text{arg} \, \text{max}}_{R_T \text{ a probability measure on } \mathbb{R}_{\geq 0}} \mathbb{E}^{P_T \circ S_T^{-1}} \left[ \log \left( \frac{dR_T}{dm} \right) \: \middle\vert \: S_T(F_t) \right], \] the usual definition of cross-entropy from \(R_T\) to \(P_T\), by properties of the KL divergence. Here, composing \(P_T\) with \(S_T^{-1}\) is used to view \(S_T\)'s possible returns, rather than \(\Omega\), as the underlying probability space, so that we can use the Lebesgue measure (as \(\Omega\) need not have the same sort of "canonical" measure). So we have shown that finding the optimal position is equivalent to a cross-entropy minimization problem.

The results of the last section – in particular, the fact that we can optimize expected utility without even considering \(Q\) – implies that we don't even need to know the prices of derivatives contracts at all in order to optimize our position. If we could apply this to real-world market conditions, it would greatly simplify mathematical models informing trades, and high-frequency traders wouldn't even have to wait for (or even pay for) price quotes anymore!

Now, we'll consider the case where bid-ask spreads are significant. We'll start by assuming there is infinite liquidity at bid and ask prices, but that they are not equal. Under these assumptions, we now need two "option payoff vectors" for each vanilla option, one for buying and one for selling, assessed at the ask and bid prices respectively (assuming a purely liquidity-taking strategy). Let's call \(\hat{P}\) and \(\check{P}\) the buying and selling versions respectively, similar for calls, and we'll drop the superscript \((0, 1)\) (there's no ambiguity with the probability measure \(P\) anymore). Then the problem becomes

It seems like we should be able to use our discovery that we don't need to sell calls to get rid of \(\check{C}\) and \(\check{\alpha}\), which would at minimum reduce our computations required per iteration by a quarter. In general, we might not be able to: even if we can construct the same position using only calls for cheaper than only puts (or vice versa), bid-ask spreads could be large enough as to deny us an arbitrage because it might be too expensive to short the more expensive one. However, if we assume that \(\hat{C}_0 = \check{C}_0\) (i.e. that the futures contract and/or underlying instrument has negligible bid-ask spread, which is typical of options markets on individual equities), put-call parity dictates that if market makers behave rationally, the spread should have the same effect for calls and puts at the same strike (otherwise they would be offering a synthetic call for cheaper than an ordinary call, or vice versa). So in such markets, we can get away with solving

Either way, now that we've seen how to proceed when there is infinite liquidity at the bid-ask spread, it becomes trivial to adapt the problem to a deep order book with finite liquidity at a number of different bid and ask ticks – we simply add more position vectors and constraints (including constraints on the maximum number of transactions at a particular tick) until we run out of orders, or until we don't think we feasibly have the liquidity to exhaust them.

This does increase complexity as the number of ticks we need increases for most numerical optimization methods, but not inordinately; the gradient remains quite tractable, and the ~20 useful strikes and ~5 useful ticks each way characteristic of even the most liquid options chains should only result in a \(400 \times 400\) Hessian at most. Moreover, we'll see later that on real-world data, the bottleneck is usually numerically approximating the Hessian rather than inverting it, so the complexity for using the Hessian is quadratic in the number of ticks rather than being governed by the complexity of matrix inversion (somewhere between quadratic and cubic).

We also observe that while the objective function has been modified when we write it in this way, this can still just be viewed as another constraint on the set of measures encoding our "wagers" with the same objective function as in the complete market case, with the caveat that some may now integrate to less than 1. Even with the more practical way to write the objective function, the incomplete markets case could simply be viewed as a special case of the objective function when \(\hat{C}_0 = \check{C}_0\) and \(\hat{P}_0 = \check{P}_0\). This supports the idea that imperfect markets can be viewed as constraining the space of possible strategies rather than changing market dynamics altogether.

In light of this, we may once again attempt to compare our optimal expected geometric growth rate in the case of significant bid-ask spreads compared to \(D_{KL}(P \: \Vert \: Q)\) to get some idea of how "imperfect" the market is with respect to our preferences. However, one should note that the assumption of significant bid-ask spreads again adds further freedom to how we choose \(Q\).

Now, suppose that in addition to having a bid-ask spread, options contracts are only available in discrete increments. Then the above constrained optimization problem with arbitrary increments is a relaxation of the discrete-increment problem, so a numerical optimization of this relaxation can be run until convergence and the optimum converted to discrete increments via randomized rounding. We won't go into implementation details as code can be found in the implementation section, but vanilla randomized rounding (based on relative coordinate-wise distances from the nearest grid points) will work fine here.

Obviously, smaller increments (meaning larger initial wealth, since contract multipliers for options are generally fixed) result in returns can be to the relaxed optimum. Also, care must be taken to ensure the rounding doesn't result in a position that could possibly lose more than one's initial wealth, which is theoretically possible for large increments. But assuming the randomized rounding gives a position with objective function not too far from the optimal one for this integer programming problem, our subjective heuristic of comparing the result to our theoretically-optimal geometric growth rate \(D_{KL}(P \: \Vert \: Q)\) is still meaningful, it just now also depends somewhat on our initial wealth.

This is essentially just some downloads (if needed) and imports.

using Pkg
Pkg.add("YFinance")
Pkg.add("DataFrames")
Pkg.add("Distributions")
Pkg.add("Plots")
Pkg.add("CSV")
Pkg.add("Integrals")
Pkg.add("Random")
Pkg.add("Optim")
Pkg.add("JuMP")
Pkg.add("JSON3")
Pkg.add("BlackBoxOptim")

using YFinance, DataFrames, Distributions, Plots, CSV, LinearAlgebra, Integrals, Random, Optim, JuMP, JSON3, Interpolations, Dates, BlackBoxOptim

There's a package for Julia called YFinance (similar to, but distinct from, the python package) that easily allows us to get options data from Yahoo Finance. This won't give us historical data, only 15-minute delayed quotes, but that should be all we nee for demonstration purposes. We'll use options on the S&P500 index. We can get this data like so:

esputs=DataFrame();
escalls=DataFrame();

if (isfile("esputs.csv") && isfile("escalls.csv") && isfile("esquotesummary.json"))
    esputs = DataFrame(CSV.File("esputs.csv"));
    escalls = DataFrame(CSV.File("escalls.csv"));
    esquotesummary = JSON3.read("esquotesummary.json");
else
    esopts = get_Options("^SPX");
    escalls = DataFrame(esopts["calls"]);
    CSV.write("escalls.csv", escalls);
    esputs = DataFrame(esopts["puts"]);
    CSV.write("esputs.csv", esputs);
    esquotesummary = get_quoteSummary("^SPX");
    open("esquotesummary.json", "w") do f
        JSON3.pretty(f, JSON3.write(esquotesummary));
    end;
end;
return (esputs, escalls)

This data contains a lot of strikes where the bid or ask is very small or the difference between the two is very big – we can remove these as they're not a good reflection of the actual price:

#Filter out strikes where the bid-ask spread is more than 15% of the ask price, or where the bid price is 0, or too far in the money
ctemp = select(escalls[map((x) -> !ismissing(x) .& x, ((escalls.ask .- escalls.bid)./escalls.ask .< 0.15) .& (escalls.bid .> 0) .& (escalls.strike .> 5300)), :], :strike, :bid, :ask);
ptemp = select(esputs[map((x) -> !ismissing(x) .& x, ((esputs.ask .- esputs.bid)./esputs.ask .< 0.15) .& (esputs.bid .> 0) .& (esputs.strike .< 6500)), :], :strike, :bid, :ask);
#const ess0 = esquotesummary["price"].regularMarketPrice #This didn't work well because the price quote was 15 minutes off the options data. So I manually added a rough price and bid-ask spread from that time.
const ess0 = 5662.8
const ess0a = 5663.2
const ess0b = 5662.4

return nothing

Finally, we also need to assume some parameters about the market, particularly interest rates. We'll just use SOFR for simplicity – at the time of writing, that's around \(5%\) in archaic annualized daily simple interest terms, and these options have 5 days to expiry, which means our discounting factor is \((1 + \frac{0.05}{365})^5 \approx 1.00027 \). We'll also use this to remove any static arbitrages stemming from put-call parity (Yahoo finance data somehow actually has some)

const disc = 1.00027
#Remove (some, not guaranteed to be all) static arbitrages
badstrikes = intersect(ctemp.strike, ptemp.strike)[.!vec(stack([(ctemp[ctemp.strike .== s, :ask] .- ptemp[ptemp.strike .== s, :bid] .+ s .>= ess0b*disc) .&
    (ctemp[ctemp.strike .== s, :bid] .- ptemp[ptemp.strike .== s, :ask] .+ s .<= ess0a*disc) for s in intersect(ctemp.strike, ptemp.strike)]))]
pre_escss = ctemp[.!in.(ctemp.strike, Ref(badstrikes)), :]
pre_espss = ptemp[.!in.(ptemp.strike, Ref(badstrikes)), :]

#Remove butterfly arbitrages
badputs = []
for i=1:size(pre_espss)[1];
    for j=i+1:size(pre_espss)[1];
        if 2*j - i <= size(pre_espss)[1] && 2*pre_espss.bid[j] - pre_espss.ask[i] - pre_espss.ask[2*j - i] > 0;
            push!(badputs, pre_espss.strike[j])
        end;
    end;
end;
const espss = pre_espss[.!in.(pre_espss.strike, Ref(badputs)), :]

badcalls = []
for i=1:size(pre_escss)[1];
    for j=i+1:size(pre_escss)[1];
        if 2*j - i <= size(pre_escss)[1] && 2*pre_escss.bid[j] - pre_escss.ask[i] - pre_escss.ask[2*j - i] > 0;
            push!(badcalls, pre_escss.strike[j])
        end;
    end;
end;
const escss = pre_escss[.!in.(pre_escss.strike, Ref(badcalls)), :]

return nothing

While not strictly necessary to simply optimize the position, computing a risk-neutral density is necessary to use the heuristic we devised in the theory section to see just how much better the optimal strategy would be in perfect markets. We don't need anything special, just a quick-and-dirty heuristic, so we'll compute it using typical numerical second-differentiation on bid-ask midpoints, last trade prices for recent enough trades, and last trade prices capped at the bid-ask spread, and compare these to the PDF we'll use with mean shifted to the right place (to give an idea of the shape)

Here's our implementation:

crecent = subset(escalls, :lastTradeDate => d -> d .> DateTime(2024, 9, 17, 15, 17, 30));
clast = max.(min.(crecent.lastPrice, crecent.bid), crecent.ask);
cmid = (crecent.bid .+ crecent.ask)./2
cstr = crecent[:, "strike"];

precent = subset(esputs, :lastTradeDate => d -> d .> DateTime(2024, 9, 17, 15, 17, 30));
plast = max.(min.(precent.lastPrice, precent.bid), precent.ask);
pmid = (precent.bid .+ precent.ask)./2
pstr = precent[:, "strike"];

numdiff(y, x) = [(y[i] - y[i - 1])/(2*(x[i] - x[i - 1])) + (y[i + 1] - y[i])/(2*(x[i + 1] - x[i])) for i=2:length(x)-1]
convolve(y::Vector{Float64}, x::Union{Vector{Int64}, Vector{Float64}}, f::Function) = [sum(y[j]*f(x[j] - x[i]) for j=1:length(y)) for i=1:length(x)]

bothstr = sort(vcat(cstr, setdiff(pstr, cstr)));
bothl = [K in cstr ? (K in pstr ?  clast[findfirst(isequal(K), cstr)]/2 + (ess0 - K + plast[findfirst(isequal(K), pstr)])/2 : clast[findfirst(isequal(K), cstr)]) : ess0 - K + plast[findfirst(isequal(K), pstr)] for K in bothstr];
bothl2 = [K < ess0 ? ess0 - K + plast[findfirst(isequal(K), pstr)] : clast[findfirst(isequal(K), cstr)] for K in bothstr];
bothl3 = [K < ess0 ? ess0 - K + pmid[findfirst(isequal(K), pstr)] : cmid[findfirst(isequal(K), cstr)] for K in bothstr];

bfirst = numdiff(bothl, bothstr);
bsecond = numdiff(bfirst, bothstr[2:end-1]);
blinterp = linear_interpolation(bothstr[3:end-2], bsecond);

bfirst2 = numdiff(bothl2, bothstr);
bsecond2 = numdiff(bfirst2, bothstr[2:end-1]);
blinterp2 = linear_interpolation(bothstr[3:end-2], bsecond2);

bfirst3 = numdiff(bothl3, bothstr);
bsecond3 = numdiff(bfirst3, bothstr[2:end-1]);
blinterp3 = linear_interpolation(bothstr[3:end-2], bsecond3);

x = range(0.99, 1.01, length=1000)
plot(x, [blinterp.(x.*ess0).*ess0, blinterp2.(x.*ess0).*ess0, blinterp3.(x.*ess0).*ess0, pdf.(LogNormal(log((disc)^2/sqrt((disc)^2 + 0.00004)), sqrt(log(1 + 0.00004/(disc^2)))), x)]; labels=["capped" "recent" "midpoints" "lognormal"])
savefig("../images/kelly-options/bl1.png")
return nothing

Here are the Breeden-Litzenberger implied PDFs.

Note that they're bumpy right now, but that's OK, we'll apply Gaussian smoothing later.

dist(s) = pdf(LogNormal(log(disc^2/sqrt(disc + 0.001)), sqrt(log(1 + 0.001/(disc^2)))), s)
gauss(s) = pdf(Normal(0, 0.001), s)/10000
testrange = (bothstr[3]/ess0):0.0001:(bothstr[end-2]/ess0)
plot(testrange, vcat((f -> convolve(dist.(testrange)./(f.(testrange.*ess0).*ess0), collect(testrange), gauss)).([blinterp, blinterp2, blinterp3]), (f -> dist.(testrange)./(f.(testrange.*ess0).*ess0)).([blinterp, blinterp2, blinterp3])); labels=["capped (smoothed)" "recent (smoothed)" "midpoints (smoothed)" "capped" "recent" "midpoints"])
savefig("../images/kelly-options/bl2.png")
return nothing

And here are the solutions to the unconstrained optimization problem based on these risk-neutral densities, smoothed with a bit of Gaussian noise:

All of these look pretty similar for most of the domain, though the one computed from "capped" last trades looks pretty bad in a low region. We choose to use the bid-ask midpoints because the original Breeden-Litzenberger density looked best, but either of the latter two would have worked fine. Finally, since our example belief PDF gives significant probability that the underlying falls outside the area we could compute by Breeden-Litzenberger, well apply tails from a generalized extreme value (here, Frechet due to fat tails) distribution:

function rtailloss3(params::Vector{Float64})
    xi = params[1]
    mu = params[2]
    sigma = params[3]
    return 3*abs(exp(-(1 + xi*(bothstr[end-4] - mu)/sigma)^(-1/xi)) + bfirst3[end-3])^2 +
        abs((1/sigma)*(1 + xi*(bothstr[end-4] - mu)/sigma)^(-(xi + 1)/xi)*exp(-(1 + xi*(bothstr[end-4] - mu)/sigma)^(-1/xi)) - bsecond3[end-2])^2 +
        abs((1/sigma)*(1 + xi*(bothstr[end-3] - mu)/sigma)^(-(xi + 1)/xi)*exp(-(1 + xi*(bothstr[end-3] - mu)/sigma)^(-1/xi)) - bsecond3[end-1])^2;
        #abs((1/sigma)*(1 + xi*(bothstr[end-2] - mu)/sigma)^(-(xi + 1)/xi)*exp(-(1 + xi*(bothstr[end-2] - mu)/sigma)^(-1/xi)) - bsecond3[end])^2;
end
function ltailloss3(params::Vector{Float64})
    xi = params[1]
    mu = params[2]
    sigma = params[3]
    return 3*abs(exp(-(1 + xi*(bothstr[end-4] - mu)/sigma)^(-1/xi)) + bfirst3[end-3])^2 +
        abs((1/sigma)*(1 + xi*(bothstr[end-4] - mu)/sigma)^(-(xi + 1)/xi)*exp(-(1 + xi*(bothstr[end-4] - mu)/sigma)^(-1/xi)) - bsecond3[end-2])^2 +
        abs((1/sigma)*(1 + xi*(bothstr[end-3] - mu)/sigma)^(-(xi + 1)/xi)*exp(-(1 + xi*(bothstr[end-3] - mu)/sigma)^(-1/xi)) - bsecond3[end-1])^2;
end
res = best_candidate(bboptimize(rtailloss3; SearchRange = [(0, 100), (bothstr[1], bothstr[end-4]), (0, 6000)], MaxTime=3, Method = :de_rand_1_bin_radiuslimited));
testrange2 = 5720:0.001:6000
testrange = 1.008:0.0001:(bothstr[end-2]/ess0)

combined3(x) = x > bothstr[end - 3]/ess0 ? (1/res[3])*(1 + res[1]*(x*ess0 - res[2])/res[3])^(-(res[1] + 1)/res[1])*exp(-(1 + res[1]*(x*ess0 - res[2])/res[3])^(-1/res[1])) : blinterp3(x*ess0)
testrange = (bothstr[3]/ess0):0.0001:1.03
plot(testrange, dist.(testrange)./(combined3.(testrange)*ess0))
savefig("../images/kelly-options/bl3.png")
return nothing

And here is the resulting right-extended solution to the unconstrained optimization problem:

Now, we'll implement an options position optimizer using what we found in the theory section. We'll use numerical integration, which is slow but allows us to use one implementation no matter the choice of \(P\). However, we do need to assume a distribution for testing purposes, so we'll say a lognormal distribution with mean \(1.01\) and variance \(0.001\) (the VIX, which roughly reflects 1-year implied volatility, was about 17 at the time, and these options have only two days until expiry, meaning a lognormal variance of \(\frac{0.17^2}{\sqrt{52*2.5}} \approx 0.0025\), but due to fat tails, if we did not lower the variance further, the directional part would be completely dwarfed by a bet on high variance once we optimized). We "hardcode" this distribution to potentially help the compiler out a bit, as the numerical integration will be pretty slow and we want all the speed we can get. Julia has just-in-time compilation, so this doesn't preclude all flexibility, but it would require a second or two to recompile if we wanted to update our distribution. If we wanted to do this in production, it would be better to pass it (or its parameters) as an argument.

put(k, pr) = (s) -> max(0, k/ess0 - s)/disc - pr/ess0
call(k, pr) = (s) -> max(0, s - k/ess0)/disc - pr/ess0
long = (s) -> s/disc - ess0a/ess0
short = (s) -> - s/disc + ess0b/ess0

const lcalls = call.(escss[:, "strike"], escss[:, "ask"]);
const lputs = put.(espss[:, "strike"], espss[:, "ask"]);
const sputs = (-).∘put.(espss[:, "strike"], espss[:, "bid"]);

const payoff = vcat(lcalls, lputs, sputs, long, short);
return nothing

#Proper conversions to log-mean and log-standard deviation
dist(s) = pdf(LogNormal(log(1.01^2/sqrt(1.01^2 + 0.001)), sqrt(log(1 + 0.001/(1.01^2)))), s)

objective(pos::Vector{Float64}) = solve(IntegralProblem.((x, p) -> log(abs(max(2^-1000, (x .|> payoff) ⋅ pos + 1)))*dist(x), 0.2, 3), QuadGKJL()).u

gradient(pos, tightness) = solve(IntegralProblem.((x, p) -> (x .|> payoff).*dist(x)./max(-2^-1000, (x .|> payoff) ⋅ pos + 1), 0.2, 3), QuadGKJL()).u .+ (tightness.*min.(0, pos)).^2

lgrad(pos) = -pos

gradstep = 1600
lagrangestep = 0.1

function optimize(pos, l, firstidx, lastidx)
    last = pos;
    for j=firstidx:lastidx
        #Gradient descent with Nesterov momentum
        grad = gradient(pos .+ 0.5*(pos .- last), sqrt(j))
        last = pos;
        pos += (gradstep/sqrt(j + 1600)).*grad
        if j % 50 == 0
            @show j
            @show objective(pos)
            @show minimum(grad)
        end
    end
    return (pos, l)
end

pos = zeros(length(payoff));
l = zeros(length(payoff));

Base.exit_on_sigint(false)
try
    (pos, l) = optimize(pos, l, 1, 100)
catch e
    if typeof(e) <: InterruptException
        println("caught Interrupt")
        return 0
    end
end

return nothing

(pos, l) = optimize(pos, l, 100, 50000)

x = range(5300, 5800, length=10000)
testrange = (bothstr[3]/ess0):0.0001:1.03
plot([(x, [(x[i]/ess0 .|> payoff) ⋅ pos + 1 for i=1:10000]), (testrange.*ess0, dist.(testrange)./(combined3.(testrange)*ess0))]; ylims=(0, 5), labels=["Real-world position" "\"Perfect market\"-optimal"])
savefig("../images/kelly-options/payoff.png");

return objective(pos)

2.5037119444049063

And here is a graph of the returns of the resulting position (which currently allows arbitrary increments) as a function of the price of the S&P500 index at expiry, compared with the theoretical optimal payoff function based on smoothed bid-ask midpoints (see the information theory part earlier):

Next, we'll implement a randomized rounding routine to find a solution using real-life contract multiples. To do this, we'll have to assume an initial wealth. Since the contract multiple on these options is \($100\) and the index price is roughly \($5650\), one options contract has a notional value of around \($565,000\), and thus the lowest feasible initial wealth for a trader in options like these is something like \($100,000\), which we'll take as our initial wealth. Note that S&P 500 "mini" options are also offered and have a multiplier of \($10\), meaning an investor with \($10,000\) could apply the same strategy to essentially the same effect (though the bid-ask spreads and/or liquidity might be worse, which would change the optimal strategy slightly. Mini options data doesn't seem to be available on Yahoo finance, though). Also, note that the multiplier is 1 for long/short positions as these can effectively be replicated by a long or short position in an ETF tracking the index.

wealth = 10^5
multiplier = 100

scaledpos(pos) = vcat(pos[1:end-1] .* wealth ./ (ess0 * multiplier), pos[end]*wealth/ess0)
invscaledpos(scaledpos) = vcat(scaledpos[1:end-1] .*(ess0 * multiplier) ./ wealth, scaledpos[end]*ess0/wealth)

best = zeros(length(pos));
bestobj = 0
for i=1:100
    candidate = trunc.(scaledpos(pos)) .+ (rand() .<= scaledpos(pos) .% 1)
    candidateobj = objective(invscaledpos(candidate))
    if candidateobj >= bestobj
        best = candidate
        bestobj = candidateobj
    end
end

x = range(5550, 5750, length=10000)
plot(x, [[(x[i]/ess0 .|> payoff) ⋅ pos + 1 for i=1:10000], [(x[i]/ess0 .|> payoff) ⋅ invscaledpos(best) + 1 for i=1:10000]].*wealth; label=["Arbitrary increments" "Discrete increments"])
savefig("../images/kelly-options/payoff2.png");

annotationvec = string.("Calls at strike " , escss[:, "strike"] , ": " , best[1:length(escss[:, "strike"])])
append!(annotationvec, string.("Puts at strike " , espss[:, "strike"] , ": " , best[length(escss[:, "strike"]) .+ (1:length(espss[:, "strike"]))] .- best[length(escss[:, "strike"]) .+ length(espss[:, "strike"]) .+ (1:length(espss[:, "strike"]))]))
push!(annotationvec, "Long/short position: " * string(best[end]))
return annotationvec

Finally, here's a chart comparing the returns of the rounded position and the original position (observe that the rounded position is extremely close):

In practice, this exact algorithm is likely to be unusably slow in production because of the required numerical integration, requiring a few minutes to converge on a Ryzen 5900H. Newton's method also won't help: our bottleneck here is computing numerical integrals in each entry of the gradient or Hessian, not inverting the Hessian, and we'd have to compute the same number of those numerical integrals for the same theoretical convergence guarantees. Quasi-Newton methods might get us into the territory of seconds, but at the cost of possibly failing to converge depending on \(\frac{dP}{dm}\) – that would be bad in production, too! However, if we had an analytical gradient or Hessian, the algorithm would be orders of magnitude faster, likely converging in milliseconds.

Unfortunately, for most parametric families we would want to use, computing the gradient or Hessian analytically is a fairly difficult task, and some families may not even admit a closed-form representation. Let's go into some details of how one might proceed.

First, we can deal with the piecewise nature of the payoff functions by breaking up the domain into intervals between strike prices. Within such intervals the payoff is a linear function of the position. In particular, if \(K_1, K_2\) are strikes with no intermediate strikes between them, then on \((K_1, K_2)\) we have

which means the inside of the log (in the objective function), the denominator (in the gradient), or the square root of the denominator (in the Hessian) can be written as \(ax + b\) on such intervals, greatly simplifying them. Thus for arbitrary \(P\), the extent of the calculus we need to do is finding the antiderivatives of \[ \frac{dP}{dm} (x) \log( \vert ax + b \vert ) \] for the objective function, \[ x \frac{dP}{dm} (x) \frac{1}{ax + b} \] and \[ \frac{dP}{dm} (x) \frac{1}{ax + b} \] for the gradient, and

for the Hessian. Any of these groups of antiderivatives would give a major speedup (even an analytical Hessian with a numerical gradient would be a big improvement), and it may turn out that some are easier than others. But even if we only had a closed-form antiderivative for the objective function, we could probably find closed-form gradients and Hessians with symbolic differentiation software packages, and approximating them numerically would also be feasible as the objective function is fairly smooth in the feasible region. So our main focus will be on finding \[ \int_0^x \frac{dP}{dm} (y) \log( \vert ay + b \vert ) \mathop{dm(y)}.\]

This seems fairly straightforward when written this way, but let's consider what \(\frac{dP}{dm}\) actually is. Most plausible distributions that could model stock returns are the exponent of some roughly symmetric distribution with mean 0. In the case of our toy example with lognormal returns with mean 1.04 and variance 0.004, the density function is something like \[\frac{dP}{dm}(x) = \frac{\exp\left(- \frac{(\log(x) - 0.026036)^2}{0.0263689}\right)}{\left(0.162385x \sqrt{2\pi}\right)}.\]

This results in integrals like

even for our "toy" example, which are highly nontrivial. SymPy and several other symbolic and symbolic-numeric integrators seem to agree, and either error out or return an expression that still contains an integral. Thus, a direct analytical solution and implementation is out of scope for this article.

However, astute readers may recognize that even if a direct solution is impossible, polynomial approximations of \(\frac{dP}{dm}\) should result in integrals that are more or less trivial if used in place of \(\frac{dP}{dm}\), and thus we can (analytically) approximate the objective function, gradient, etc by their Taylor series assuming these are analytic in the appropriate interval. In our toy example, aside from perhaps cases where a denominator or parameter of log is zero in an interval (which will never occur in a feasible region by proofs in earlier sections, though might occur in practice depending on the optimization routine), the integrand is clearly smooth with bounded derivatives over intervals between strikes, so its integral, regardless of whether it has a closed-form expression, also has this property (by the fundamental theorem of calculus) and is thus analytic over such intervals between strikes.

This wouldn't actually change the results from our earlier numerical integration: QuadGKJL, the Gauss-Kronrod quadrature we used, is approximating the integrals based on a slightly more complicated version of the same idea under the hood. However, it could save time because QuadGKJL doesn't know that the terms of the Taylor series will be of the same parametric form each time, just with different coefficients, but we do know this and thus could hard-code it into our own quadrature method. We won't give a second implementation using this idea as it would be somewhat involved, but anyone who might want to implement this for a real-world purpose might benefit from it.