GigaProjects

← Back to WealthDist

simulation.js

export class SimulationUtils {
    static normalQuantile(p) {
        const a = [-39.69683028665376, 220.9460984245205, -275.9285104469687, 138.357751867269, -30.66479806614716, 2.506628277459239];
        const b = [-54.47609879822406, 161.5858368580409, -155.6989798598866, 66.80131188771972, -13.28068155288572];
        const c = [-0.007784894002430293, -0.3223964580411365, -2.400758277161838, -2.549732539343734, 4.374664141464968, 2.938163982698783];
        const d = [0.007784695709041462, 0.3224671290700398, 2.445134137142996, 3.754408661907416];
        const pLow = 0.02425;
        const pHigh = 1 - pLow;

        if (p <= 0) return -Infinity;
        if (p >= 1) return Infinity;

        let q;
        let r;
        if (p < pLow) {
            q = Math.sqrt(-2 * Math.log(p));
            return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
                ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1);
        }
        if (p > pHigh) {
            q = Math.sqrt(-2 * Math.log(1 - p));
            return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
                ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1);
        }

        q = p - 0.5;
        r = q * q;
        return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q /
            (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1);
    }

    static generateLogNormal(n, mean, sigma) {
        const mu = Math.log(mean) - 0.5 * Math.pow(sigma, 2);
        const results = new Float32Array(n);
        for (let i = 0; i < n; i++) {
            const p = (i + 0.5) / n;
            const z = SimulationUtils.normalQuantile(p);
            results[i] = Math.exp(mu + sigma * z);
        }
        return results;
    }

    static generatePareto(n, shape, scale) {
        const results = new Float32Array(n);
        for (let i = 0; i < n; i++) {
            const p = (i + 0.5) / n;
            results[i] = scale / Math.pow(1 - p, 1 / shape);
        }
        return results;
    }

    static gini(array) {
        if (array.length === 0) return 0;
        let sum = 0;
        let n = array.length;
        // Gini formula for sorted array
        // G = (2 * sum(i * xi) - (n + 1) * sum(xi)) / (n * sum(xi))
        // Note: i is 1-based index
        let numerator = 0;
        let total = 0;
        for (let i = 0; i < n; i++) {
            total += array[i];
            numerator += (i + 1) * array[i];
        }
        if (total === 0) return 0;
        return (2 * numerator) / (n * total) - (n + 1) / n;
    }
}

export class TaxSystem {
    constructor(rate = 0.2, deduction = 10000, type = 'flat', brackets = []) {
        this.rate = rate;
        this.deduction = deduction;
        this.type = type; // 'flat', 'custom', 'us2024', 'eu25'

        // Hardcoded Presets for simplicity in JS port
        if (type === 'us2024') {
            // 2024 Single Filer Brackets
            this.brackets = [
                { threshold: 11600, rate: 0.10 },
                { threshold: 47150, rate: 0.12 },
                { threshold: 100525, rate: 0.22 },
                { threshold: 191950, rate: 0.24 },
                { threshold: 243725, rate: 0.32 },
                { threshold: 609350, rate: 0.35 },
                { threshold: Infinity, rate: 0.37 }
            ];
            this.deduction = 14600; // Standard Deduction 2024
        } else if (type === 'eu25') {
            this.rate = 0.25;
            this.deduction = 0; // Flat tax often has no deduction in simple models
            this.type = 'flat';
        } else if (type === 'custom-progressive') {
            this.brackets = brackets.sort((a, b) => a.threshold - b.threshold);
            this.deduction = deduction; // Allow deduction for custom brackets too? Or assume 0?
            // Let's assume standard deduction applies first, then brackets.
        } else {
            this.brackets = brackets;
        }
    }

    calculateTax(incomes) {
        const taxes = new Float32Array(incomes.length);

        if (this.type === 'us2024' || this.type === 'custom-progressive') {
            // Progressive Logic
            for (let i = 0; i < incomes.length; i++) {
                let income = Math.max(0, incomes[i] - this.deduction);
                let tax = 0;
                let previous_threshold = 0;

                for (let b of this.brackets) {
                    if (income > previous_threshold) {
                        let taxable_in_bracket = Math.min(income, b.threshold) - previous_threshold;
                        tax += taxable_in_bracket * b.rate;
                        previous_threshold = b.threshold;
                    } else {
                        break;
                    }
                }
                taxes[i] = tax;
            }
        } else {
            // Flat Logic
            for (let i = 0; i < incomes.length; i++) {
                let taxable = Math.max(0, incomes[i] - this.deduction);
                taxes[i] = taxable * this.rate;
            }
        }
        return taxes;
    }
}

