Hi @ignacio!
I’m not the great @rtrangucci, but let me give it a try.
So, unfortunately the notation changed a bit over the course of this thread—but, basically what’s u[i] in the code is \mu_i in equation, right? So far so good.
In the answers above you can find that \mathbb{E}[Y_{i,t}]=\mu_i, so in the code that is the u[i]—the mean of y[1,i].
For the variance we have
\begin{align}
\text{Var}[Y_{i,t}]&=\text{Var}[\delta Y_{i,t-1}]+\text{Var}[\epsilon_{i,t}] \\
\text{Var}[Y_{i,t}]&=\delta^2\text{Var}[Y_{i,t-1}] +\sigma^2_\epsilon\\
\text{Var}[Y_{i,t}] - \delta^2\text{Var}[Y_{i,t-1}]&=\sigma^2_\epsilon.
\end{align}
Now, I think we need to assume \text{Var}[Y_{i,t}] = \text{Var}[Y_{i,t-1}], which is reasonable (assume iid residuals / a stationary process). Then,
\begin{align}
\text{Var}[Y_{i,t}] - \delta^2\text{Var}[Y_{i,t}]&=\sigma^2_\epsilon\\
\text{Var}[Y_{i,t}](1 - \delta^2)&=\sigma^2_\epsilon\\
\text{Var}[Y_{i,t}]&=\frac{\sigma^2_\epsilon}{(1 - \delta^2)}\\
\text{Sd}[Y_{i,t}]&=\frac{\sigma_\epsilon}{\sqrt{1 - \delta^2}},
\end{align}
which is the sigma_e / sqrt(one_minus_delta_sq) part of the code.
I learned this at the Helsinki StanCon Tutorial with Jonah. Have a look—he also discusses a GP formulation of this AR(1) process.
Cheers! :)