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[[1]])
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)
