From dd237de06d19e00a42fd7f4a331184a36b7315ee Mon Sep 17 00:00:00 2001 From: davidanthoff Date: Fri, 26 Aug 2022 06:11:38 +0000 Subject: [PATCH] Format files using DocumentFormat --- src/MimiRICE2010.jl | 2 +- src/components/climatedynamics_component.jl | 2 +- src/components/co2cycle_component.jl | 2 +- src/components/damages_component.jl | 6 +- src/components/emissions_component.jl | 14 +- src/components/grosseconomy_component.jl | 8 +- src/components/neteconomy_component.jl | 14 +- src/components/radiativeforcing_component.jl | 4 +- src/components/slr_component.jl | 152 +++++++------- src/components/slrdamages_component.jl | 4 +- src/components/welfare_component.jl | 16 +- src/helpers.jl | 20 +- src/marginaldamage.jl | 44 ++-- src/parameters.jl | 36 ++-- test/runtests.jl | 200 +++++++++---------- 15 files changed, 262 insertions(+), 262 deletions(-) diff --git a/src/MimiRICE2010.jl b/src/MimiRICE2010.jl index 90e6be0..cd655dc 100644 --- a/src/MimiRICE2010.jl +++ b/src/MimiRICE2010.jl @@ -92,7 +92,7 @@ function constructrice(p) return m end #function -function get_model(;datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")) +function get_model(; datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")) params = getrice2010parameters(datafile) m = constructrice(params) diff --git a/src/components/climatedynamics_component.jl b/src/components/climatedynamics_component.jl index 52b42cc..c4161cf 100644 --- a/src/components/climatedynamics_component.jl +++ b/src/components/climatedynamics_component.jl @@ -21,7 +21,7 @@ elseif t.t == 2 v.TATM[t] = p.tatm1 else - v.TATM[t] = v.TATM[t-1] + p.c1 * ((p.FORC[t] - (p.fco22x/p.t2xco2) * v.TATM[t-1]) - (p.c3 * (v.TATM[t-1] - v.TOCEAN[t-1]))) + v.TATM[t] = v.TATM[t-1] + p.c1 * ((p.FORC[t] - (p.fco22x / p.t2xco2) * v.TATM[t-1]) - (p.c3 * (v.TATM[t-1] - v.TOCEAN[t-1]))) end #Define function for TOCEAN diff --git a/src/components/co2cycle_component.jl b/src/components/co2cycle_component.jl index a050c89..358d0b8 100644 --- a/src/components/co2cycle_component.jl +++ b/src/components/co2cycle_component.jl @@ -49,7 +49,7 @@ if is_first(t) v.MATSUM[t] = 0 else - v.MATSUM[t] = v.MAT[t] + (v.MAT[t] * p.b11 + v.MU[t] * p.b21 + (p.E[t] * 10)) + v.MATSUM[t] = v.MAT[t] + (v.MAT[t] * p.b11 + v.MU[t] * p.b21 + (p.E[t] * 10)) end end end diff --git a/src/components/damages_component.jl b/src/components/damages_component.jl index 04c6878..07e76bc 100644 --- a/src/components/damages_component.jl +++ b/src/components/damages_component.jl @@ -15,15 +15,15 @@ #Define function for DAMFRAC for r in d.regions - v.DAMFRAC[t,r] = (((p.a1[r] * p.TATM[t]) + (p.a2[r] * p.TATM[t]^p.a3[r])) / 100) + (p.SLRDAMAGES[t,r] / 100) + v.DAMFRAC[t, r] = (((p.a1[r] * p.TATM[t]) + (p.a2[r] * p.TATM[t]^p.a3[r])) / 100) + (p.SLRDAMAGES[t, r] / 100) end #Define function for DAMAGES for r in d.regions if is_first(t) - v.DAMAGES[t,r] = p.YGROSS[t,r] * (1 - 1 / (1+v.DAMFRAC[t,r])) + v.DAMAGES[t, r] = p.YGROSS[t, r] * (1 - 1 / (1 + v.DAMFRAC[t, r])) else - v.DAMAGES[t,r] = (p.YGROSS[t,r] * v.DAMFRAC[t,r]) / (1. + v.DAMFRAC[t,r]^10) + v.DAMAGES[t, r] = (p.YGROSS[t, r] * v.DAMFRAC[t, r]) / (1.0 + v.DAMFRAC[t, r]^10) end end end diff --git a/src/components/emissions_component.jl b/src/components/emissions_component.jl index 3c7782b..91f4ad5 100644 --- a/src/components/emissions_component.jl +++ b/src/components/emissions_component.jl @@ -21,32 +21,32 @@ #Define function for EIND for r in d.regions - v.EIND[t,r] = p.sigma[t,r] * p.YGROSS[t,r] * (1-p.MIU[t,r]) + v.EIND[t, r] = p.sigma[t, r] * p.YGROSS[t, r] * (1 - p.MIU[t, r]) end #Define function for E - v.E[t] = sum(v.EIND[t,:]) + p.etree[t] + v.E[t] = sum(v.EIND[t, :]) + p.etree[t] #Define function for CCA if is_first(t) - v.CCA[t] = sum(v.EIND[t,:]) * 10. + v.CCA[t] = sum(v.EIND[t, :]) * 10.0 else - v.CCA[t] = v.CCA[t-1] + (sum(v.EIND[t,:]) * 10.) + v.CCA[t] = v.CCA[t-1] + (sum(v.EIND[t, :]) * 10.0) end #Define function for ABATECOST for r in d.regions - v.ABATECOST[t,r] = p.YGROSS[t,r] * p.cost1[t,r] * (p.MIU[t,r]^p.expcost2[r]) * (p.partfract[t,r]^(1 - p.expcost2[r])) + v.ABATECOST[t, r] = p.YGROSS[t, r] * p.cost1[t, r] * (p.MIU[t, r]^p.expcost2[r]) * (p.partfract[t, r]^(1 - p.expcost2[r])) end #Define function for MCABATE for r in d.regions - v.MCABATE[t,r] = p.pbacktime[t,r] * p.MIU[t,r]^(p.expcost2[r] - 1) + v.MCABATE[t, r] = p.pbacktime[t, r] * p.MIU[t, r]^(p.expcost2[r] - 1) end #Define function for CPRICE for r in d.regions - v.CPRICE[t,r] = p.pbacktime[t,r] * 1000 * p.MIU[t,r]^(p.expcost2[r] - 1) + v.CPRICE[t, r] = p.pbacktime[t, r] * 1000 * p.MIU[t, r]^(p.expcost2[r] - 1) end end diff --git a/src/components/grosseconomy_component.jl b/src/components/grosseconomy_component.jl index f028349..883beff 100644 --- a/src/components/grosseconomy_component.jl +++ b/src/components/grosseconomy_component.jl @@ -18,20 +18,20 @@ #Define function for K for r in d.regions if is_first(t) - v.K[t,r] = p.k0[r] + v.K[t, r] = p.k0[r] else - v.K[t,r] = (1 - p.dk[r])^10 * v.K[t-1,r] + 10 * p.I[t-1,r] + v.K[t, r] = (1 - p.dk[r])^10 * v.K[t-1, r] + 10 * p.I[t-1, r] end end #Define function for YGROSS for r in d.regions - v.YGROSS[t,r] = (p.al[t,r] * (p.l[t,r]/1000)^(1-p.gama)) * (v.K[t,r]^p.gama) + v.YGROSS[t, r] = (p.al[t, r] * (p.l[t, r] / 1000)^(1 - p.gama)) * (v.K[t, r]^p.gama) end # TODO remove this, just a temporary output trick for r in d.regions - v.L[t,r] = p.l[t,r] + v.L[t, r] = p.l[t, r] end end end diff --git a/src/components/neteconomy_component.jl b/src/components/neteconomy_component.jl index ce47b6c..8b11c5d 100644 --- a/src/components/neteconomy_component.jl +++ b/src/components/neteconomy_component.jl @@ -19,34 +19,34 @@ #Define function for YNET for r in d.regions if is_first(t) - v.YNET[t,r] = p.YGROSS[t,r]/(1+p.DAMFRAC[t,r]) + v.YNET[t, r] = p.YGROSS[t, r] / (1 + p.DAMFRAC[t, r]) else - v.YNET[t,r] = p.YGROSS[t,r] - p.DAMAGES[t,r] + v.YNET[t, r] = p.YGROSS[t, r] - p.DAMAGES[t, r] end end #Define function for Y for r in d.regions - v.Y[t,r] = v.YNET[t,r] - p.ABATECOST[t,r] + v.Y[t, r] = v.YNET[t, r] - p.ABATECOST[t, r] end #Define function for I for r in d.regions - v.I[t,r] = p.S[t,r] * v.Y[t,r] + v.I[t, r] = p.S[t, r] * v.Y[t, r] end #Define function for C for r in d.regions if t.t != 60 - v.C[t,r] = v.Y[t,r] - v.I[t,r] + v.C[t, r] = v.Y[t, r] - v.I[t, r] else - v.C[t,r] = v.C[t-1, r] + v.C[t, r] = v.C[t-1, r] end end #Define function for CPC for r in d.regions - v.CPC[t,r] = 1000 * v.C[t,r] / p.l[t,r] + v.CPC[t, r] = 1000 * v.C[t, r] / p.l[t, r] end end end diff --git a/src/components/radiativeforcing_component.jl b/src/components/radiativeforcing_component.jl index c1ed626..47b3614 100644 --- a/src/components/radiativeforcing_component.jl +++ b/src/components/radiativeforcing_component.jl @@ -10,9 +10,9 @@ function run_timestep(p, v, d, t) #Define function for FORC if is_first(t) - v.FORC[t] = p.fco22x * ((log10((((p.MAT[t] + p.mat1)/2)+0.000001)/596.4)/log10(2))) + p.forcoth[t] #TEMP NOTE: Uses mat1 because it's given as a parameter...not calculated so couldn't use MATSUM + v.FORC[t] = p.fco22x * ((log10((((p.MAT[t] + p.mat1) / 2) + 0.000001) / 596.4) / log10(2))) + p.forcoth[t] #TEMP NOTE: Uses mat1 because it's given as a parameter...not calculated so couldn't use MATSUM else - v.FORC[t] = p.fco22x * ((log10((((p.MATSUM[t])/2)+0.000001)/596.4)/log10(2))) + p.forcoth[t] + v.FORC[t] = p.fco22x * ((log10((((p.MATSUM[t]) / 2) + 0.000001) / 596.4) / log10(2))) + p.forcoth[t] end end end diff --git a/src/components/slr_component.jl b/src/components/slr_component.jl index 9dee7eb..066c9a0 100644 --- a/src/components/slr_component.jl +++ b/src/components/slr_component.jl @@ -39,98 +39,98 @@ function run_timestep(p, v, d, t) - #THERMAL EXPANSION + #THERMAL EXPANSION #Define function for THERMEQUIL v.THERMEQUIL[t] = p.TATM[t] * p.thermeq - #Define function for SLRTHERM - if is_first(t) - v.SLRTHERM[t] = p.therm0 + p.thermadj * (v.THERMEQUIL[t] - p.therm0) - else - v.SLRTHERM[t] = v.SLRTHERM[t-1] + p.thermadj * (v.THERMEQUIL[t] - v.SLRTHERM[t-1]) - end + #Define function for SLRTHERM + if is_first(t) + v.SLRTHERM[t] = p.therm0 + p.thermadj * (v.THERMEQUIL[t] - p.therm0) + else + v.SLRTHERM[t] = v.SLRTHERM[t-1] + p.thermadj * (v.THERMEQUIL[t] - v.SLRTHERM[t-1]) + end #GLACIERS AND SMALL ICE CAPS (GSIC) - #Define function for GSICREMAIN - if is_first(t) - v.GSICREMAIN[t] = p.gsictotal - else - v.GSICREMAIN[t] = p.gsictotal - v.GSICCUM[t-1] - end - - #Define function for GSICMELTRATE - if is_first(t) - v.GSICMELTRATE[t] = 0.01464 # #Bug in Excel Model? gsiceq parameter drops out after first period (this is replicated here by just setting period 1 and exluding gsiceq from equation) - else - v.GSICMELTRATE[t] = p.gsicmelt * 10 * (v.GSICREMAIN[t] / p.gsictotal)^(p.gsicexp) * p.TATM[t] - end - - #Define function for GSICCUM - if is_first(t) - v.GSICCUM[t] = v.GSICMELTRATE[t] - else - v.GSICCUM[t] = v.GSICCUM[t-1] + v.GSICMELTRATE[t] - end + #Define function for GSICREMAIN + if is_first(t) + v.GSICREMAIN[t] = p.gsictotal + else + v.GSICREMAIN[t] = p.gsictotal - v.GSICCUM[t-1] + end + + #Define function for GSICMELTRATE + if is_first(t) + v.GSICMELTRATE[t] = 0.01464 # #Bug in Excel Model? gsiceq parameter drops out after first period (this is replicated here by just setting period 1 and exluding gsiceq from equation) + else + v.GSICMELTRATE[t] = p.gsicmelt * 10 * (v.GSICREMAIN[t] / p.gsictotal)^(p.gsicexp) * p.TATM[t] + end + + #Define function for GSICCUM + if is_first(t) + v.GSICCUM[t] = v.GSICMELTRATE[t] + else + v.GSICCUM[t] = v.GSICCUM[t-1] + v.GSICMELTRATE[t] + end #GREENLAND ICE SHEETS (GIS) - #Define function for GISREMAIN - if is_first(t) - v.GISREMAIN[t] = p.gis0 - else - v.GISREMAIN[t] = v.GISREMAIN[t-1] - v.GISMELTRATE[t-1] / 100 - end - - #Define function for GISMELTRATE - if is_first(t) || t.t == 2 - v.GISMELTRATE[t] = p.gismelt0 - else - v.GISMELTRATE[t] = (p.gismeltabove * (p.TATM[t] - p.gismineq) + p.gismelt0) * v.GISEXPONENT[t-1] - end - - #Define function for GISCUM - if is_first(t) - v.GISCUM[t] = p.gismelt0 / 100 - else - v.GISCUM[t] = v.GISCUM[t-1] + v.GISMELTRATE[t] / 100 - end - - #Define function for GISEXPONENT - if is_first(t) || t.t == 2 - v.GISEXPONENT[t] = 1. - else - v.GISEXPONENT[t] = 1. - (v.GISCUM[t] / p.gis0)^p.gisexp - end + #Define function for GISREMAIN + if is_first(t) + v.GISREMAIN[t] = p.gis0 + else + v.GISREMAIN[t] = v.GISREMAIN[t-1] - v.GISMELTRATE[t-1] / 100 + end + + #Define function for GISMELTRATE + if is_first(t) || t.t == 2 + v.GISMELTRATE[t] = p.gismelt0 + else + v.GISMELTRATE[t] = (p.gismeltabove * (p.TATM[t] - p.gismineq) + p.gismelt0) * v.GISEXPONENT[t-1] + end + + #Define function for GISCUM + if is_first(t) + v.GISCUM[t] = p.gismelt0 / 100 + else + v.GISCUM[t] = v.GISCUM[t-1] + v.GISMELTRATE[t] / 100 + end + + #Define function for GISEXPONENT + if is_first(t) || t.t == 2 + v.GISEXPONENT[t] = 1.0 + else + v.GISEXPONENT[t] = 1.0 - (v.GISCUM[t] / p.gis0)^p.gisexp + end #ANTARCTIC ICE SHEET (AIS) - #Define function for AISMELTRATE - if t.t <= 11 - v.AISMELTRATE[t] = ifelse(p.TATM[t] < 3., (p.aismeltlow * p.TATM[t] * p.aisratio + p.aisintercept), (p.aisinflection * p.aismeltlow + p.aismeltup * (p.TATM[t] - 3.) + p.aisintercept)) - else - v.AISMELTRATE[t] = ifelse(p.TATM[t] < 3., (p.aismeltlow * p.TATM[t] * p.aisratio + p.aismelt0), (p.aisinflection * p.aismeltlow + p.aismeltup * (p.TATM[t] - 3.) + p.aismelt0)) - end - - #Define function for AISCUM - if is_first(t) - v.AISCUM[t] = v.AISMELTRATE[t] / 100 - else - v.AISCUM[t] = v.AISCUM[t-1] + v.AISMELTRATE[t] / 100 - end - - #Define function for AISREMAIN - if is_first(t) - v.AISREMAIN[t] = p.aiswais + p.aisother - else - v.AISREMAIN[t] = v.AISREMAIN[TimestepIndex(1)] - v.AISCUM[t] - end + #Define function for AISMELTRATE + if t.t <= 11 + v.AISMELTRATE[t] = ifelse(p.TATM[t] < 3.0, (p.aismeltlow * p.TATM[t] * p.aisratio + p.aisintercept), (p.aisinflection * p.aismeltlow + p.aismeltup * (p.TATM[t] - 3.0) + p.aisintercept)) + else + v.AISMELTRATE[t] = ifelse(p.TATM[t] < 3.0, (p.aismeltlow * p.TATM[t] * p.aisratio + p.aismelt0), (p.aisinflection * p.aismeltlow + p.aismeltup * (p.TATM[t] - 3.0) + p.aismelt0)) + end + + #Define function for AISCUM + if is_first(t) + v.AISCUM[t] = v.AISMELTRATE[t] / 100 + else + v.AISCUM[t] = v.AISCUM[t-1] + v.AISMELTRATE[t] / 100 + end + + #Define function for AISREMAIN + if is_first(t) + v.AISREMAIN[t] = p.aiswais + p.aisother + else + v.AISREMAIN[t] = v.AISREMAIN[TimestepIndex(1)] - v.AISCUM[t] + end #TOTAL SEA LEVEL RISE AND DAMAGES - #Define function for TOTALSLR - v.TOTALSLR[t] = v.SLRTHERM[t] + v.GSICCUM[t] + v.GISCUM[t] + v.AISCUM[t] + #Define function for TOTALSLR + v.TOTALSLR[t] = v.SLRTHERM[t] + v.GSICCUM[t] + v.GISCUM[t] + v.AISCUM[t] end end diff --git a/src/components/slrdamages_component.jl b/src/components/slrdamages_component.jl index 042911c..2d9b039 100644 --- a/src/components/slrdamages_component.jl +++ b/src/components/slrdamages_component.jl @@ -14,9 +14,9 @@ #Define function for SLRDAMAGES for r in d.regions if is_first(t) - v.SLRDAMAGES[t,r] = 0. + v.SLRDAMAGES[t, r] = 0.0 else - v.SLRDAMAGES[t,r] = 100. * p.slrmultiplier[r] * (p.TOTALSLR[t-1] * p.slrdamlinear[r] + p.TOTALSLR[t-1]^2 * p.slrdamquadratic[r]) * (p.YGROSS[t-1,r] / p.YGROSS[TimestepIndex(1),r])^(1/p.slrelasticity[r]) + v.SLRDAMAGES[t, r] = 100.0 * p.slrmultiplier[r] * (p.TOTALSLR[t-1] * p.slrdamlinear[r] + p.TOTALSLR[t-1]^2 * p.slrdamquadratic[r]) * (p.YGROSS[t-1, r] / p.YGROSS[TimestepIndex(1), r])^(1 / p.slrelasticity[r]) end end end diff --git a/src/components/welfare_component.jl b/src/components/welfare_component.jl index a828069..f7da19b 100644 --- a/src/components/welfare_component.jl +++ b/src/components/welfare_component.jl @@ -19,35 +19,35 @@ #Define function for PERIODU #NEED TO ADD IF STATEMENT LIKE IN JUMP MODEL OR IS THAT ONLY ISSUES WHEN ELASMU = 1.0? for r in d.regions - if p.elasmu[r]==1. - v.PERIODU[t,r] = log(p.CPC[t,r]) * p.alpha[t,r] + if p.elasmu[r] == 1.0 + v.PERIODU[t, r] = log(p.CPC[t, r]) * p.alpha[t, r] else - v.PERIODU[t,r] = ((1. / (1. - p.elasmu[r])) * (p.CPC[t,r])^(1. - p.elasmu[r]) + 1.) * p.alpha[t,r] + v.PERIODU[t, r] = ((1.0 / (1.0 - p.elasmu[r])) * (p.CPC[t, r])^(1.0 - p.elasmu[r]) + 1.0) * p.alpha[t, r] end end #Define function for CEMUTOTPER for r in d.regions if t.t != 60 - v.CEMUTOTPER[t,r] = v.PERIODU[t,r] * p.l[t,r] * p.rr[t,r] + v.CEMUTOTPER[t, r] = v.PERIODU[t, r] * p.l[t, r] * p.rr[t, r] else - v.CEMUTOTPER[t,r] = v.PERIODU[t,r] * p.l[t,r] * p.rr[t,r] / (1. - ((p.rr[t-1,r] / (1. + 0.015)^10) / p.rr[t-1,r])) + v.CEMUTOTPER[t, r] = v.PERIODU[t, r] * p.l[t, r] * p.rr[t, r] / (1.0 - ((p.rr[t-1, r] / (1.0 + 0.015)^10) / p.rr[t-1, r])) end end #Define function for REGCUMCEMUTOTPER for r in d.regions if is_first(t) - v.REGCUMCEMUTOTPER[t,r] = v.CEMUTOTPER[t,r] + v.REGCUMCEMUTOTPER[t, r] = v.CEMUTOTPER[t, r] else - v.REGCUMCEMUTOTPER[t,r] = v.REGCUMCEMUTOTPER[t-1, r] + v.CEMUTOTPER[t,r] + v.REGCUMCEMUTOTPER[t, r] = v.REGCUMCEMUTOTPER[t-1, r] + v.CEMUTOTPER[t, r] end end if t.t == 60 #Define function for REGUTILITY for r in d.regions - v.REGUTILITY[r] = 10 * p.scale1[r] * v.REGCUMCEMUTOTPER[t,r] + p.scale2[r] + v.REGUTILITY[r] = 10 * p.scale1[r] * v.REGCUMCEMUTOTPER[t, r] + p.scale2[r] end #Define function for UTILITY v.UTILITY = sum(v.REGUTILITY[:]) diff --git a/src/helpers.jl b/src/helpers.jl index 754d2fb..9218bec 100644 --- a/src/helpers.jl +++ b/src/helpers.jl @@ -10,22 +10,22 @@ end #Function to read a single parameter value from original RICE 2010 model. function getparam_single(f, range::AbstractString, regions) - vals= Array{Float64}(undef, length(regions)) - for (i,r) = enumerate(regions) - data=f[r][range] - vals[i]=data[1] + vals = Array{Float64}(undef, length(regions)) + for (i, r) = enumerate(regions) + data = f[r][range] + vals[i] = data[1] end return vals end #Function to read a time series of parameter values from original RICE 2010 model. function getparam_timeseries(f, range::AbstractString, regions, T) - vals= Array{Float64}(undef, T, length(regions)) - for (i,r) = enumerate(regions) - data=f[r][range] - for n=1:T - vals[n,i] = data[n] + vals = Array{Float64}(undef, T, length(regions)) + for (i, r) = enumerate(regions) + data = f[r][range] + for n = 1:T + vals[n, i] = data[n] end end return vals -end \ No newline at end of file +end diff --git a/src/marginaldamage.jl b/src/marginaldamage.jl index 7f61c4a..e6d9ca4 100644 --- a/src/marginaldamage.jl +++ b/src/marginaldamage.jl @@ -5,12 +5,12 @@ Computes the social cost of CO2 for an emissions pulse in `year` for the provide If no model is provided, the default model from MimiRICE2010.get_model() is used. Constant discounting is used from the specified pure rate of time preference `prtp`. """ -function compute_scc(m::Model=get_model(); year::Union{Int, Nothing} = nothing, last_year::Int = model_years[end], prtp::Float64 = 0.015, eta::Float64=1.5) +function compute_scc(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta::Float64=1.5) year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2015)`.") : nothing !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing !(year in model_years[1]:10:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing - mm = get_marginal_model(m; year = year) + mm = get_marginal_model(m; year=year) return _compute_scc(mm, year=year, last_year=last_year, prtp=prtp, eta=eta) end @@ -23,15 +23,15 @@ Computes the social cost of CO2 for an emissions pulse in `year` for the provide If no model is provided, the default model from MimiRICE2010.get_model() is used. Constant discounting is used from the specified pure rate of time preference `prtp`. """ -function compute_scc_mm(m::Model=get_model(); year::Union{Int, Nothing} = nothing, last_year::Int = model_years[end], prtp::Float64 = 0.015, eta::Float64=1.5) +function compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, last_year::Int=model_years[end], prtp::Float64=0.015, eta::Float64=1.5) year === nothing ? error("Must specify an emission year. Try `compute_scc_mm(m, year=2015)`.") : nothing !(last_year in model_years) ? error("Invlaid value of $last_year for last_year. last_year must be within the model's time index $model_years.") : nothing !(year in model_years[1]:10:last_year) ? error("Cannot compute the scc for year $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing - mm = get_marginal_model(m; year = year) + mm = get_marginal_model(m; year=year) scc = _compute_scc(mm; year=year, last_year=last_year, prtp=prtp, eta=eta) - - return (scc = scc, mm = mm) + + return (scc=scc, mm=mm) end # helper function for computing SCC from a MarginalModel, not to be exported or advertised to users @@ -39,16 +39,16 @@ function _compute_scc(mm::MarginalModel; year::Int, last_year::Int, prtp::Float6 ntimesteps = findfirst(isequal(last_year), model_years) # Will run through the timestep of the specified last_year run(mm, ntimesteps=ntimesteps) - marginal_damages = -1 * mm[:neteconomy, :C][1:ntimesteps, :] * 10.0^12 * 12/44 # convert from trillion $/ton C to $/ton CO2; multiply by -1 to get positive value for damages - global_marginal_damages = dropdims(sum(marginal_damages, dims = 2), dims=2) + marginal_damages = -1 * mm[:neteconomy, :C][1:ntimesteps, :] * 10.0^12 * 12 / 44 # convert from trillion $/ton C to $/ton CO2; multiply by -1 to get positive value for damages + global_marginal_damages = dropdims(sum(marginal_damages, dims=2), dims=2) - global_c = dropdims(sum(mm.base[:neteconomy, :C], dims = 2), dims=2) - global_pop = dropdims(sum(mm.base[:neteconomy, :l], dims = 2), dims=2) + global_c = dropdims(sum(mm.base[:neteconomy, :C], dims=2), dims=2) + global_pop = dropdims(sum(mm.base[:neteconomy, :l], dims=2), dims=2) global_cpc = 1000 .* global_c ./ global_pop year_index = findfirst(isequal(year), model_years) - df = [zeros(year_index-1)..., ((global_cpc[year_index]/global_cpc[i])^eta * 1/(1+prtp)^(t-year) for (i,t) in enumerate(model_years) if year<=t<=last_year)...] + df = [zeros(year_index - 1)..., ((global_cpc[year_index] / global_cpc[i])^eta * 1 / (1 + prtp)^(t - year) for (i, t) in enumerate(model_years) if year <= t <= last_year)...] scc = sum(df .* marginal_damages * 10) # currently implemented as a 10year step function; so each timestep of discounted marginal damages is multiplied by 10 return scc end @@ -59,7 +59,7 @@ get_marginal_model(m::Model = get_model(); year::Int = nothing) Creates a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`. If no Model m is provided, the default model from MimiDICE2010.get_model() is used as the base model. """ -function get_marginal_model(m::Model=get_model(); year::Union{Int, Nothing} = nothing) +function get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing) year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2015)`.") : nothing !(year in model_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $(model_years[1]):10:$last_year.") : nothing @@ -72,7 +72,7 @@ end """ Adds a marginal emission component to year m which adds 1Gt of additional C emissions per year for ten years starting in the specified `year`. """ -function add_marginal_emissions!(m::Model, year::Int) +function add_marginal_emissions!(m::Model, year::Int) add_comp!(m, Mimi.adder, :marginalemission, before=:co2cycle) time = Mimi.dimension(m, :time) @@ -87,8 +87,8 @@ end # OLD previously used marginal damages functions -function getmarginal_rice_models(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")) - +function getmarginal_rice_models(; emissionyear=2005, datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")) + RICE = get_model() run(RICE) @@ -102,9 +102,9 @@ function getmarginal_rice_models(;emissionyear=2005,datafile=joinpath(@__DIR__, addem = zeros(length(time)) addem[time[emissionyear]] = 1.0 - set_param!(m2,:marginalemission,:add,addem) - connect_param!(m2,:marginalemission,:input,:emissions,:E) - connect_param!(m2, :co2cycle,:E,:marginalemission,:output) + set_param!(m2, :marginalemission, :add, addem) + connect_param!(m2, :marginalemission, :input, :emissions, :E) + connect_param!(m2, :co2cycle, :E, :marginalemission, :output) run(m1) run(m2) @@ -114,18 +114,18 @@ end # This function returns a matrix of marginal damages per one t of carbon emission in the # emissionyear parameter year. -function getmarginaldamages_rice(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")) +function getmarginaldamages_rice(; emissionyear=2005, datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")) m1, m2 = getmarginal_rice_models(emissionyear=emissionyear, datafile=datafile) run(m1) run(m2) - damage1 = m1[:damages,:DAMAGES] - damage2 = m2[:damages,:DAMAGES] + damage1 = m1[:damages, :DAMAGES] + damage2 = m2[:damages, :DAMAGES] # Calculate the marginal damage between run 1 and 2 for each # year/region - marginaldamage = (damage2.-damage1) .* 1e12 / 1e9 / 10. + marginaldamage = (damage2 .- damage1) .* 1e12 / 1e9 / 10.0 return marginaldamage end diff --git a/src/parameters.jl b/src/parameters.jl index a8cf292..55a6e58 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -16,14 +16,14 @@ function getrice2010parameters(filename) ifopt = true # Indicator where optimized is 1 and base is 0 # Preferences - p[:elasmu] = getparam_single(f, "B18:B18", regions) # Elasticity of MU of consumption + p[:elasmu] = getparam_single(f, "B18:B18", regions) # Elasticity of MU of consumption # p[:prstp] = getparam_single(f, "B15:B15", regions) # Rate of Social Time Preference; not currently used in constructrice # Population and technology # Capital elasticity in production function p[:gama] = 0.300 - p[:dk] = getparam_single(f, "B8:B8", regions) # Depreciation rate on capital (per year) + p[:dk] = getparam_single(f, "B8:B8", regions) # Depreciation rate on capital (per year) p[:k0] = getparam_single(f, "B11:B11", regions) #Initial capital # p[:miu0] = getparam_single(f, "B103:B103", regions) # Initial emissions control rate for base case 2010; not currently used in constructrice # p[:miubase] = getparam_timeseries(f, "B103:BI103", regions, T) # Optimized emission control rate results from RICE2010 (base case); duplicate line to p[:MIU] @@ -31,29 +31,29 @@ function getrice2010parameters(filename) # Carbon cycle # Initial Conditions - p[:mat0] = 787.0 # Initial Concentration in atmosphere 2010 (GtC) - p[:mat1] = 829.0 - p[:mu0] = 1600. # Initial Concentration in upper strata 2010 (GtC) - p[:ml0] = 10010. # Initial Concentration in lower strata 2010 (GtC) + p[:mat0] = 787.0 # Initial Concentration in atmosphere 2010 (GtC) + p[:mat1] = 829.0 + p[:mu0] = 1600.0 # Initial Concentration in upper strata 2010 (GtC) + p[:ml0] = 10010.0 # Initial Concentration in lower strata 2010 (GtC) # Carbon cycle transition matrix # Flow paramaters - p[:b12] = 12.0/100 # Carbon cycle transition matrix atmosphere to shallow ocean - p[:b23] = 0.5/100 # Carbon cycle transition matrix shallow to deep ocean + p[:b12] = 12.0 / 100 # Carbon cycle transition matrix atmosphere to shallow ocean + p[:b23] = 0.5 / 100 # Carbon cycle transition matrix shallow to deep ocean # Parameters for long-run consistency of carbon cycle - p[:b11] = 88.0/100 # Carbon cycle transition matrix atmosphere to atmosphere - p[:b21] = 4.704/100 # Carbon cycle transition matrix biosphere/shallow oceans to atmosphere - p[:b22] = 94.796/100 # Carbon cycle transition matrix shallow ocean to shallow oceans - p[:b32] = 0.075/100 # Carbon cycle transition matrix deep ocean to shallow ocean - p[:b33] = 99.925/100 # Carbon cycle transition matrix deep ocean to deep oceans + p[:b11] = 88.0 / 100 # Carbon cycle transition matrix atmosphere to atmosphere + p[:b21] = 4.704 / 100 # Carbon cycle transition matrix biosphere/shallow oceans to atmosphere + p[:b22] = 94.796 / 100 # Carbon cycle transition matrix shallow ocean to shallow oceans + p[:b32] = 0.075 / 100 # Carbon cycle transition matrix deep ocean to shallow ocean + p[:b33] = 99.925 / 100 # Carbon cycle transition matrix deep ocean to deep oceans # Climate model parameters p[:t2xco2] = 3.2 # Equilibrium temp impact (oC per doubling CO2) fex0 = 0.83 # 2010 forcings of non-CO2 GHG (Wm-2) fex1 = 0.3 # 2100 forcings of non-CO2 GHG (Wm-2) - p[:tocean0] = .0068 # Initial lower stratum temp change (C from 1900) + p[:tocean0] = 0.0068 # Initial lower stratum temp change (C from 1900) p[:tatm0] = 0.83 # Initial atmospheric temp change 2005 (C from 1900) p[:tatm1] = 0.98 # Initial atmospheric temp change 2015 (C from 1900) @@ -87,7 +87,7 @@ function getrice2010parameters(filename) # Additive scaling coefficient (combines two additive scaling coefficients from RICE for calculating utility with welfare weights) scale2 = Array{Float64}(undef, length(regions)) - for (i,r) in enumerate(regions) + for (i, r) in enumerate(regions) data = f[r]["B53:C53"] scale2[i] = data[1] - data[2] end @@ -106,14 +106,14 @@ function getrice2010parameters(filename) # Global Emissions from Land Use Change (Sum of regional emissions for land use change in RICE model) etree = Array{Float64}(undef, T) for i = 1:T - etree[i] = sum(regtree[i,:]) + etree[i] = sum(regtree[i, :]) end p[:etree] = etree # Exogenous forcing for other greenhouse gases - forcoth = Array{Float64}(undef, 60) + forcoth = Array{Float64}(undef, 60) data = f["Global"]["B21:BI21"] - for i=1:T + for i = 1:T forcoth[i] = data[i] end p[:forcoth] = forcoth diff --git a/test/runtests.jl b/test/runtests.jl index 0725e31..6442ec7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,16 +10,16 @@ run(m) parameter_filename = joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm") -f=readxlsx(parameter_filename) +f = readxlsx(parameter_filename) regions = ["US", "EU", "Japan", "Russia", "Eurasia", "China", "India", "MidEast", "Africa", "LatAm", "OHI", "OthAsia"] #Function to get true values form Rice2010 Excel function Truth(range::AbstractString) - true_vals=Array{Float64}(undef, 60, length(regions)) - for (i,r) = enumerate(regions) - data=f[r][range] - for n=1:60 - true_vals[n,i] = data[n] + true_vals = Array{Float64}(undef, 60, length(regions)) + for (i, r) = enumerate(regions) + data = f[r][range] + for n = 1:60 + true_vals[n, i] = data[n] end end return true_vals @@ -30,134 +30,134 @@ Precision = 1.0e-10 @testset "MimiRICE2010" begin -#------------------------------------------------------------------------------ -# 1. Run tests on the whole model -#------------------------------------------------------------------------------ + #------------------------------------------------------------------------------ + # 1. Run tests on the whole model + #------------------------------------------------------------------------------ -@testset "mimi-rice-2010-model" begin + @testset "mimi-rice-2010-model" begin -# TATM Test (temperature increase) -True_TATM = vec(f["Global"]["B70:BI70"]) -@test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0. atol=Precision + # TATM Test (temperature increase) + True_TATM = vec(f["Global"]["B70:BI70"]) + @test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0.0 atol = Precision -# MAT Test (carbon concentration atmosphere) -True_MAT = vec(f["Global"]["B59:BI59"]) -@test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0. atol=Precision + # MAT Test (carbon concentration atmosphere) + True_MAT = vec(f["Global"]["B59:BI59"]) + @test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0.0 atol = Precision -# DAMFRAC Test (damages fraction) -True_DAMFRAC = Truth("B63:BI63") -@test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0. atol=Precision + # DAMFRAC Test (damages fraction) + True_DAMFRAC = Truth("B63:BI63") + @test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0.0 atol = Precision -# DAMAGES Test (damages $) -True_DAMAGES = Truth("B64:BI64") -@test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0. atol=Precision + # DAMAGES Test (damages $) + True_DAMAGES = Truth("B64:BI64") + @test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0.0 atol = Precision -# E Test (emissions) -True_E = vec(f["Global"]["B55:BI55"]) -@test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0. atol=Precision + # E Test (emissions) + True_E = vec(f["Global"]["B55:BI55"]) + @test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0.0 atol = Precision -# YGROSS Test (gross output) -True_YGROSS = Truth("B61:BI61") -@test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0. atol=Precision + # YGROSS Test (gross output) + True_YGROSS = Truth("B61:BI61") + @test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0.0 atol = Precision -# CPC Test (per capita consumption) -True_CPC = Truth("B88:BI88") -@test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0. atol=Precision + # CPC Test (per capita consumption) + True_CPC = Truth("B88:BI88") + @test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0.0 atol = Precision -# FORC Test (radiative forcing) -True_FORC = vec(f["Global"]["B71:BI71"]) -@test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0. atol=Precision + # FORC Test (radiative forcing) + True_FORC = vec(f["Global"]["B71:BI71"]) + @test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0.0 atol = Precision -# TOTALSLR Test (total sea level rise) -True_TOTALSLR = vec(f["SLR"]["B62:BI62"]) -@test maximum(abs, m[:sealevelrise, :TOTALSLR] .- True_TOTALSLR) ≈ 0. atol=Precision + # TOTALSLR Test (total sea level rise) + True_TOTALSLR = vec(f["SLR"]["B62:BI62"]) + @test maximum(abs, m[:sealevelrise, :TOTALSLR] .- True_TOTALSLR) ≈ 0.0 atol = Precision -# SLRDAMAGES Test (damages from sea level rise) -True_SLRDAMAGES = Truth("B50:BI50") -@test maximum(abs, m[:sealeveldamages, :SLRDAMAGES] .- True_SLRDAMAGES) ≈ 0. atol=Precision + # SLRDAMAGES Test (damages from sea level rise) + True_SLRDAMAGES = Truth("B50:BI50") + @test maximum(abs, m[:sealeveldamages, :SLRDAMAGES] .- True_SLRDAMAGES) ≈ 0.0 atol = Precision -True_PERIODUTILITY = Truth("B94:BI94") -@test maximum(abs, m[:welfare, :CEMUTOTPER] .- True_PERIODUTILITY) ≈ 0. atol=Precision + True_PERIODUTILITY = Truth("B94:BI94") + @test maximum(abs, m[:welfare, :CEMUTOTPER] .- True_PERIODUTILITY) ≈ 0.0 atol = Precision -True_UTILITY = f["Global"]["B77:B77"] -@test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0. atol=Precision + True_UTILITY = f["Global"]["B77:B77"] + @test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0.0 atol = Precision -end #mimi-rice-2010-model testset + end #mimi-rice-2010-model testset -#------------------------------------------------------------------------------ -# 2. Run tests to make sure integration version (Mimi v0.5.0) -# values match Mimi 0.4.0 values -#------------------------------------------------------------------------------ + #------------------------------------------------------------------------------ + # 2. Run tests to make sure integration version (Mimi v0.5.0) + # values match Mimi 0.4.0 values + #------------------------------------------------------------------------------ -@testset "mimi-rice-2010-integration" begin + @testset "mimi-rice-2010-integration" begin -nullvalue = -999.999 + nullvalue = -999.999 -for c in map(nameof, Mimi.compdefs(m)), v in Mimi.variable_names(m, c) + for c in map(nameof, Mimi.compdefs(m)), v in Mimi.variable_names(m, c) - #load data for comparison - filepath = joinpath(@__DIR__, "..", "data", "validation_data_v040", "$c-$v.csv") - results = m[c, v] + #load data for comparison + filepath = joinpath(@__DIR__, "..", "data", "validation_data_v040", "$c-$v.csv") + results = m[c, v] - df = load(filepath) |> DataFrame - if typeof(results) <: Number - validation_results = df[1,1] + df = load(filepath) |> DataFrame + if typeof(results) <: Number + validation_results = df[1, 1] - else - validation_results = Matrix(df) + else + validation_results = Matrix(df) - #remove NaNs and Missings - results[ismissing.(results)] .= nullvalue - results[isnan.(results)] .= nullvalue - validation_results[ismissing.(validation_results)] .= nullvalue - validation_results[isnan.(validation_results)] .= nullvalue + #remove NaNs and Missings + results[ismissing.(results)] .= nullvalue + results[isnan.(results)] .= nullvalue + validation_results[ismissing.(validation_results)] .= nullvalue + validation_results[isnan.(validation_results)] .= nullvalue - #match dimensions - if size(validation_results,1) == 1 - validation_results = validation_results' - end - end + #match dimensions + if size(validation_results, 1) == 1 + validation_results = validation_results' + end + end - @test results ≈ validation_results atol = Precision + @test results ≈ validation_results atol = Precision -end #for loop + end #for loop -end #mimi-rice-2010-integration testset + end #mimi-rice-2010-integration testset -@testset "Standard API functions" begin + @testset "Standard API functions" begin -m = MimiRICE2010.get_model() -run(m) + m = MimiRICE2010.get_model() + run(m) -# Test the errors -@test_throws ErrorException MimiRICE2010.compute_scc() # test that it errors if you don't specify a year -@test_throws ErrorException MimiRICE2010.compute_scc(year=2020) # test that it errors if the year isn't in the time index -@test_throws ErrorException MimiRICE2010.compute_scc(last_year=2300) # test that it errors if the last_year isn't in the time index -@test_throws ErrorException MimiRICE2010.compute_scc(year=2105, last_year=2100) # test that it errors if the year is after last_year + # Test the errors + @test_throws ErrorException MimiRICE2010.compute_scc() # test that it errors if you don't specify a year + @test_throws ErrorException MimiRICE2010.compute_scc(year=2020) # test that it errors if the year isn't in the time index + @test_throws ErrorException MimiRICE2010.compute_scc(last_year=2300) # test that it errors if the last_year isn't in the time index + @test_throws ErrorException MimiRICE2010.compute_scc(year=2105, last_year=2100) # test that it errors if the year is after last_year -# Test the SCC -scc1 = MimiRICE2010.compute_scc(year=2015) -@test scc1 isa Float64 + # Test the SCC + scc1 = MimiRICE2010.compute_scc(year=2015) + @test scc1 isa Float64 -# Test that it's smaller with a shorter horizon -scc2 = MimiRICE2010.compute_scc(year=2015, last_year=2295) -@test scc2 < scc1 + # Test that it's smaller with a shorter horizon + scc2 = MimiRICE2010.compute_scc(year=2015, last_year=2295) + @test scc2 < scc1 -# Test that it's smaller with a smaller prtp -scc3 = MimiRICE2010.compute_scc(year=2015, last_year=2295, prtp=0.02) -@test scc3 < scc2 + # Test that it's smaller with a smaller prtp + scc3 = MimiRICE2010.compute_scc(year=2015, last_year=2295, prtp=0.02) + @test scc3 < scc2 -# Test with a modified model -m = MimiRICE2010.get_model() -update_param!(m, :t2xco2, 5) -scc4 = MimiRICE2010.compute_scc(m, year=2015) -@test scc4 > scc1 # Test that a higher value of climate sensitivty makes the SCC bigger + # Test with a modified model + m = MimiRICE2010.get_model() + update_param!(m, :t2xco2, 5) + scc4 = MimiRICE2010.compute_scc(m, year=2015) + @test scc4 > scc1 # Test that a higher value of climate sensitivty makes the SCC bigger -# Test compute_scc_mm -result = MimiRICE2010.compute_scc_mm(year=2035) -@test result.scc isa Float64 -@test result.mm isa Mimi.MarginalModel + # Test compute_scc_mm + result = MimiRICE2010.compute_scc_mm(year=2035) + @test result.scc isa Float64 + @test result.mm isa Mimi.MarginalModel -end + end end #mimi-rice-2010 testset