Emgo decision rules.Rmd

---
title: "DRAFT:  Emperor Goose Decision Rules"
author: "Erik Osnas and Charles Frost"
date: "June 27, 2016"
output: html_document
---

This is an Rmarkdown document.  Code and associated files that produced this document can be found at https://github.com/eosnas/Emperor-Goose-Harvest-Strategy.git

```{r glob_opts, include=FALSE}
knitr::opts_chunk$set(cache=TRUE, echo=FALSE, eval=TRUE, warning=FALSE, message=FALSE)
```

# Introduction
We report on a population model and analysis to derive optimal harvest policy thresholds for emperor geese, given a statement of stakeholder values and uncertainty in goose population dynamics and future harvest.  The methods used are similar to those of the adaptive harvest management (AHM) program of the United States Fish and Wildlife Service (USFWS; USFWS 2015).  Uncertainty was characterized by sampling from probability distributions; using those from a complimentary analysis of emperor goose harvest potential (Dooley et al. 2016) or estimating distributions from existing data when possible. The analysis has four important components:  (1) a population dynamics model used to estimate population dynamic parameters from survey and harvest data and to *predict* future emperor goose population size, given a harvest policy and the estimated population demographic parameters, (2) a model of the realized harvest under the harvest policy, (3) a statement of stakeholder values expressed as a function of population size, harvest, and harvest season regulatory burden, and (4) an optimization method to find the "best" harvest policy decision rule, given stakeholder values and the population model. We also conducted a sensitivity analysis where we evaluated the effect of different models of realized harvest (2, above) on the outcome of the optimal harvest policy.  Finally, we derived the optimal harvest policy under a desire to manage for maximum sustained yield.   

# Population Model and Parameter Estimation
We used the theta-logistic model to predict emperor goose populations under different harvest rules.  The theta-logistic model is  
$$E[N_{t+1}] = N_t + N_t r \left[1 - \left(\frac{N_t}{K}\right)^\theta \right] - \frac{H_t}{(1 - c)}$$
$$N_{t+1} \sim LogNormal(log(E[N_{t+1}]), \sigma)$$
$$Y_t \sim Normal(q N_t, \sigma_{obs,t})$$
Where $N_t$ is the total population size in year $t$, $E[\cdot{}]$ is the expectation of the quantity inside the brackets, $r$ is the maximum growth rate, $K$ is the carrying capacity determined by the strength of density-dependence, $\theta$ is a parameter that controls the non-linearity in density-dependence, $H_t$ is the observed harvest in year $t$, and $c$ is the proportion of the kill that is unharvested (due to crippling) or unreported, $\sigma$ is the standard deviation of random year effects that describe the departure from the deterministic expectation, $E[N_t]$. Observations of a population index $Y_t$ are some constant proportion $q$ of the total population and observed with error that is normally distributed with know standard deviation $\sigma_{obs,t}$. Harvest in year $t$ was distributed as $H_t \sim Normal(\bar{H}, \sigma_H)$.    

To obtain parameter estimates, we used the observed YKD Coastal Zone Goose Survey and harvest data as reported in Dooley et al. (2016) and fit the model using Bayesian Markov chain Monte Carlo methods implemented in program JAGS (Plummer 2013). Bayesian analysis requires specifying prior distributions on parameters; therefore, we used priors that were either identical to or consistent with the distributions used in Dooley et al. (2016).  These priors can be found in the Supplemental code.  It is important to note that much of the information for the parameter estimates of the model come from the priors and not the data; therefore, these prior distribution are very important to the analysis.  

#Derivation of Harvest Decision Thresholds
```{r functions}
## Define important functions
HarvestPolicy <- function(n, plist = list()){
  hpol <- ifelse(n < plist[[1]], "Red", ifelse(n < plist[[2]],"Yellow","Green"))
  return(hpol)
}

utility <- function(n=NA, h=NA, policy=NULL, par=list(npar=c(1,1), hpar=c(1,1)), weight=c(0.5, 0.5), fn=c("logistic", "linear"), policy.parameters=c("Red"=1, "Yellow"=1, "Green"=1)){
  if(fn[1]=="exp"){util.n <- 1 - exp(-par$npar[1]*n)}
  if(fn[2]=="exp"){util.h <- 1 - exp(-par$hpar[1]*h)}
  if(fn[1]=="logistic"){util.n <- plogis(q=n, par$npar[1], par$npar[2])}
  if(fn[2]=="logistic"){util.h <- plogis(q=h, par$hpar[1], par$hpar[2])}
  if(fn[1]=="linear"){util.n <- ifelse(n <= par$npar[1], n/par$npar[1], 1)}
  if(fn[2]=="linear"){util.h <- ifelse(h <= par$hpar[1], h/par$hpar[1], 1)}
  if(fn[1]=="msy"){util.n <- rep(0,length(h))}
  if(fn[2]=="msy"){util.h <- h}
  if(length(policy)>0){
    if(sum(!policy%in%names(policy.parameters))==0){ util.h <- util.h*policy.parameters[policy]}
    else(stop("Names of policy parameters do not match policy"))
  }
  
  if(length(n)==1){
    util <- matrix(c(rep(util.n, length(h)), util.h),length(h),2)%*%matrix(weight,2,1)
  }
  if(length(h)==1){
    util <- matrix(c(util.n, rep(util.h, length(n))),length(n),2)%*%matrix(weight,2,1)
  }
  if(length(n)==length(h)){
    util <- matrix(c(util.n, util.h),length(h),2)%*%matrix(weight,2,1)
  }
  return(util)
}
```
### Elicitation of Population Utility
We elicited a utility function for emperor goose population size from Eric Taylor on 27 April 2016 and ADFG representatives on 20 May 2016 using two standard methods (probability and certainty equivalence). A utility function measures how value of a good or service (here population size) changes over a range of outcomes when it is uncertain which outcome will be realized.  Here the key uncertainty is the future population size of emperor geese after a subsistence harvest season is opened. In the first method (probability equivalence), we gave Eric Taylor a choice between a certain outcome of goose population size and entering a "bet" where the could be better or worse than the certain outcome.  He was then asked to state the minimum probability for the better outcome where he would switch to favoring the bet over the certain outcome.  The second method ("certainty equivalence") consisted of the same decision between a "bet" and a certain outcome but this time the probabilities for the better or worse outcome where equal (0.5), and the value of the certain outcome that gave indifference to the bet was elicited.  The range of goose population size ranged from 50,000 to 200,000. These values represented an approximate range of the emperor goose population sizes that were found under simulation from the theta-logistic model using estimated parameters and historical harvest.  It is assumed that Eric Taylor's utility function for population size represents that of the USFWS. We then pooled the data from the two replicate elicitations and fit a logistic function with an intercept at zero using maximum likelihood.  ADFG's utility curve was very similar to that of USFWS (not shown).  The resulting estimates were used as the population size utility function, shown here:


