Randomized Model of Markets populated by Gode and Sunder (1993) ZI robots + Todd Kaplan's Snipers with Sniper Failure -- for section 3 of Brewer and Ratan (2018)

RandomizedSniperFailureModel.py

#!/usr/bin/env python3
# File:  RandomizedSniperFailureModel.py
# Copyright 2018 Paul Brewer, Economic and Financial Technology Consulting LLC
# This code approximates allocation efficiency, source of efficiency losses, and
# Gini Coefficient for profits accumulated over many periods in a market populated
# by a mix of Gode and Sunder (1993) ZI robots and Todd Kaplan's Sniper robots.
# In this code, we only model the "failsafe" case where the snipers fail to find or
# act on exceptional prices or low bid-ask spread and so execute the end-of-period failsafe
# strategy of accepting any existing profitable bid/ask.  Taking an average
# over random shuffling of Snipers imperfectly approximates the interactions over multiple
# periods of the double auction, which is considerably more complex.
# Calculations produced by this code are mentioned as a "software test"   
# in Brewer and Ratan (2019) with the full output cut during peer review.
# This file is open source software.  You may use it according to the LICENSE.
# LICENSE: The MIT LICENSE
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
#  of the Software, and to permit persons to whom the Software is furnished to
# do so, subject to the following conditions:
#  The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
# INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
# PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
# COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
# IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# IMPORTANT:  This code was written for a single use and is fragile to changes.  
# You can safely change the number of repetitions (e.g. 100, 1000, 1000000) 
# Any other changes will probably break mathematical relationships that are hard coded.
# For example, CEP is hard coded to 250, which is appropriate to the research.
# If you change the MV or MC curves, the correct CEP for the model will no longer be 250.
#
from random import shuffle

def MarketMV(i,market):
    if i<1 or i>20 or market<1 or market>2:
        raise IndexError
    D = 1 if market==1 else 0
    unit1 = 450-5*(i-1)
    unit2 = D*(255+5*(i-1))+(1-D)*(350-5*(i-1))
    unit3 = 250-5*(i-1)
    return [unit1,unit2,unit3]

def MarketMC(j,market):
    if j<1 or j>20 or market<1 or market>2:
        raise IndexError
    D = 1 if market==1 else 0
    unit1 = 50+5*(j-1)
    unit2 = D*(245-5*(j-1))+(1-D)*(150+5*(j-1))
    unit3 = 250+5*(j-1)
    return [unit1,unit2,unit3]

CEP = 250
# replications = 1000  5 sec on my machine
# replications = 10000 51 sec on my machine
# 1 million should take about 90 minutes or so
replications = 1000000
maxprofit = 8200.0*replications

for snipers in [0,1,2,4,6,8,10,12,14,16,18,19]:
    out = []
    for market in [1,2]:
        buyerProfits = [0.0]*20
        sellerProfits = [0.0]*20
        includedExtra = 0.0
        omittedInfra = 0.0
        for replication in range(replications):
            ZIBuyers = range(1+snipers,21)
            ZISellers = range(1+snipers,21)
            SniperBuyers = list(range(1,snipers+1))
            SniperSellers = list(range(1,snipers+1))
            for i in ZIBuyers:
                for mv in MarketMV(i,market):
                    if mv>CEP:
                        buyerProfits[i-1] += mv-CEP
            for j in ZISellers:
                for mc in MarketMC(j,market):
                    if CEP>mc:
                        sellerProfits[j-1] += CEP-mc
            remainingZIBuyUnits = [(MarketMV(i,market)[2],i) for i in ZIBuyers]
            remainingZISellUnits = [(MarketMC(j,market)[2],j) for j in ZISellers]
            firstSniper = True
            for k in range(3):
                shuffle(SniperBuyers)
                for i in SniperBuyers:
                    mv = MarketMV(i,market)[k]
                    if len(remainingZISellUnits)>0:
                        (mc,j) = remainingZISellUnits[0]
                        if (mv>mc):
                            remainingZISellUnits.pop(0)
                            includedExtra += CEP-mc
                            profit = mv-mc
                            if firstSniper:
                                buyerProfits[i-1] += profit
                            else:
                                buyerProfits[i-1] += 0.8*profit
                                sellerProfits[j-1] += 0.2*profit
                            firstSniper = False
                        else:
                            if mv>CEP:
                                omittedInfra += -(mv-CEP)
                    else:
                        if mv>CEP:
                            omittedInfra += -(mv-CEP)
            firstSniper = True
            for k in range(3):
                shuffle(SniperSellers)
                for j in SniperSellers:
                    mc = MarketMC(j,market)[k]
                    if len(remainingZIBuyUnits)>0:
                        (mv,i) = remainingZIBuyUnits[0]
                        if (mv>mc):
                            remainingZIBuyUnits.pop(0)
                            includedExtra += mv-CEP
                            profit = mv-mc
                            if firstSniper:
                                sellerProfits[j-1] += profit
                            else:
                                buyerProfits[i-1] += 0.2*profit
                                sellerProfits[j-1] += 0.8*profit
                            firstSniper = False
                        else:
                            if mc<CEP:
                                omittedInfra += -(CEP-mc)
                    else:
                        if mc<CEP:
                            omittedInfra += -(CEP-mc)
        allprofits = []
        allprofits.extend(buyerProfits)
        allprofits.extend(sellerProfits)
        assert(len(allprofits)==40)
        sumprofits = sum(allprofits)
        sumdiff = sum((abs(allprofits[i]-allprofits[j]) for i in range(len(allprofits)) for j in range(i)))
        GiniCoefficient = (0.0+sumdiff)/((len(allprofits)-1.0)*sumprofits)
        out.extend((5*snipers,market,includedExtra/maxprofit,omittedInfra/maxprofit,(sum(buyerProfits)+sum(sellerProfits))/maxprofit, GiniCoefficient))
    print(",".join(map(str,out)))