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?
 A: 
I think $\rho$ can be identified as follows:
We have $v_{it}=y_{it}-\epsilon_{it}$ from the first equation [...]

No, you can't. For that, you would need to know $\epsilon_{it}$, while you don't as it is not observed. The rest of the equations also doesn't work as they are based on variables that are not known. To identify $\nu_{it}$ you need to know $\nu_{it-1}$ and $\rho$, both unknown. All those variables are unknown, so there is an infinite number of solutions that "match" them. For it to be identifiable you would need additional sources of information to identify at least some of the parameters.
A: 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.
A: 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))

