library("methods")
library("knitr")
basename <- "supplement"
opts_chunk$set(fig.path = paste("components/figure/", basename, "-", sep=""),
               cache.path = paste("components/cache/", basename, "/", sep=""))
opts_chunk$set(cache = 2)
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, 
               comment = NA, verbose = TRUE, echo=FALSE)
# PDF-based figures
opts_chunk$set(dev='pdf')
fig.cache <- TRUE

library("rmarkdown")
library("pdgControl")
library("reshape2")
library("plyr")
library("ggplot2")
library("data.table")
library("pander")
library("cboettigR")
library("ggthemes")
library("snowfall")
library("tidyr")
library("dplyr")
options(digits=2)

theme_set(theme_tufte(base_size=12) + theme(axis.line = element_line()))

\makeatletter \renewcommand\appendix{\par \setcounter{section}{0}% \setcounter{subsection}{0}% \setcounter{equation}{0}% \setcounter{table}{0}%------------ << add \setcounter{figure}{0}%----------- << add \gdef\theequation{A\arabic{equation}}% \gdef\thefigure{A\arabic{figure}}% \gdef\thesection{A@section}% \@addtoreset{equation}{section}% \@addtoreset{table}{section}%----- << add \@addtoreset{figure}{section}%---- << add } \makeatother

\appendix

Appendix A

Code and data for methods and analysis ##

All code and data used in this analysis can be found under the project's GitHub repository: https://github.com/cboettig/pdg_control. A snapshot of this repository representing the precise versions of the code used in this analysis has been permanently archived in the scientific software and data repository Zenodo, managed by CERN; and should be referenced or cited using its assigned digital object identifier, doi:10.5281/zenodo.21416 (This DOI is also clearly linked from the GitHub repository).

Organization of code

The project repository is organized as an R package. The intent is not to provide a general-purpose software tool for this kind of analysis, but rather to rely on a standardized directory structure for representing complex projects. It is our hope that this facilities understanding, reproducibility, and further exploration of these results sensu Gentleman & Temple Lang 2004, but caution that additional software engineering would be required before any of this could be applied in a different context. The functions provided are largely specific to the analyses presented here and may undocumented or incompletely documented.

The source documents for both the manuscript and supplement are also found in the code repository using the dynamic document .Rmd format which embeds all of the code required for the analysis; see http://rmarkdown.rstudio.com/.

Within the .Rmd files themselves are found only the blocks of code required to generate each of the figures shown in the manuscript. The analysis is sourced in from an external R script, analysis.R by a code block at the beginning of the .Rmd file. This script in turn relies on more modular functions defined in the R/ directory that define the stochastic dynamic programming optimization algorithms, simulation routines, and so forth.

Data

All empirical or simulated data shown in the figures can also be found as tabular .csv files in this repository.

Dependencies and machine images

The package DESCRIPTION file lists the software dependencies required to run the included scripts. To facilitate portability and preservation of these dependencies, we also provide a machine image with all the necessary dependencies already installed which can be run on most common computer platforms using Docker. Our Docker image, cboettig/pdg-control can be freely downloaded from the Docker Hub, or built using the Dockerfile included in the package repository.

Impact on the stock dynamics

Figure A1 shows the corresponding stock dynamics for the replicates illustrated in Figure 3 of the main text. The impact of the adjustment cost on the stock dynamics is less dramatic than it is on the harvest dynamics, due to the stochastic growth in stock sizes and the sequence in which the steps are applied: stock values are assessed without error, then an optimal harvest is determined (using the pre-computed optimal policy based on both the stock assessment and the previous harvest, to account for possible adjustment costs), and then the stock that escapes harvesting recruits stochastically to set the population size for the next year.

