13
$\begingroup$

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.
$\endgroup$
  • $\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 Aug 23 '13 at 18:07
  • $\begingroup$ Thanks for the input, @Metrics. I'll check out the book! $\endgroup$ – Néstor Aug 24 '13 at 18:33
1
$\begingroup$
  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$.

$\endgroup$
  • $\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 Aug 22 '13 at 3:43
  • $\begingroup$ Hi Nestor, I've edited the answer to respond to your comments. Hope that helps. $\endgroup$ – Jamie Hall Aug 23 '13 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 Aug 23 '13 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 Aug 25 '13 at 20:56
-1
$\begingroup$

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

$\endgroup$
  • 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 Jan 25 '14 at 22:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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