How do we identify parameters in this simple model? Consider the following model:
$$ y_{it}=\nu_{it}+\epsilon_{it}$$
$$\nu_{it}=\rho \nu_{it-1}+\zeta_{it}$$
Where $y_{it}$ is the income for $i$ at time $t$. $\epsilon_{it}$ is the idiosyncratic income shock. $\nu_{it}$ denotes the permanent component that follows an AR(1) process, like unobserved productivity. The only observed data is $\{y_{it}\}$. We assume that $\epsilon_{it}\sim N(0,\sigma_\epsilon^2)$ and $\zeta_{it}\sim N(0,\sigma_{\zeta}^2)$, and $0<|\rho|<1$ so that the process is stationary.
I think $\rho$ can be identified as follows:
We have $v_{it}=y_{it}-\epsilon_{it}$ from the first equation, and $v_{it}=\rho v_{it-1}+\zeta_{it}$ from the second. Therefore, we have
\begin{equation*}
    y_{it}-\epsilon_{it}=\rho v_{it-1}+\zeta_{it}
\end{equation*}
We also have $v_{it-1}=y_{it-1}-\epsilon_{it-1}$. Thus,
\begin{equation*}
    y_{it}-\epsilon_{it}=\rho(y_{it-1}-\epsilon_{it-1})+\zeta_{it}
\end{equation*}
Then we have
\begin{equation*}
    \rho=\frac{ y_{it}-\epsilon_{it}-\zeta_{it}}{y_{it-1}-\epsilon_{it-1}}
\end{equation*}
Thus, $\rho$ is identified.
But how can we identify $\sigma_{\epsilon}^2$ and $\sigma_{\zeta}^2$, the two variances, in this case?

The problem you are trying to solve (i.e. inferring the parameters of a state space system) is called system identification. To the best of my knowledge, there is no closed-form, analytical solution to find an estimate of $\rho$. Given the relative simplicity of your model (univariate, linear trend), I would recommend to use the Expectation-Maximization algorithm (EM), which is an iterative method to find (local) maximum likelihood or maximum a posteriori (MAP) estimates of parameters in statistical models, where the model depends on unobserved latent variables (cf here).
This method is reviewed in the BRML textbook (Chapter 24.5.3) accessible online.