Though this results in an effect that is less pronounced, the stock dynamics of Figure A1 reflect the same general trends observed in Main text Figure 3; e.g. the fixed fee $\Pi_3$ results in the most exaggerated fluctuations in the stock sizes. It is worth noting how changes in policy costs alone, such as introducing or rising fixed costs, can have important knock-on effects to population dynamics such as increased variance, a recent phenomenon observed in important fisheries that has received much attention elsewhere [@Anderson2008].

labeller <- function(variable,value){
    return(relabel[paste(value)])
}

dt_s <- dt %>% 
  filter(replicate=="rep_17" & time < OptTime) %>%
  select(time, fishstock, alternate, penalty_fn) %>%
  gather(variable, value, -time, -penalty_fn) 
dt_s$penalty_fn <- factor(dt3$penalty_fn, levels=c("L1", "L2", "fixed"))


ggplot(dt_s, aes(time, value, col=variable, lty=variable)) +
  geom_line(lwd=1) +
  facet_grid(penalty_fn~., labeller = labeller) + 
  labs(x="time", y="stock size", title = "Example Stock Dynamics")  +
  scale_color_discrete(labels=c("With Penalty", "Reed Solution"), name="") +
  scale_linetype_discrete(labels=c("With Penalty", "Reed Solution"), name="")  

write.csv(dt_s, "components/data/figureS1.csv")
stats_df %>% 
  filter(variable != "cross.correlation") %>%
  separate(variable, c("measurement", "statistic"), sep="\\.") %>%
  filter(measurement == 'stock') %>% 
  spread(statistic, value) ->
  s2df

s2df$penalty_fn <- factor(f4df$penalty_fn, levels = c("L1", "L2", "fixed"))

p4 <- ggplot(s2df, aes(x = penalty_fraction, fill=penalty_fn, col=penalty_fn)) + 
  coord_cartesian(xlim = c(0, .3)) + 
  scale_color_discrete(labels=relabel, name="cost model") + 
  scale_fill_discrete(labels=relabel, name="cost model") 

p4a <- p4 + 
  stat_summary(aes(y=autocorrelation),
               fun.y = mean, 
               fun.ymin = function(x) mean(x) - sd(x), 
               fun.ymax = function(x) mean(x) + sd(x), 
               geom = "ribbon", alpha=0.3, colour=NA) +
  stat_summary(aes(y=autocorrelation), fun.y = mean, geom = "line") + 
  xlab("")

p4b <- p4 + 
  stat_summary(aes(y=variance),
               fun.y = mean, 
               fun.ymin = function(x) mean(x) - sd(x), 
               fun.ymax = function(x) mean(x) + sd(x), 
               geom = "ribbon", alpha=0.3, colour=NA) +
  stat_summary(aes(y=variance), fun.y = mean, geom = "line") +
  xlab(expression(paste("Penalty size (as fraction of ", NPV[0], ")"))) 


legend <- gtable_filter(ggplot_gtable(ggplot_build(p4a)), "guide-box") 
grid.arrange(do.call(arrangeGrob, 
            lapply(list(p4a, p4b),
                   function(x) x + theme_bw() + theme(legend.position="none"))),
    legend,
    ncol = 2,
    widths = c(5,1))

write.csv(s2df, "components/data/figureS2.csv")

Recall that in our model definition, stock sizes are measured before harvest ($h_t$ being determined after observing stock size $N_t$) and thus show less influence of the policy cost than the harvest itself.

Figure A2 shows the same summary statistics illustrating the trends of increasing adjustment costs on the stock dynamics, instead of the harvest dynamics shown in Figure 4 of the main text. Note that the variance in stock sizes under the quadratic penalty ($\Pi_2$) increases, while in Figure 4 of the main text the corresponding variance in the harvest decreases. Under costly adjustment, the undershooting behavior of the $\Pi_2$ harvest policy cannot absorb as much of the natural fluctuation, resulting in increased variation in the stock sizes.

Additional stochastic realizations ##

Figure 3 of the main text shows only a single stochastic realization. Here in Figure A3 we provide a sample of further realizations (including both stock and harvest dynamics).

