I suppose you mean $P_t$ and $P_{t-1}$ are i.i.d. Note that we may express
$$ P_t = e^{\mu + \sigma Z_t}, P_{t-1} = e^{\mu + \sigma Z_{t-1}}$$
where $Z_t, Z_{t-1}$ are i.i.d. standard normal. Then
$$ R_t = \frac {P_t - P_{t-1}} {P_{t-1}} = \frac {P_t} {P_{t-1}} -1
= \frac {e^{\mu + \sigma Z_t}} {e^{\mu + \sigma Z_{t-1}}} -1
= e^{\sigma (Z_t - Z_{t-1})} - 1$$
Since $Z_t - Z_{t-1} \sim \mathcal{N}(0, 2)$, we have $e^{\sigma (Z_t - Z_{t-1})} \sim \text{lognormal}(0,2\sigma^2)$ and thus the resulting $R_t$ is a shifted lognormal.
I had a go at calculating it with the Delta method like they state and also don't get exactly the same result the authors get. For some transformation $g(X)$ of a random variable $X$, the (univariate ) Delta method states that
$$\sqrt{n}[g(X) - g(\mu)] \xrightarrow{D} \mathcal{N}(0, \sigma_X^2[g'(\mu)]^2) $$
In other words, the variance of $g(X)$ should be given by $\sigma_X^2[g'(\mu)]^2$. In this case call $X = ln(\alpha)$ such that $g(X) = e^X = \alpha$. Then also $g'(X) = e^X$. So
$$\sigma_{g(X)} \approx \sqrt{0.806}*e^{-1.327} \approx 0.2381545$$
Then $S.E. = \sigma_{g(X)} / \sqrt{171} \approx 0.0182121$. Slightly under the quoted 0.020.
Maybe they used the multivariate Delta method and the discrepancy is due to correlations with a different parameters being estimated, for example?
Edit: with the additional information, the result of the multivariate Delta method should be the same given that the variables are transformed independently. Let's check. In the multivariate case the approximate covariance matrix of the parameters after transformation $g(X, Y)$ should be given by
$$\bf{J(\mu)} \Sigma \bf{J(\mu)^T}$$
With $\bf{J}$ the Jacobian of the transformation, this is very similar in form to the univariate case. We call $X = ln(\alpha)$ and $Y = ln(\gamma)$ so that $g_1(X, Y) = e^X = \alpha$ and $g_2(X, Y) = e^Y = \gamma$
Now the Jacobian is
$$ \bf{J} = \begin{pmatrix}
\frac{\partial g_1}{\partial X} & \frac{\partial g_1}{\partial Y} \\
\frac{\partial g_2}{\partial X} & \frac{\partial g_2}{\partial Y}
\end{pmatrix}
=
\begin{pmatrix}
e^X & 0 \\
0 & e^Y
\end{pmatrix}$$
And the original covariance is
$$\Sigma =
\begin{pmatrix}
0.806 & 0.306 * \sqrt{0.806*0.803} \\
0.306 * \sqrt{0.806*0.803} & 0.803
\end{pmatrix}
=
\begin{pmatrix}
0.806 & 0.246 \\
0.246 & 0.803
\end{pmatrix}$$
So that the covariance matrix of $\alpha$ and $\gamma$ is approximated by
$$
\Sigma_{\alpha, \gamma}
\approx
\begin{pmatrix}
e^{-1.327} & 0 \\
0 & e^{-0.735}
\end{pmatrix}
\begin{pmatrix}
0.806 & 0.246 \\
0.246 & 0.803
\end{pmatrix}
\begin{pmatrix}
e^{-1.327} & 0 \\
0 & e^{-0.735}
\end{pmatrix}
$$
$$
=
\begin{pmatrix}
0.806 (e^{-1.327})^2 & 0.246 e^{-1.327} e^{-0.735} \\
0.246 e^{-1.327} e^{-0.735} & 0.803 (e^{-0.735})^2
\end{pmatrix}
=
\begin{pmatrix}
0.0567 & 0.0313 \\
0.0313 & 0.1846
\end{pmatrix}
$$
So the formula for the univariate case reappears on the diagonal. This gives standard errors of $0.0182$ and $0.0329$ respectively for $\alpha$ and $\gamma$. Still off from what the authors did, perhaps my calculation is off somewhere, or the software the authors used doesn't calculate it exactly this way (maybe higher order Delta method, for example).
Best Answer
The actual model here is the one described by the equation with the tilde (~). In that model:
The equation with $\Delta$'s is less a definition of the model, and more a guideline to its numerical implementation.