Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como simular uma rede de firmas que participam de licitações públicas lançadas por entidades emissoras por intermédio de simulações Monte Carlo?

0 votos
16 visitas
perguntada Nov 1 em Sistemas Complexos por bonfim_tiago (6 pontos)  

O presente tópico possui como objetivo o detalhamento do código elaborado para simular o algoritmo descrito no paper intitulado "A network approach to cartel
detection in public auction markets", dos autores Johannes Wachs & János Kertész, por intermédio de simulações Monte Carlo.

Compartilhe

2 Respostas

0 votos
respondida Nov 2 por bonfim_tiago (6 pontos)  

O código elaborado necessita das seguintes bibliotecas para funcionar corretamente:

import numpy as np
import random
import matplotlib.pyplot as plt
from itertools import combinations

A biblioteca numpy possui funcionalidades matemáticas bastante úteis, como seria o caso de gerar uma distribuição normal. Isso foi necessário para distribuir as firmas e os contratantes no terreno de forma aleatória. A seguir, tem-se a biblioteca random que foi utilizada para gerar números aleatórios no intervalo entre 0 - 1, por diversas razões.

A seguir, observa-se o uso da biblioteca matplotlib que foi empregada para gerar figuras ao término das simulações. Por fim, observa-se a utilização da biblioteca itertools que foi utilizada para gerar combinações distintas de conjuntos. Isso foi relevante na hora de se observar quais combinações de firmas participaram de cada contrato no passado.

Todas os parâmetros que foram empregados no código foram adicionados a um dicionário com nomes definidos em strings, de forma a facilitar seu reconhecimento, modificação e manutenção futura do código, conforme fragmento:

# initial string parameters
STR_FIRMS_NUMBER = "firms_number"
FIRMS_INIT_VECTOR = "firms_init_vector"
ISSUERS_NUMBER = "issuers_number"
ROUNDS_NUMBER = "rounds_number"
TOTAL_ITERATIONS = "total_iterations"
PARTICIPANTS_RADIUS = "participants_radius"
POSITIONING_ISSUER_STD = "positioning_issuer_std"
TERRAIN_SIZE = "terrain_size"
FAMILIAR_FIRM_THRESHOLD = "familiar_firm_threshold"
COLLUSION_FIRM_THRESHOLD = "collusion_firm_threshold"
F_FREQUENCY_K = "f_frequency_k"
# initial parameters value
firms_number = 50
parameter_map = {STR_FIRMS_NUMBER: firms_number, ISSUERS_NUMBER: 75, ROUNDS_NUMBER: 2000, PARTICIPANTS_RADIUS: 0.1,
                 POSITIONING_ISSUER_STD: 0.3, TERRAIN_SIZE: 1, FAMILIAR_FIRM_THRESHOLD: float(2 / 3),
                 COLLUSION_FIRM_THRESHOLD: 0.9, FIRMS_INIT_VECTOR: [0] * firms_number, F_FREQUENCY_K: 10,
                 TOTAL_ITERATIONS: 250}

Depois, passa-se a declaração de algumas classes. A modelagem empregada fez uso de uma classe Firms para modelar as empresas que fazem propostas pelos contratos:

class Firms:

    def __init__(self, id: int, x_coord: float, y_coord: float):
        self.id = id
        self.x_coord = x_coord
        self.y_coord = y_coord

Uma classe Issuers para modelar as entidades contratantes que lançam as licitações no terreno, conforme:

class Issuers:

    def __init__(self, id: int, x_coord: float, y_coord: float):
        self.id = id
        self.x_coord = x_coord
        self.y_coord = y_coord

E uma classe Contract para modelar os contratos propriamente ditos. Essa última classe vai ser detalhada em etapas. Começaremos pela declaração da classe e seu construtor:

class Contract:

    def __init__(self, id: int, issuer_id: int):
        self.id = id
        self.issuer_id = issuer_id
        self.x_coord = 0
        self.y_coord = 0
        self.participants = list()

Esse construtor possui como parâmetro o id do contrato para referência/localização futura, o id da entidade contratante responsável pelo contrato, elemento da classe Issuers, as coordenadas em que o contrato está sendo lançado no terreno e, por fim, uma lista que contém todas as firmas que vão participar dessa licitação/contrato público.

A seguir, são definidas duas funções "getters" que seriam para pegar as coordenadas do contrato de uma distribuição aleatória centralizada na entidade emissora da licitação, mas segundo uma distribuição normal com desvio padrão correspondente ao parâmetro POSITIONINGISSUERSTD, conforme:

