function evaluate_tax_reform()


% Assigning parameters
crra              = 1/2  ;
wage              = 1.0  ;
interest_rate     = 1.04 ;
capital           = 0.0  ;


% Creating a structure of parameters
parameters.crra              = crra              ;
parameters.wage              = wage              ;
parameters.interest_rate     = interest_rate     ;
parameters.capital           = capital           ;


% Solving for benchmark utility
lump_sum_tax          = 0.1 ;
proportional_tax_rate = 0   ;
[benchmark_government_revenue, benchmark_utility] = evaluate_tax_regime( lump_sum_tax, proportional_tax_rate, parameters ) ;


% Testing to ensure that benchmark government revenue is equal to the lump sum tax
if ~eq( benchmark_government_revenue, lump_sum_tax )
    error('Government tax revenue does not equal lump-sum tax in benchmark case.')
end


% Initializing the bisection boundaries
tax_rate_ub  = 0.6 ;
tax_rate_lb  = 0   ;

% Solving for proportional tax reform
lump_sum_tax = 0        ;
tolerance    = 0.000001 ;
iteration    = 0        ;
while true
    iteration = iteration + 1;
    
    % Calculating government revenue and utility
    proportional_tax_rate = (tax_rate_ub + tax_rate_lb)/2                                                  ;
    [reform_government_revenue, reform_utility] = evaluate_tax_regime( lump_sum_tax, proportional_tax_rate, parameters ) ;
    
    % Testing proximity to benchmark revenue
    deviation = abs(reform_government_revenue - benchmark_government_revenue);
    if deviation<tolerance
        break
    end
    
    % Updating guess using bisection
    if reform_government_revenue<benchmark_government_revenue
        
        % Tax revenue too low: raise tax rate
        tax_rate_lb = proportional_tax_rate;
        
    elseif reform_government_revenue>benchmark_government_revenue
        
        % Tax revenue too high: reduce tax rate
        tax_rate_ub = proportional_tax_rate;
    end
    
end

% Printing results
fprintf( 'Government budget equalized after %d iterations.\n'             , iteration                        )
fprintf( 'Government revenue equalized at proportional tax rate %0.4f.\n' , proportional_tax_rate            )
fprintf( 'Percentage change in utility from the reform is %0.4f.\n'       , reform_utility/benchmark_utility )


end


% Calculating total government revenue and utility
function [total_government_revenue, optimal_utility] = evaluate_tax_regime( lump_sum_tax, proportional_tax_rate, parameters )

    initial_guess                    = 0.15                                                                                       ;
    [optimal_labor, optimal_utility] = fminsearch(@(labor) calculate_utility( labor, proportional_tax_rate, lump_sum_tax, parameters), initial_guess )  ;
    optimal_utility                  = -1*optimal_utility                                                                         ;
    proportional_tax_revenue         = proportional_tax_rate * parameters.wage * optimal_labor                                    ;
    lump_sum_tax_revenue             = lump_sum_tax                                                                               ;
    total_government_revenue         = proportional_tax_revenue + lump_sum_tax_revenue                                            ;

end




% Creating objective function

function utility = calculate_utility( labor, proportional_tax_rate, lump_sum_tax, parameters )
    
    consumption = parameters.wage*(1-proportional_tax_rate) * labor + parameters.interest_rate * parameters.capital  - lump_sum_tax             ;
    
    utility     = ( consumption^parameters.crra ) / parameters.crra  - labor                                                                    ;
    
    utility     = -1 * utility                                                                                                                  ;    % Converting into minimization problem.
    
end