Merge pull request #27 from fund-model/pulse-size

pulse_size option and other changes to compute_scco2

Project.toml

@@ -7,18 +7,19 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Mimi = "e4e893b0-ee5e-52ea-8111-44b3bdec128c"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
-[extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-CSVFiles = "5d742f6a-9f54-50ce-8119-2520741973ca"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-
 [compat]
 julia = "1"
-Mimi = "0.9.3"
+Mimi = "0.9.4"
 Distributions = "0.20, 0.21"
 StatsBase = "0.30, 0.31, 0.32"
 
+[extras]
+CSVFiles = "5d742f6a-9f54-50ce-8119-2520741973ca"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 [targets]
 test = ["Test", "DataFrames", "CSVFiles"]

README.md

@@ -43,7 +43,7 @@ pkg> add MimiFUND
 
 You probably also want to install the Mimi package into your julia environment, so that you can use some of the tools in there:
 
-```
+```julia
 pkg> add Mimi
 ```
 
@@ -59,49 +59,87 @@ pkg> up
 The model uses the Mimi framework and it is highly recommended to read the Mimi documentation first to understand the code structure. For starter code on running the model just once, see the code in the file `examples/main.jl`.
 
 The basic way of accessing a copy of the default MimiFUND model is the following:
-```
+```julia
 using MimiFUND
 
 m = MimiFUND.get_model()
 run(m)
 ```
-## Calculating the Social Cost of Carbon
+## Calculating the Social Cost of CO2 and other gases
 
-Here is an example of computing the social cost of carbon with MimiFUND. Note that the units of the returned value are $/ton CO2.
+Here is an example of computing the Social Cost of CO2 with MimiFUND. Note that the units of the returned value are 1995$ per metric tonne of CO2.
 
-```
+```julia
 using Mimi
 using MimiFUND
 
-# Get the social cost of carbon in year 2020 from the default MimiFUND model:
-scc = MimiFUND.compute_scc(year = 2020)
+# Get the Social Cost of CO2 in year 2020 from the default MimiFUND model:
+scc = MimiFUND.compute_scco2(year = 2020, eta = 0., prtp = 0.03, equity_weights = false)
 
-# Or, you can also compute the SCC from a modified version of a MimiFUND model:
-m = MimiFUND.get_model() # Get the default version of the FUND model
-update_param!(m, :climatesensitivity, 5) # make any modifications to your model
-scc = MimiFUND.compute_scc(m, year = 2020) # Compute the SCC from your model
+# Or, you can also compute the SC-CO2 from a modified version of a MimiFUND model:
+m = MimiFUND.get_model()                        # Get the default version of the FUND model
+update_param!(m, :climatesensitivity, 5)        # make any modifications to your model using Mimi
+scc = MimiFUND.compute_scco2(m, year = 2020)    # Compute the SC-CO2 from your model
 ```
+There are also functions for computing the Social Cost of CH4, N2O, and SF6: `compute_scch4`, `compute_scn2o`, and `compute_scsf6`.
+These functions are all wrappers for the generic social cost function `compute_sc`, which takes a keyword `gas` with default value `:CO2`.
 
-There are several keyword arguments available to `compute_scc`. Note that the user must specify a `year` for the SCC calculation, but the rest of the keyword arguments have default values.
-```
-compute_scc(m = get_model(),  # if no model provided, will use the default MimiFUND model
-    year = nothing,  # user must specify an emission year for the SCC calculation
-    gas = :CO2,  # which greenhouse gas to use. Other options are :CH4, :N2O, and :SF6.
-    last_year = 3000,  # the last year to run and use for the SCC calculation. Default is the last year of the time dimension, 3000.
-    eta = 1.45,  # eta parameter for ramsey discounting representing the elasticity of marginal utility
-    prtp = 0.015,  # pure rate of time preference parameter for discounting
-    equity_weights = false  # whether or not to use regional equity weighting
-)
+There are several other keyword arguments available to `compute_sc`. Note that the user must specify a `year` for the SC calculation, 
+but the rest of the keyword arguments have default values.
+```julia
+MimiFUND.compute_sc(m::Model=get_model();
+        gas::Symbol = :CO2,                     
+        year::Union{Int, Nothing} = nothing,    
+        eta::Float64 = 1.45,                    
+        prtp::Float64 = 0.015,                  
+        equity_weights::Bool = false,           
+        last_year::Int = 3000,                  
+        pulse_size::Float64 = 1e7,              
+        return_mm::Bool = false,
+        n::Union{Int, Nothing} = nothing,
+        trials_output_filename::Union{String, Nothing} = nothing,
+        seed::Union{Int, Nothing} = nothing)
 ```
+Description of keyword arguments:
+- **`m`**: a MimiFUND model from which to calculate the social cost. If no model is provided, the default MimiFUND model will be used. Note that the provided model `m` can be a highly modified MimiFUND model, but certain internal structures of the model need to remain in order for the `compute_sc` function to work. They are:
+    - The original parameter connection between the `emissions` component and the climate cycling component for the specified `gas` must still be intact (this is where the pulse of emissions is added).
+    - There must still be a `:socioeconomic` component with fields `:ypc` and `:globalypc` (used for discounting).
+    - There must still be an `:impactaggregation` component with field `:loss`, which is the total damages value from which the social cost is calculated.
+- **`gas`**: which greenhouse gas to calculate the social cost for. The default is `:CO2`. Other options are `:CH4`, `:N2O`, and `:SF6`.
+- **`year`**: the user must specify an emission year for the SC calculation. Valid years are 1951 to 2990.
+- **`eta`**: the elasticity of marginal utility to be used in ramsey discounting. Setting `eta = 0` is equivalent to constant discounting with rate `prtp`.
+- **`prtp`**: pure rate of time preference parameter for discounting
+- **`equity_weights`**: whether or not to use regional equity weighting in discounting
+- **`last_year`**: the last year to run and use for the SC calculation. Default is the last year of FUND's time index, 3000.
+- **`pulse_size`**: the size of the marginal emissions pulse, in metric tonnes of the specified `gas`. Changing this value will not change the units of the returned value, which are always "1995$ per metric tonne of `gas`". The returned value is always normalized by the size of `pulse_size` that is used.
+- **`return_mm`**: whether or not to also return the MarginalModel used in the social cost calculation. If set to `true`, then a NamedTuple `(sc = sc, mm = mm)` of the social cost value and the MarginalModel used to compute it is returned.
+- **`n`**: By default, `n = nothing`, and a single value for the "best guess" social cost is returned. If a positive value for keyword `n` is specified, then a Monte Carlo simulation with sample size `n` will run, sampling from all of FUND's random variables, and a vector of `n` social cost values will be returned. Note that if the user has provided a modified model `m`, the user modifications may be overwritten by the Monte Carlo simulation. If the user has modified certain parameter values, but they are parameters that have assigned random variable distributions in FUND, then they will be overwritten by the sampled values. For a list of which parameters have assigned random variable definitions, see "src/montecarlo/defmcs.jl"
+- **`trials_output_filename`**: an optional CSV file path to save all of the sampled trial data.
+- **`seed`**: the user can optionally provide a seed value, which will set the random seed before the simulation is run. This allows results to be replicated. 
+
+
+Example Monte Carlo simulation:
+```julia
+using Mimi
+using MimiFUND
 
-There is an additional function for computing the SCC that also returns the MarginalModel that was used to compute it. It returns these two values as a NamedTuple of the form (scc=scc, mm=mm). The same keyword arguments from the `compute_scc` function are available for the `compute_scc_mm` function. Example:
+scco2_values = MimiFUND.compute_sc(year = 2020, gas = :CO2, eta = 1.0, prtp = 0.01, n = 1000)
+mean(scco2_values)
+median(scco2_values)
+
+# Experiment with the same set of trial data by setting the seed (any Integer value)
+values_lo_discounting = MimiFUND.compute_sc(year = 2020, gas = :CO2, eta = 1., prtp = 0.015, n = 1000, seed = 999)
+values_hi_discounting = MimiFUND.compute_sc(year = 2020, gas = :CO2, eta = 1., prtp = 0.05, n = 1000, seed = 999)
 ```
+
+Example of working with the MarginalModel from setting `return_mm = true`:
+```julia
 using Mimi
 using MimiFUND
 
-result = compute_scc_mm(year=2020, last_year=2300, eta=0, prtp=0.03)
+result = MimiFUND.compute_sc(year = 2020, gas = :CH4, last_year = 2300, eta = 0, prtp = 0.03, return_mm = true)
 
-result.scc  # returns the computed SCC value
+result.sc  # returns the computed SC-CH4 value
 
 result.mm   # returns the Mimi MarginalModel
 

src/MimiFUND.jl

@@ -1,7 +1,8 @@
 module MimiFUND
 
 using Mimi
-using DelimitedFiles #base.DelimitedFiles
+using DelimitedFiles
+using Random
 
 include("helper.jl")
 

src/montecarlo/run_fund_scc_mcs.jl

@@ -29,8 +29,7 @@ function run_fund_scc_mcs(trials = 10000; years = [2020], rates = [0.03], ntimes
 
     # get models and sim
     mcs = getmcs()
-    mm = create_marginal_FUND_model()
-    models = [mm.base, mm.marginal] #fixing bug in Mimi so this can go back to models = mm
+    mm = get_marginal_model()
 
     # Define scenario function
     function _scenario_func(mcs::SimulationInstance, tup::Tuple)
@@ -39,10 +38,10 @@ function run_fund_scc_mcs(trials = 10000; years = [2020], rates = [0.03], ntimes
         (rate, emissionyear) = tup
 
         # Get models
-        base, marginal = mcs.models
+        mm = mcs.models[1]
 
         # Perturb emissons in the marginal model
-        perturb_marginal_emissions!(marginal, emissionyear)
+        perturb_marginal_emissions!(mm.marginal, emissionyear)
 
     end
 
@@ -53,11 +52,11 @@ function run_fund_scc_mcs(trials = 10000; years = [2020], rates = [0.03], ntimes
         (rate, emissionyear) = tup
 
         # Get marginal damages
-        base, marginal = mcs.models
-        marginaldamages = (marginal[:impactaggregation, :loss] - base[:impactaggregation, :loss]) / 10000000.0
+        mm = mcs.models[1]
+        marginaldamages = mm[:impactaggregation, :loss]
 
         # Calculate discount factor
-        T = ntimesteps == typemax(Int) ? length(Mimi.dimension(base, :time)) : ntimesteps
+        T = ntimesteps == typemax(Int) ? length(Mimi.dimension(mm.base, :time)) : ntimesteps
         discount_factor = zeros(T)
         idx = getindexfromyear(emissionyear)
         discount_factor[idx:T] = [1 / ((1 + rate) ^ t) for t in 0:T-idx]
@@ -74,7 +73,7 @@ function run_fund_scc_mcs(trials = 10000; years = [2020], rates = [0.03], ntimes
     end
 
     # Run monte carlo trials
-    res = run(mcs, models, trials;
+    res = run(mcs, mm, trials;
         ntimesteps = ntimesteps,
         trials_output_filename = trials_output_filename,
         results_output_dir = "$output_dir/results",