def gen_contract_coord(self):

for issuer in issuers:
    if issuer.id == self.issuer_id:
        self.x_coord = min(np.random.normal(issuer.x_coord, parameter_map[POSITIONING_ISSUER_STD] ** 2), 1)
        self.y_coord = min(np.random.normal(issuer.y_coord, parameter_map[POSITIONING_ISSUER_STD] ** 2), 1)
        break

O próximo método, devo confessar, acabou ficando bastante longo. Idealmente, esse método deveria ser quebrado em várias funcionalidades menores. Contudo, acabei não tendo o tempo necessário para efetuar essa implementação da melhor maneira possível, correndo o risco de pagar o preço por isso no futuro quando de possíveis etapas de manutenção do código. Vamos analisar esse método por etapas:

O método funciona para buscar os participantes de determinado contrato, até aí, tudo certo, mas ele também acaba calculando os parâmetros relevantes da simulação para essas firmas, conforme:

def get_participants(self):

for firm in firms:
    if ((firm.x_coord - self.x_coord) ** 2 + (firm.y_coord - self.y_coord) ** 2) ** 0.5 <= \
            parameter_map[PARTICIPANTS_RADIUS]:
        # This firm is inside the range of the contract
        if parameter_map[FIRMS_INIT_VECTOR][firm.id] != 0:
            self.participants.append({"firm_id": firm.id})
        else:
            if random.random() > 1 / 2:
                self.participants.append({"firm_id": firm.id, "collude": True})
            else:
                self.participants.append({"firm_id": firm.id, "collude": False})

            # firm is initialized with random choice
            parameter_map[FIRMS_INIT_VECTOR][firm.id] = 1

Esse fragmento avalia quais firmas podem participar do referido contrato, segundo um critério de distância entre os elementos com base na distância definida por PARTICIPANTS_RADIUS. A seguir, essa firma vai ser adicionada como se fosse participante do leilão e vai ter seu parâmetro de propensão a participar de um cartel inicializado aleatoriamente.

O próximo bloco é responsável por calcular o f_memory de um contrato, que consiste em um levantamento, para todos os nós que estão participante de determinada licitação se, relativo aos demais participantes, houve cooperação deles em coluir na última vez em que eles participaram conjuntamente de uma licitação, podendo ter ocorrido em qualquer contrato anterior.

# Eval f_memory
participant_id_list = list()
for participant in self.participants:
    participant_id_list.append(participant["firm_id"])

for participant_a in participant_id_list:
    f_memory = 0
    for participant_b in participant_id_list:
        if participant_a != participant_b:
            # check if participant_a and participant_b have colluded in the past
            contract_iterator = self.id
            while contract_iterator:
                contract_iterator -= 1
                previous_participants = contracts[contract_iterator].participants
                if len(previous_participants):
                    previous_participant_id_list = list()
                    for participant in previous_participants:
                        previous_participant_id_list.append(participant["firm_id"])
                    if participant_a in previous_participant_id_list:
                        if participant_b in previous_participant_id_list:
                            # last run where both participants met
                            # check if b colluded with a
                            for participant in previous_participants:
                                if participant["firm_id"] == participant_b and "collude" in participant:
                                    if participant["collude"]:
                                        # participant b have colluded with participant a in the last run
                                        f_memory += 1 / (len(self.participants) - 1)

                                        # update contract object
                                        index = 0
                                        for firm_dict in self.participants:
                                            # add f_memory parameter to participant_a
                                            if firm_dict["firm_id"] == participant_a:
                                                firm_dict["f_memory"] = f_memory
                                                # no need to keep searching for the last time this group met
                                                contract_iterator = 0
                                                break
                                            else:
                                                index += 1

O próximo bloco é o responsável pelo cálculo do outro parâmetro relevante no algoritmo, que seria o f_frequency. Ele é responsável por estimar se os demais nós participantes do leilão, relativo ao nó remanescente, cooperaram no passado um número expressivo de vezes, analisando os últimos k procedimentos concorrenciais.

   # Eval f_frequency
    for participant_a in participant_id_list:
        other_participants = participant_id_list
        other_participants.remove(participant_a)
    # check if the list is not empty
    if len(other_participants) == 0:
        continue

    f_frequency = 0
    participant_a_counter = 0

    #  get up to k last times participant_a entered an auction
    contract_iterator = self.id
    while contract_iterator:
        contract_iterator -= 1
        previous_participants = contracts[contract_iterator].participants
        # make a list of ids with previous participants
        previous_participant_id_list = list()
        for participant in previous_participants:
            previous_participant_id_list.append(participant["firm_id"])

        # check if pacticipant_a was in this auction
        if participant_a in previous_participant_id_list:
            participant_a_counter += 1
            # check if other_participants were in this auction
            if set(other_participants).issubset(set(previous_participant_id_list)):
                f_frequency += 1
            # check if we can iterate another time
            if contract_iterator:
                if (self.id - contract_iterator) == parameter_map[F_FREQUENCY_K]:
                    # finished the k trials, lets register f_frequency
                    for firm_dict in self.participants:
                        # add f_frequency parameter to participant_a
                        if firm_dict["firm_id"] == participant_a:
                            firm_dict["f_frequency"] = 1 if (f_frequency/participant_a_counter) >= 2/3 else 0
                    break
                else:
                    continue
            else:
                # this was the last loop
                # finished the n possible trials, lets register f_frequency
                for firm_dict in self.participants:
                    # add f_frequency parameter to participant_a
                    if firm_dict["firm_id"] == participant_a:
                        firm_dict["f_frequency"] = 1 if (f_frequency / participant_a_counter) >= 2 / 3 else 0
                break
        else:
            # participant_a was not in the last auction
            if contract_iterator == 0:
                for firm_dict in self.participants:
                    # add f_frequency parameter to participant_a
                    if firm_dict["firm_id"] == participant_a and participant_a_counter > 0:
                        firm_dict["f_frequency"] = 1 if (f_frequency / participant_a_counter) >= 2 / 3 else 0
            continue

Com base no cálculo dos últimos parâmetros, checamos agora se os nós participantes do contrato estão dispostos a conluiar no presente contrato e, com base nisso, atualizados os metadados de cada objeto.

# Check if Nodes will collude in this run
for firm_dict in self.participants:
    # check if firm has both f_frequency and f_memory defined
    if "f_frequency" in firm_dict and "f_memory" in firm_dict:
        if firm_dict["f_frequency"] * firm_dict["f_memory"] >= 0.9:
            firm_dict["collude"] = True

    # Now, adds the random possibility of collusion
    if random.random() <= 0.01:
        firm_dict["collude"] = True

# Check if a Cartel was formed this run
# Eager to collude participants counter
collusion_counter = 0
for firm_dict in self.participants:
    # Check if the firm has collude parameter defined
    if "collude" in firm_dict and len(self.participants) > 1:
        if firm_dict["collude"]:
            collusion_counter += 1
            if collusion_counter == len(self.participants):
                # A Cartel was formed
                cartel_list.append({"Contract_run": self.id, "Cartel": True})

Externamente às classes do programa, temos um método responsável por gerar os gráficos necessários para visualização do terreno e dos elementos da rede.

def plot():
    # Plot firms on the map
    all_firms_x_coord_list = list()
    all_firms_y_coord_list = list()
    for firm in firms:
        all_firms_x_coord_list.append(firm.x_coord)
        all_firms_y_coord_list.append(firm.y_coord)

    plt.scatter(all_firms_x_coord_list, all_firms_y_coord_list, color='g')
    # Loop to write the number of the firm by its side
    for i in range(len(firms)):
        plt.annotate(str(i), (all_firms_x_coord_list[i], all_firms_y_coord_list[i]))

    # Plot issuers on the map
    all_issuers_x_coord_list = list()
    all_issuers_y_coord_list = list()
    for issuer in issuers:
        all_issuers_x_coord_list.append(issuer.x_coord)
        all_issuers_y_coord_list.append(issuer.y_coord)

    plt.scatter(all_issuers_x_coord_list, all_issuers_y_coord_list, color='b')
    # Loop to write the number of the issuer by its side
    for i in range(len(issuers)):
        plt.annotate(str(i), (all_issuers_x_coord_list[i], all_issuers_y_coord_list[i]))

    # Replot Cartel Nodes on the map in red
    all_cartel_nodes_x_coord_list = list()
    all_cartel_nodes_y_coord_list = list()
    for cartel in cartel_list:
        nodes_involved = contracts[cartel["Contract_run"]].participants

        for node in nodes_involved:
            firm_id = node["firm_id"]
            all_cartel_nodes_x_coord_list.append(firms[firm_id].x_coord)
            all_cartel_nodes_y_coord_list.append(firms[firm_id].y_coord)

    plt.scatter(all_cartel_nodes_x_coord_list, all_cartel_nodes_y_coord_list, color='r')

    # Display Graph
    plt.show(block=False)
    plt.show()

