#Exercises Yield and economic model III - capture fishery
#Multiple projections for risk analysis. Projection matrix method.
#as is Per Recruit model

#Author: Jorge Santos, Tromsų, 2013, 2014

#Monte Carlo analysis; here it explores the uncertainty in Recruitment and prices of fish.
#Projects population for t years, n times, to calculate risks in Yield, Stock biomass, 
#   Spawning stock and Profit, given a selected fishing strategy (an f (or F) and tc combination).
#Install necessary libraries (packages).
#Adjust 1ary and 2ary input according to species biology and fishery.
#Check suitability of random variates (uncertain input variables).
#If available, indicate acceptable risk levels (minimum acceptable levels)
#To run script highlight the whole script and Run (or ctrl-enter).
#Aggregated output in Results.mat BiomassTraj.mat files, 2 tables of output statistics,
# one figure and several histograms and cumulative histograms of distributions in the final year.

#needs
library(Hmisc)

#age structured growth - INPUT
Rec=1                         #Recruitment mean, subject to random error (see below)
M=?                           #Inst. rate natural mortality (year-1)
f=?                           #desired fishing effort f (e.g. days)
Linf=?                        #Linf von Bertalanffy (cm)
K=?                           #K von Bertalanffy
t0=?                          #t0 von Bertalanffy (y) (set to 0)
a=?                           #Constant weight length (g/cm-3)
b=?                           #Power weight length relation (unitless)
tr=?                          #Age recruitment to fishery (y)
tc=?                          #Age 1st capture (knife-edge) (y)
ts=?                          #Age maturity (knife-edge) (y)
price1=?                      #Sales price constant ($/kg), price=price1+price2*ln(Weight,kg)
price2=?                      #Sales price -size multiplier (if =0, then price constant)
#pricevariation               #price has now a random component (see below)
cost=?                        #Production (capture) cost ($/f unit); variable cost only
q=?                           #catchability (=F/(f^elast)); stock and fishery specific
elast=?                       #elasticity (exponent); denotes aggregation of stock (elast<1) or search time (>1); specific

#cost=20000                   #OLD Production (capture) cost ($/F unit)

#Random factors
Rlnormmean=log(Rec)           #Number of Recruits (R) (mean in log scale) lognormal variate
Rlnormsd=?                    #sd of R (in log scale) lognormal variate
pricenormsd=?                 #random normal component of price (sd in price estimate)

#Limit values and risk   (unit= tons and US$ millions, in Krill fishery example)
TotBiomin=NA                  #Minimum acceptable total exploitable biomass (or NA, not declared)
SpawnBiomin=NA                #Minimum acceptable spawning biomass (or NA, not declared)
Yieldmin=NA                   #Minimum acceptable yield (or NA, not declared)
Profitmin=N                   #Minimum acceptable profit (or NA, not declared) 


#time steps (2ary input)
sims=100                      #Number of simulations
ages=30                       #Maximum age of species (y), be generous (large ages)
years=30                      #Projection horizon (y)
step=1                        #Time step of analysis [0,1] (y), = number of cohorts/year
                              #     =age groups/year; R=1 each step

###############
#Simulation control 

#Recruitment (variable, lognormal distributed with mean R=1)
R.vec=rlnorm(sims,mean=Rlnormmean,sd=Rlnormsd)
#Price (variable by sim, normal distributed, variable by age if price2!=0)
#calculates price index= price of a fish 1kg in size in that simulation
PriceIndx.vec=rnorm(sims,mean=price1,sd=pricenormsd)

Results.mat=matrix(NA,nrow=10,ncol=sims)
rownames(Results.mat)=c("R","Survivors", "TotBio", "SpawnBio",
                "CatchN","Yield", "Costs","PriceIndx","Revenue","Profit")
BiomassTraj.mat=matrix(NA,nrow=years/step+1,ncol=sims)

for (j in 1:sims) {
R=R.vec[j]
Results.mat["R",j]=R
PriceIndx=PriceIndx.vec[j]
Results.mat["PriceIndx",j]=PriceIndx

#CALCULATIONS

F=q*f^elast                                         #NEW F calculated from f

#time steps
time=years/step
Age=c(seq(tr,(tr+ages),by=step))

#Survival and mortality (instantaneous and fractional)
Fage=c(rep(0,times=((tc-tr)/step)),rep(F*step,times=(tr+ages-tc)/step+1))
Mage=c(rep(M*step,times=(ages/step)+1))
Zage=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) Zage[i]= Mage[i]+Fage[i]
SurvAge=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) SurvAge[i]=exp(-Zage[i])

#Growth and weight (von Bertalanffy)
Wage=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) Wage[i]=a*(Linf*(1-exp(-K*(Age[i]-t0))))^b