labels <- c(harvest = substitute(paste("Harvest Policy under ", paste(Pi[1-3]))), 
            harvest_alt = substitute(paste("Harvest under ", paste(Pi[0]), " (Reed)")),
            fishstock = substitute(paste("Stock dynamics under under ", paste(Pi[1-3]))), 
            alternate = substitute(paste("Stock under ", paste(Pi[0]), " (Reed)")))
labeller <- function(variable,value){
if (variable=='penalty_fn') {
    return(relabel[paste(value)])
  } else {
    return(value)
  }
}

dt$penalty_fn <- factor(dt$penalty_fn, levels = c("L1", "L2", "fixed"))


ex <- dt %>% 
  filter(replicate %in% paste0("rep_", 16:20) & time < OptTime) %>%
  select(time, fishstock, harvest, alternate, harvest_alt, replicate, penalty_fn) %>%
  gather(variable, stock, -time, -replicate, -penalty_fn)


ggplot(ex) + 
  geom_line(aes(time, stock, col=variable), lwd=1) + 
  scale_color_discrete(labels=labels, name="") +  
  facet_grid(replicate ~ penalty_fn, labeller=labeller) 



write.csv(ex, "components/data/figureS3.csv")

Economic penalties using quadratic costs

Figure A4 shows that adding quadratic terms to the cost of fishing effort directly, as in @Lewis1981 or @McGough2009, results in a similar smoothing effect as we observe in quadratic adjustment costs ($\Pi_2$), which contrasts sharply with the patterns we observe in other functional forms. This case corresponds to a profit function defined as:

$$\Pi_4(N_t,h_t) = p h_t - c_0 E_t - c_4 E_t^2$$

Where $c_4$ is calibrated as in the case of the other coefficients to induce the fixed reduction of net present value relative to $NPV_0$.

source("components/quad_costs.R")
labels <- c(harvest = substitute(paste("Harvest Policy under ", paste(Pi[4]))), 
                                harvest_alt = substitute(paste("Harvest under ", paste(Pi[0]), " (Reed)")))
ex <- dt_quad %>% 
  filter(replicate=="rep_17" & time < OptTime) %>%
  select(time, harvest, harvest_alt) %>%
  gather(variable, stock, -time) 
ggplot(ex) + 
  geom_line(aes(time, stock, col=variable), lwd=1) + 
  scale_color_discrete(labels=labels, 
                       name="") 
shifts <- dt[, fraction_no_shift(harvest), by=c("replicate")]
ave <- shifts[, mean(V1)]

Like the $\Pi_2$ adjustment cost based penalty, this quadratic cost based penalty shows periods of constant, strictly positive policy only r ave * 100% of the time, and exhibits the same kind of smoothing of more small adjustments that is seen in $\Pi_2$ scenario.

Ignoring adjustment costs is worst than wrongly assuming adjustment costs

Figure A5 extends Figure 5 of the main text across all three functional forms of adjustment costs, $\Pi_{1-3}$. To reduce over-plotting of Figure 5 we have also separated the distributions of dockside revenue and adjustment fees by scenario, such that the top row shows the distributions when the policy accounts for adjustment costs, and the bottom shows the policy when it does not (that is, the original @Reed1979 model applied to the scenario of costly adjustment).

assume <- levels(as.factor(hist_dat$Assumption))
labeller <- function(variable,value){
if (variable=='penalty_fn') {
    return(relabel[paste(value)])
  } else {
    return(paste(value))
    #return(assume[paste(value)])
  }
}

ggplot(hist_dat) + 
  geom_density(aes(value, fill=variable, color=variable), alpha=0.8)+
  facet_grid(Assumption~penalty_fn, labeller = labeller)

write.csv(hist_dat, "components/data/figureS5.csv")

Figure A6 extends the results shown in Figure 6 across different severities of adjustment cost, totaling from 10% up to 25% of the $NPV$ of the stock. While the same qualitative pattern is observed in all cases, it is strongest in the higher penalties.

