%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{An introduction to the metafolio R package}

% tw=68 for code
\documentclass[12pt]{article}
\usepackage{geometry}
\usepackage{url}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}

\setlength{\parskip}{10pt}
\setlength{\parindent}{0cm}

\textheight 21.5cm

\begin{document}

<<setup, include=FALSE, cache=FALSE>>=
library(knitr)
# set global chunk options
opts_chunk$set(fig.path='figure/', fig.align='center', fig.show='hold')
options(replace.assign=TRUE, width=90)
opts_knit$set(out.format = "latex")
opts_chunk$set(warning=FALSE, message=FALSE, comment=NA, tidy=FALSE,
refresh=TRUE, cache=FALSE, autodep=TRUE)
@

\title{An introduction to the \texttt{metafolio} \texttt{R} package}
\author{
Sean C. Anderson\textsuperscript{1*},
Jonathan W. Moore\textsuperscript{1,2},
Michelle M. McClure\textsuperscript{3},\\
Nicholas K. Dulvy\textsuperscript{1}
Andrew B. Cooper\textsuperscript{2}
}

\date{}
\maketitle

\textsuperscript{1}Department of Biological Sciences, Simon Fraser University, Burnaby BC, V5A 1S6, Canada

\textsuperscript{2}School of Resource and Environmental Management, Simon Fraser University, Burnaby, BC, V5A 1S6, Canada

\textsuperscript{3}Northwest Fisheries Science Center, National Marine Fisheries Service, Seattle, WA 98110, USA

\bigskip