#(Lesley) A-matrix, N-matrix
A.mat=matrix(0,nrow=((ages/step)+1),ncol=((ages/step)+1))
A.mat[1,1]=1
for (i in 1:(ages/step)) A.mat[i+1,i]=SurvAge[i]
#Initial abundance vector
N0=c(R,rep(0,times=ages/step))
#Projection matrices
N.proj=matrix(0,nrow=nrow(A.mat),ncol=time+1)
N.proj[,1]=t(N0)
for (i in 1:time)N.proj[,i+1]=A.mat%*%N.proj[,i]
#matplot(0:time+1,t(N.proj),type="l",lty=1:3,
#        col=1,ylab="Age abundance N", xlab="Year")
N.totals=apply(N.proj,2,sum,na.rm=TRUE)
N=N.totals[time+1]
Results.mat["Survivors",j]=N

#Biomass matrix
B.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
B.proj=N.proj*Wage
B.totals=apply(B.proj,2,sum,na.rm=TRUE)
B=B.totals[time+1]
Results.mat["TotBio",j]=B
BiomassTraj.mat[,j]=B.totals

#Spawning biomass matrix
SB.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in ((ts-tr)/step+1):((ages/step)+1)) SB.proj[i,]=N.proj[i,]*Wage[i]
SB.totals=apply(SB.proj,2,sum,na.rm=TRUE)
SB=SB.totals[time+1]
Results.mat["SpawnBio",j]=SB

#Catch matrix
C.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in 1:(ages/step+1))C.proj[i,]=N.proj[i,]*(1-exp(-Fage[i]))
C.totals=apply(C.proj,2,sum,na.rm=TRUE)
Catch=C.totals[time+1]
Results.mat["CatchN",j]=Catch
  
#Yield matrix
Y.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in 1:(ages/step+1))Y.proj[i,]=C.proj[i,]*(Wage[i])              #/1000  #NEW (/1000) yield in kgs - not useful in per Recruit analysis
Y.totals=apply(Y.proj,2,sum,na.rm=TRUE)
Yield=Y.totals[time+1]
Results.mat["Yield",j]=Yield

#Economics (Price at size ($/kg), Profit=total cost-revenue)
Pricsize=numeric((ages/step)+1)
for(i in 1:(ages/step+1)) Pricsize[i]=PriceIndx+price2*log(Wage[i]/1000)

# Results.mat["Costs",j]=F*cost                           #OLD costs function of F
Results.mat["Costs",j]=f*cost                             #NEW costs function of f (=effort)

Revenue.proj=matrix(0,nrow=nrow(N.proj),ncol=ncol(N.proj))
for (i in 1:(ages/step+1))Revenue.proj[i,]=Y.proj[i,]*(Pricsize[i])     #/1000      #NEW /1000 to make $/g - not useful in per Recruit
R.totals=apply(Revenue.proj,2,sum,na.rm=TRUE)
Revenue=R.totals[time+1]
Results.mat["Revenue",j]=Revenue

Results.mat["Profit",j]=
  Results.mat["Revenue",j]-Results.mat["Costs",j]      
Profit=max(Results.mat["Profit",j])


}

#transpose results matrix
Results.mat=t(Results.mat)

#TABLES (summary statistics, including percentiles)
#needs
require(Hmisc)
summary(Results.mat)
describe(Results.mat)

#PLOTS

#INPUT PARAMETERS IN RISK ANALYSIS

#Plot - Distribution of (uncertain) simulated Recruitment (R) - a sample

hist(Results.mat[,"R"],probability=TRUE,col=gray(.9),
     xlab="Simulated R", main=c(paste("Recruitment (mean= ",
     round(mean(Results.mat[,"R"]),0),") and 95% PI")),cex.main=0.8)
abline(v=mean(Results.mat[,"R"]),lty=1)
abline(v=quantile(Results.mat[,"R"],p=c(0.0275,0.975)),lty=3)
abline(v=mean(Results.mat[,"R"]),lty=1)

#Plot - Distribution of (uncertain) simulated Price index

hist(Results.mat[,"PriceIndx"],probability=TRUE,col=gray(.9),
     xlab="Simulated Price Index (of 1 kg fish product)", main=c(paste("Price (mean= ",
     round(mean(Results.mat[,"PriceIndx"]),2),") and 95% PI")),cex.main=0.8)
abline(v=mean(Results.mat[,"PriceIndx"]),lty=1)
abline(v=quantile(Results.mat[,"PriceIndx"],p=c(0.0275,0.975)),lty=3)
abline(v=mean(Results.mat[,"PriceIndx"]),lty=1)

#RESULTS OF SIMULATIONS

#Trajectories of population stock biomass (sub-sample of 100 simulations)

matplot(1:(years/step+1),BiomassTraj.mat,type="l",                       
        main=c(paste("Population trajectories (n=",sims," simulations)")),
        xlab="Time (y)",ylab="Biomass / R (g)",xlim=c(5,30),cex.main=0.8)
abline(h=mean(BiomassTraj.mat[15:nrow(N.proj),]),lty=1,lwd=3)             

#Plot - Distributions of Total stock biomass in last year of series

hist(Results.mat[,"TotBio"],probability=TRUE,col=gray(.9),               
     xlab="Total stock biomass/ R  (g)", cex.main=0.8,
     main="Stock bio distribution, risk level, mean and 95% PI")
