free_decrease.md

Policy Costs

• author Carl Boettiger, cboettig@gmail.com
• license: CC0

Setup the system

rm(list = ls())
require(pdgControl)
require(reshape2)
require(ggplot2)
require(data.table)

delta <- 0.05  # economic discounting rate
OptTime <- 50  # stopping time
gridsize <- 50  # gridsize (discretized population)
sigma_g <- 0.2  # Noise in population growth
sigma_m <- 0  # noise in stock assessment measurement
sigma_i <- 0  # noise in implementation of the quota
reward <- 0  # bonus for satisfying the boundary condition


## @knitr noise_dists
z_g <- function() rlnorm(1, 0, sigma_g)  # mean 1
z_m <- function() rlnorm(1, 0, sigma_m)  # mean 1
z_i <- function() rlnorm(1, 0, sigma_i)  # mean 1



f <- BevHolt  # Select the state equation
pars <- c(1.5, 0.05)  # parameters for the state equation
K <- (pars[1] - 1)/pars[2]  # Carrying capacity (for reference
xT <- 0  # boundary conditions
x0 <- K


## @knitr profit_
profit <- profit_harvest(price = 10, c0 = 30, c1 = 10)


## @knitr create_grid_
x_grid <- seq(0.01, 1.2 * K, length = gridsize)
h_grid <- seq(0.01, 0.8 * K, length = gridsize)


## @knitr reed_sdp
SDP_Mat <- determine_SDP_matrix(f, pars, x_grid, h_grid, sigma_g)
opt <- find_dp_optim(SDP_Mat, x_grid, h_grid, OptTime, xT, profit, 
    delta, reward = reward)

L1 <- function(c2) function(h, h_prev) c2 * abs(h - h_prev)
free_increase <- function(c2) function(h, h_prev) c2 * abs(min(h - 
    h_prev, 0))  # increasing harvest is free
free_decrease <- function(c2) function(h, h_prev) c2 * max(h - h_prev, 
    0)  # decreasing harvest is free
fixed <- function(c2) function(h, h_prev) c2
L2 <- function(c2) function(h, h_prev) c2 * (h - h_prev)^2

Solve the policy cost for the specified penalty function

c2 <- 4
penalty <- free_decrease(c2)
policycost <- optim_policy(SDP_Mat, x_grid, h_grid, OptTime, xT, 
    profit, delta, reward, penalty = penalty)
cache = FALSE

Simulate

Now we'll simulate 100 replicates of this stochastic process under the optimal harvest policy determined above. We use a modified simulation function that can simulate an alternate policy (the Reed optimum, where policy costs are zero, opt$D ) and a focal policy, policycost$D

sims <- lapply(1:100, function(i) simulate_optim(f, pars, x_grid, 
    h_grid, x0, policycost$D, z_g, z_m, z_i, opt$D, profit = profit, penalty = penalty))

Make data tidy (melt), fast (data.tables), and nicely labeled.

dat <- melt(sims, id = names(sims[[1]]))
dt <- data.table(dat)
setnames(dt, "L1", "reps")  # names are nice

Plots

A single replicate, alternate dynamics show the Reed optimum, while harvest/fishstock show the impact of having policy costs.

ggplot(subset(dt, reps == 1)) + geom_line(aes(time, alternate)) + 
    geom_line(aes(time, fishstock), col = "darkblue") + geom_line(aes(time, 
    harvest), col = "purple") + geom_line(aes(time, harvest_alt), col = "darkgreen")

A second replicate

ggplot(subset(dt, reps == 2)) + geom_line(aes(time, alternate)) + 
    geom_line(aes(time, fishstock), col = "darkblue") + geom_line(aes(time, 
    harvest), col = "purple") + geom_line(aes(time, harvest_alt), col = "darkgreen")

Profits

ggplot(subset(dt, reps == 1)) + geom_line(aes(time, profit_fishing)) + 
    geom_line(aes(time, policy_cost), col = "darkblue")

costs <- dt[, sum(policy_cost), by = reps]
profits <- dt[, sum(profit_fishing), by = reps]

qplot(costs$V1)

qplot(profits$V1)

qplot(profits - costs$V1)

Error: stat_bin requires the following missing aesthetics: x

We can visualize the equilibrium policy for each possible harvest:

policy <- sapply(1:length(h_grid), function(i) policycost$D[[i]][, 
    1])
ggplot(melt(policy)) + geom_point(aes(h_grid[Var2], (x_grid[Var1]), 
    col = h_grid[value] - h_grid[Var2])) + labs(x = "prev harvest", y = "fishstock") + 
    scale_colour_gradientn(colours = rainbow(4))

Here we plot previous harvest against the recommended harvest, coloring by stocksize. Note this swaps the y axis from above with the color density. Hence each x-axis value has all possible colors, but they map down onto a subset of optimal harvest values (depending on their stock).

policy <- sapply(1:length(h_grid), function(i) policycost$D[[i]][, 
    1])
ggplot(melt(policy)) + geom_point(aes(h_grid[Var2], (h_grid[value]), 
    col = x_grid[Var1]), position = position_jitter(w = 0.005, h = 0.005), alpha = 0.5) + 
    labs(x = "prev harvest", y = "harvest") + scale_colour_gradientn(colours = rainbow(4))