value_of_information_allee.md

Calculating the value of information with highly nonlinear value function

Implements a numerical version of the SDP described in (Sethi et. al. 2005). Compute the optimal solution under different forms of uncertainty, but under strong allee dynamics

We consider the state dynamics

f <- RickerAllee
K <- 4
xT <- 2  # final value, also allee threshold
pars <- c(1.5, K, xT)
f

function (x, h, p) 
{
    sapply(x, function(x) {
        x <- max(0, x - h)
        x * exp(p[1] * (1 - x/p[2]) * (x - p[3])/p[2])
    })
}
<environment: namespace:pdgControl>

We consider a profits from fishing to be a function of harvest h and stock size x, \( \Pi(x,h) = h - \left( c_0 + c_1 \frac{h}{x} \right) \frac{h}{x} \), conditioned on \( h > x \) and \(x > 0 \),

price <- 1
c0 <- 0.01
c1 <- 0
profit <- profit_harvest(price = price, c0 = c0, 
    c1 = c1)

with price 1, c0 0.01 and c1 0.

xmin <- 0
xmax <- 1.5 * K
grid_n <- 100

We seek a harvest policy which maximizes the discounted profit from the fishery using a stochastic dynamic programming approach over a discrete grid of stock sizes from 0 to 6 on a grid of 100 points, and over an identical discrete grid of possible harvest values.

x_grid <- seq(xmin, xmax, length = grid_n)
h_grid <- x_grid

delta <- 0.05
xT <- 0
OptTime <- 25

We will determine the optimal solution over a 25 time step window with boundary condition for stock at 0 and discounting rate of 0.05.

Scenarios:

We use Monte Carlo integration over the noise processes to determine the transition matrix.

require(snowfall)
sfInit(parallel = TRUE, cpu = 16)

R Version:  R version 2.14.1 (2011-12-22) 

scenario <- function(policy_g, policy_m, policy_i) {
    
    z_g <- function() rlnorm(1, 0, policy_g)
    z_m <- function() rlnorm(1, 0, policy_m)
    z_i <- function() rlnorm(1, 0, policy_i)
    
    SDP_Mat <- SDP_by_simulation(f, pars, x_grid, h_grid, 
        z_g, z_m, z_i, reps = 20000)
    opt <- find_dp_optim(SDP_Mat, x_grid, h_grid, OptTime, 
        xT, profit, delta, reward = 0)
}

det <- scenario(0, 0, 0)

Library ggplot2 loaded.

g <- scenario(0.1, 0, 0)

Library ggplot2 loaded.

m <- scenario(0, 0.1, 0)

Library ggplot2 loaded.

i <- scenario(0, 0, 0.1)

Library ggplot2 loaded.

gm <- scenario(0.1, 0.1, 0)

Library ggplot2 loaded.

gi <- scenario(0.1, 0, 0.1)

Library ggplot2 loaded.

gmi <- scenario(0.1, 0.1, 0.1)

Library ggplot2 loaded.

det <- scenario(0.001, 0, 0)

Library ggplot2 loaded.

plots

require(reshape2)
policy <- melt(data.frame(stock = x_grid, det = det$D[, 
    1], g = g$D[, 1], m = m$D[, 1], gm = gm$D[, 1], gi = gi$D[, 
    1], gmi = gmi$D[, 1]), id = "stock")
ggplot(policy) + geom_point(aes(stock, stock - 
    x_grid[value], color = variable)) + ylab("escapement")

ggplot(policy) + geom_smooth(aes(stock, stock - 
    x_grid[value], color = variable)) + ylab("escapement")

value <- melt(data.frame(stock = x_grid, det = det$V, 
    g = g$V, m = m$V, gm = gm$V, gi = gi$V, gmi = gmi$V), 
    id = "stock")
ggplot(value) + geom_point(aes(stock, value, color = variable)) + 
    ylab("Net Present Value")

ggplot(value) + geom_smooth(aes(stock, value, 
    color = variable)) + ylab("Net Present Value")

Simulations

