# Does this Monte Carlo Technique Have a Name?

I sketched this algorithm out the other night. I am sure it has a name, I just do not know what it is yet. It would be helpful if someone could point me in the right direction for research. I provide motivation and then define the technique.

## Motivation

I am in a situation where I have to recursively estimate Bayesian time series models. That is, I have to estimate a posterior for my data over a time periods $[1,2,...,t]$, then $[1,2,...,t,t+1]$, then $[1,2,...,t,t+1,t+2]$, going all the way toward the end of my sample.

I am using STAN to estimate posterior samples. Recursively estimating the posteriors through STAN is excessively burdensome. The motivation for the below MC technique is to speed up computation time by taking weighted re-samples from the posterior sample estimated over $[1,...,t]$ to approximate the posteriors for $[1,..,t+1]$,$[1,..,t+2]$, up to some time period $[1,..,t+h]$ where the weights are functions of the predictive likelihoods generated at time $t$.

The general logic is that the posterior for $[1,...,t]$ and $[1,...,t+h]$ should be "similar" if $t$ is large and $h$ is relativity small.

## The MC Algorithm

Consider two conditional pdfs, $\pi(\theta|M_1)$ and $\pi(\theta|M_2)$. Assuming we have a sample $\theta_1, \theta_2,...,\theta_G$ where each $\theta_g$ was generated by one of the above distributions. The posterior odds ratio can be calculated as $$\frac{\pi(\theta|M_1)}{\pi(\theta|M_2)} \times r$$ where $r=\frac{\pi(M_1)}{\pi(M_2)}$ is the prior odds ratio. Using this we know $$Pr(\theta_g \sim \pi(\theta|M_1)) = w_g = 1-\frac{1}{\frac{\pi(\theta|M_1)}{\pi(\theta|M_2)}r +1}$$ We can then take a weighted sample of $\theta_1,...,\theta_G$, re-sampling each $\theta_g$ with probability $w_g$ to approximate $\pi(\theta|M_1)$.

Now consider the case where we new before hand that all the $\theta_g$ where sampled from $\pi(\theta|M_2)$. Then $w_g=r=0$. However, this also means that as $r \rightarrow 0$ the pdf of the weighted sample approaches $\pi(\theta|M_1)$. Thus, one could approximate $P(\theta|M_1)$ using a reasonably small $r$ (I have ran simulations for simple examples that seem to support this claim).

## Application

The above example is trivial, the following application better exemplifies my motives for using this technique. Consider the time series data $y_1,...,y_T$. Let $Y_{n:m}=y_n,y_{n+1},...,y_{m}$, for $m>n$. Denote the non-normalized posterior distribution as $P(\theta|Y_{n:m})$, the likelihood as $f(Y_{n:m}|\theta)$, and the normalization constant (or marginal likelihood) as $m(Y_{n:m})$.

In this example I have a sample $\theta_1,...,\theta_G$ which is drawn exclusively from $P(\theta|Y_{1:t})$ and which I would like to re-sample to approximate $P(\theta|Y_{1:(t+h)})$. Following the procedure in the above section and using Corollary 1 and 2 given below I calculate the posterior odds ratio :

$$\frac{P(\theta|Y_{1:(t+h)})/m(Y_{1:(t+h)})}{P(\theta|Y_{1:t})/m(Y_{1:t})} \times r=\frac{P(\theta|Y_{1:t})f(Y_{(t+1):(t+h)}|\theta)/m(Y_{1:(t+h)})}{P(\theta|Y_{1:t})/m(Y_{1:t})} \times r =$$ $$\frac{f(Y_{(t+1):(t+h)}|\theta)}{E_{1:t}[ f(Y_{(t+1):(t+h)}|\theta)]} \times r$$

Which makes the re-sampling weights $$w_g = 1-\frac{1}{\frac{f(Y_{(t+1):(t+h)}|\theta_g)}{E_{1:t}[ f(Y_{(t+1):(t+h)}|\theta)]}r + 1}$$

I am hoping that with a small enough choice of $r$, I can estimate the majority of posteriors in this fashion rather than re-estimating all the posteriors through STAN.

Corollary 1: $P(\theta|Y_{1:(t+h)}) = P(\theta|Y_{1:t})f(Y_{(t+1):(t+h)}|\theta)$.

