Complex exponentiation error in SCC Monte Carlo calculation

1

Hi all,

When running a 10,000-trial Monte Carlo SCC simulation with a modified version of MimiGIVE, I run into an error part of the way through the trials: “DomainError with -0.988…[more decimals]: Exponentiation yielding a complex result requires a complex argument.”

There is a stacktrace, but it’s not particularly informative as far as I can tell. In one of the MC iterations, the SCC calculation is requiring an imaginary/complex number (as a result of taking a square root of a negative, maybe?) and breaking. One of the elements of the stacktrace mentions a variable that includes :prtp and :eta, so it’s possible that this is occurring during some calculation involving a discount rate.

Changing the number of iterations to 1,000 avoids the problem, potentially since it doesn’t reach the particular trial that triggers it. Going back to 10,000 but changing the seed yields the same error, but on a different negative number and on a different trial.

Has anyone run into this kind of error before? Are there ways to avoid it?

Thank you!

2

Hi @christophercallahan your interpretation here is likely correct. That DomainError about complex result requires a complex argument generally comes from something like taking the square root of a of a negative, and to be more specific as you point out it often occurs in the discounting calculation in the post_trial_function of the Monte Carlo Simulation, likely in the Ramsey discounting formula here:

df = [((cpc[year_index]/cpc[i])^dr.eta * 1/(1+dr.prtp)^(t-year) for (i,t) in enumerate(_model_years) if year<=t<=last_year)...]

Looking at that formula, and the fact you are adding a damage function and modifying the base GIVE model, I would presume something in your damage function might be causing (in extreme outlier cases, given as you say it only happens a few times 10k trials) consumption per capita to go negative, and thus ((cpc[year_index]/cpc[i]) becomes negative. The stack trace of your error should point to a line of code.

First off, you should feel free to use a seed (Random.seed(1234) to run the same mcs every time, I recommend setting the seed right before compute_scc. This is not completely fool proof over time, it can change in Julia versions or OS in some cases etc., but for a local machine it’s quite dependable.

In terms of debugging, one thing that’s worked for me in the past is as follows:

(1) Add a domain_error_trialnums empty vector to the payload (anywhere you create, load, or unload the payload you add this in).

(2) Add a try ... catch block to the df computation that does something like the following

        df = try
                [((cpc[year_index]/cpc[i])^dr.eta * 1/(1+dr.prtp)^(t-year) for (i,t) in enumerate(MimiGIVE._model_years) if year<=t<=last_year)...]
            catch err 
                if isa(err, DomainError)
                    [(1. for (i,t) in enumerate(MimiGIVE._model_years) if year<=t<=last_year)...]
                    append!(domain_error_trialnums, trialnum)
                end
            end

which basically just sets the df to all ones (you could use anything you just want something that won’t break) and then adds the trial number to your list of errors.

(3) Additionally add some things to the save_list like the Socioeconomic components id, and perhaps consumption per capita, damages from sectors, etc. anything that will help you figure out what is causing consumption to drop. You can also exogenously set the fair_parameter_set_ids to some random set, so that when you are done you will be able to recover the fair parameter set and the RFF-SPs id, which is often enough to recreate a model with that error.

(4) run!

For the sake of simplicity, I am going to move any further conversation about behavior of, and debugging of, MimiGIVE to an issue here: Complex Exponentiation Error in SCC Monte Carlo Calculation · Issue #21 · rffscghg/MimiGIVE.jl · GitHub since further discussion will likely be very model-specific. That said don’t hesitate to reach out about this I’m happy to help further, I’ve debugged this before and have some good insights on how to do it so you don’t have to reinvent the wheel :slight_smile: