basic_income_monte_carlo.py

from pylab import *
from scipy.stats import *

num_adults = 227e6
basic_income = 7.25*40*50
labor_force = 154e6
disabled_adults = 21e6
current_wealth_transfers = 3369e9

def jk_rowling(num_non_workers):
    num_of_jk_rowlings = binom(num_non_workers, 1e-7).rvs()
    return num_of_jk_rowlings * 1e9

def basic_income_cost_benefit():
    direct_costs = num_adults * basic_income / 2
    administrative_cost_per_person = norm(250,75)
    non_worker_multiplier = uniform(-0.10, 0.15).rvs()
    non_workers = (num_adults-labor_force-disabled_adults) * (1+non_worker_multiplier)
    marginal_worker_hourly_productivity = norm(10,1)

    administrative_costs = num_adults * administrative_cost_per_person.rvs()

    undisabled = 2e6
    undisabled_hourly_wage = uniform(0, 10).rvs()
    undisabled_cost_benefit = -1 * undisabled * (40*52*undisabled_hourly_wage)

    avg_hourly_wage = uniform(10, 30).rvs()
    productivity_multiplier = uniform(0.0, 0.2).rvs()
    productivity_cost_benefit = (-1 * labor_force * (40*52*avg_hourly_wage) *
                                 productivity_multiplier)

    return direct_costs + administrative_costs + undisabled_cost_benefit + productivity_cost_benefit - jk_rowling(non_workers)

def basic_job_cost_benefit():
    administrative_cost_per_disabled_person = norm(500,150).rvs()
    administrative_cost_per_worker = norm(5000, 1500).rvs()
    non_worker_multiplier = uniform(-0.20, 0.25).rvs()
    basic_job_hourly_productivity = uniform(-7.25, 7.25).rvs()

    disabled_cost = disabled_adults * (basic_income + administrative_cost_per_disabled_person)
    num_basic_workers = ((num_adults - disabled_adults - labor_force) *
                         (1+non_worker_multiplier)
                         )

    basic_worker_cost_benefit = num_basic_workers * (
        basic_income +
        administrative_cost_per_worker -
        40*50*basic_job_hourly_productivity
        )
    return disabled_cost + basic_worker_cost_benefit


N = 1024*32
bi = zeros(shape=(N,), dtype=float)
bj = zeros(shape=(N,), dtype=float)

for k in range(N):
    bi[k] = basic_income_cost_benefit()
    bj[k] = basic_job_cost_benefit()

subplot(211)
height=50*N/1024

title("Basic Income")
hist(bi, bins=50)
axis([-1e12,2.5e12,0,height])

subplot(212)
title("Basic Job")
hist(bj, bins=50)

axis([-1e12,2.5e12,0,height])

show()

See http://www.jeremyscheff.com/2013/11/basic-income-vs-basic-job/