The “Simulation Modification Functions” section of Tutorial 4 at https://www.mimiframework.org/Mimi.jl/dev/tutorials/tutorial_4/ says the function “replaceRV!” is “available to modify an existing SimulationDef.” Can this function be used to update a random variable for use in a Monte Carlo simulation? I want to update the climatesensitivity random variable which is currently “Truncated(Gamma(6.47815626,0.547629469), 1.0, Inf)”. How can I do this?
Hi Ken, yes that is what that function is for! Sorry the documentation is sparse, we might change some of these functions slightly in the future, but in Mimi v0.9.5 (the most recent version) you should be able to do something like the following:
using Distributions
using Mimi
Mimi.replaceRV!(sim_def, :climatesensitivity, Uniform(1, 10))
Where sim_def is your Mimi SimulationDefinition object. The second argument is a Symbol for the name of the random variable you want to update, and the third argument is the distribution you want to use instead.
Let me know if this doesn’t work for some reason!
Thanks for you reply. Unfortunately the suggested function format didn’t work. I wrote:
using Distributions
using Mimi
using MimiFUND
m = MimiFUND.get_model()
Mimi.replaceRV!(m, :climatesensitivity, Uniform(1,3))
This gives:
MethodError: no method matching replaceRV!(::Model, ::Symbol, ::Uniform{Float64})
Closest candidates are:
replaceRV!(!Matched::SimulationDef, ::Symbol, ::Distribution) at C:\Users\Ken\.julia\packages\Mimi\rkRPr\src\mcs\defmcs.jl:236
It doesn't like the MimiFUND model as the 'sim_def'.
Further suggestions would be most appreciated!Yep, get_model returns a Mimi Model type, not a Mimi SimulationDefinition type. To access MimiFUND's default Monte Carlo simulation definition, try the following:
using Distributions
using Mimi
using MimiFUND
simdef = MimiFUND.getmcs()
Mimi.replaceRV!(simdef, :climatesensitivity, Uniform(1,3))
And then to use this modified simulation definition, try the following:
m = MimiFUND.get_model()
sim_results = run(simdef, m, N) # N is the number of trials you want to run
There are other keyword arguments to the run function that you may be interested in for saving different output.
An additional function you can use to save specific variables from the simulation is the addSave! function, which you can do before running the simulation:
Mimi.addSave!(simdef, :climatedynamics, :temp) # examples
Mimi.addSave!(simdef, :impactaggregation, :loss)
# can save others with: Mimi.addSave!(simdef, :component_name, :variable_name)
# And then after running the simulation, you can access these results
sim_results = run(simdef, m, N)
sim_results[:climatedynamics, :temp] # gives the results for this variable as a DataFrame
explore(sim_results) # Displays the results graphically
HI Cora,
Thanks for your reply. You said “There are other keyword arguments to the run function”, but I looked at the Mimi User Guide and at ‘Running a Model’ it only shows ‘run(mymodel)’. It seems I am missing a listing of functions with their descriptions and arguments. Where can I find the other keyword arguments to the ‘run’ function?
I input into Jupyter notebook as per your instructions (I think);
simdef = MimiFUND.getmcs()
Mimi.replaceRV!(simdef, :climatesensitivity, Normal(1.5,0.003))
m = MimiFUND.get_model()
Mimi.addSave!(simdef, :climatedynamics, :climatesensitivity)
sim_results = run(simdef, m, 100);
sim_results[:climatedynamics, :climatesensitivity]
I used a very small std deviation for the normal distribution of climate sensitivity (CS) so I expected the 100 results to be all close to 1.5. The results returned 30 values, averaging 3.52, which is near the FUND default, not my modified CS. The other 70 values were not displayed, only …
How do I get results using my modified CS random variable? Does the command ‘m = MimiFUND.get_model()’ get the default model including the default random variables and overwrite my CS variable?
How do I display all the values? Notebook usually presents a long list of values, such as 300 years of temperatures, in a window with a slider.
I also want to use the the function:
scc = MimiFUND.compute_scco2(m, year = 2020, eta = 0., prtp = 0.05, equity_weights = false, n = 100) which is described at https://github.com/fund-model/MimiFUND.jl/
But when preceding this with Mimi.replaceRV!(simdef, :climatesensitivity, Normal(1.5,0.003))
the compute_scc ignores my CS and uses the FUND default values.
I would like to run the scco2 function in Monte Carlo mode while using my modified CS values. How can I do this?
Hi Ken, hopefully I’ve addressed each of your questions below! Let me know if there is anything else that still isn’t clear.
replacing the RV correctly:
I just ran the code I had suggested, and you are right, it does not correctly replace the values as I had thought. My apologies for leading you astray, the functionality is a bit finicky because Mimi currently appends unique numbers to the names of the random variables in the simulation’s definition under the hood, so the random variable that you are trying to replace isn’t actually called :climatesensitivity, it’s called :climatesensitivity!265. I had forgotten this, and it is obviously not ideal behavior, given how misleading that is, and is something we are hoping to change about these functions in the future.
For now though, you can use the replaceRV! function with the name :climatesensitivity!265 in the following way:
using Distributions
using Mimi
using MimiFUND
simdef = MimiFUND.getmcs();
Mimi.replaceRV!(simdef, :climatesensitivity!265, Normal(1.5, 0.003));
Mimi.addSave!(simdef, :climatedynamics, :climatesensitivity);
m = MimiFUND.get_model();
sim_results = run(simdef, m, 100);
It should be understood though that this solution is a bit brittle, because if something changes in MimiFUND’s SimulationDefinition, the number appended to the climatesensitivity’s name could also change. I more robust (but lengthier) solution would be to simply copy the code from MimiFUND for defining the whole Simulation (including all random variables) and simply modify the line for climatesensitivity (here: https://github.com/fund-model/MimiFUND.jl/blob/master/src/montecarlo/defmcs.jl)
viewing the full list of values:
the returned value from the following is a DataFrame object, but if you extract the values from the desired column as a vector it might display the way you are hoping in notebook:
df = sim_results[:climatedynamics, :climatesensitivity];
values = df[!, :climatesensitivity]
Computing the social cost of carbon:
The compute_scco2 function in MimiFUND does not currently support passing in a modified SimulationDefinition (something we may add later). So in order to calculate the SC-CO2 with your modified simulation, you will have to dig in to MimiFUND’s code a bit more, copy out the relevant functions yourself, and just replace your modified SimulationDefinition where necessary. That code is all available on MimiFUND’s github repository, most of the relevant code should be in the “src/new_marginaldamages.jl” file.
other keyword arguments to the run function:
So there are two run methods implemented in Mimi: one which runs a single model over one set of it’s default values, which looks like:
run(m)
and the second which runs a Simulation for a given model and a SimulationDefinition (which defines all of the random variables to sample from) and looks like:
run(simdef, m, N)
This second one is the one I mentioned that has other optional keyword arguments. In a Julia REPL, if you type a single question mark ? that should open the “help” bar, and if you then type a function name it should display the docstrings, which will describe the other optional arguments.
julia> ?
help?> run
If you are primarily working in notebook, I’m not sure what the best way to view the docstring is. If you can’t find a way to view it, let me know and I can just paste it here instead.
Alternatively, this page in the Docs under “Step 4” shows the various arguments, but doesn’t go into detail about all of them: https://www.mimiframework.org/Mimi.jl/stable/tutorials/tutorial_4/#Advanced-Post-trial-Functions-1 So let me know if you have specific questions about any of them.
I may want to change other RVs. How do you know the number that you used to append to ‘climatesensitivity’ to become ‘climatesensitivity!265’? It would be easy if the number was the line number in the file defmcs.jl but climatesensitivity is on line 28. I will want to change the RV agcbm. How do I find the agcbm unique number and numbers for other RVs?
The command ‘values = df[!, :climatesensitivity]’ only displays some of the numbers. But println(values[i]) for i = 1:100 prints all of them.
For computing the social cost of carbon with modified RVs, you are recommending to create two models with ‘Mimi.replaceRV!’, m1 and m2. Then copy the scc code into notebook to add a CO2 pulse to m2, create a marginal model of losses from m2 - m1, then discount and calculate the scc. You do not recommend me modifying the RVs is the original source code of the file defmcs.jl which is located at C:\Users\Ken.julia\packages\MimiFUND\3NP0Z\src\montecarlo\defmcs.jl. Correct? I’ll work on this.
I used the following code to access the full name of the random variable:
simdef = MimiFUND.getmcs()
filter(x->startswith(String(x),"climatesensitivity"), collect(keys(simdef.rvdict)))
The reason the number isn’t just 28 for climatesensitivity is that some of the previous lines in the defmcs.jl file are actually defining more than one random variable, because some of them are defining different distributions for different regions of the model, so Mimi has to store multiple random variables for the different regions (hence the unique number appended to the name). If you check this out for the “agcbm” random variable, you’ll see that there are 16 of them, one for each region:
filter(x->startswith(String(x),"agcbm"), collect(keys(simdef.rvdict)))
For this reason, I would probably just recommend you copy the whole simulation definition and change the distributions you want to modify directly. However, you could still modify the ones you want with replaceRV! if you find the correct full name.
For calculating the SC-CO2:
you don’t actually call replaceRV! on the two models, you just call it on the SimulationDefinition object obtained by MimiFUND.getmcs() which stores all of the information about the random variables to use in the monte carlo simulation. The simulation is then run for two models (one with a pulse of emissions) to calculate the SCC. What I’m recommending is that you copy the MimiFUND code for the compute_scco2 function, and just replace the part in that code where it calls simdef = MimiFUND.getmcs() and replace that with your modified SimulationDefinition. Everything else about that function (setting up the two models, using a post-trial function for calcualting the scc, etc.) should remain the same.
Thanks for your patience!
I copied the function ‘sc_compute’ from file ‘new_marginaldamages.jl’ into notebook. (The compute_scco2 function just calls the ‘sc_compute’ function.) The ‘sc_compute’ function calls other functions, some of which then calls other functions. So I also copied functions:
function getindexfromyear
function add_marginal_emissions!
function get_marginal_model
function _compute_sc_from_mm
function _fund_sc_post_trial
@defcomp emissionspulse begin
I copied all of file defmcs.jl into notebook. The first two lines of the file are:
import Mimi.compinstance
fund_default_mcs = @defsim begin
I am not sure what the first line does but I left it in.
I revised the second line to fund_rev_mcs = @defsim begin. The last 3 lines of the file are:
function getmcs()
return deepcopy(fund_default_mcs)
end
I changed the first 2 of these lines to:
function getmcsrev()
return deepcopy(fund_rev_mcs).
I don’t know what 'return deepcopy means, but I thought it best to change ‘default’ to ‘rev’.
The only change to the RV I made was “climatesensitivity = Normal(1.5,0.00000000000000001)”.
In function compute_sc I changed “simdef = getmcs()” to “simdef = getmcsrev()” which corresponds to the function line of the copied and revised dfmcs.fl file.
To calculate the SCCO2 with Monte Carlo with a changed climate sensitivity of 1.5 I ran:
SCCO2FileName = “scCO2_ecs1_5.csv”
file1 = open(SCCO2FileName, “w”)
m = MimiFUND.get_model()
SCCO2 = MimiFUND.compute_sc(m; gas = :CO2, year = 2020, eta = 0.0, prtp = 0.03, equity_weights = false, equity_weights_normalization_region = 0, last_year = 2400, pulse_size = 1e7, return_mm = true, n = 100, trials_output_filename = “/Users/Ken/Documents/JuliaFUND/scCO2_ecs1_5.csv”);
This runs without errors. This print command “print(mean(SCCO2.sc),” ",median(SCCO2.sc), " “,mode(SCCO2.sc))” gave:
10.687570241976385 9.034459619378786 2.975231737754152
However, the file scCO2_ecs1_5.csv should have displays climate sensitivities of 1.5 as per my modified defmcs. It actualy displays 100 larger CS, of median 3.28 and average 3.48 !! So my attempt to calculate SCC with a modified CS failed. Can you please tell me where I went wrong?
I also investigated the marginal model returned by ;
marg_mco2 = SCCO2.mm[:emissions, :mco2]
for i = 71:451
print(marg_mco2[i],", ")
end
The incremental emissions of CO2 should be a million tonnes CO2 from 10 years starting 2020, but the maximum emissions returned in any year was 3.22e-12 !
My notebook file is at: https://www.dropbox.com/s/wuh5kupkgptepbb/FUND%20Random%20Variables%20for%20SCCO2.ipynb?dl=0
I think your problem might be that you shouldn’t call MimiFUND.compute_sc, which will still use MimiFUND’s definition of that function.
I think instead you should just call compute_sc which should call the version of the function that you modified in your notebook. (You could also rename the one that you changed to compute_sc_rev to make it a little more clear. Either way though, it’s defined within your environment, and not within the MimiFUND module, so you don’t prefix it with a MimiFUND.)
Regarding the small emissions in the returned marginal model though, they return values normalized by the pulse size so that the units should just be the marginal emissions from 1 ton of CO2 (even though the pulse was a million tons of C for ten years). Given that though, I think the units on :mco2 are MtC, so I would expect to see a value of at least 1e-7 * 12/44 in 2020. I’m curious why that’s not the case, so I’ll look into that…
Ah I figured it out: the “pulse” of emissions gets added as an intermediate component between the emissions component and the CO2 cycle component. If you try:
marg_mco2 = SCCO2.mm[:climateco2cycle, :mco2]
marg_mco2[71] * 44/12
you should see the correct numbers.
That is wonderful. It now works fine! I graphed the inc. temperature, inc. loss and the CO2 pulse from the marginal model. The inc. temperature goes up to 5.4 E-13 °C by 2029 then declines. Emissions are 0.1 tCO2/y for 10 years = 1 tCO2. The loss is positive for 4 years, then mostly negative to 2100 then positive, but with big up and down spike.
I am unclear in what the marginal loss number represents. I ran the model with n = 5000 trials, so there should be 5000 values for the marginal loss in each year. From “marg_loss = SCCO2.mm[:impactaggregation, :loss]”, the loss in 2020 was 0.021. Is this the loss in 2020 of the n = 5000th trial, or the mean of all 5000 trials, or the median of all 5000 trials?
The calculated “median(SCCO2.sc)” is 0.63 with eta = 0.0, prtp = 0.03. But when I apply 3% discounting to the loss for years 2020 to 2400 in Excel, the discounted loss is -0.74. If the loss returned by the SCCO2.mm is the median of 5000 trials, I expected these two number to be the same. Can you explain why these are different?
You wrote “even though the pulse was a million tons of C for ten years” but I think the pulse size is 10 MtCO2, or 1.0 MtCO2/yr for 10 years. That is the pulse size in the compute_scco2 is in tCO2. It is converted in the function add_marginal_emissions! by “pulse_size / 1e7 * (gas == :CO2 ? 12/44 : 1)” to 2.727 MtC, then we convert it back to 10 MtCO2 in notebook via * 44/12. Do you agree?
The values stored in the returned marginal model are just from the final run of the simulation (not very useful). To save results from all 5000 trials, you need to modify the simulation definition with the Mimi.addSave! function described above, which will save the results from all trials to a file (and then you will have to read the files in and perform desired averages/medians yourself).
And yes, that’s what I meant by the pulse size is one million for ten years-- one million in each year for ten years.