Nesse momento, temos então o laço principal do programa, que é responsável pela geração dos objetos no número necessário de vezes. Aqui reside a essência da modelagem de Monte Carlo. Vamos ter um duplo laço de repetição for, pois externamente temos o número de iterações a serem realizadas, com mapas, firmas e issuers distintos e internamente, temos o número de licitações a serem realizados no mapa atualmente gerado.

final_cartel_list = list()
# Outer loop with the map, firms and issuers iterations
for i in range(parameter_map[TOTAL_ITERATIONS]):
    print("Iteration Number: ", i, "/", parameter_map[TOTAL_ITERATIONS])
    # Create the initial firms and issuers
    firms = [Firms(id=i, x_coord=random.random(), y_coord=random.random()) for i in range(parameter_map[STR_FIRMS_NUMBER])]
    issuers = [Issuers(id=i, x_coord=random.random(), y_coord=random.random()) for i in
               range(parameter_map[ISSUERS_NUMBER])]

    cartel_list = list()

    # Main Routine to create contracts
    contracts = list()
    for j in range(parameter_map[ROUNDS_NUMBER]):
        # print("run number: ", j, "/", parameter_map[ROUNDS_NUMBER])
        # create object contract
        contracts.append(Contract(id=j, issuer_id=random.randint(0, parameter_map[ISSUERS_NUMBER] + 1)))
        # generate contract coordinates
        contracts[j].gen_contract_coord()
        contracts[j].get_participants()

    plot()

Nesse ponto, vamos calcular os pesos (para todas as combinações de pares de nós) dentre os elementos que estão participando do atual contrato.

  # Time to calculate weights of the connections
    weights = list()
    neighbor_weights = list()

    for cartel in cartel_list:
        participants = contracts[cartel["Contract_run"]].participants

        # convert participants_id into a list
        nodes_list = list()
        for participant in participants:
            nodes_list.append(participant["firm_id"])

        # get the distinct combinations of 2 nodes in the list of participants [(node_a, node_b), ...]
        nodes_comb = [tuple(comb) for comb in combinations(nodes_list, 2)]

        # the weight is calculated for every combination of two nodes
        for nodes_tuple in nodes_comb:
            # first, lets run through all previous contracts
            if cartel["Contract_run"]:
                contract_counter = cartel["Contract_run"]
            else:
                continue

            contract_with_a_counter = 0
            contract_with_b_counter = 0
            contract_with_a_and_b_counter = 0

            while contract_counter:
                contract_participants = contracts[contract_counter].participants

                # add all contract participants in a list
                contract_participants_list = list()
                for firm in contract_participants:
                    contract_participants_list.append(firm["firm_id"])

                # check if node a participated
                if nodes_tuple[0] in contract_participants_list:
                    contract_with_a_counter += 1

                # check if node b participated
                if nodes_tuple[1] in contract_participants_list:
                    contract_with_b_counter += 1

                #  check if both nodes participated
                if nodes_tuple[0] in contract_participants_list:
                    if nodes_tuple[1] in contract_participants_list:
                        contract_with_a_and_b_counter += 1

                contract_counter -= 1

            tuple_union = contract_with_a_counter + contract_with_b_counter - contract_with_a_and_b_counter
            tuple_intersection = contract_with_a_and_b_counter

            # Register in a vector the results for the weight parameter
            tuple_weight = tuple_intersection/tuple_union
            weight_dict = {"Cartel": cartel["Contract_run"], "node_A": nodes_tuple[0], "node_B": nodes_tuple[1],
                           "weight": tuple_weight}
            weights.append(weight_dict)
0 votos
respondida Nov 2 por bonfim_tiago (6 pontos)  

Devido ao limite de 20.000 caracteres do Prorum, vou continuar com a resposta da referida questão no presente post.

Podemos calcular então o fator de coerência da rede:

 # Calculate de Coherence Factor
    cartel_weight_sum = 0
    cartel_weight_product = 1
    tuples_size = len(nodes_comb)
for weight in weights:
    if weight["Cartel"] == cartel["Contract_run"]:
        cartel_weight_sum += weight["weight"]
        cartel_weight_product = cartel_weight_product*weight["weight"]