The \texttt{metafolio} \texttt{R} package is a tool to simulate
metapopulations and apply financial portfolio optimization concepts to those
metapopulations. The package is primarily intended for salmon simulations, but
could be adapted for other taxonomic groups. This vignette accompanies the
paper \textit{Portfolio conservation of metapopulations under climate change},
which is In Press at Ecological Applications. In this vignette we will
describe the main functions in the \texttt{metafolio} package, show how to
re-create the main analyses in our paper, and demonstrate some of the included
plotting functions for exploring the output. You can get more detailed
installation instructions, read the source code, and report bugs on GitHub:
\url{https://github.com/seananderson/metafolio}

\section{Loading the package and getting help}

<<my-load-package, include=FALSE>>=
library(metafolio)
set.seed(1)
@

%\texttt{metafolio} will be submitted to CRAN and available through GitHub. Until then, you can install the package as follows:

%In an \texttt{R} console (version $\ge 3.1.0$), you can install \texttt{metafolio} by first setting your working directory with \texttt{setwd()} to wherever you saved the package and then running:

% <<install-source, eval=FALSE>>=
% install.packages("devtools")
%
% # OS X:
% devtools::install("metafolio_0.3.0.tgz")
%
% # Windows:
% devtools::install("metafolio_0.3.0.zip")
%
% # Or, install from source if you have a C++ compiler installed:
% devtools::install("metafolio_0.3.0.tar.gz")
% @

%The development version of \texttt{metafolio} can be installed via \texttt{devtools}:
%\textit{This website is a private repository while the paper is in review and so the following code will not work.}

%<<github-package, eval=FALSE>>=
%install.packages("devtools") # if needed
%devtools::install_github("metafolio", "seananderson",
  %dependencies = TRUE)
%@


You can load the package with:

<<load-package, eval=FALSE>>=
library(metafolio)
@

You can find a copy of this vignette and access the help files with:

<<help, eval=FALSE>>=
vignette("metafolio")
?metafolio
help(package = "metafolio")
@

\section{An example simulation}

The main simulation function in \texttt{metafolio} is \texttt{meta\_sim()}.
We'll start by running a simulation using the base-case parameter values from
the paper. First, we'll set up an \texttt{R} list object that contains the
argument values for the stationary environmental stochasticity scenario.

<<setup-arma>>=
arma_env_params <- list(mean_value = 16, ar = 0.1, sigma_env = 2,
  ma = 0)
@

The arguments \texttt{mean\_value}, \texttt{ar}, and \texttt{ma} refer to the
mean, autoregressive (AR), and moving-average (MA) components of an ARIMA
model. See \texttt{?arima} for further explanation of ARIMA models in
\texttt{R}. The argument \texttt{sigma\_env} refers to the standard deviation
of the environmental signal.

Next we'll run the simulation for one iteration. We'll simulate ten
populations and re-assess the fishery every five years.  Re-assessing the
fishery means that the simulation fits a Ricker curve to the observed
spawner-return data and updates the harvest-level targets. See the
accompanying paper for details.

<<base-case-run, cache=FALSE>>=
base1 <- meta_sim(n_pop = 10, env_params = arma_env_params,
  env_type = "arma", assess_freq = 5)
@

We can plot the output with the function \texttt{plot\_sim\_ts()}. The output is
shown in Figure~\ref{fig:plot-base-case-ts}.

<<plot-base-case-ts, fig.width=5, fig.height=6.5, out.width="4in", fig.cap="An example simulation with stationary environmental stochasticity and the base-case parameter values.">>=
plot_sim_ts(base1, years_to_show = 70, burn = 30)
@

\clearpage

\section{Exploring prioritization strategies}

One of the key elements to the analysis in our paper is the choice of
which populations to prioritize for conservation. We can try different
prioritization strategies by manipulating the ``investment weights'' in each
stream of salmon. In the case of salmon, we represent these with the Ricker
$b_i$ parameters, which indicate the maximum population capacity of streams
$i$ 1 through $n$.

In this example, we'll create one scenario in which we conserve response
diversity from across the spectrum of possible responses and another scenario
in which we conserve one half of the response diversity. We've set up the
weights carefully so that each metapopulation contains the same number
of populations (10) and the same total capacity. We set the Ricker $b_i$
parameters equal to the nominal level of five salmon for the streams that we
choose not to prioritize.

<<spatial-plans>>=
w_plans <- list()
w_plans[["balanced"]] <- c(5, 1000, 5, 1000, 5, 5, 1000, 5,
    1000, 5)
w_plans[["one_half"]] <- c(rep(1000, 4), rep(5, 6))
w <- list()
for(i in 1:2) { # loop over plans
 w[[i]] <- list()
 for(j in 1:80) { # loop over iterations
   w[[i]][[j]] <- matrix(w_plans[[i]], nrow = 1)
 }
}
@

We've now created a nested list of stream capacities (\texttt{w}). (\texttt{w}
refers to the mnemonic ``weights'' for investment weights.) The first level of the
list contains the two different scenarios The second level of the list contains
the $b_i$ values for each iteration. Each iteration will have re-sampled process
noise and observation error when run with \texttt{meta\_sim()}. Here, we're
keeping the $b_i$ values the same between iterations, but that might not be the
case if we wanted to stochastically simulate the $b_i$ values. We're only going
to run 80 iterations to reduce the runtime of the example, but in reality you
would likely want to run many more iterations.

We can now stochastically simulate with these strategies using the function\\
\texttt{run\_cons\_plans()} and plot the output with \texttt{plot\_cons\_plans()}
(Figure~\ref{fig:spatial-runs}).

<<spatial-runs, fig.width=5, fig.height=5, out.width="4in", cache=FALSE, fig.cap="Two spatial conservation strategies shown in risk-return space with stationary environmental stochasticity. The dots show simulated metapopulations and the contours show 25\\% and 75\\% quantiles across 80 simulations per scenario. The grey line indicates the efficient frontier across all simulated metapopulations. The efficient frontier represents the minimum expected mean growth rate for a given expected variance in growth rate.">>=
set.seed(1)
arma_sp <- run_cons_plans(w, env_type = "arma", env_params =
  arma_env_params)
plot_cons_plans(arma_sp$plans_mv,
  plans_name = c("Balanced", "One half"),
  cols = c("#E41A1C", "#377EB8"), xlab = "Variance of growth rate",
  ylab = "Mean growth rate")
@

\section{Generating alternative environmental time series}

We can use the function \texttt{generate\_env\_ts()} to create alternative
environmental time series. The function can generate five kinds of time
series: sine waves, regime shifts, linear changes, and constant values. Each
type has its own set of parameter arguments that are passed in a list
format. See \texttt{?generate\_env\_ts} for examples of all of these.

