Universal Basic Income Comparison Script for USA/AUS

income_distribution.py

# TODO: do the conversion from the data directly

# Data Format is a list of [weekly income, percentage of population from 0-100]
wealth_distributions = {
    # Graph 1 from https://www.abs.gov.au/statistics/economy/finance/government-benefits-taxes-and-household-income-australia/latest-release#income-redistribution-and-inequality
    "AUS": [[50.0, 1.9], [100.0, 1.8], [150.0, 2.2], [200.0, 2.1], [250.0, 2.1], [300.0, 1.8], [350.0, 2.1], [400.0, 2.3], [450.0, 2.3], [500.0, 2.9], [550.0, 2.2], [600.0, 2.6], [650.0, 2.4], [700.0, 3.3], [750.0, 2.3], [800.0, 2.4], [850.0, 2.6], [900.0, 2.7], [950.0, 2.6], [1000.0, 2.5], [1050.0, 2.8], [1100.0, 2.3], [1150.0, 2.3], [1200.0, 2.2], [1250.0, 2.1], [1300.0, 2.0], [1350.0, 1.7], [1400.0, 2.5], [1450.0, 2.0], [1500.0, 1.6], [1550.0, 1.6], [1600.0, 1.5], [1650.0, 1.4], [1700.0, 1.2], [1750.0, 1.3], [1800.0, 1.2], [1850.0, 1.1], [1900.0, 0.9], [1950.0, 0.9], [2000.0, 1.2], [2050.0, 0.8], [2100.0, 0.9], [2150.0, 0.9], [2200.0, 0.8], [2250.0, 0.8], [2300.0, 0.8], [2350.0, 0.3], [2400.0, 0.6], [2450.0, 0.7], [2500.0, 0.3], [2550.0, 0.4], [2600.0, 0.3], [2650.0, 0.4], [2700.0, 0.3], [2750.0, 0.5], [2800.0, 0.1]],
    # https://www.census.gov/data/tables/time-series/demo/income-poverty/cps-hinc/hinc-06.html, using the "mean" and "number" columns
    "USA": [[21.057692307692307, 2.945870409728223], [152.25, 2.1035258581093177], [238.71153846153845, 4.002304380658773], [333.8076923076923, 4.025659590038225], [428.3076923076923, 3.943916357210142], [520.2692307692307, 4.08482612046617], [618.5576923076923, 4.248312586122335], [713.9807692307693, 3.660539816739457], [807.3269230769231, 4.186032027777129], [905.7307692307693, 3.8707367011545255], [1000.0, 3.8030065939541147], [1099.0, 3.3974044577309632], [1192.0192307692307, 3.5274151232765805], [1291.2884615384614, 2.8454430093965795], [1386.5192307692307, 2.96377607025247], [1481.0576923076924, 2.762921269589182], [1576.7884615384614, 2.6484807436298667], [1674.4615384615386, 2.3456415286763046], [1771.076923076923, 2.3534265984694556], [1870.8076923076924, 2.1930541607305507], [1963.2115384615386, 2.220301905006578], [2060.6923076923076, 1.9392608854738382], [2153.4423076923076, 1.8785373410872628], [2251.9423076923076, 1.5297662143541118], [2346.0576923076924, 1.8629672015009615], [2444.2884615384614, 1.3896349580773992], [2540.2884615384614, 1.3841854092221937], [2638.903846153846, 1.23159804127644], [2732.4423076923076, 1.0899097710410972], [2829.75, 1.020622649882056], [2922.3653846153848, 1.3398105114012346], [3024.4615384615386, 0.940436431012604], [3114.596153846154, 0.9809187939369876], [3215.096153846154, 0.7847350351495901], [3311.403846153846, 0.9739122311231521], [3406.7884615384614, 0.7761714583771243], [3499.3846153846152, 0.7068843372180831], [3600.3076923076924, 0.638375723038357], [3693.5384615384614, 0.6446037788728776], [3791.846153846154, 0.50369401561685], [4250.961538461538, 4.168126367252882], [8016.75, 6.0848105503265835]],
}

ubi_cost.py

import matplotlib.pyplot as plt
import numpy as np

from income_distribution import wealth_distributions


# AUS values in AUD, USA in USD
tax_brackets = {
    "AUS": {
        18_200: .19,
        45_000: .325,
        120_000: .37,
        180_000: .45,
    },
    "USA": {
        0: 0.1,
        9_950: 0.12,
        40_525: 0.22,
        86_375: 0.24,
        164_925: 0.32,
        209_425: 0.35,
        523_600: 0.37,
    },
}

max_tax_brackets = {country: max(brackets.keys()) for (country, brackets) in tax_brackets.items()}

poverty_levels = {
    "AUS": 457 * 52,
    "USA": 13300,
}

bar_thickness = {
    "AUS": 0.8,
    "USA": 0.4,
}

assert poverty_levels.keys() == tax_brackets.keys() == wealth_distributions.keys() == bar_thickness.keys()


