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?


3 Answers 3


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.

  • $\begingroup$ No. $\nu$ is not observed. $\endgroup$ Commented Sep 25, 2022 at 11:34
  • $\begingroup$ Shouldn't the AR model be $y_t=\rho y_{t-1}+\epsilon_t$? $\endgroup$ Commented Sep 25, 2022 at 11:43
  • $\begingroup$ To be honest, I checked this with my professor, and he said that we should consider $y_{t-1}$. I think $\rho$ can be identified, as shown in my question. $\endgroup$ Commented Sep 25, 2022 at 11:57
  • $\begingroup$ Then you mean $\rho$ and $\sigma_v$ and $\sigma_\epsilon$ cannot be identified? Then what went wrong in my analysis? $\endgroup$ Commented Sep 25, 2022 at 12:05
  • $\begingroup$ I think we can identify $\rho$, as in $E\left(\frac{y_t}{y_{t-1}}\right)$. $\endgroup$ Commented Sep 25, 2022 at 12:51

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.

  • $\begingroup$ It's a bit vague. How would you define the E and M steps? $\endgroup$
    – Tim
    Commented Sep 25, 2022 at 17:35

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


# 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)

# 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',

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

print("Estimated rho, sigma_eps, sigma_zeta")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.