export class MultiYearSimulation {
    constructor({
        n_people = 1000,
        n_years = 50,
        distribution_type = 'lognormal',
        dist_params = { mean: 50000, sigma: 0.75 },
        initial_wealth_multiplier = 6.0,
        net_wealth_yield = 0.03,
        tax_capital_gains = false,
        policyA = { type: 'custom', rate: 0.5, deduction: 10000, brackets: [], ubi_pct: 1.0 },
        policyB = { type: 'custom', rate: 0.2, deduction: 10000, brackets: [], ubi_pct: 1.0, capital_tax_rate: 0.0 },
        labor_displacement = 0.0,
        capital_efficiency_boost = 0.0
    }) {
        this.n_people = n_people;
        this.n_years = n_years;
        this.net_wealth_yield = net_wealth_yield;
        this.initial_wealth_multiplier = initial_wealth_multiplier;
        this.tax_capital_gains = tax_capital_gains;
        this.labor_displacement = labor_displacement;
        this.capital_efficiency_boost = capital_efficiency_boost;

        // Generate Distribution
        if (distribution_type === 'pareto') {
            this.base_labor_income = SimulationUtils.generatePareto(n_people, dist_params.shape, dist_params.scale);
        } else {
            this.base_labor_income = SimulationUtils.generateLogNormal(n_people, dist_params.mean, dist_params.sigma);
        }

        this.tax_system_a = new TaxSystem(policyA.rate, policyA.deduction, policyA.type, policyA.brackets);
        this.tax_system_b = new TaxSystem(policyB.rate, policyB.deduction, policyB.type, policyB.brackets);
        this.ubi_pct_a = policyA.ubi_pct;
        this.ubi_pct_b = policyB.ubi_pct;
        this.capital_tax_rate_b = policyB.capital_tax_rate ?? 0.0;
    }

