# Variance of Stochastic Growth Rate for Average Population Size Over Time

I'm interested in calculating the variance of the stochastic growth rate for the average population size over time,
$$\ln \beta = \lim_{t\to\infty}\dfrac{1}{t}\ln V[N(t)]$$ which shows that the variance of $N(t)$ grows at a rate $\beta$. The equation to compute this is given in the book "Matrix Population Models" by Hal Caswell. The variance can be calculated as, $$\beta = \lambda_1^{(\bf{B}_2)}$$ where $\lambda_1$ is the dominant eigenvalue, \begin{align*} \bf{B_2} &= \bf{F_2}[\bf{P}\otimes(\bf{I}\otimes\bf{I})]\mathrm{,} \\[1em] \bf{F}_2 &= \mathrm{diag}({\bf A_1 \otimes A_1, \ldots, A_k \otimes A_k})\mathrm{,} \end{align*} $\bf A$ are the population projection matrices of dimension $s$ from $k$ independent environments, $\bf P$ is the column-stochastic transition matrix of dimension $k$, and $\bf I$ is the identity matrix of dimension $s$. I tried doing this in R using the goldenheather (Hudsonia) example from the Caswell book that can be loaded into R as follows,

library(popbio); library(Matrix); data(hudsonia)


I ended up with a value for $\beta$ of 0.9366351. This cannot be right, the variance is way too high. I must be doing it wrong, or perhaps there is a mistake in one of the equations (very, very likely the former).

Here's my code (if it helps),

s = NROW(hudsonia[])
k = length(hudsonia)
P = Matrix(1 / k, k, k)
I = Diagonal(s, 1)
F2 = bdiag(lapply(hudsonia, function(i) kronecker(i, i)))
I2 = kronecker(I, I)
P2 = kronecker(P, I2)
B2 = F2 %*% P2
lambda(B2)