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?

Given the Gaussian distribution of the error terms $\varepsilon$ and $\zeta$, this is a Kalman filtering problem. See my answer here for a worked-out example.
Generally speaking, $\rho$ and both $\sigma$ terms can be estimated by maximum likelihood even if we only observe $y$.
The intuition of how it works is that once we imposed the structure of the two equations for $\nu$ and $y$ (so called state space representation) in addition to the Guassianity of errors, we can iteratively construct the likelihood. That is, given what we know in period $t$, the distribution for the realization of $y_{t+1}$ and $\nu_{t+1}$ is Gaussian with a particular, time-dependent, mean and variance. We can construct the whole likelihood this way and then search for $\rho$ and $\sigma$s that maximize it.
For an illustration, let's consider an example from this this book.
# this based on the codes acompanying 
# Time Series Analysis and Its Applications, Ed 4 
# by Shumway & Stoffer
# https://github.com/nickpoison/tsa4/blob/master/chap6.md
library("astsa")

set.seed(5)

# ground truth
rho_true <- 0.8
sigma_eps_true  <- .5
sigma_zeta_true <- .25

# create simmulated data
N <- 100
nu_t <- sarima.sim(ar=rho_true, n=N+1,sd=sigma_zeta_true)
y_t <- ts(nu_t[-1] + rnorm(N,0,sigma_eps_true))    

# Function to Calculate the Likelihood, to be supplied to a solver
Linn = function(para){
  rho <- exp(para[1]) # exp to make sure solver does not explore negative vales
  sigma_eps <- exp(para[2])
  sigma_zeta <- exp(para[3])

  kf <- Kfilter0(N,y_t,1,mu0=0,1e-4,rho,sigma_zeta,sigma_eps)
  return(kf$like)
}

# Estimation via BFGS solver on the likelihood
init.par <- log(c(.5,.5,.5)) # initial values of parameters
(est <- optim(init.par, Linn, NULL, method='BFGS',
  hessian=TRUE,control=list(trace=1,REPORT=1)))

rho <- exp(est$par[1])
sigma_eps <- exp(est$par[2])
sigma_zeta <- exp(est$par[3])

print("Estimated rho, sigma_eps, sigma_zeta")
print(c(rho,sigma_eps,sigma_zeta))

