Gambler's Ruin when the payout isn't 1:1 (open problem in math?) - Gambling and Probability

Two Plus Two Forums Other Topics Probability

Gambler's Ruin when the payout isn't 1:1 (open problem in math?)

Post Reply Subscribe

...

08-14-2019 , 11:13 AM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

A gambler starts with a bankroll of b dollars and will repeatedly bet on a coinflip, wagering L dollars to win W with probability p of winning and q=1-p of losing. Given W≠L and p≠.5, what is the probability of busting if the gambler will play until bust?

For now, let's ignore short-stack situations. Maybe we say the minimum bet is L and the game is over if b<L. Or maybe we allow the gambler to go all-in when b<L. But let's deal with that later because the problem is hard enough without that wrinkle. We'll pretend the gambler can still win W even if b<L.

After taking a few unsuccessful stabs at this, I skimmed lots of books, journals and preprints for the answer. The closest thing I could find is this paper by Cong & Sato, but their solution only works when giving/receiving integer odds (n:1). What I seek is a formula that works for any parameters. Given how many sources I've checked (including the few that cite the linked paper), I'm led to believe this is an open problem in mathematics.

So far I've attempted a basic combinatorial approach and I chose the parameters W=2, L=3 and b=2. The gambler is embarking on a one-dimensional random walk (aka a drunkard's walk), but it's more convenient to think of it as a restricted lattice path which I'll describe. (I also encourage you to try thinking of a better way to visualize the problem, since the lattice approach hasn't much helped me in my attempts at this.)

On a piece of graph paper, draw the line y = 1.5x.

The lines of the graph paper are the lattice. Hero starts 2 units left of the line. Hero moves one unit up on a winning flip or one unit right on a losing flip. To touch or cross the line is to go broke.

Another way to depict it is with disjoint line segments:
(0,0) to (2,2)
(2,3) to (4,5)
(4,6) to (6,8)
and so on.

That 2nd representation shows why this problem isn't like the 1:1 case (whose diagram is just y=x). There is an easy pattern for a 45° line, but in this problem the line keeps shifting and changing the pattern.

Suppose Hero starts at (-2,0).

Hero can reach (0,0) with probability q².

Hero can reach (3,1) via 4 different paths, but 2 of the paths cross through (0,0) and thus were already counted. So there are only 2 new paths to count.

There are then 5 new paths to (2,2) and then 14 new paths to (2,3). So far the coefficients have been Catalan numbers, but next the pattern breaks: there are 28 new paths to (3,4).

Here's a list of the coefs: 1, 2, 5, 14, 28, 76, 227, 488, 1391, 4308, 9647, 28301, 89617, 205893, 615444, 1978069, 4624876, 14007082, ...

(After the first several coefs I reverse-engineered the rest by coding a transition matrix.)

No hits on OEIS.

For demonstration, the probability of going broke within the first 17 or 18 flips is as follows (letting v=pq):

q² + 2q²v + 5q²v² + 14qv³ + 28qv⁴ + 76qv⁵ + 227v⁶ + 488v⁷ + 1391v⁸ + 4308v⁸p

If you can find a good way to get those coefficients, or find another combinatorial solution altogether, then all power to you.

A better route might be to express the problem recursively: R(b) = p*R(b+W) + q*R(b-L)

That's a variable-degree difference equation, which I know nothing about and IDK if it's solvable. Specifying W and L might make it easier, though I'm not even sure, because when the odds are irrational it's still a fractional-order equation.

Even with rational odds, we probably need a bit of sophisticated machinery to solve the recurrence because there doesn't seem to be an algebraic shortcut (unlike the 1:1 case). For W=2 and L=3 there is a way to form a closed system of 6 equations, but I think it would only give solutions to specific b's, not a formula for arbitrary b.

Specifying W and L sacrifices generality, but maybe after separately solving for multiple values of W and L we can combine the formulas into one or at least gain clues for how to think about the problem.

Closing remarks:

- I get the impression that problems like these are where analytic combinatorics (generating functions etc) really comes in handy and that without understanding that stuff, one can only have scratched the surface of combinatorics.

- I wouldn't be surprised if the pro move here is to find tight bounds instead of an exact solution.

I'm interested in people's thoughts on the problem.

Quote

08-14-2019 , 04:20 PM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

Quote:

I chose the parameters W=2, L=3 and b=2

Sorry, that should say b=6.

I forgot to mention that on the lattice I described, each eastbound unit represents -3 and each upward unit represents +2.

Just now, on a whim I flipped open The Mathematics of Poker and found a very good approximation on page 285.

Chen & Ankenman start with the 1:1 case and show the familiar algebraic solution, but then with some more algebra they derive an exponential version of the formula. That's neat, but then they extend the exponential solution to apply to any payout distribution, including non-binary outcomes. Furthermore, the shortstack scenario I mentioned isn't an issue.

The formula is e^-ab where a is the real nonzero solution to p⋅e^-a⋅W + q⋅e^a⋅L = 1

For W=2 and L=3, WolframAlpha tells me a = .174212, so with b=6 and p=.7 the the probability of ruin is about 35.16%.

They don't mention this, but that's slightly high because the formula rests on the assumption that R(b) = R(1)^b, which isn't quite true when the payout isn't 1:1.

R(1) = q + p*R(3)

But the chance of losing $1 when you have $6 is smaller than R(1) because unlike when you have $1, losing a flip would lose you $3.

I wrote a simulation in Julia and the actual answer is closer to 34.2%.

Solving a different exponential equation for each parameter set isn't the most practical, but this is better than nothing!

Here's my sim code for those who want to try the problem and check their answers:

Code:

function rorsim(w::Float64, capital::Float64, lessp::Int64, s::Int64, sims::Int64, games::Int64, strict=true)
    busts = 0
    for n=1:sims
        b = capital
        t = 0
        while t<games && b>0
            r = rand(1:s)
            if r<lessp
                if strict
                    b≥1 ? b+=w : b += w*b
                else b+=w end
            else b-=1 end
            t += 1
        end
        b≤0 && begin busts+=1 end
    end
    return busts/sims
end

I may be the only one here with Julia installed, but there's a site where you can run julia code in your browser: https://tio.run/#julia1x

Edit: oh, if you want to run it on that site, you'll have to replace

Code:

return busts/sims

with

Code:

print(busts/sims)

Quote

08-14-2019 , 10:34 PM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

I've coded Chen & Ankenman's solution:

Code:

function rorexp(p::Vector{Float64}, w::Vector{T}, b::Number, tol=2^-10) where T<:Number
    len = length(p)
    bound = 1+tol  # Technically (1-tol, 1+tol), but >1 gives a more accurate answer.
    step = -1/8
    x = step
    while true
        sum = 0.0
        for k in 1:len begin sum += p[k]*exp(w[k]*x) end end
        if (sum>bound && step<0) || (sum<1 && step>0)
            step /= -2
        elseif 1 < sum < bound
            return exp(b*x)
        end
        x += step
    end
end

Code:

julia> rorexp([.7, .3], [2, -3], 6)
0.34830125444308524

It handles an arbitrary number of outcomes. It solves for x, getting the equation within desired proximity (tol) to 1, defaulted to 1/1024 if you don't enter a 4th parameter. It always errs on the side >1 (between 1 and 1+tol) because that gives a lower probability estimate than if the equation is <1. If you don't input a tight tolerance, this makes the code's answer slightly smaller than the formula's, which partly compensates (in an uncontrolled way) for the formula's answer being high.

As for the two-outcome sim I shared earlier, I don't remember why I did the inputs and RNG like that but it seems goofy. Here's an improvement:

Code:

function rorsim(p::Float64, w::Number, capital::Number, sims::Int64, games::Int64)
    busts = 0
    for n=1:sims
        b = capital
        t = 0
        while t<games && b>0
            r = rand()
            if r<=p
                b≥1 ? b+=w : b += w*b
            else b-=1 end
            t += 1
        end
        b≤0 && begin busts+=1 end
    end
    return busts/sims
end

Quote

08-15-2019 , 11:35 AM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

In the off chance that anyone cares, here's a revision to the MOP formula code:

Code:

function rorexp(p::Vector{Float64}, w::Vector{T}, b::Number, tol=2^-10) where T<:Number
    len = length(p)
    bound = 1+tol  # Technically (1-tol, 1+tol), but >1 gives a more accurate answer.
    step = -1/8
    x = step
    while true
        sum = 0.0
        for k in 1:len begin sum += p[k]*exp(w[k]*x) end end
        if (sum>bound && step<0) || (sum<1 && step>0)
            step /= -2
        elseif 1 < sum < bound
            return exp(b*x)
        end
        x += step
        if x>=0
            x += (step /= -2)
            while x>=0
                x += (step/=2)
            end
        end
    end
end

If anyone is looking at MOP and wondering why their table on 287 has lower percentages than my code (for the SNG example), it's because their value for alpha is a little off (it should be .0006056065). When you plug that in for b=500 you get 73.87% and my code gets that or a little lower depending on the tolerance parameter, as expected.

My guess is they used Excel's solver.

I'll stop rambling to myself ITT until I come across / come up with another solution worth mentioning.

Quote

08-21-2019 , 01:36 AM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

My research prior to this thread was shoddy.

The first source I had checked was Feller (vol 1), but not thoroughly enough, because he addresses this very problem. He uses a generating function to arrive at an upper and lower bound, though the bounds are a bit loose.

Feller's bounds aren't very computable for humans, so others after him came up with approximations to them (e.g. Ethier/Khoshnevisan and Hurlimann 2004). These days, I think Feller's method would be practical for a computer, but I haven't tried.

Others came up with point estimates. There's a simple one by Peter Griffin (not the Peter Griffin) in his gambling books and in the Ethier & Khoshnevisan paper. There is Sileo 1992, which isn't so good according to Canjar 2000.

Not sure if relevant, but there is a trove of actuarial science literature discussing ruin problems, ruin theory and risk theory. After stumbling on that, I noticed that Feller vol 2 covers some of that stuff.

The above is all well and good, but I may have hit the jackpot finding Katriel 2018, who has what I interpret to be an exact formula applying to rational payouts. I'll try to coding it when I get a chance.

But meanwhile, here's an algorithm that iterates the recursive formula until it converges to the answer (for rational payouts):

Code:

function seidelgauss(p::Vector{Float64}, w::Vector{Int64}, B::Int64, T::Int64, tol=10^-9)
    dig = ceil(Int64, abs(log(10,tol))-2)
    pwlen = length(p)
    maxloss = abs(minimum(w))
    r = zeros(Float64, maximum(w)+T)  # T is the dollar target at which we assume RoR ≈ 0
    r[1 : maxloss] .= 1.0
    while true
        prev = r[B+1]
        for j=(maxloss+1):T
            r[j] = 0.0
            for k=1:pwlen begin r[j] += p[k]*r[j+w[k]] end end
        end
        abs(r[B+1]-prev)<tol && break
    end
    return round(r[B+1]; digits=dig)
end

In addition, I made multiple changes to the monte carlo:

- I expanded it to allow for more than two outcomes.

- The house's bankroll can now be finite, which is also useful for modeling the infinite case because you can set a high profit target (aka a high house bankroll) that, once reached, makes RoR virtually impossible. This is more efficient than setting a firm # of games per sim, because your target doesn't need to be high to get very close to the answer. For this reason, the game limit is now infinite by default.

- There are now two ways to handle the shortstack case, described in the comments.

Code:

# For better accuracy, use integer W & L when practical.
# For no grossTarget (ie to play against an infinite bankroll)), put 0.
# Shortstacking false == consider it a bust when you have less capital than the max possible loss.
# Shortstacking true == allow the gambler to go all-in when their stack is less than the normal bet.

function rorsim(p::Vector{Float64}, w::Vector{T}, initCapital::Number, grossTarget::Int64,
         sims::Int64, gameLimit=0, shortstacking=false) where T<:Number
    greedy = (grossTarget <= initCapital)
    eternal = (gameLimit < 1)
    greedy && eternal && error("need a grossTarget or a gameLimit")
    len = length(p)
    maxloss = abs(minimum(w))
    busts = 0
    for n=1:sims
        !eternal && begin g=0 end
        b = initCapital
        while (b<grossTarget || greedy) && ( b ≥ maxloss || (shortstacking && b>0)) && (eternal || g<gameLimit)
            r = rand()
            sum = p[1]
            for k=1:len
                if r<=sum || k==len
                    (b≥maxloss || w[k]≤0) ? b+=w[k] : b += w[k]*b/maxloss
                    break
                else sum += p[k+1] end
            end
            !eternal && begin g+=1 end
        end
        (b≤0 || (b<maxloss && !shortstacking)) && begin busts+=1 end
    end
    return busts/sims
end

As for the combinatorial version of the problem, it sounds to me like that has been tackled too, thanks to a "kernel method".

Quote

08-21-2019 , 01:54 AM

nickthegeek

journeyman

Join Date: Sep 2011 Posts: 273

I used the Chen & Ankenman formula 10/15 years ago when I was trying to determine the risk of ruin when playing 6-max sit and go's.

Given your results, within the bayesian approach you could estimate p1, p2 and pL (the probabilities of arriving first, second or losing) and, based on them, the relative risk of ruin.

Quote

08-21-2019 , 10:05 AM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

Worth noting is that the Chen & Ankenman approach seems to correspond to the "shortstacking=true" mode (now renamed to "movedown") of my monte carlo algo, which I wouldn't have guessed but it's cool if true. It could be a coincidence given the few parameters I've tried, but hopefully not, otherwise the formula's error is larger than I thought.

That's the correct ruin definition to use whenever it's possible to move down in stakes. Seidelgauss() can only handle the other definition.

I'm not sure if Katriel's formula will be faster than seidelgauss(); it depends on how many iterations it will take to numerically solve for multiple complex roots. The formula will at least be more memory-efficient and potentially the better choice for payout distributions that have very large jackpots.

For fun I added an optional progress indicator to rorsim() and seidelgauss(). Also some minor tweaks to rorsim().

rorsim()

Code:

# For better accuracy, use integer W & L when practical.
# For no grossTarget (ie to play against an infinite bankroll)), put Inf.
# Movedown == whether or not you have the ability to move down in stakes if your bankroll falls below the max possible loss.
function rorsim(p::Vector{Float64}, w::Vector{T}, initCapital::Number, grossTarget::Number,
         sims::Int64, movedown=true, gameLimit=Inf, showprogress=false) where T<:Number
    eternal = (gameLimit==Inf)
    grossTarget==Inf && eternal && error("need a grossTarget or a gameLimit")
    len = length(p)
    maxloss = abs(minimum(w))
    busts = 0
    showprogress && begin qtr1, qtr2, qtr3 = ceil(Int64,sims/4), ceil(Int64,sims/2), ceil(Int64, 3*sims/4) end
    for n=1:sims
        !eternal && begin g=0 end
        b = initCapital
        while b<grossTarget && (b ≥ maxloss || (movedown && b>0)) && (eternal || g<gameLimit)
            r = rand()
            sum = p[1]
            for k=1:len
                if r<=sum || k==len
                    (b≥maxloss || w[k]≤0) ? b+=w[k] : b += w[k]*b/maxloss
                    break
                else sum += p[k+1] end
            end
            !eternal && begin g+=1 end
        end
        (b≤0 || (b<maxloss && !movedown)) && begin busts+=1 end
        showprogress && (n==qtr1 || n==qtr2 || n==qtr3) && begin println("n = ",n) end
    end
    return busts/sims
end

seidelgauss()

Code:

function seidelgauss(p::Vector{Float64}, w::Vector{Int64}, b::Int64, T::Int64, tol=10^-9, showprogress=false)
    dig = ceil(Int64, abs(log(10,tol))-2)  # This many digits will be accurate.
    pwlen = length(p)
    maxloss = abs(minimum(w))
    r = zeros(Float64, maximum(w)+T)  # T is the dollar target at which we assume RoR ≈ 0
    r[1 : maxloss] .= 1.0
    showprogress && begin tolprog=[.001, .000001, .000000001, .000000000001, .0000000000000001] end
    while true
        prev = r[b+1]
        for j=(maxloss+1):T
            r[j] = 0.0
            for k=1:pwlen begin r[j] += p[k]*r[j+w[k]] end end
        end
        diff = abs(r[b+1] - prev)
        diff<tol && break
        if showprogress
            for tp in tolprog
                if diff < tp
                    println("Current: ",r[b+1])
                    popfirst!(tolprog)
                    break
                end
            end
        end
    end
    return round(r[b+1]; digits=dig)
end

Quote

08-24-2019 , 07:00 PM

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

I've sped up my MOP implementation (not that it was slow) by making my root-finding algorithm smarter. Now it requires about 2/3 the iterations compared to before.

Code:

function rorexp(p::Vector{Float64}, w::Vector{T}, b::Number, tol=2^-10) where T<:Number
    minimum(w)≥0 && return 0
    len = length(p)
    nearone = 1+tol  # Technically (1-tol, 1+tol), but >1 gives a more accurate answer.
    goal = 1+tol/2
    step = -1/8
    x = step
    prevx = step
    prevdiff = Inf
    zerohi = true
    foundlow = false
    while true
        iter += 1
        sum = 0.0
        for k=1:len begin sum += p[k]*exp(w[k]*x) end end
        diff = abs(sum-goal)
        if diff > prevdiff
            if zerohi && step>0
                zerohi = false
            elseif !foundlow && step<0
                foundlow = true
                !zerohi && begin step/=2 end
            end
            x = prevx + (step /= -2)
        elseif 1 < sum < nearone
            return exp(b*x)
        else
            prevx = x
            prevdiff = diff
            x += (step/=2)
        end
    end
end

Apologies for the lack of commenting. If you wanna see how it works, insert a println(x) at an appropriate line. To pause after each iteration until pressing enter, you can do:

Code:

print(x)
readline()

Quote

08-28-2019 , 07:16 PM

antialias

grinder

Join Date: May 2018 Posts: 590

Quote:

100%

(unless p = 1 and/or L = 0 in which case it's indeterminate because the final condition cannot be met...end with the exception of b=0 in which case the game ends immediately)

Quote

08-29-2019 , 09:44 AM

#10

heehaww

Pooh-Bah

Join Date: Aug 2011 Posts: 5,081

No, when the game is +EV the RoR is not 100%. I worded it "play until bust", but that could mean you're playing forever.

When the payout is 1:1 the exact formula can be derived via algebra and comes to (q/p)^b

It can also be expressed as an infinite (and convergent) series involving binomial coefficients. The lattice path approach gets you that series, whose advantage is that it also gets the answer for finite time.

Btw, the codes I posted earlier aren't quite right (e.g. in one of them I forgot to delete a line and it doesn't even run), so I'll paste again if someone wants to use them. I might be the only one here who uses Julia though.

Quote

Post Reply Subscribe

...