NASA HANDY Model in R

handy.R

library(deSolve)

handyMod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    # First calculate threshold below which famine begins
    # rho = minimum required consumption
    wThresh = rho * popC + rho * kappa * popE
    # Check to avoid dividing by zero then estimate consumption of
    # commons = cC and elites = cE
    if(wThresh > 0) {
      cC <- min(1, (w/wThresh)) * s * popC
      cE <- min(1, (w/wThresh)) * s * kappa * popE
    }
    else {
      cC <- s * popC
      cE <- s * kappa * popE
    }
    # Estimate death rates for commons and elites
    if(popC > 0) {
      alphaC <- alphaNormal + max(0,(1-(cC/(s*popC)))) * (alphaFamine - alphaNormal)
    }
    else {
      alphaC <- alphaNormal
    }
    if(popE > 0) {
      alphaE <- alphaNormal + max(0,(1-(cE/(s*popE)))) * (alphaFamine - alphaNormal)
    }
    else {
      alphaE <- alphaNormal
    }
    # Run solver on derivatives of population of commons, elites, natural resources, and
    # wealth.
    xpC <- betaC * popC - alphaC * popC
    xpE <- betaE * popE - alphaE * popE
    yp  <- gamma * y * (lambda - y) - delta * popC * y
    w  <- delta * popC * y - cC - cE
    return(list(c(xpC, xpE, yp, w)))
  })
}

# Death rates
alphaNormal = 0.01
alphaFamine = 0.07
# Birth rates
betaC = 0.03
betaE = 0.03

filename <- "/tmp/civ01"
# Rate of extraction of renewables
gamma = 0.01
# Environmental carrying capacity
lambda = 100.0
# Subsistence wage
#s = 0.0005
s = 0.004
# Differential between elites and commons
#kappa = 1.0
kappa = 2.0
# Minimum survival consumption below which famine begins
rho = 0.005
# Initial conditions
popC0 = 100
#popE0 = 0
popE0 = 10
y0  = lambda
w0  = 0
# Calculating delta (environmental extraction rate)
eta <- (alphaFamine - betaC)/(alphaFamine - alphaNormal)
phi <- popE0/popC0
# Maximized carrying capacity is achieved at extraction rate of lambda/2
#delta <- (2 * eta * s)/lambda
# A delta which supposedly will produce a stable state with an elite
# population
delta <- ((2 * eta * s)/lambda) * (1 + phi)
print(delta)

Pars <- c(alphaNormal=alphaNormal, alphaFamine=alphaFamine, betaC=betaC,
          betaE=betaE, gamma=gamma, lambda=lambda, kappa=kappa, rho=rho,
          delta=delta, s=s)
State <- c(popC=popC0, popE=popE0, y=y0, w=w0)
Time <- seq(0, 1000, by=1)

out <- as.data.frame(ode(func=handyMod, y=State, parms=Pars, times=Time))
t <- out[,1]

f <- paste(filename, ".png", sep="")
png(f, width=648, height=432)
par(mfrow=c(1,2))
matplot(out[,-1], type = "n", xlab = "Time(yrs)", ylab = "Pop.")
matlines(out[,2], type='l', lty=1, col=1)
matlines(out[,3], type='l', lty=1, col=2)
legend("right", legend=c("Pop Commons", "Pop Elites"), lty = 1, col = c(1,2))

matplot(out[,4],type='n', xlab="Time(yrs)", ylab="Eco-Dollars", ylim=c(0,150))
matlines(t,out[,4], type='l', lty=1, col=3)
matlines(t,out[,5], type='l', lty=1, col=4)
legend("topright", legend=c("Nature", "Wealth"), lty = 1, col = c(3,4))
dev.off()