# 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)