We can see demonstrations of all environmental time series types with the
default argument values with the following code (Figure~\ref{fig:env-eg}).

<<env-eg, fig.width=5, fig.height=6, out.width="3in", fig.cap="Example environmental time series.", fig.pos="hpb">>=
types <- c("sine", "arma", "regime", "linear", "constant")
x <- list()
for(i in 1:5) x[[i]] <- generate_env_ts(n_t = 100, type = types[i])
par(mfrow = c(5, 1), mar = c(3,3,1,0), cex = 0.7)
for(i in 1:5) plot(x[[i]], type = "o", main = types[i])
@

\clearpage

\section{Additional plotting functions}

We can visualize the variability in the Ricker $a$ parameters using the
function \texttt{plot\_rickers()} (Figure~\ref{fig:plot-rickers}).

<<plot-rickers, fig.width=7, fig.height=4.5, out.width="5in", fig.pos="hpb", fig.cap="Ricker curves from a sample of 40 years in the example simulation. Each panel represents a different stream population. Population 1 is more productive in cool conditions and population 10 is more producitive in warm conditions. The colour of the Ricker curves represents the relative temperatue in that year (warm: red; cool: blue). The grey shaded area represents the variation in spawners observed within the simulation.">>=
plot_rickers(base1, pal = rep("black", 10))
@

We can look at the correlation between salmon returns in the various
streams using the function \texttt{plot\_correlation\_between\_returns()}
(Figure~\ref{fig:return-correlations}).

<<return-correlations, fig.width=5, fig.height=5, out.width="5in", fig.cap="A plot comparing the log(returns) between each population. The population IDs are coloured from warm tolerant (warm colours) to cool tolerant (cool colours). Note how populations 1 and 10 have asynchronous returns whereas populations with more similar thermal-tolerance curves (say populations 9 and 10) have more synchronous dynamics. Populations with thermal tolerance curves in the middle (e.g.~population 6) are less correlated with other populations. Their population dynamics end up primarily driven by demographic stochasticity and less so by temperature-induced systematic changes in productivity.">>=
plot_correlation_between_returns(base1)

@

\section{Optimizing metapopulation portfolios}

A standard procedure in financial portfolio management is to determine
optimal investment weights of the assets in a portfolio. The portfolios made
up of these optimal investment weights are referred to as the efficient
frontier. This efficient frontier describes a set of portfolios which
have minimal risk for a specified level of return or maximum return for a
specified level of risk. While it would be complicated to determine the
optimal metapopulation portfolios via algebra we can do so by letting
\texttt{metafolio} sample from possible investment weights --- a form of Monte
Carlo sampling (Figure~\ref{fig:monte-carlo-plot}).

<<monte-carlo-eg, eval=FALSE>>=
set.seed(1)
weights_matrix <- create_asset_weights(n_pop = 6, n_sims = 3000,
  weight_lower_limit = 0.001)
mc_ports <- monte_carlo_portfolios(weights_matrix = weights_matrix,
  n_sims = 3000, mean_b = 1000)
@

<<monte-carlo-eg-load, eval=TRUE, echo=FALSE>>=
# To make the vignette compile more quickly:
# port_vals <- mc_ports$port_vals
# save(port_vals, file = "port_vals.rda")
weights_matrix <- create_asset_weights(n_pop = 6, n_sims = 3000,
  weight_lower_limit = 0.001)
load("port_vals.rda")
mc_ports <- list()
mc_ports$port_vals <- port_vals
@

<<monte-carlo-plot, fig.width=7, fig.height=4, out.width="5in", fig.cap="Efficient frontier of metapopulation portfolios (red dots). Each dot represents a different set of weights of the Ricker $b$ parameters. The colours in the right panel correspond to the five populations with warm tolerant populations in warmer colours and cool tolerant populations in cooler colours.">>=
col_pal <- rev(gg_color_hue(6))
ef_dat <- plot_efficient_portfolios(port_vals = mc_ports$port_vals,
  pal = col_pal, weights_matrix = weights_matrix)
@

\end{document}