1. The problem

I have some measurements of a variable $y_t$, where $t=1,2,..,n$, for which I have a distribution $f_{y_t}(y_t)$ obtained via MCMC, which for simplicity I'll assume is a gaussian of mean $\mu_t$ and variance $\sigma_t^2$.

I have a physical model for those observations, say $g(t)$, but the residuals $r_t = \mu_t-g(t)$ appear to be correlated; in particular, I have physical reasons to think that an $AR(1)$ process will suffice to take into account the correlation, and I plan on obtaining the coefficients of the fit via MCMC, for which I need the likelihood. I think the solution is rather simple, but I'm not pretty sure (it seems so simple, that I think I'm missing something).

2. Deriving the likelihood

A zero-mean $AR(1)$ process can be written as: $$X_t = \phi X_{t-1}+\varepsilon_t,\ \ \ (1)$$ where I'll assume $\varepsilon_t\sim N(0,\sigma_w^2)$. The parameters to be estimated are, therefore, $\theta = \{\phi,\sigma_w^2\}$ (in my case, I also have to add the parameters of the model $g(t)$, but that's not the problem). What I observe, however, is the variable $$R_t = X_t+\eta_t,\ \ \ (2)$$ where I'm assuming $\eta_t\sim N(0,\sigma_t^2)$, and the $\sigma_t^2$ are known (the measurement errors). Because $X_t$ is a gaussian process, $R_t$ is also. In particular, I know that $$X_1 \sim N(0,\sigma_w^2/[1-\phi^2]),$$ therefore, $$R_1 \sim N(0,\sigma_w^2/[1-\phi^2]+\sigma_t^2).$$ The next challenge is to obtain $R_t|R_{t-1}$ for $t\neq 1$. To derive the distribution of this random variable, note that, using eq. $(2)$ I can write $$X_{t-1} = R_{t-1}-\eta_{t-1}.\ \ \ (3)$$ Using eq. $(2)$, and using the definition of eq. $(1)$, I can write, $$R_{t} = X_t+\eta_t = \phi X_{t-1}+\varepsilon_{t}+\eta_t.$$ Using eq. $(3)$ in this last expression, then, I obtain, $$R_{t} = \phi (R_{t-1}-\eta_{t-1})+\varepsilon_{t}+\eta_t,$$ thus, $$R_t|R_{t-1} = \phi (r_{t-1}-\eta_{t-1})+\varepsilon_{t}+\eta_t,$$ and, therefore, $$R_t|R_{t-1} \sim N(\phi r_{t-1},\sigma_w^2+\sigma_t^2-\phi^2\sigma^2_{t-1}).$$ Finally, I can write the likelihood function as $$L(\theta) = f_{R_1}(R_1=r_1) \prod_{t=2}^{n} f_{R_{t}|R_{t-1}}(R_t=r_t|R_{t-1}=r_{t-1}),$$ where the $f(\cdot)$ are the distributions of the variables that I just defined, .i.e., defining $\sigma'^2 = \sigma_w^2/[1-\phi^2]+\sigma_t^2,$ $$f_{R_1}(R_1=r_1) = \frac{1}{\sqrt{2\pi \sigma'^2}}\text{exp}\left(-\frac{r_1^2}{2\sigma'^2}\right),$$ and defining $\sigma^2(t) = \sigma_w^2+\sigma_t^2-\phi^2\sigma^2_{t-1}$, $$f_{R_{t}|R_{t-1}}(R_t=r_t|R_{t-1}=r_{t-1})=\frac{1}{\sqrt{2\pi \sigma^2(t)}}\text{exp}\left(-\frac{(r_t-\phi r_{t-1})^2}{2\sigma^2(t)}\right)$$

