# Set-up parameters for stochastic exponential growth
N0 = 20
yrs = 300

mu.Lambda=1.06
sd.Lambda=sqrt(0.05)
sd.Lambda

# calculate variance and mean for stochastic 'r' or 'r.bar'
r.bar.var = (sd.Lambda^2)/(mu.Lambda^2)
r.bar = log(mu.Lambda)- r.bar.var/2
r.bar
sd.r.bar = sqrt(r.bar.var)
sd.r.bar

# plot the expected trajectory if there were no temporal variation
# working on log scale for population size
plot(1:(yrs+1),log(N0*mu.Lambda^(0:yrs)),'l',col='blue', ylim=c(0,20),
  xlab="Time", ylab="ln N")

# add a line for the expected stochastic trajectory
# again working on log scale for population size
lines(1:(yrs+1),log(N0*exp(r.bar)^(0:yrs)),'l',col='green')

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# Work with simulated stochastic trajectories #
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

# choose how many stochastic trajectories to run
sims=5

# create object to use for storing random growth rates
rand.r = matrix(data=NA, nrow=yrs, ncol=sims)

# create object to use for storing population sizes
N = matrix(data=NA, nrow=yrs+1, ncol=sims)

# store starting value in 1st row of each column
N[1,]=N0

# generate dynamics and store values
for (j in 1:sims) {
  for (i in 1:yrs) {
     # generate a random growth rate
     rand.r[i,j] = rnorm(1, r.bar, sd.r.bar)
     # apply growth rate to population for 1 time step
     N[i+1,j]=N[i,j]*exp(rand.r[i,j])
  }
  # add simulated trajectory for each simulation
  # We'll slow things down to better see each simulation trajectory
  Sys.sleep(0.25)  
  lines(1:(yrs+1), log(N[,j]), 'l', col='red')
}

# Re-plot the expected lines in thicker lines over red lines
lines(1:(yrs+1),log(N0*mu.Lambda^(0:yrs)),'l',col='blue',lwd=3)
lines(1:(yrs+1),log(N0*exp(r.bar)^(0:yrs)),'l',col='green',lwd=3)

# calculate realized growth rates from simulations
RGR=(N[(yrs+1),]/N[1,])^(1/yrs)
RGR
mean(RGR)
exp(r.bar)

# calculate the sum of 'r' values for each simulation (i.e., by column)
# and check that these can be used along with N0 
# to get final population size
r.sums=apply(rand.r,2,sum)
alt.N=exp(log(N0)+r.sums)
diffs=alt.N - N[(yrs+1),] # should be ~0
summary(diffs) # check that all are ~0