```{r P_utility}
# results from Elicitation of utiity function from Eric Tayor
# on 4-27-2016
u1 <- list(N = c(0,50,60,80,100, 150, 175, 200), 
           U = c(0,0, 0, 0.6, 0.75, 0.9, 1, 1))
u2 <- list(N = c(0,60,80,100,125,155, 200), 
           U = c(0,0.25,0.5,0.625, 0.75, 0.875, 1))
#pdf("plot1.pdf")
plot(u1$N, u1$U, pch=16, col=1, xlab="Population Size (x1000)", ylab="Utility", main = "Utility function from Eric Taylor")
points(u2$N, u2$U, pch=16, col=2)
legend("topleft", legend=c("First replicate", "Second replicate", "Best fit logistic function"), pch=c(16,16,NA), col=c(1,2,1), lty=c(NA, NA, 1))

u <- rbind(as.data.frame(u1),as.data.frame(u2))
u$N2 <- u$N*u$N
u$N3 <- u$N*u$N*u$N 

library(bbmle)
nloglike <- function(m = 80, s=20, sigma=0.1){
  -sum(dnorm(U, plogis(N, location=m, scale=s), sigma, log=TRUE))
}

fit2 <- mle2(nloglike, start=list(m=80, s=20, sigma=0.1), data=u)
#summary(fit2)
lines(1:200, plogis(1:200, location=coef(fit2)[1], scale=coef(fit2)[2]))
#dev.off()
```