table1_long$penalty.fn <- factor(table1_long$penalty.fn, levels = c("L1", "L2", "fixed"))

ggplot(table1_long, aes(penalty.fn, value, fill = variable)) + 
  geom_bar(stat="identity", position="dodge") + 
  facet_wrap(~reduction, ncol=2) + 
  scale_x_discrete(labels = c('L1' = substitute(paste(Pi[1])),
                              'L2' = substitute(paste(Pi[2])),
                              'fixed' = substitute(paste(Pi[3])))) + 
  xlab("Penalty function")
write.csv(table1_long, "components/data/figureS6.csv")

Mismatched penalty functions

When policy costs are present, assuming the wrong functional form may or may not be worse than ignoring the presence of policy costs all together. Figure A7 considers several such cases. The X axis labels indicate first which policy is being used to drive the harvest decisions, and then which penalty is actually being applied. Thus $\Pi_1_\Pi_2$ indicates a policy calculated under the linear ($\Pi_1$) costs, when the reality involves quadratic ($\Pi_2$) costs. Though the coefficients of the ($\Pi_2$) costs have been scaled such that the net effect is the same, the inferred policies differ. As a result, the value is slightly more than half of the optimal free solution when the expected reduction should only be 10% (panel A) less than optimal free based on the $c_i$ calibration. In this case, ignoring the costs all together is still worse than using $\Pi_1$ costs -- the red bar is lower than the blue.

For fixed costs ($\Pi_3$), the reverse is true. Because fixed costs increase the volatility of the optimal solution while $\Pi_2$ costs decrease it, simply ignoring the cost of adjustment all together results in a higher value than either assuming costs are quadratic ($\Pi_2$) when they are fixed ($\Pi_3$), or assuming they are fixed when they are quadratic. Scenarios in which true costs are quadratic (regardless of what cost structure is assumed by the policy) do significantly worse than those that are not (i.e. all cases that assume quadratic but actually have linear or fixed costs).

who <- c("penalty_fn", "ignore_fraction", "mismatched_fraction", "reduction")
table2 <- arrange(mismatches_df[who], reduction) 
names(table2) = c("penalty.fn", "ignoring", "mismatched", "reduction")
table2_long <- melt(table2, id = c('penalty.fn', 'reduction'))
table2_long$reduction <- as.factor(table2_long$reduction)

ggplot(filter(table2_long, reduction==reduction), 
  aes(penalty.fn, value, fill = variable)) + 
  geom_bar(stat="identity", position="dodge") + 
  facet_wrap(~reduction, scales="free_y") +
  scale_x_discrete(labels = c('L1_L2' = substitute(paste(Pi[1],"_", Pi[2])), 
                              'L2_L1' =  substitute(paste(Pi[2],"_", Pi[1])),
                              'L1_fixed' = substitute(paste(Pi[1],"_", Pi[3])), 
                              'fixed_L1' =  substitute(paste(Pi[3],"_", Pi[1])),
                              'fixed_L2' = substitute(paste(Pi[3],"_", Pi[2])), 
                              'L2_fixed' =  substitute(paste(Pi[2],"_", Pi[3])))) + 
  xlab("Penalty function")
write.csv(table2_long, "components/data/figureS7.csv")

These patterns hold across the different size reductions shown (from 10% to 25%, Figure A8 panels a-d respectively), though net present value becomes negative under the quadratic form for the larger penalties as costs paid for adjustments exceed profits derived from harvest. (Though optimal solutions will always avoid negative net present values, these scenarios in which policies are applied under conditions other than assumed in the optimality calculation have no such guarantee.) For each scenario, results are averaged over 100 replicates.

Robustness to reduction level

