bioeconomic-milner.R

# Bioeconomic model from Milner-Gulland & Rowcliffe, p. 163-...

# define constants
K = 1000 # carrying capacity
r = 0.2 # intrinsic rate of increase
P = 105 # price
a = 200 # constant for cost calculation
b = 0.2 # constant for cost calculation
s = 10 # SD of the distribution for the cost
N1 = 500 # initial pop size
Years = 30 # run the model for 30 years
Hunters = 200 # nb of hunters

Pop <- rep(0,Years+1) # pop size
Harvest <- rep(0,Years) # harvest
Prod <- rep(0,Years) # productivity

Pop[1] <- N1

for (t in 1:Years){
	
	PrF <- r * Pop[t] * (1 - Pop[t] / K) # logistic productivity
	Prod[t] <- round(PrF) # nb of births must be a whole number
	Cost <- a - b * Pop[t] # mean hunter cost
	
	for (i in 1:Hunters){
		
		zval <- rnorm(1,0,1)
		IndC <- Cost + zval * s # hunter cost
		if (IndC < 0) IndC = 0 # ensure no negative costs
		B = P - IndC # profitability
		if (B > 0) Harvest[t] <- Harvest[t] + 1 # hunting decision
		
	}
	
	if (Harvest[t] > Pop[t] + Prod[t]) Harvest[t] <- Pop[t] + Prod[t] # can't harvest more than is there
	Pop[t+1] = Pop[t] + Prod[t] - Harvest[t]
	
}
 
Pop
Prod
Harvest 

par(mfrow=(c(3,1)))
plot(1:Years,Pop[-1],xlab='Years',ylab='Pop size',type='o')
plot(1:Years,Prod,xlab='Years',ylab='Productivity',type='o')
plot(1:Years,Harvest,xlab='Years',ylab='Harvest',type='o')