World-Earth System Model

#xppaut version of world-earth model.  To run the model, download and install xppaut
# from www.math.pitt.edu/~bard/xpp/xpp.html.  Otherwise, porting to Julia or Python
# should be straightforward.  Note: because there are both quite fast and quite slow
# dynamics in the model, a stiff ode integrator is advised (required).

# ===================Human population subsystem=======================================

# Human population is in billions.  The parameters were fitted
# using UN population scenarios.
# Note: inter-regional migration is set to zero in this model version.

par Kp1=1.5,Kp2=9.7,r1=0.038,r2=0.042

dP1/dt = r1*P1*(1-P1/Kp1)
dP2/dt = r2*P2*(1-P2/Kp2)

#===============End of Human population subsystem===============================


#================Economic Subsystem================================================

# First some parameter settings

# total factor productivities for backstop technology

number a1b=1, a2b=1

#Total factor productivities for industrial technology

par a1=2.7,a2=1.7

# capital elasticities of output for each region:

number alpha1=0.5, alpha2=0.5

# Then by definition, for Cobb-Douglass, constant returns to scale

beta1=1-alpha1
beta2=1-alpha2   

# social foundation for consumption. Note, backstop technology percapita output
# is 1. If min_C is bigger than one, the population can't survive on the backstop
# technology.  

par min_C=0.7

# Capital depreciation rates (delta)

number delta1=0.05, delta2=0.05

# Now we can write the output for industrial production
# note the max function just prevents K from going negative
# purely a computational issue.
 
Y1i = a1*(max(0,K1)**alpha1)*(P1**beta1)
Y2i = a2*(max(0,K2)**alpha2)*(P2**beta2)

# Backstop technology.  This allows K to be completely destroyed,
# but people can mobilize old, simple technology. We assume this
# is fully labor based (same as setting alpha =1 in CB technology.

Y1b = a1b*P1 
Y2b = a2b*P2

#total domestic output: society chooses most productive technology
# that is, "industrializes".  First we need a threshold function.

# Threshold function thr(x) = 1, x>0, 0, x<0 Use a rational power
# function instead and force to be zero for negative argument
# This acts like a 'smooth' max(x) function.
# first define general form

thr2(x,b,s) = if(x>0)then(x**s/(b**s + x**s))else(0)

# Then set for actual parameters, b and S

thr(x) = thr2(x,0.0001,5)

#Switching functions (this generates a "smooth" ramp function):

Y1ir = thr(Y1i - Y1b)*(Y1i-Y1b)
Y2ir = thr(Y2i - Y2b)*(Y2i-Y2b)

#Then total output can be written as

Y1 = Y1b +  Y1ir
Y2 = Y2b +  Y2ir

# Note - output is measured in trillions of $$.

# Carbon intensity of economic activites. 

# baseline emssions. eb < e1, e2
# express eb as a ratio of e1 for easy interpretation
# that is, emission from background tehcnology are
# for example, an order of magnitude smaller than
# industrial tech.

par eb1r = 0.1

#then

eb = e1*eb1r

# then "realized" emissions depend on which technology is
# employed.

e1r = thr(Y1b-Y1i)*eb + thr(Y1i - Y1b)*e1
e2r = thr(Y2b-Y2i)*eb + thr(Y2i - Y2b)*e2

#trade (transfers)
# It is likely that trade depends on local information
# and some desire for some bundle of goods from the other region.... Not sure
# how to model this easily.
#next step - make trade a function of min consumption.  Also, do you make trade
#decisions before or after investment decisions?


trade12 = tau1*thr(Y1 - P1*min_C)*(Y1 - P1*min_C)
trade21 = tau2*thr(Y2 - P2*min_C)*(Y2 - P2*min_C)

# parameters for trade/transfer rate (0 <= tau <= 1)
# set defaults to zero, and vary in model analysis.
par tau1=0,tau2=0

# post-trade budgets: 
Z1 = Y1 - trade12 + trade21
Z2 = Y2 + trade12 - trade21

# savings fractions as a function of disposable income

# disposable income  (proportion above subsistence) note, min_C is per capita:

DispI1=thr(Z1 - P1*min_C)*(Z1 - P1*min_C)
DispI2= thr(Z2 - P2*min_C)*(Z2 - P2*min_C)

#savings rates
par s1=0.25, s2=0.21

#private investment proportion of disposable income

I1 = s1*DispI1
I2 = s2*DispI2

#==================End of Economic Subsystem model========================



#===============Earth System model ========================================

# The main issue for the Earth system is to estimate the emmision
# rates, e1 and e2 and uptake rate, u, and how the global externality
# generates shocks to infrastructure. 

# First estimate emmision rates per gdp based on the current situation
# e = average of 0.2 kg/$gdp and there is 1 ppm atmoshperic carbon 
# concentration per 2.1 Pg of carbon.

# then e = 0.2 kg/$ =  0.2 kg x 10^12/$ x 10^12  = 0.2 Pg/TR$ * 1 ppm/2 Pg  = 0.1 ppm/TR$.
# So if  we measure GDP in units of trillions of dollars (i.e. the model output of 5 = 5 TR$)
# and G in hundereds of ppm, e.g. 4 = 400 ppm,  then e = 0.001 100 ppm/TR$.   That is,
# e1 and e2 are on the order of 0.001.