## Reset parameter values
parallel = TRUE
ncpu = parallel::detectCores()
price = 10
c0 = 30
profit <- profit_harvest(price = price, c0 = c0, c1 = 0)
c2 <- seq(0, 40, length.out=100)
reduction <- 0.25 
seed <- 123                 # Random seed (replicable results)
delta <- 0.05               # economic discounting rate
OptTime <- 50               # stopping time
gridsize <- 50              # grid size for fish stock and harvest rate (discretized population)
sigma_g <- 0.2              # Noise in population growth
reward <- 0                 # bonus for satisfying the boundary condition
z_g <- function() rlnorm(1,  0, sigma_g) # mean 1
z_m <- function() 1         # No measurement noise, 
z_i <- function() 1         # No implemenation noise
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
x_grid <- seq(0.01, 1.2 * K, length = gridsize)  
h_grid <- seq(0.01, 0.8 * K, length = gridsize)

## Recompute
source("components/vary_reduction.R")

Figures A8-A11 illustrate that the these qualitative patterns still hold at when the models are calibrated to a reduction in NPV of 10% or 30%, though naturally the magnitude of the pattern is greater when the penalty is larger (e.g. 30% vs 10%).

p3 <- ggplot(low$dt3, aes(time, value, col=variable, lty=variable)) +
  geom_line(lwd=1) +
  facet_grid(penalty_fn~., labeller = labeller) + 
  labs(x="time", y="stock size", title = "Example Harvest Dynamics")  +
  scale_color_discrete(labels=c("Adjustment penalty", "No penalty"), name="") +
  scale_linetype_discrete(labels=c("Adjustment penalty", "No penalty"), name="") + 
  theme_bw() + labs(title="example trajectories at 10% reduction in NPV")

p3
p3 <- ggplot(high$dt3, aes(time, value, col=variable, lty=variable)) +
  geom_line(lwd=1) +
  facet_grid(penalty_fn~., labeller = labeller) + 
  labs(x="time", y="stock size", title = "Example Harvest Dynamics")  +
  scale_color_discrete(labels=c("Adjustment penalty", "No penalty"), name="") +
  scale_linetype_discrete(labels=c("Adjustment penalty", "No penalty"), name="") + 
  theme_bw() + labs(title="example trajectories at 30% reduction in NPV")
p3
fig5_dat <- low$dt5 %>% 
  filter(penalty_fn=="L1") %>%
  select(-penalty_fn) %>%
  unite(dist, Assumption, variable, sep = " : ")
labels <- c(expression(paste("Costs under ",  Pi[1])), 
            expression(paste("Revenue under ", Pi[1])),
            expression(paste("Cost under ", Pi[0])), 
            expression(paste("Revenue under ", Pi[0])))
ggplot(fig5_dat) + 
  geom_density(aes(value, fill=dist, color=dist, lty=dist), alpha=0.5, size=1) +
  scale_fill_discrete(name="", labels=labels) + 
  scale_color_discrete(labels=labels, name="") + 
   scale_linetype_discrete(labels=labels, name="") +
  xlab("Realized NPV") + 
  ylab("Probability density of revenue/cost")
fig5_dat <- high$dt5 %>% 
  filter(penalty_fn=="L1") %>%
  select(-penalty_fn) %>%
  unite(dist, Assumption, variable, sep = " : ")
labels <- c(expression(paste("Costs under ",  Pi[1])), 
            expression(paste("Revenue under ", Pi[1])),
            expression(paste("Cost under ", Pi[0])), 
            expression(paste("Revenue under ", Pi[0])))
ggplot(fig5_dat) + 
  geom_density(aes(value, fill=dist, color=dist, lty=dist), alpha=0.5, size=1) +
  scale_fill_discrete(name="", labels=labels) + 
  scale_color_discrete(labels=labels, name="") + 
   scale_linetype_discrete(labels=labels, name="") +
  xlab("Realized NPV") + 
  ylab("Probability density of revenue/cost")


write.csv(hist_dat, "components/data/figure5.csv")

\newpage

References



cboettig/pdg_control documentation built on May 13, 2019, 2:10 p.m.