abline(v=mean(Results.mat[,"TotBio"]),lty=1)                             
abline(v=quantile(Results.mat[,"TotBio"],p=c(0.0275,0.975)),lty=3)       
abline(v=mean(Results.mat[,"TotBio"]),lty=1)                             
abline(v=TotBiomin,lty=8, col=2, lwd=2)                             #input in grams

plot(ecdf(Results.mat[,"TotBio"]), verticals = TRUE, do.points = FALSE,  
     main="Stock biomass / R cum distrib, mean and risk level", xlab="Biomass (g)", 
     ylab="Cumulative probability", cex.main=0.8)
minor.tick(nx=5, ny=5, tick.ratio=0.5)
abline(v=mean(Results.mat[,"TotBio"]),lty=1)                              
abline(v=TotBiomin,lty=8, col=2, lwd=2)                               #input in g

#Plot - Distributions of Spawning stock biomass in last year of series

hist(Results.mat[,"SpawnBio"],probability=TRUE,col=gray(.9),               
     xlab="Spawning stock biomass/ R  (g)", cex.main=0.8,                          
     main="Spawn biomass distribution, risk level, mean and 95% PI")
abline(v=mean(Results.mat[,"SpawnBio"]),lty=1)                            
abline(v=quantile(Results.mat[,"SpawnBio"],p=c(0.0275,0.975)),lty=3)      
abline(v=mean(Results.mat[,"SpawnBio"]),lty=1)                            
abline(v=SpawnBiomin,lty=8, col=2, lwd=2)                              #input in g                                   

plot(ecdf(Results.mat[,"SpawnBio"]), verticals = TRUE, do.points = FALSE, 
     main="Spawn bio cum distrib, mean and risk level", cex.main=0.8, 
     xlab="Spawning biomass (g)", ylab="Cumulative probability")               
minor.tick(nx=5, ny=5, tick.ratio=0.5)
abline(v=mean(Results.mat[,"SpawnBio"]),lty=1)                            
abline(v=SpawnBiomin,lty=8, col=2,lwd=2)                               #input in g


#Plot - Distributions of Yield in last year of series

hist(Results.mat[,"Yield"],probability=TRUE,col=gray(.9),                     
     xlab="Yield / R (g)", cex.main=0.8, 
     main="Yield: annual distribution, risk level, mean and 95% PI")
abline(v=mean(Results.mat[,"Yield"]),lty=1)                                  
abline(v=quantile(Results.mat[,"Yield"],p=c(0.0275,0.975)),lty=3)            
abline(v=Yieldmin,lty=8, col=2, lwd=2)                                #input in g

plot(ecdf(Results.mat[,"Yield"]), verticals = TRUE, do.points = FALSE,        
     main="Yield / R annual cum distrib, mean and risk level", xlab="Yield / R (g)",    
     ylab="Cumulative probability", cex.main=0.8)
minor.tick(nx=5, ny=5, tick.ratio=0.5)
abline(v=mean(Results.mat[,"Yield"]),lty=1)                                 
abline(v=Yieldmin,lty=8, col=2, lwd=2)                               #input in g


#Plot - Distributions of Profit in last year of series

hist(Results.mat[,"Profit"],probability=TRUE,col=gray(.9),                
     xlab="Profit/ R  (US$)", cex.main=0.8,                                   
     main="Profit/ R : annual distribution, risk level, mean and 95% PI")
abline(v=mean(Results.mat[,"Profit"]),lty=1)                              
abline(v=quantile(Results.mat[,"Profit"],p=c(0.0275,0.975)),lty=3)        
abline(v=Profitmin,lty=8, col=2, lwd=2)                            #input in US$ 

plot(ecdf(Results.mat[,"Profit"]), verticals = TRUE, do.points = FALSE,   
     main="Profit/ R : annual cum dist, mean and risk level", cex.main=0.8,
     xlab="Profit / R (US$ )", ylab="Cumulative probability")                   
minor.tick(nx=5, ny=5, tick.ratio=0.5)
abline(v=mean(Results.mat[,"Profit"]),lty=1)                             
abline(v=Profitmin,lty=8, col=2, lwd=2)                            #input in US$ 


#END 1

# keep only necessary variables in workspace
#rm(list=(ls()[ls()!="Results.mat"]))
rm(list= ls()[!(ls() %in% c('Results.mat','BiomassTraj.mat'))])

#to pass results files (e.g. Results.mat) to another software as text (.txt) 
#         or Excel (.xls) file do:
#set working directory (Files in bottom right pane of RStudio)
#write.table(Results.mat,file = "Results.txt")
#write.table(Results.mat,file = "Results.xls")
# obs: in Excel may have to import as space delimited file

#END 2 (not run)


#LICENSE & REFERENCE

#Santos, J. 2015. FI??H IT 1.0 - Student Manual: A Training System for Aquatic Resource Managers. Septentrio Educational 2015(3).
#DOI:http://dx.doi.org/10.7557/se.2015.3

#Chapter 7a R-script YR I PR  by projection contour F tc
#DOI: http://dx.doi.org/10.7557/8.3602


#This work is licensed under a Creative Commons Attribution 4.0 International License.