# A recalculation based on 40GT at global np = 80 trillion -> 0.5 gigatons/trillion $.
#Nasa data -> 127 ppm/GT (easy calculation) 1ppm = 2.12 GT c  x 3.66 t co2/c = 7.76 gt co2.
# = 0.127 ppm/gt.  So we have 0.5 GT/TR$ * 0.127 ppm/gt = 0.0635 ppm/TR$.
#So we are measuring  ppm in hecto-ppm or Hppm.
#Thus, e1 = 0.0635 ppm/TR$ * 1 Hppm/100 ppm = 0.000635 Hppm/TR$

# these would be the parameter settings if the model is run without decarbonization
# low par e1=0.0001, e2=0.0001, medium e1=e2=0.0004.
#Note with decarbonization these parameter setting will become initial conditions
# i.e. the baseline from which the decarbonization process begins.


### Decarbonization

# e1 and e2 become functions of time or are differential
#  equations instead of being parameters.
# de1/dt = -r_d1(function of income?)*e1
# de2/dt = -r_d2()*e2

# So, e1 and e2 obey simple odes with exponential decay. z is a Markov
# variable that turns on decarobonization (represents a policy change
# to begin decarbonization) at a certian level of the global
# externality. It is zero until G>Tg (Tg = threshold to begin
# decarbonization.

de1/dt = -z*re1*e1
de2/dt = -z*re2*e2

markov z 2
{0} {if(G<Tg)then(0)else(1)}
{0} {0}

init z=0

# where re1 is the decarbonization rate (10% per year).

par  re1=0.1

#and the decarbonization rate in region 2 depends on inequality
# turn on and off inequality with 'ineq'


#Turn off and on whehter inequality affects decarbonization

par ineqon=0

ineq=Z2*P1/(P2*Z1)

re2= re1*(1-ineqon + ineqon*ineq)

# set the level of G at which decarbonization begins.  

par Tg=20

# note initial conditions are the baseline estimates of carbon
# intensity of the present industrial complex.

init e1=0.0004,e2=0.0004

#What is u=?.  Similar caluculation - from
#https://science.sciencemag.org/content/363/6432/1193/tab-figures-data, ocean and
# land have absored 280 Pg over 200 years or around 140 ppm over that period, so
# u = 1.4/200*G. Take average G just to ballpark it of say 3, then u = 1.4/600 = 0.0025

par u=0.0025

# shock severity coefficients (note sigma has to be much larger for
#noisy shocks.

par sigma1=0.03, sigma21r=1

sigma2=sigma1*sigma21r


# global irreversibility, note threshold at Go=4.5 seems empirically
# grounded.  Vary in analysis.

par alpha=0.1, Go=20

#weiner processes for global externality

weiner w1,w2

# set levels below which shocks have no impact

par s1min=10,s2min=10

# turn on and off stochastic shocks. If discon=0,
# damages increase smoothly with G as it surpasses
# G1.  If discon=1, damages are stochasitic.

par discOn=0

par G1=5
 
#Damage functions

dmg1 =  sigma1*((1-discOn)*exp(G-G1) + discon*thr(G*(w1-s1min)))
dmg2 =  sigma2*((1-discOn)*exp(G-G1) + discon*thr(G*(w1-s2min)))


#================End of Earth system model===================


#=================Full ODE Model============================

# Finally, the differential equations...... for 5 stocks of
# built and natural infrastrcutures

dK1/dt = I1 - delta1*K1 - dmg1*K1
dK2/dt  = I2 - delta2*K2 - dmg2*K2
dG/dt  = e1r*Y1 + e2r*Y2 -u* (G-2.8) + alpha*thr(G-Go)

# set default initial conditions:

init K1=0, K2=0,P1=0.24, P2=0.24

#Note, G = 2.8 at start of industrial revolution....

init G=2.8

#=============End ODE Model=========================



#set a bunch of auxiliary variables for plotting.....

#set scale on plot for atmospheric carbon...

par carbscale=10

aux atmcarb = G*carbscale
aux gdp1 = Y1
aux gdp1pc = Y1/P1
aux gdp2 = Y2
aux gdp2pc = Y2/P2
aux gni1 = Z1
aux gni1pc = Z1/P1
aux gni2 = Z2
aux gni2pc = Z2/P2
aux gnitot = Z1 + Z2
aux gdp1i = Y1i
aux gdp1ir = Y1ir
aux gdp1b = Y1b
aux gdp2i = Y2i
aux gdp2ir = Y2ir 
aux gdp2b = Y2b
aux t12=trade12
aux t21-trade21
aux incineq = ineq
aux em1 = e1r*Y1
aux em2 = e2r*Y2 
aux e1rr = e1r
aux e2rr = e2r

#these are settings for xppaut - total time of integration
#what to plot, and integration method (stiff).

@ total=500, bounds=10000,yp=gni1,xp=t
@ nplot=3 
@ xp2=t,yp2=gni2
@ xp3=t,yp3=atmcarb
@ xlo=0,xhi=500,ylo=0,yhi=100
@ meth=stiff

done