simulatereps <- function(opt, true_g, true_m, 
    true_i) {
    z_g <- function() rlnorm(1, 0, true_g)
    z_m <- function() rlnorm(1, 0, true_m)
    z_i <- function() rlnorm(1, 0, true_i)
    
    sims <- lapply(1:100, function(i) {
        ForwardSimulate(f, pars, x_grid, h_grid, x0 = K, 
            opt$D, z_g, z_m, z_i, profit)
    })
    
    # list(sims = sims, opt = opt, true_stochasticity =
    #   c(true_g, true_m, true_i))
    sims
}

All cases

policyfn <- list(det = det, g = g, m = m, i = i, 
    gm = gm, gi = gi, gmi = gmi)
noise <- list(s0 = c(0, 0, 0), sg = c(0.2, 0, 
    0), sm = c(0, 0.2, 0), si = c(0, 0, 0.2), sgm = c(0.2, 
    0.2, 0), sgi = c(0.2, 0, 0.2), sgmi = c(0.2, 0.2, 0.2))
allcases <- lapply(policyfn, function(policyfn_i) {
    lapply(noise, function(noise_i) {
        simulatereps(policyfn_i, noise_i[1], noise_i[2], 
            noise_i[3])
    })
})

sims <- unlist(allcases, recursive = FALSE)
dat <- melt(sims, id = names(sims[[1]][[1]]))
dt <- data.table(dat)
setnames(dt, c("L2", "L1"), c("reps", "uncertainty"))  # names are nice

Plots

ggplot(subset(dt, reps == 1)) + geom_line(aes(time, 
    fishstock)) + geom_line(aes(time, harvest), col = "darkgreen") + 
    facet_wrap(~uncertainty)

This plot summarizes the stock dynamics by visualizing the replicates.

p1 <- ggplot(dt)
p1 + geom_line(aes(time, fishstock, group = reps), 
    alpha = 0.1) + facet_wrap(~uncertainty)

ggplot(subset(dt, reps == 1)) + geom_line(aes(time, 
    profit)) + facet_wrap(~uncertainty)

profits <- dt[, sum(profit), by = c("reps", "uncertainty")]
ggplot(profits) + geom_histogram(aes(V1)) + facet_wrap(~uncertainty)

Summary statistics

means <- profits[, mean(V1), by = uncertainty]
sds <- profits[, sd(V1), by = uncertainty]

require(xtable)
uncertainties <- c("deter", "growth", "measure", 
    "implement", "growth+measure", "growth+implement", "all")
print(xtable(matrix(means$V1, nrow = 7, dimnames = list(uncertainties, 
    uncertainties)), caption = "Mean realized net present value over simulations"), 
    type = "html")

Mean realized net present value over simulations

	deter	growth	measure	implement	growth+measure	growth+implement	all
deter	10.92	10.91	10.92	10.86	10.88	10.64	10.68
growth	9.73	9.43	9.94	10.68	9.69	10.47	11.02
measure	5.41	6.10	5.81	5.75	5.84	6.33	6.52
implement	10.63	10.62	10.62	10.51	10.55	10.37	10.36
growth+measure	5.27	5.64	4.68	5.73	6.07	5.98	6.76
growth+implement	9.53	9.33	8.83	9.64	9.23	10.37	10.47
all	4.94	5.63	5.39	5.29	5.98	6.25	5.09

print(xtable(matrix(sds$V1, nrow = 7, dimnames = list(uncertainties, 
    uncertainties)), caption = "Standard deviation in realized net present value over simulations"), 
    type = "html")

Standard deviation in realized net present value over simulations

	deter	growth	measure	implement	growth+measure	growth+implement	all
deter	0.00	0.00	0.00	0.00	0.00	0.00	0.00
growth	5.08	4.06	4.46	4.04	4.83	4.58	4.69
measure	1.52	1.95	1.84	1.83	1.72	2.06	2.04
implement	0.35	0.33	0.36	0.39	0.36	0.37	0.34
growth+measure	2.30	2.87	1.80	3.00	2.86	2.86	3.40
growth+implement	4.45	4.33	4.26	4.41	4.33	4.39	4.41
all	1.94	2.84	2.32	2.31	3.04	3.35	2.46

References

Sethi G, Costello C, Fisher A, Hanemann M and Karp L (2005). “Fishery Management Under Multiple Uncertainty.” Journal of Environmental Economics And Management, 50. ISSN 00950696, http://dx.doi.org/10.1016/j.jeem.2004.11.005.