### Harvest Utility
We did not use the same elicitation process for harvest utility.  Instead, we specified a function based on listening to stakeholders at various AMBCC meetings. From this, we determined that two important characteristics of harvest utility existed:  (1) utility increases positively with harvest until all demand is met and (2) harvest under a regulation-free ("traditional") subsistence season is more valuable than when regulations are imposed.  The decision rule will consider thresholds where additional regulations are imposed to restrict harvest from no additional regulations beyond the standard subsistence regulations, ("traditional"" or "Green") to mild restriction that reduce harvest but are yet to be determined ("Restricted" or "Yellow") to  full season closure for emperor geese (the current status quo with a "closed" but with an illegal harvest of approximately 4000).  Therefore, we used a utility function that increased linearly up to a maximum that depended on the harvest restrictions regime.  Under the "traditional" harvest season, the utility function reached a maximum of 1.0 at 15000 harvested geese, after which it was constant.  Under a "Redistricted (Yellow)" season and "Red (Closed)" season full utility is 1/2 and 1/4 that of the "traditional" season, respectively, as seen in the figure. Until we can get direct feedback from stakeholders, we feel this utility function approximates that of the subsistence community.  

```{r H_utility}
har=seq(1,2e4,length=100)
u <- utility(n=2e4, h=har, par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                      policy="Green",
                      policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
plot(har, u, type="l", col="green", lwd=2, xlab="Harvest", ylab="Utility", main="Utilty as a function of harvest and season type")
u <- utility(n=2e4, h=har, par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                      policy="Yellow",
                      policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
lines(har, u, col="orange", lwd=2)
u <- utility(n=2e4, h=har, par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                      policy="Red",
                      policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
lines(har, u, col="red", lwd=2)
legend("topleft", legend=c("Traditional Season", "Restricted/Managed Season", "Very Limited Season"), lty=1, lwd=2, col=c("Green", "orange", "Red"))
```
 
### Total Utility
We assumed that total utility for a year was determined by weighting population and harvest utility equally and adding
$$U_t(H,N) = 0.5 U_t (H) + 0.5 U_t(N)$$
Under a "traditional" subsistence season and the population utility described above, the utility function looks like this:

```{r tot_utility}
num=seq(1,2e5, length=100)
har=seq(1,2e4,length=100)
reward=matrix(NA, length(num), length(har))
for(i in 1:length(har)){
  reward[,i] = utility(n=num, h=har[i], par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                      policy="Green",
                      policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}

image(reward, x=num, y=har, col=gray(1:10/10), xlab="Population Size", ylab="Harvest") #terrain.colors(10)))
text(x=175000, y=14000, labels="+", font = 2, cex = 2)
text(x=25000, y=1000, labels="-", font = 2, cex = 2)
```

In the above figure , white is high utility and dark is low utility.  

### Harvest Under Each Policy
The consequence of different harvest policies depend on the actual number harvested under these policies. We only have harvest data from the past under a "closed" season (1987-present) and only two years of data before 1987 when the season was open in a "traditional" manner. During the period when the was no legal take of emperor geese, the *reported* harvest ranged between about 2600 and 5700 (Dooley et al. 2016). Therefore, we used log-normal distributions to represent our uncertainty in harvest under each harvest policy.  For the "Closed" or "Very Limited" or "Red" policy we assumed harvest was distributed as $H_{Red} \sim LogNormal(log(4000), 0.1)$, for the policy under "Restricted" or "Yellow" we assumed harvest was distributed as $H_{Yellow} \sim LogNormal(log(7000), 0.1)$, and under the "Traditional" or "Green" policy we assumed harvest was distributed as $H_{Green} \sim LogNormal(log(15000), 0.2)$. This gives probability densities:
```{r harvest_plot}
#plot harvest distribution under each policy
plot(density(rlnorm(10000, log(4000), 0.1)), xlim=c(0,25000), col="red", xlab="Harvest", main="Harvest Distribution Under Each Policy", lwd=2)
lines(density(rlnorm(10000, log(7000), 0.1)), col="orange", lwd=2)
lines(density(rlnorm(10000, log(15000), 0.2)), col="green", lwd=2)
legend("topright", legend=c("Traditional Season", "Restricted/Managed Season", "Very Limited Season"), lty=1, lwd=2, col=c("Green", "orange", "Red"))
```

### Optimization of Policy Thresholds
The task is to find the decision threshold that maximizes expected utility over a long period of time.  There are two decision thresholds.  $T_1$ is the threshold between a "Very Restricted" or "Red" season and a "Limited/Managed" or "Yellow" season. $T_2$ is the threshold between the "Yellow" and "Traditional" or "Green" season; thus, $T_1 \leq T_2$.  For each combination of threshold values, we calculated expected utility over (1) a sample of the posterior distribution of model parameters and (2) a sample of the harvest consequence for each policy ("Red", "Yellow" or "Green") and then projected the stochastic population dynamics forward for 200 years for each sample. At each point in time, the harvest policy was determined from the model-predicted observed value of the summer survey and then the utility contribution for that year was calculated. We repeated this for a large sample of the posterior ($S$ = 500) and then calculated cumulative utility and averaged this over the sample of parameters

$$E[U(T_1,T_2)] = \frac{1}{S}\sum_{i=1}^{S} \sum_{t=1}^{200} U_t(T_1, T_2, \phi_i) \delta^{t-1}$$

where $\phi_i$ is a parameter sample from the posterior and harvest distribution and $\delta$ is a decay parameter that discounts future utility.  We used a small discounting rate, $\delta = 0.99$. We then found the set $\{T_1,T_2\}$ that gave the maximum expected utility.  

We chose to use the observed summer survey value, $Y_t$, as the decision metric; thus, when $Y_t < T_1$, the hares season policy would be "Red" and "Very Restrictive" regulations would be in place.  This rule can be translated to the total population by $T'_1 = \hat{q} T_1$, where $\hat{q}$ is the mean or median (or other quantile) of the posterior distribution of $q$.  The rule then becomes "Red" when $\hat{N}_t < T'_1$, where $\hat{}$ is taken as the same quantile used for $q$ to project $T$ to $T'$.  

# Results:  Optimal Thresholds
```{r main_results}
posterior <- TRUE
filename <- "summer_theta_logistic_output.2.txt"
Reps <- 2000   #1500 works well
T <- 200
a <- 0.0001   #functional response parameter for harvest, ensures that har --> 0 as N --> 0
d <- 0.99     #utility discounting, 1 - discount rate
N.last <- 26235
if(posterior==FALSE){  #need to modify to match below at posterior == TRUE
NULL
}
if(posterior==TRUE){
  post <- dget(file=filename)
  dims <- dim(post$N.tot)
  spots <- sample(1:dims[1], Reps)
  N0 <- post$N.tot[spots,dims[2]]
  Rmax <- post$r.max[spots,1]
  CC <- post$CC[spots,1]
  Theta <- post$theta[spots,1]
  q <- post$q[spots,1]
  sigma.proc <- post$sigma.proc[spots,1]
  cv.obs <- rnorm(Reps, 0.07, 0.003)  #from lm(dat$SUM_TOT_SE~dat$SUM_TOT-1)
}
t.min <- 10000
t.max <- 50000
t.delta <- 1000

Har = list(Red = rlnorm(Reps, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(5000), 0.1), 
              Green = rlnorm(Reps, log(15000), 0.2))

thresholds = seq(t.min, t.max, by=t.delta)
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))

#Iterate through threshold combinations
for(ii in 1:length(thresholds)){for(jj in ii:length(thresholds)){
  
#define harvest thresholds
H.pol = list(Red = thresholds[ii], Yellow = thresholds[jj])

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )

for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <- utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}
discount <- matrix(d^c(0:(T-1)), nrow=Reps, ncol=T, byrow=TRUE)
e.util[ii,jj] <- mean(apply(Results$Utility*discount,1,sum))
}}

smart.spot = which(e.util==max(e.util, na.rm=TRUE), arr.ind = TRUE)
harvest.rule=list(Red=thresholds[smart.spot[1]], Yellow=thresholds[smart.spot[2]])
```
Given the utility functions and assumptions about harvest under each policy, described above, the management thresholds that maximize expected utility were found to be:  

Harvest Policy   | Summer Index Value  | Total Population (Median)
---------------------  | ------------ | -------------------
Very Restricted (Red) | $N_t <$ `r as.integer(thresholds[smart.spot[1]])` | $Y_t <$ `r as.integer(round(thresholds[smart.spot[1]]/median(post$q),0))`  
Restricted (Yellow) | `r as.integer(thresholds[smart.spot[1]])` $\leq N_t <$ `r as.integer(thresholds[smart.spot[2]])` | `r as.integer(round(thresholds[smart.spot[1]]/median(post$q),0))` $\leq Y_t <$ `r as.integer(round(thresholds[smart.spot[2]]/median(post$q),0))` 
Traditional (Green) | $N_t \geq$ `r as.integer(thresholds[smart.spot[2]])` | $Y_t \geq$ `r as.integer(round(thresholds[smart.spot[2]]/median(post$q),0))`  


```{r plot_eutil_main}
image(e.util, axes=FALSE, col=terrain.colors(10))
contour(e.util, add=TRUE)
axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```


Here is a plot of these thresholds applied to historical data:
```{r plot_on_data}
require(scales)
dat <- list(y = c(19805, 12430, 13035, 16392, 16855, 17347, 14888, 15416, 17147, 18733, 18764, 24413, 
                  23287, 21741, 21406, 18667, 27297, 19504, 21378, 21396, 19798, 26562, 24362, 22100, 
                  20684, 20167, 21223, 20388, 29840, 32550, 26235),
            YR = c(1985:2015),
            har.obs = c(7152, 5378, 3638, NA, 3900, 5726, 4679, 4955, 4888, 3780, 4487, 4661, 3796, 
                        4185, 3246, 3637, 3207, 3373, NA, 3414, 2712, 4147, 3574, 3132, 4274, 3810, 
                        2629, NA, 3725, NA, NA)
            )
plot(1,1, type="n", pch=16, xlim=c(1985,2015), ylim=c(0, 50000),  ylab="Summer Index", xlab="Year")
polygon(x=c(0,0,2100,2100),y=c(-1000,thresholds[smart.spot[1]],thresholds[smart.spot[1]],-1000), col=alpha("red", 0.5))
polygon(x=c(0,0,2100,2100),y=c(thresholds[smart.spot[1]],thresholds[smart.spot[2]], thresholds[smart.spot[2]],thresholds[smart.spot[1]]),col=alpha("yellow", 0.5))
polygon(x=c(0,0,2100,2100),y=c(thresholds[smart.spot[2]],1000000,1000000,thresholds[smart.spot[2]]), col=alpha("green", 0.5))
#abline(h=harvest.rule, col=c("yellow", "red"))
points(dat$YR, dat$y, pch=16)
ci = apply(post$N.est, 2, quantile, probs=c(0.025, 0.5, 0.975))
arrows(x0=dat$YR, x1=dat$YR, y0=ci[1,],y1=ci[3,], length=0)
points(dat$YR, ci[2,], pch=21, bg="gray50")
points(dat$YR, dat$har.obs, bg=1, pch=22)
ci = apply(post$har, 2, quantile, probs=c(0.025, 0.5, 0.975))
arrows(x0=dat$YR, x1=dat$YR, y0=ci[1,],y1=ci[3,], length=0)
points(dat$YR, ci[2,], pch=22, bg="gray50")
legend("topleft", legend=c("Observed Summer Index", "Estimated Summer Total", "Reported Harvest", "Estimated Harvest"), pch=c(16,16,22,22), col=c(1,"gray50", 1,1), pt.bg=c(NA,NA,1,"gray50"))
```

Second, we can translate the decision rule to the total population estimate:  

```{r tot_plot}
plot(1,1, type="n", pch=16, xlim=c(1985,2015), ylim=c(0, 200000),  ylab="Total Population Estimate",  xlab="Year")
polygon(x=c(0, 0, 210000, 210000),y=c(-1000, thresholds[smart.spot[1]]/median(post$q), thresholds[smart.spot[1]]/median(post$q), -1000), col=alpha("red", 0.5))

polygon(x=c(0, 0, 210000, 210000),y=c(thresholds[smart.spot[1]]/median(post$q), thresholds[smart.spot[2]]/median(post$q), thresholds[smart.spot[2]]/median(post$q), thresholds[smart.spot[1]]/median(post$q)), col=alpha("yellow", 0.5))

polygon(x=c(0, 0, 210000, 210000), y=c(thresholds[smart.spot[2]]/median(post$q), 1000000, 1000000, thresholds[smart.spot[2]]/median(post$q)), col=alpha("green", 0.5))

points(dat$YR, dat$y/median(post$q), pch=21, bg=1)
```

```{r sim_main}
#Simulate dynamics at smart threshold
H.pol = harvest.rule

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}

p.regs <- as.matrix(table(Results$HarPol)/(T*Reps))
regs <- matrix(0, 3, 1)
row.names(regs) <- c("Red","Yellow","Green")

#n.plot = 100
#plot(1,1, type="n",xlim=c(0, T), ylim=c(0, 1e5), xlab="Year", ylab="Obs. Population")
#for(i in 1:min(100,Reps)){
#  lines(1:T, Results$N.obs[i,], col=1)
#}
#abline(h=H.pol, col=c("Red", "Yellow"), lwd=3)
#text(x=T*(1:dim(p.regs)[1])/10, y=9000, labels=row.names(p.regs))
#text(x=T*(1:dim(p.regs)[1])/10, y=5000, labels=round(p.regs,2))

#image(e.util, axes=FALSE, col=terrain.colors(10))
#axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
#mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
#mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
#mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```

Mean total population size under this policy is `r as.integer(round(mean(Results$N),0))` and the frequency of "Red", "Yellow", and "Green" seasons is 
```{r freq_main} 
t(p.regs)
```

# Sensitivity of Thresholds to Harvest Assumptions 
A key uncertainty is the harvest realized under each policy, as shown in the figure above ("Harvest Distribution Under Each Policy").  The distribution used under a "Green" or "traditional" subsistence season is especially wide. In order to evaluate the effect of resolving this uncertainty, we found the optimum thresholds when harvest under a "Green" and "Yellow" policy was on the low and high end of their respective distributions (shown in the figure "Harvest Distribution Under Each Policy", above). Harvest under a "Yellow" policy was simulated as $H_{Yellow} \sim LogNormal(log(5000), 0.1)$ and harvest under a "Green" policy was $H_{Green} \sim LogNormal(log(7000), 0.1)$.  

Another key harvest issue is the magnitude of harvest under a "Red" or "very Restrictive" season. The purpose of this harvest policy is to allow the population to grow allowing larger future population sizes and harvest. Yet, historically, harvest has occurred even though no legal take was allowed. Because the harvest rate under each policy affects the derived thresholds, we wanted to show the effect of different harvest levels under a "Red" policy.  Therefore, we simulated harvest under "Red" as above and as $H_{Red} \sim LogNormal(log(300), 0.1)$, which is a very low harvest that might be realized in the most optimistic scenario. Thus, in total we simulated four combinations of actual harvest under difference harvest policies:  (1) high harvest under "Red" and low harvest under "Yellow and "Green", (2) high harvest under "Red" and high harvest under "Yellow" and "Green" (presented above), (3) low harvest under "Red" and low harvest under "Yellow" and "Green", and (4) low harvest under "Red" and high harvest under "Yellow" and "Green".  

```{r har_cases_plot}
par(mfrow=c(2,2))
#Case 1
Har = list(Red = rlnorm(10000, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(10000, log(5000), 0.1), 
              Green = rlnorm(10000, log(7000), 0.1))
plot(density(Har$Red), xlim=c(0,25000), col="red", xlab="Harvest", main="Case (1)", lwd=2)
lines(density(Har$Yellow), col="orange", lwd=2)
lines(density(Har$Green), col="green", lwd=2)
#Case 2
Har = list(Red = rlnorm(10000, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(10000, log(7000), 0.1), 
              Green = rlnorm(10000, log(15000), 0.1))
plot(density(Har$Red), xlim=c(0,25000), col="red", xlab="Harvest", main="Case (2)", lwd=2)
lines(density(Har$Yellow), col="orange", lwd=2)
lines(density(Har$Green), col="green", lwd=2)
#Case 3
Har = list(Red = rlnorm(10000, log(300), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(10000, log(5000), 0.1), 
              Green = rlnorm(10000, log(7000), 0.1))
plot(density(Har$Red), xlim=c(0,25000), ylim=c(0,0.001), col="red", xlab="Harvest", main="Case (3)", lwd=2)
lines(density(Har$Yellow), col="orange", lwd=2)
lines(density(Har$Green), col="green", lwd=2)
#Case 4
Har = list(Red = rlnorm(10000, log(300), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(10000, log(7000), 0.1), 
              Green = rlnorm(10000, log(15000), 0.1))
plot(density(Har$Red), xlim=c(0,25000), ylim=c(0,0.001), col="red", xlab="Harvest", main="Case (4)", lwd=2)
lines(density(Har$Yellow), col="orange", lwd=2)
lines(density(Har$Green), col="green", lwd=2)
par(mfrow=c(1,1))
```

### Case (1): high harvest under Red, low harvest under Green
```{r low_har}
Reps <- 1500
Har = list(Red = rlnorm(Reps, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(5000), 0.1), 
              Green = rlnorm(Reps, log(7000), 0.1))

thresholds = seq(t.min, t.max, by=t.delta)
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
```

```{r opt_low}
#Iterate through threshold combinations
for(ii in 1:length(thresholds)){for(jj in ii:length(thresholds)){
  
#define harvest thresholds
H.pol = list(Red = thresholds[ii], Yellow = thresholds[jj])

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )

for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <- utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}
discount <- matrix(d^c(0:(T-1)), nrow=Reps, ncol=T, byrow=TRUE)
e.util[ii,jj] <- mean(apply(Results$Utility*discount,1,sum))
}}

smart.spot = which(e.util==max(e.util, na.rm=TRUE), arr.ind = TRUE)
harvest.rule=list(Red=thresholds[smart.spot[1]], Yellow=thresholds[smart.spot[2]])
```

The optimal thresholds under the above harvest scenario 

```{r rule_low, echo=FALSE}
harvest.rule
```

```{r sim_low}
#Simulate dynamics at smart threshold
H.pol = harvest.rule

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}

p.regs <- as.matrix(table(Results$HarPol)/(T*Reps))
regs <- matrix(0, 3, 1)
row.names(regs) <- c("Red","Yellow","Green")

#n.plot = 100
#plot(1,1, type="n",xlim=c(0, T), ylim=c(0, 1e5), xlab="Year", ylab="Obs. Population")
#for(i in 1:min(100,Reps)){
#  lines(1:T, Results$N.obs[i,], col=1)
#}
#abline(h=H.pol, col=c("Red", "Yellow"), lwd=3)
#text(x=T*(1:dim(p.regs)[1])/10, y=9000, labels=row.names(p.regs))
#text(x=T*(1:dim(p.regs)[1])/10, y=5000, labels=round(p.regs,2))

#image(e.util, axes=FALSE, col=terrain.colors(10))
#axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
#mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
#mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
#mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```

Mean total population size under this policy is `r as.integer(round(mean(Results$N),0))` and the frequency of "Red", "Yellow", and "Green" seasons is 
```{r freq_low} 
t(p.regs)
```

### Case (2): high harvest under Red, high harvest under Green
```{r har_high}

Har = list(Red = rlnorm(Reps, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(7000), 0.1), 
              Green = rlnorm(Reps, log(15000), 0.1))

thresholds = seq(t.min, t.max, by=t.delta)
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
```
```{r opt_high}
#Iterate through threshold combinations
Reps <- 1000
for(ii in 1:length(thresholds)){for(jj in ii:length(thresholds)){
  
#define harvest thresholds
H.pol = list(Red = thresholds[ii], Yellow = thresholds[jj])

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )

for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <- utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}
discount <- matrix(d^c(0:(T-1)), nrow=Reps, ncol=T, byrow=TRUE)
e.util[ii,jj] <- mean(apply(Results$Utility*discount,1,sum))
}}

smart.spot = which(e.util==max(e.util, na.rm=TRUE), arr.ind = TRUE)
harvest.rule=list(Red=thresholds[smart.spot[1]], Yellow=thresholds[smart.spot[2]])
```

The optimal thresholds under the above harvest scenario 

```{r rule_high, echo=FALSE}
harvest.rule
```


```{r sim_high}
#Simulate dynamics at smart threshold
H.pol = harvest.rule

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}

p.regs <- as.matrix(table(Results$HarPol)/(T*Reps))
regs <- matrix(0, 3, 1)
row.names(regs) <- c("Red","Yellow","Green")

#n.plot = 100
#plot(1,1, type="n",xlim=c(0, T), ylim=c(0, 1e5), xlab="Year", ylab="Obs. Population")
#for(i in 1:min(100,Reps)){
#  lines(1:T, Results$N.obs[i,], col=1)
#}
#abline(h=H.pol, col=c("Red", "Yellow"), lwd=3)
#text(x=T*(1:dim(p.regs)[1])/10, y=9000, labels=row.names(p.regs))
#text(x=T*(1:dim(p.regs)[1])/10, y=5000, labels=round(p.regs,2))

#image(e.util, axes=FALSE, col=terrain.colors(10))
#axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
#mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
#mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
#mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```

Mean total population size under this policy is `r as.integer(round(mean(Results$N),0))` and the frequency of "Red", "Yellow", and "Green" seasons is 
```{r freq_high} 
t(p.regs)
```

### Case (3): low harvest under Red, low harvest under Green
```{r har_no_low}
Har = list(Red = rlnorm(Reps, log(300), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(5000), 0.1), 
              Green = rlnorm(Reps, log(7000), 0.1))

thresholds = seq(t.min, t.max, by=t.delta)
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
```

```{r opt_no_low}
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
Reps <- 1000
for(ii in 1:length(thresholds)){for(jj in ii:length(thresholds)){
  
#define harvest thresholds
H.pol = list(Red = thresholds[ii], Yellow = thresholds[jj])

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )

for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <- utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}
discount <- matrix(d^c(0:(T-1)), nrow=Reps, ncol=T, byrow=TRUE)
e.util[ii,jj] <- mean(apply(Results$Utility*discount,1,sum))
}}

smart.spot = which(e.util==max(e.util, na.rm=TRUE), arr.ind = TRUE)
harvest.rule=list(Red=thresholds[smart.spot[1]], Yellow=thresholds[smart.spot[2]])
```

The optimal thresholds are 

```{r rule_no_low}
harvest.rule
```

```{r sim_no_low}
#Simulate dynamics at smart threshold
H.pol = harvest.rule

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}

p.regs <- as.matrix(table(Results$HarPol)/(T*Reps))
regs <- matrix(0, 3, 1)
row.names(regs) <- c("Red","Yellow","Green")

#n.plot = 100
#plot(1,1, type="n",xlim=c(0, T), ylim=c(0, 1e5), xlab="Year", ylab="Obs. Population")
#for(i in 1:min(100,Reps)){
#  lines(1:T, Results$N.obs[i,], col=1)
#}
#abline(h=H.pol, col=c("Red", "Yellow"), lwd=3)
#text(x=T*(1:dim(p.regs)[1])/10, y=9000, labels=row.names(p.regs))
#text(x=T*(1:dim(p.regs)[1])/10, y=5000, labels=round(p.regs,2))

#image(e.util, axes=FALSE, col=terrain.colors(10))
#axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
#mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
#mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
#mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```

Mean total population size under this policy is `r as.integer(round(mean(Results$N),0))` and the frequency of "Red", "Yellow", and "Green" seasons is 
```{r freq_no_low} 
t(p.regs)
```

### Case (4): low harvest under Red, high harvest under Green
```{r har_no_high}
Har = list(Red = rlnorm(Reps, log(300), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(7000), 0.1), 
              Green = rlnorm(Reps, log(15000), 0.1))

thresholds = seq(t.min, t.max, by=t.delta)
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
```
```{r opt_no_high}
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
Reps <- 1000
for(ii in 1:length(thresholds)){for(jj in ii:length(thresholds)){
  
#define harvest thresholds
H.pol = list(Red = thresholds[ii], Yellow = thresholds[jj])

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )

for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <- utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                  policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}
discount <- matrix(d^c(0:(T-1)), nrow=Reps, ncol=T, byrow=TRUE)
e.util[ii,jj] <- mean(apply(Results$Utility*discount,1,sum))
}}

smart.spot = which(e.util==max(e.util, na.rm=TRUE), arr.ind = TRUE)
harvest.rule=list(Red=thresholds[smart.spot[1]], Yellow=thresholds[smart.spot[2]])
```

The optimal thresholds are 

```{r rule_no_high}
harvest.rule
```

```{r sim_no_high}
#Simulate dynamics at smart threshold
#Simulate dynamics at smart threshold
H.pol = harvest.rule

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], 
                                 par=list(npar=coef(fit2)[1:2]*1000, hpar=c(15000, 1)),
                                 policy=Results$HarPol[i,],
                                 policy.parameters=c("Red"=0.25, "Yellow"=0.5, "Green"=1))
}

p.regs <- as.matrix(table(Results$HarPol)/(T*Reps))
regs <- matrix(0, 3, 1)
row.names(regs) <- c("Red","Yellow","Green")

#n.plot = 100
#plot(1,1, type="n",xlim=c(0, T), ylim=c(0, 1e5), xlab="Year", ylab="Obs. Population")
#for(i in 1:min(100,Reps)){
#  lines(1:T, Results$N.obs[i,], col=1)
#}
#abline(h=H.pol, col=c("Red", "Yellow"), lwd=3)
#text(x=T*(1:dim(p.regs)[1])/10, y=9000, labels=row.names(p.regs))
#text(x=T*(1:dim(p.regs)[1])/10, y=5000, labels=round(p.regs,2))

#image(e.util, axes=FALSE, col=terrain.colors(10))
#axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
#mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
#mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
#mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```

Mean total population size under this policy is `r as.integer(round(mean(Results$N),0))` and the frequency of "Red", "Yellow", and "Green" seasons is 
```{r freq_no_high} 
t(p.regs)
```

# Thresholds at Maximum Sustained Yield (MSY)
```{r har_msy}
#Set parameters
Har = list(Red = rlnorm(Reps, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(7000), 0.1), 
              Green = rlnorm(Reps, log(15000), 0.2))

thresholds = seq(t.min, t.max, by=t.delta)
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
```

```{r opt_msy}
e.util = matrix(NA, nrow=length(thresholds), ncol=length(thresholds))
Reps <- 1000
for(ii in 1:length(thresholds)){for(jj in ii:length(thresholds)){
  
#define harvest thresholds
H.pol = list(Red = thresholds[ii], Yellow = thresholds[jj])

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )

for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <- utility(n = Results$N[i,], h = Results$Har[i,], fn=c("msy", "msy"))
}
discount <- matrix(d^c(0:(T-1)), nrow=Reps, ncol=T, byrow=TRUE)
e.util[ii,jj] <- mean(apply(Results$Utility*discount,1,sum))
}}

smart.spot = which(e.util==max(e.util, na.rm=TRUE), arr.ind = TRUE)
harvest.rule=list(Red=thresholds[smart.spot[1]], Yellow=thresholds[smart.spot[2]])
```

The optimal thresholds are 

```{r rule_msy}
harvest.rule
```

```{r sim_msy}
#Simulate dynamics at smart threshold
#Simulate dynamics at smart threshold
H.pol = harvest.rule

Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- HarvestPolicy(n = N.obs[t-1], plist=H.pol)
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], fn=c("msy", "msy"))
}

p.regs <- as.matrix(table(Results$HarPol)/(T*Reps))
regs <- matrix(0, 3, 1)
row.names(regs) <- c("Red","Yellow","Green")

#n.plot = 100
#plot(1,1, type="n",xlim=c(0, T), ylim=c(0, 1e5), xlab="Year", ylab="Obs. Population")
#for(i in 1:min(100,Reps)){
#  lines(1:T, Results$N.obs[i,], col=1)
#}
#abline(h=H.pol, col=c("Red", "Yellow"), lwd=3)
#text(x=T*(1:dim(p.regs)[1])/10, y=9000, labels=row.names(p.regs))
#text(x=T*(1:dim(p.regs)[1])/10, y=5000, labels=round(p.regs,2))

#image(e.util, axes=FALSE, col=terrain.colors(10))
#axis(1, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#axis(2, at=seq(0,1, by=0.2), labels=round(seq(t.min, t.max, length=6),0))
#points(smart.spot[1,1]/length(thresholds), smart.spot[1,2]/length(thresholds), pch=16)
#mtext("Closure Threshold (Red Zone)", 1, line=3, font=2)
#mtext("Restriction Threshold (Yellow Zone)", 2, line=3, font=2)
#mtext(paste("Expected Utility over ",T,"years, discount = ",1-d), 3, line=2, font=2)
```

Mean total population size under this policy is `r as.integer(round(mean(Results$N),0))` and the frequency of "Red", "Yellow", and "Green" seasons is 
```{r freq_msy} 
t(p.regs)
```

# Predicted Effect of 3-year Traditional Season
The predicted population response under an experimental 3-year "traditional" season was modeled.  We assumed the most pessimistic population scenario where harvest is was log normally distributed with log mean of 20,000 and log variance of 0.1.  All other parameter came from the posterior of the theta-logistic model.  We found that population decline was highly likely, often > 50% of the current population, see figures below.  

```{r exp_season}
T <- 3

Har = list(Red = rlnorm(Reps, log(4000), 0.1), #increased for crippling loss ~25%
              Yellow = rlnorm(Reps, log(7000), 0.1), 
              Green = rlnorm(Reps, log(20000), 0.1))

#Normal harvest (red) for three years
Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- "Red"
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], fn=c("msy", "msy"))
}
Normal <- Results
#High harvest for 3 years
Results = list(N=array(NA, dim=c(Reps,T)), 
               N.obs = array(NA, dim=c(Reps,T)), 
               HarPol = array(NA, dim=c(Reps,T)),
               Har = array(NA, dim=c(Reps,T)),
               Utility = array(NA, dim=c(Reps,T))
               )
for(i in 1:Reps){
  N <- N.obs<- numeric()
  N[1] <- N0[i]
  N.obs[1] <- N.last 
  
  for(t in 2:(T+1)){
    HarPol <- "Green"
    Results$HarPol[i,t-1] <- HarPol
    Results$Har[i,t-1] <- (1-exp(-a*N[t-1]))*Har[[HarPol]][i]  #added functional response  
    E.N.pre <- N[t-1]*exp(Rmax[i]*(1 - (N[t-1]/CC[i])^Theta[i])) - Results$Har[i,t-1]  
    N[t] <- ifelse(E.N.pre <= 0,0,rlnorm(1,log(E.N.pre/CC[i]), sigma.proc[i])*CC[i])
    N.obs[t] <- ifelse(N[t]<=0,0,rnorm(1,N[t]*q[i], cv.obs[i]*N[t]*q[i]))
  }
  Results$N[i,] <- N[-1]
  Results$N.obs[i,] <- N.obs[-1]
  Results$Utility[i,] <-  utility(n = Results$N[i,], h = Results$Har[i,], fn=c("msy", "msy"))
}
par(mfrow=c(1,2), pty="s")
hist((Results$N[,3]-Normal$N[,3])/Normal$N[,3], col="lightgray", main="", xlab="Proportional Change After 3 Years", breaks=20)
require(scales)
plot(0:3, 1:4, type="n", ylim=c(0, 5e5), ylab="Population", xlab="Year")
for(i in 1:1000){
  lines(0:3, c(N0[i],Results$N[i,]), col=alpha("gray50", 0.5))
}
par(mfrow=c(1,1), pty="m")
```



# Summary
The main conclusions of this work are summarized as:  

* The optimal policy, given the population model and utility function, is to close the season when the YKD summer index is < 26,000, have a modest harvest when the index is between 26,000 and 33,000, and have a traditional liberal harvest when the index exceeds 33,000.  
* The exact location of these thresholds does not matter much to the realized utility, as long a restrictions or season closure occurs when the index drops below 26000, approximately.  
* The exact location of harvest thresholds were highly dependent of the realized harvest model.  Therefore, the value of information on harvest is high.  This could best be addressed by observing the harvest consequence of a liberal, "traditional," harvest season during a short (1-3 years) series of seasons.  
* After gaining additional data on harvest under a liberal season, the model should be updated and new policy thresholds derived.  
* All modeled scenarios resulted in a relatively high frequency of closed seasons.  Thus, closed seasons should be expected in the future in order to maintain the emperor goose population and harvests into the future.  
* There is substantial risk of a large (~50%) population decline under the most pessimistic scenario of a 3-year experimental "traditional" season. Thus, careful monitoring and emergency closure of the experimental season may be necessary.