    run() {
        const income_sorted = Float32Array.from(this.base_labor_income).sort();
        function get_rank(val) {
            let low = 0, high = income_sorted.length;
            while (low < high) {
                let mid = (low + high) >>> 1;
                if (income_sorted[mid] < val) low = mid + 1;
                else high = mid;
            }
            return Math.min(low / income_sorted.length, 1.0);
        }
        function get_savings_rate(flow) {
            let r = get_rank(flow);
            return 0.02 + Math.pow(r, 3) * 0.48;
        }

        // Initialize Wealth
        let w_pre = new Float32Array(this.n_people);
        let w_a = new Float32Array(this.n_people);
        let w_b = new Float32Array(this.n_people);

        for (let i = 0; i < this.n_people; i++) {
            let w = this.base_labor_income[i] * this.initial_wealth_multiplier;
            w_pre[i] = w_a[i] = w_b[i] = w;
        }

        this.gini_history = { pre: [], a: [], b: [] };
        this.dist_history = { pre: [], a: [], b: [] };
        this.history_employment = [];
        this.revenue_history = { a: [], b: [] };
        this.ubi_history = { a: [], b: [] };

        // Dynamic State
        let labor_participation = new Float32Array(this.n_people);
        labor_participation.fill(1);
        let current_labor_income = new Float32Array(this.base_labor_income);
        let current_yield = this.net_wealth_yield;

        for (let year = 0; year <= this.n_years; year++) {
            this.gini_history.pre.push(SimulationUtils.gini(w_pre));
            this.gini_history.a.push(SimulationUtils.gini(w_a));
            this.gini_history.b.push(SimulationUtils.gini(w_b));

            this.dist_history.pre.push(new Float32Array(w_pre));
            this.dist_history.a.push(new Float32Array(w_a));
            this.dist_history.b.push(new Float32Array(w_b));

            // Employment is tracked as aggregate labor participation.
            let employed_sum = 0;
            for (let i = 0; i < this.n_people; i++) {
                employed_sum += labor_participation[i];
            }
            this.history_employment.push(employed_sum / this.n_people);

            let preview_yield = Math.max(0, current_yield);
            let preview_taxes_a = this.tax_system_a.calculateTax(current_labor_income);
            let preview_taxes_b = this.tax_system_b.calculateTax(current_labor_income);
            let preview_total_tax_a = preview_taxes_a.reduce((a, b) => a + b, 0);
            let preview_total_tax_b = preview_taxes_b.reduce((a, b) => a + b, 0);

            if (this.tax_capital_gains || this.capital_tax_rate_b > 0) {
                for (let i = 0; i < this.n_people; i++) {
                    if (this.tax_capital_gains) {
                        preview_total_tax_a += this.tax_system_a.calculateTax([w_a[i] * preview_yield])[0];
                    }
                    preview_total_tax_b += (w_b[i] * preview_yield) * this.capital_tax_rate_b;
                }
            }

            this.revenue_history.a.push(preview_total_tax_a);
            this.revenue_history.b.push(preview_total_tax_b);
            this.ubi_history.a.push((preview_total_tax_a / this.n_people) * this.ubi_pct_a);
            this.ubi_history.b.push((preview_total_tax_b / this.n_people) * this.ubi_pct_b);

            if (year < this.n_years) {
                // Apply deterministic AI disruption to annual labor participation.
                if (year > 0) {
                    if (this.labor_displacement > 0) {
                        for (let i = 0; i < this.n_people; i++) {
                            labor_participation[i] *= (1 - this.labor_displacement);
                            current_labor_income[i] = this.base_labor_income[i] * labor_participation[i];
                        }
                    }
                    current_yield += this.capital_efficiency_boost;
                }

                let eff_yield = Math.max(0, current_yield);
                let decay = Math.min(0, current_yield);

                // Recalculate Taxes/UBI for this year based on current labor and capital income.
                const taxes_a = this.tax_system_a.calculateTax(current_labor_income);
                const taxes_b = this.tax_system_b.calculateTax(current_labor_income);
                const capital_taxes_a = new Float32Array(this.n_people);
                const capital_taxes_b = new Float32Array(this.n_people);
                let total_tax_a = taxes_a.reduce((a, b) => a + b, 0);
                let total_tax_b = taxes_b.reduce((a, b) => a + b, 0);

                for (let i = 0; i < this.n_people; i++) {
                    const cap_yield_a = w_a[i] * eff_yield;
                    const cap_yield_b = w_b[i] * eff_yield;
                    capital_taxes_a[i] = this.tax_capital_gains ? this.tax_system_a.calculateTax([cap_yield_a])[0] : 0;
                    capital_taxes_b[i] = cap_yield_b * this.capital_tax_rate_b;
                    total_tax_a += capital_taxes_a[i];
                    total_tax_b += capital_taxes_b[i];
                }

                const ubi_a = (total_tax_a / this.n_people) * this.ubi_pct_a;
                const ubi_b = (total_tax_b / this.n_people) * this.ubi_pct_b;

                const labor_a = new Float32Array(this.n_people);
                const labor_b = new Float32Array(this.n_people);

                for (let i = 0; i < this.n_people; i++) {
                    labor_a[i] = current_labor_income[i] - taxes_a[i] + ubi_a;
                    labor_b[i] = current_labor_income[i] - taxes_b[i] + ubi_b;
                }

                for (let i = 0; i < this.n_people; i++) {
                    // Pre-Tax Baseline
                    let flow_pre = current_labor_income[i] + (w_pre[i] * eff_yield);
                    w_pre[i] = w_pre[i] + (w_pre[i] * decay) + (flow_pre * get_savings_rate(flow_pre));

                    // Policy A
                    let cap_yield_a = w_a[i] * eff_yield;
                    let cap_tax_a = capital_taxes_a[i];
                    let flow_a = labor_a[i] + (cap_yield_a - cap_tax_a);
                    w_a[i] = w_a[i] + (w_a[i] * decay) + (flow_a * get_savings_rate(flow_a));

                    // Policy B
                    let cap_yield_b = w_b[i] * eff_yield;
                    let cap_tax_b = capital_taxes_b[i];
                    let flow_b = labor_b[i] + (cap_yield_b - cap_tax_b);
                    w_b[i] = w_b[i] + (w_b[i] * decay) + (flow_b * get_savings_rate(flow_b));
                }
                w_pre.sort(); w_a.sort(); w_b.sort();
            }
        }

        return {
            gini: this.gini_history,
            distribution: this.dist_history,
            employment: this.history_employment,
            revenue: this.revenue_history,
            ubi: this.ubi_history
        };
    }
}