3. Questions

  1. Is my derivation ok? I don't have any resources to compare other than simulations (which seem to agree), and I'm not a statistician!
  2. Are there any derivation of this kind of things in the literature for $MA(1)$ proccesses or $ARMA(1,1)$ proccesses? A study for $ARMA(p,q)$ proccesses in general that could be particularized to this case would be nice.
  • $\begingroup$ I exactly don't have a solution for you. But, I think this is a kind of error-in variables problem. I have seen this stuff in Macroeconomic Theory by Thomas Sergent (1980's book). You may want to look at that one. $\endgroup$
    – Metrics
    Commented Aug 23, 2013 at 18:07
  • $\begingroup$ Thanks for the input, @Metrics. I'll check out the book! $\endgroup$
    – Néstor
    Commented Aug 24, 2013 at 18:33

2 Answers 2

  1. You're on the right track, but you've made a mistake in deriving the distribution of $R_t$ given $R_{t-1}$: the conditional mean isn't $\phi r_{t-1}$. It's $\phi \widehat{x}_{t-1}$, where $\widehat{x}_{t-1}$ is your best estimate of $X$ from the previous period. The value of $\widehat{x}_{t-1}$ includes information from previous observations as well as $r_{t-1}$. (To see this, consider a situation where $\sigma_w$ and $\phi$ are negligible, so you're effectively estimating a fixed mean. After lots of observations, your uncertainty about $X$ will be much smaller than $\sigma_{\eta}$.) This can be confusing at first, because you observe $R$ and not $X$. That just means that you're dealing with a state-space model.

  2. Yes, there is a very general framework for using linear-Gaussian models with noisy observations, called the Kalman filter. This applies to anything with an ARIMA structure and many more models too. Time-varying $\sigma_{\eta}$ is OK for the Kalman filter, provided it's not stochastic. Models with, e.g., stochastic volatility need more general methods. To see how the Kalman filter is derived, try Durbin-Koopman or chapter 3 of Harvey. In Harvey's notation, your model has $Z=1$, $d=c=0$, $H_t = \sigma_{\eta,t}^2$, $T=\phi$, $R=1$ and $Q=\sigma^2_w$.

  • $\begingroup$ Hi Jamie, thanks for your input. A couple of comments: 1. I'm not sure about that. It was, actually, my first attempt as a solution but both my intuition and simulations do not agree with that. The thing is that I actually don't observe $X_t$, I observe $R_{t}$; plus, can you prove (arithmetically) that the conditional mean of the random variable $R_{t}|R_{t-1}=r_{t-1}$ (note that it is not $R_{t}|X_{t-1}=x_{t-1}$) is actually $\phi \hat{x}_{t-1}$? 2. Can you elaborate on the application of the Kalman filter to this particular problem? $\endgroup$
    – Néstor
    Commented Aug 22, 2013 at 3:43
  • $\begingroup$ Hi Nestor, I've edited the answer to respond to your comments. Hope that helps. $\endgroup$
    – Jamie Hall
    Commented Aug 23, 2013 at 4:29
  • $\begingroup$ Hi Jamie: about the second point, that's ok, thanks :-)! However, I still can't see your first point. Can you point me to a formal derivation? In particular, I would like to know what part of my reasoning is wrong (and why)! $\endgroup$
    – Néstor
    Commented Aug 23, 2013 at 4:50
  • $\begingroup$ You skipped a step: the distribution of $X_1$ given $R_1$. It's $N(\frac{\sigma^2_{x,1}}{(\sigma^2_{x,1}+\sigma^2_{\eta,1})} r_1, \sigma^2_{x,2})$, where $\sigma_{x,1}^2$ is the variance you calculated in the first step, and $\sigma_{x,2}^2$ is twice the harmonic mean of $\sigma_{x,1}^2$ and $\sigma_{\eta,1}^2$. (This is just like Bayesian updating with two Gaussian pdfs.) Your equation (3) is formally correct, but you're throwing away information by using that instead of $p(X_{t-1} | R_{1:t-1})$. $\endgroup$
    – Jamie Hall
    Commented Aug 25, 2013 at 20:56

Honestly, you should code this in BUGs or STAN and not worry about it from there. Unless this is a theoretical question.

  • 2
    $\begingroup$ (-1) To this response; this is clearly a theoretical question ;-). Consider improving why you think I should code it in BUGs or STAN and what it does have to do with the original question? $\endgroup$
    – Néstor
    Commented Jan 25, 2014 at 22:01

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.