How to approximate (log-)likelihood from model specification using particle filters In Calvet et al., "Robust Filtering" (JASA, 2015), the authors obtain (pseudo?) maximum likelihood estimates of the parameters of a time series model of the trading volume of a futures contract. The model is specified by
$$V_t = e^{x_t+\epsilon_t}, \epsilon \sim \mathcal{N}(-\sigma^2_y/2,\sigma^2_y);$$
$$x_t = a + bx_{t-1} + \sigma_x|x_{t-1}|w_t, w \sim \mathcal{N}(0,1);$$
$$y_t = \log(V_t).$$
Here $V$ is trading volume, $x$ is the system state, and $y$ is the observation. The parameters to be estimated are $a$, $b$, $\sigma_x$, and $\sigma_y$.
I am familiar with basic closed-form maximum likelihood estimation, but the likelihood is in this case not available in closed form. How can I go about approximating the likelihood to numerically estimate the parameters?
The authors indicate that they employ the downhill simplex algorithm for the actual estimation (once the likelihood approximation has been computed), which I believe I could figure out how to use, but I struggle with the step of simulating a numerical approximation to the likelihood in the first place. Neither the paper nor the online appendix gives any more details about the simulation process itself. The paper is about particle filters, and the authors obtain different likelihoods for different particle filters. Below are their simulated likelihoods. My understanding is that these likelihoods are the result of simulations of the particle filter (state trajectories for each parameter value). But this simulation requires some parameter values to be run in the first place, so it seems circular. How these initial parameters are obtained is central to my question. Do they just start with random parameters and then compute the likelihoods in order to determine which parameters give a bitter fit? Wouldn't this introduce a dependence of the estimates on the initial conditions?
Any hints as to how to obtain these likelihood approximations and the role played in the process by the particle filter would be much appreciated.

 A: You can think of these different particle filters as different pieces of measurement equipment (like scales or rulers).  When we're outside of pure geometry, it's difficult to know exactly how big an object is, and we might get slightly different answers if we measure it a few times with different equipment.  Similarly, when we don't have the closed form of the likelihood, we might get different answers depending on exactly how we measure.
When the full likelihood is difficult to compute exactly, particle filters can give an approximate likelihood by walking through the steps of your simulated process and estimating the likelihood associated with each step. Then you can get the full likelihood by multiplying the steps together (or adding up the log-likelihoods), as they suggest just after "step 3" of their description in section 3.5.2.
These are Monte Carlo estimates, so there will be some sampling error (e.g. the jaggedness in the red lines of your figure), but they should asymptotically converge to the true value as the number of particles increases.
The authors were concerned about how well this algorithm will perform when the number of samples is small/finite, so they give it a tweak that makes it less noisy.  This tweak changes the estimates of the likelihood (hopefully by bringing them closer to the true value). If we think of the particle filter as a scale that weighs objects, their robust version is like a scale that tries to minimize random noise (e.g. from air currents occasionally pushing down on the scale), but which might introduce other artifacts.
A: Differing from your model, I can give you some ideas based on my experience. Let's say you have a state space model:
$y_t = ax_t + \alpha_t$
$x_{t+1} = bx_t + e_t, \quad e_t \sim N(0, 1)$
A regular assumption on a state space model is that $\alpha_t$ is also a random variable from Gaussian distribution, e.g., $N(0, 2)$. We can use the maximum likelihood method to estimate parameters $a$ and $b$ under a Gaussian distributed on $\alpha_t$. However, if we consider $\alpha_t$ is a random number from an $\alpha$ stable distribution, then the pdf is not analytically expressible. 
A possible way to handle this kind of problem is to use Approximate Bayesian Computation (ABC) method.  Please read the following link to learn this method if you like.
https://en.wikipedia.org/wiki/Approximate_Bayesian_computation