cartel_weight_sum = cartel_weight_sum/tuples_size
cartel_weight_product = cartel_weight_product**(1/tuples_size)

if cartel_weight_sum:
    cartel_coherence = cartel_weight_product/cartel_weight_sum
else:
    cartel_coherence = 1

for group in cartel_list:
    if group["Contract_run"] == cartel["Contract_run"]:
        group["Coherence"] = cartel_coherence

# Time to calculate the weights with the neighborhood
for node in nodes_list:

    node_neighborhood = list()
    # First, gather the neighborhood of the node
    for firm in firms:
        if ((firm.x_coord - firms[node].x_coord)**2 + (firm.y_coord - firms[node].y_coord)**2)**0.5 <= 2*parameter_map[PARTICIPANTS_RADIUS]:
            # firm is near this node
            node_neighborhood.append(firm.id)

    # Now, we must exclude the neighbors from the cartel
    static_node_neighborhood = node_neighborhood.copy()
    for neighbor in static_node_neighborhood:
        if neighbor in nodes_list:
            # this neighbor is part of the cartel
            node_neighborhood.remove(neighbor)

    # tuples should be with node, neighbor for each neighbor...
    # Now, we have to evaluate the weights of tuples in node_neighborhood

    neighbor_comb = [(node, neighbor) for neighbor in node_neighborhood]

    # the weight is calculated for every combination of two nodes
    for neighbor_tuple in neighbor_comb:
        # first, lets run through all previous contracts
        contract_counter = cartel["Contract_run"]

        contract_with_na_counter = 0
        contract_with_nb_counter = 0
        contract_with_na_and_nb_counter = 0

        while contract_counter:
            neighbor_contract_participants = contracts[contract_counter].participants

            # add all contract participants in a list
            neighbor_contract_participants_list = list()
            for neighbor in neighbor_contract_participants:
                neighbor_contract_participants_list.append(neighbor["firm_id"])

            # check if node a participated
            if neighbor_tuple[0] in neighbor_contract_participants_list:
                contract_with_na_counter += 1

            # check if node b participated
            if neighbor_tuple[1] in neighbor_contract_participants_list:
                contract_with_nb_counter += 1

            #  check if both nodes participated
            if neighbor_tuple[0] in neighbor_contract_participants_list:
                if neighbor_tuple[1] in neighbor_contract_participants_list:
                    contract_with_na_and_nb_counter += 1

            contract_counter -= 1

        neighbor_tuple_union = contract_with_na_counter + contract_with_nb_counter - contract_with_na_and_nb_counter
        neighbor_tuple_intersection = contract_with_na_and_nb_counter

        # Register in a vector the results for the weight parameter
        tuple_weight = neighbor_tuple_intersection / neighbor_tuple_union
        neighbor_weight_dict = {"Cartel": cartel["Contract_run"], "node": neighbor_tuple[0],
                                "neighbor": neighbor_tuple[1], "weight": tuple_weight}
        neighbor_weights.append(neighbor_weight_dict)

Seguimos então com o cálculo do fator de exclusividade para os participantes do atual contrato.

# Calculate de Exclusivity Factor
    s_in = 0
    s_out = 0
    s_in_size = 0
    s_out_size = 0

for weight in weights:
    if weight["Cartel"] == cartel["Contract_run"]:
        s_in += weight["weight"]
        s_in_size += 1

for weight in neighbor_weights:
    if weight["Cartel"] == cartel["Contract_run"]:
        s_out += weight["weight"]
        s_out_size += 1

s_in = s_in/s_in_size
if s_out:
    s_out = s_out/s_out_size
cartel_exclusivity = s_in/(s_in + s_out)

for group in cartel_list:
    if group["Contract_run"] == cartel["Contract_run"]:
        group["Exclusivity"] = cartel_exclusivity

for cartel in cartel_list:
    final_cartel_list.append(cartel)

Finalmente, fora dos laços de repetição principais, vamos plotar mais alguns gráficos de interesse.

# finally, lets plot cartel list by its two parameters
coord_x = list()
coord_y = list()
for cartel in final_cartel_list:
    coord_x.append(cartel["Coherence"])
    coord_y.append(cartel["Exclusivity"])

plt.scatter(coord_x, coord_y, color='b')
plt.title("Iterations: 250; Rounds: 2000")
plt.xlabel("Coherence")
plt.ylabel("Exclusivity")
plt.xlim(0, 1.05)
plt.ylim(0, 1.05)
plt.show()
...