Deriving method of moments estimator for AR(1) process

The method of moments estimator for AR processes can be had with the Yule-Walker equations. But how is it derived?

The equation for AR(1):

$$Y_t =aY_{t-1}+\epsilon_t$$

Where $\epsilon$ ~ $N(0,\sigma^2 )$.

So the moment conditions are $E(\epsilon)=0, E(\epsilon^2)=\sigma^2$

But the "standard" solution can't be used as $E(\epsilon)=E(Y_t-aY_{t-1})=0$ is true for any $a$.

Multiply the equation by $Y_{t-1}$

$Y_tY_{t-1}=aY_{t-1}Y_{t-1}+\epsilon_tY_{t-1}$

$E(Y_tY_{t-1})=E(aY_{t-1}Y_{t-1}+\epsilon_tY_{t-1})$

$E(Y_tY_{t-1})=aE(Y_{t-1}^2)$

$E(Y_tY_{t-1})/E(Y_{t-1})=a$

$E(Y_tY_{t-1})/E(Y_{t})=a$

And for $\sigma^2$:

$E((Y_t+aY_{t-1})^2)=\sigma^2$

$E(Y_t^2+a^2Y_{t-1}^2-2aY_tY_{t-1})=\sigma^2$

$E(Y_t^2+a^2Y_{t-1}^2-2aY_tY_{t-1})=\sigma^2$

$E(Y_t^2+a^2Y_{t-1}^2-2aY_{t-1}(aY_{t-1}+\epsilon_t))=\sigma^2$

$E(Y_t^2+a^2Y_{t-1}^2-2a^2Y_{t-1}^2)=\sigma^2$

$E(Y_t^2-a^2Y_{t-1}^2)=\sigma^2$

$E(Y_t^2)-a^2E(Y_{t-1}^2)=\sigma^2$

$(1-a^2)E(Y_{t}^2)=\sigma^2$