# from https://stackoverflow.com/a/29677616
def weighted_quantile(values, quantiles, sample_weight=None, 
                      values_sorted=False, old_style=False):
    """ Very close to numpy.percentile, but supports weights.
    NOTE: quantiles should be in [0, 1]!
    :param values: numpy.array with data
    :param quantiles: array-like with many quantiles needed
    :param sample_weight: array-like of the same length as `array`
    :param values_sorted: bool, if True, then will avoid sorting of
        initial array
    :param old_style: if True, will correct output to be consistent
        with numpy.percentile.
    :return: numpy.array with computed quantiles.
    """
    values = np.array(values)
    quantiles = np.array(quantiles)
    if sample_weight is None:
        sample_weight = np.ones(len(values))
    sample_weight = np.array(sample_weight)
    assert np.all(quantiles >= 0) and np.all(quantiles <= 1), \
        'quantiles should be in [0, 1]'

    if not values_sorted:
        sorter = np.argsort(values)
        values = values[sorter]
        sample_weight = sample_weight[sorter]

    weighted_quantiles = np.cumsum(sample_weight) - 0.5 * sample_weight
    if old_style:
        # To be convenient with numpy.percentile
        weighted_quantiles -= weighted_quantiles[0]
        weighted_quantiles /= weighted_quantiles[-1]
    else:
        weighted_quantiles /= np.sum(sample_weight)
    return np.interp(quantiles, weighted_quantiles, values)


def get_tax(income: int, country: str) -> int:
    tax = 0
    previous_bracket = max_tax_brackets[country]

    for (bracket, rate) in sorted(tax_brackets[country].items(), reverse=True):
        if income < bracket: # e.g. 46K < 120K or 180K brackets
            continue # move to next lowest bracket

        elif bracket <= income < previous_bracket:
            # e.g. 45K < 46K <= 120K ( between )
            tax += (income - bracket) * rate

        else:
            if bracket == max_tax_brackets[country]:
                tax += (income - bracket) * rate
            else:
                # i.e. income > bracket, e.g. after the above, 18K < 46K
                # add full tax from this lower bracket
                tax += (previous_bracket - bracket) * rate

        previous_bracket = bracket

    return tax


def get_pc_income(percent: int, country: str) -> float:
    # data is a list of (weekly, % people as [0, 100])
    # we want the sum incomes
    # percent arg is [0, 1]

    weekly_sum = 0
    for (weekly, pc) in wealth_distributions[country]:
        weekly_sum += weekly * (pc/100)

    yearly_sum = 52 * weekly_sum

    return percent * yearly_sum


def adjust_tax_for_ubi(taxes: list, incomes: list, percent: int, country: str) -> list:
    # takes in a list of income
    # applies UBI tax and adds UBI result

    # everyone gets this "40% of everyone" but removes "40% of theirs"
    income_pc = get_pc_income(percent, country)

    # "negative" tax = get, positive tax = remove
    return [ tax - income_pc + percent*inc for (tax, inc) in zip(taxes, incomes)]


def adjust_income_for_ubi(taxes: list, incomes: list, percent: int, country: str) -> list:
    return [inc - ubi_and_tax for (inc, ubi_and_tax) in zip(incomes, adjust_tax_for_ubi(taxes, incomes, percent, country))]


def test_ubi(country: str):
    # Base Income
    incomes = range(0, max_tax_brackets[country], 10) # generic "range", not real distribution
    base_taxes = [get_tax(i, country) for i in incomes]

    fig, ((tax_v_income, net_v_gross), (wealth_distribution, _)) = plt.subplots(nrows=2, ncols=2)

    fig.tight_layout() # prevent overlap

    wealth_distribution.set_title("Income Distribution")
    wealth_distribution.set_xlabel("Private Income ($)")
    wealth_distribution.set_ylabel("% People with Income")
    dis_incomes, dis_percentages = zip(*wealth_distributions[country])
    dis_incomes = 52 * np.array(dis_incomes) # convert weekly to yearly
    
    wealth_distribution.bar(dis_incomes, dis_percentages, width=bar_thickness[country] * max_tax_brackets[country] / len(dis_incomes))
    wealth_distribution.set_xlim(right=max_tax_brackets[country])

    # Display the wealth distribution percentiles
    wanted_pciles = [0.25, 0.5, 0.75, 0.9, 0.95, .99]
    c = iter(["r", "pink", "chartreuse", "cyan", "yellow", "orange"])
    percentiles = weighted_quantile(dis_incomes, wanted_pciles, sample_weight=dis_percentages)
    for pciles, pcile_value in zip(wanted_pciles, percentiles):
        wealth_distribution.axvline(x=pcile_value, color=next(c), label=f"{pciles*100:.0f}% PC")

    wealth_distribution.legend(loc="upper right")

    tax_v_income.set_title("Tax vs Income")
    
    tax_v_income.set_xlabel("Private Income ($)")
    tax_v_income.set_ylabel("Tax ($)")
    
    net_v_gross.set_title("Private Income vs Post-Tax Income")
    
    net_v_gross.set_xlabel("Private Income ($)")
    net_v_gross.set_ylabel("Post-Tax Income ($)")

    net_v_gross.axhline(y=poverty_levels[country], color="black", label="Poverty Line = ${}".format(poverty_levels[country]))

    tax_v_income.plot(incomes, base_taxes, label="Base Tax")
    net_v_gross.plot(incomes, adjust_income_for_ubi(base_taxes, incomes, 0, country), label="Base Tax")

    for percent in [0.2, 0.3, 0.4, 0.5, 0.6]:
        tax_v_income.plot(incomes, adjust_tax_for_ubi(base_taxes, incomes, percent, country), label=f"{percent*100}% UBI")
        net_v_gross.plot(incomes, adjust_income_for_ubi(base_taxes, incomes, percent, country), label=f"{percent*100}% UBI")

    tax_v_income.legend(loc="upper left")
    net_v_gross.legend(loc="upper left")

    tax_v_income.grid()
    net_v_gross.grid()
    plt.show()

if __name__ == "__main__":
    print("Choose country from:", ", ".join(tax_brackets.keys()), end=": ")
    test_ubi(input().upper())

# TODO: Display value of meeting points