Mdptoolbox Allen Model

library("MDPtoolbox", quietly = TRUE)
library("ggplot2", quietly = TRUE)

K <- 150 # state space limit
states <- 0:K # Vector of all possible states
actions <- states # Vector of actions: harvest


sigma_g = 0.1
p <- c(2, 100, 50)

f <- function (x, h){
  sapply(x, function(x) {
    x <- pmax(0, x - h)
    x * exp(p[1] * (1 - x/p[2]) * (x - p[3])/p[2])
  })
}

pdfn <- function(x, mu, sigma = sigma_g){
  dlnorm(x, log(mu), sdlog = sigma)
}

# Utility function
discount = 0.95
get_utility <- function(x,h) {
    pmin(x,h)
}

R <- outer(states, actions, get_utility)

transition_matrix <- function(states, actions, f, pdfn){
  # Initialize
  transition <- array(0, dim = c(length(states), length(states), length(actions)))
  
  K <- length(states)
  
  for (k in 1:length(states)) {
    for (i in 1:length(actions)) {
  
  # Calculate the transition state at the next step, given the 
  # current state k and action i (harvest H[i])
        nextpop <- f(states[k], actions[i])
        
        ## Population always extinct if this is negative. since multiplicitive shock z_t * f(n) < 0 for all f(n) < 0
        if(nextpop <= 0)
          transition[k, , i] <- c(1, rep(0, length(states) - 1))
    # Implement demographic stochasticity 
        else {
  
        # Cts distributions need long-tailed denominator as normalizing factor:
          fine_states <- seq(min(states), 10 * max(states), by = states[2] - states[1])
        N <- sum(pdfn(fine_states, nextpop))  
          transition[k, , i] <-pdfn(states, nextpop) / N
          
            # We need to correct this density for the final capping state ("Pile on boundary") (discrete or cts case)
          # this can be a tiny but negative value due to floating-point errors. so we take max(v,0) to avoid
        transition[k, K, i] <- max(1 - sum(transition[k, -K, i]), 0)
        }
    } 
  }
  transition
}

P <- transition_matrix(states, actions, f, pdfn)

Using toolbox

mdp_check(P = P, R = R)

[1] ""

mdp <- mdp_value_iteration(P, R, discount = discount, epsilon = 0.001, max_iter = 5e3, V0 = numeric(length(states)))

[1] "MDP Toolbox: iterations stopped, epsilon-optimal policy found"

plot(states, states - actions[mdp$policy],  xlab="Population size", ylab="Escapement")

Compare to Reed

From Reed (1979) we know that the optimal solution is a constant-escapement rule when the growth function in convex. Note that this condition is violated by the growth function with alternative stable states (Allen/Ricker-Allee model), resulting in a very different optimal policy:

\[f'(s^*) = 1/\alpha\]

For growth-rate function \(f\), where \(\alpha\) is the discount factor and \(s^*\) the stock size for the constant escapement. Analytic solutions are clearly possible for certain growth functions, but here I’ve just implemented a generic numerical solution.

fun <- function(x) - f(x,0) + x / discount
out <- optimize(f = fun, interval = c(0,K))
S_star <- out$minimum

exact_policy <- sapply(states, 
                       function(x) 
                        if(x < S_star) 0
                        else x - S_star)

plot(states, states - actions[mdp$policy],  xlab="Population size", ylab="Escapement")

# The difference between Bellman and the analytical solution is small:
lines(states, states - exact_policy)