Corollary 2: It can be shown that $\frac{m(Y_{(t+1):(t+h)})}{m(Y_{1:t})}=E_{1:t}[ f(Y_{(t+1):(t+h)}|\theta)]$, which is a joint predictive likelihood.

• Drawing Bayesian inference for a sequential problem calls for a sequential method like particle filters, SMC, and variants. Running a full MCMC run at each time t is sort of hopeless... Jan 23 '16 at 12:34
• @xian thank you, I do not think this is a sequential problem. The model is being backtested. I track real-time performance by collecting out-of sample predictive likelihoods for forecasts that would have been estimated at each time $t$. This is different then estimating dynamic coeficients over the full sample which is where I have learned about filtering problems. Jan 23 '16 at 12:51
• When I read "estimate a posterior for my data over a time periods [1,2,...,t], then [1,2,...,t,t+1][1,2,...,t,t+1]", I fail to see how the issue is not sequential: you consider sequentially a sequence of posteriors... Jan 23 '16 at 12:54
• @Xian Sorry, my fault. I think I am still a little lost in the vocabulary used to describe SMC. I will have a look at the Wikipedia entry. Jan 23 '16 at 13:08

This is a sequential Monte Carlo problem where one considers a sequence of posteriors $\pi_1,\pi_2,\ldots,\pi_T$ that are "reasonably" close from one another for the simulations from $\pi_t$ to be interesting proposals for $\pi_{t+1}$.

The solution you propose is not valid as such because of the reliance on this proportion $r$ that is actually equal to zero when running the sequential Monte Carlo scheme. (Or, in other words, there is no prior ratio for the marginal distributions of $y_{1:t}$ and of $y_{1:t+1}$.) A very similar resolution is sequential importance sampling where the current sample from $\pi_t$ $\theta_1^t,...,\theta_G^t$ is reweighted by $$\{\pi_{t+1}/\pi_t\}(\theta_g^t)$$ to become a sample from $\pi_t$. In the case of a state-space model where $$\pi_{t+1}(\theta)\propto\pi_t(\theta) f(y_{t+1}|\theta,y_{1:t})$$ the importance weight reduces to $$\{\pi_{t+1}/\pi_t\}(\theta)\propto f(y_{t+1}|\theta,y_{1:t})$$

While valid per se (in the importance sampling sense that the expectation of the weighted average produces the correct expectation under $\pi_{t+1}$), this naïve type of sequential importance sampling is dominated by more advanced particle filters where the sample itself is modified before being reweighted. This avoids in particular the degeneracy phenomenon, namely that the weights$$\prod_{j=1}^h f(y_{t+j}|\theta,y_{1:(t+j-1)})$$degenerate to zero almost surely.

I suggest you read the Wikipedia entry on particle filters, as it is quite detailed and definitely rigorous, given that one of the leaders in the field, Pierre Del Moral, has contributed to its improvement.

• +1. And as the wikipedia entry says in the first sentence, but I'll repeat here for future searches -- does the search look at comments? -- this technique is also known as a "Particle Filter". There is a: Stan State Space github project that I just found in a search, that might or might not be helpful: github.com/sinhrks/stan-statespace . Jan 23 '16 at 14:12
• Thank you for your answer. The one part I still do not understand is that $\pi_{t+1}(\theta)\propto\pi_t(\theta) f(y_{t+1}|\theta,y_{1:t})$ implies $\pi_{t+1}(\theta) = \frac{\pi_t(\theta) f(y_{t+1}|\theta,y_{1:t})}{M_{t+1}}$, $M_{t+1}>0$ which means $\{\pi_{t+1}/\pi_t\}(\theta)=\frac{f(y_{t+1}|\theta,y_{1:t})}{M_{t+1}}$. Jan 24 '16 at 21:21
• Wouldn't the importance weight reduce to $\frac{f(y_{t+1}|\theta,y_{1:t})}{M_{t+1}}$ as opposed to $f(y_{t+1}|\theta,y_{1:t})$? Jan 24 '16 at 21:31
• Yes, on principle, the ratio would involve the normalising constant. However, in most realistic settings, $M_t$ is not available and is thus estimated by a second importance sampling step. This is sometimes called the self-normalised importance sampler. Jan 25 '16 at 5:31
• In this specific case would't $M_t$ be the predictive likelihood $f(y_{t+1}|y_{1:t})$ which could be approximated very well by taking the sample average of $f(f(y_{t+1}|y_{1:t},\theta_t)$? Jan 25 '16 at 7:33