This is a question that's been bugging me for some time. The problem is this: I'm modelling the residuals of a model $f(t,\vec{\theta})$ with (what I think is) an AR process plus a white noise process via MCMC using a multivariate gaussian likelihood, where I model the covariance matrix elements equal to the elements of the autocovariance of an $AR(1)$ process plus i.i.d. white noise $\varepsilon(t)$, i.e., I model

$$r(t)=AR(1)+\varepsilon(t),$$
where $r(t)$ are the residuals of my model, $r(t)=d(t)-f(t,\vec{\theta})$, where $d(t)$ is the data. The thing is that I wanted a visual check to see if my fit was ok (because I'm pretty sure via some previous analysis that my residuals can be modelled efficiently by an $AR(1)$ plus white noise: I just want to show this result to other people "quickly").

What I did then was to analize the residuals via the `arima` and `arma` functions in the `tseries` library in R, and I saw that they actually fitted the residuals with the same coefficient as the MCMC result (within the error bounds and not taking in consideration the additive white noise, which has small variance anyways). However, I also saw that the `arma` function has some nice plots: it plots the "residuals" of the AR fit.

This is my question: **how do you plot the residuals of an AR fit if an AR model is stochastic by nature?** I think my question goes down to: **how do you generate a realization of your model that matches your particular realization of the AR process (i.e. your data) in order to substract this realization from your data?** The only way I could think of is to simulate various realizations of the $AR(1)$ model until one fits your data, but this seems like stacking the deck.