\section{Introduction}

We consider the problem of estimating parameters $\boldsymbol{\theta}$ from models characterized by missing (incomplete) data, in the sense discussed in \cite{Dempster1977} when introducing the celebrated EM algorithm. Observable quantities are denoted with $\mathbf{Y}$, while $\mathbf{X}$ denotes latent (unobservable) components. However, since our work deals with parameter inference, we explicitly distinguish between $\boldsymbol{\theta}$ and $\mathbf{X}$, i.e. even though $\boldsymbol{\theta}$ is unknown we do not consider it as embedded in $\mathbf{X}$.
While here we deal with a modification of an EM-type algorithm, for the moment our interest is to discuss the inference problem for models having so-called intractable likelihoods. These are models for which the likelihood function is unavailable in closed form and for which computation of an approximation (or the evaluation of said approximation) is computationally prohibitive. Two of the discussed examples are state-space models (SSM), and for SSM recent advancements in sequential Monte Carlo methods (also known as particle filters) have revolutionised the practical application of statistical inference, especially the Bayesian kind, see \cite{andrieu2009pseudo}, \cite{andrieu2010particle} and the review in \cite{kantas2015particle}. For more general models than SSM, approximate Bayesian computation (ABC) is often the only available solution to perform statistical inference for the parameters of complex models with intractable likelihoods. ABC (see \citealp{sisson-fan(2011)} and \citealp{marin-et-al(2011)} for reviews) is an ensemble of algorithms that only requires the ability to generate synthetic observations from the assumed data generating model: therefore, the incredible rising of ABC methods in the past 15 years is due to its generality and plug-and-play features. The essence of ABC methods is briefly discussed in sections \ref{sec:synlik} and \ref{sec:abc-mcmc}. While ABC algorithms have been developed since the '90s (\citealp{tavare1997inferring}, \citealp{pritchard1999population}), the most important issues for a successful implementation of ABC are still as relevant today as they were twenty years ago. In particular, the most typical usage of ABC requires the analyst to specify summary statistics that are ``informative'' regarding the unknown $\boldsymbol{\theta}$. Moreover, a threshold parameter is introduced to compare summary statistics computed on the available data with summaries computed on simulations from the assumed data generating model. The problem of selecting appropriate summaries is the most serious of the two, and developments of methods for their automatized construction is crucial (see \citealp{fearnhead-prangle(2011)}). The determination of the threshold for summaries comparison is also very important and has a significant impact on the computational budget (the smaller the threshold, the better the inference accuracy and the higher the number of rejected simulations). Setting an appropriate weigthing scheme for comparing summaries having different variabilities can be another crucial factor, see \cite{prangle-adaptiveABC} and our application in section \ref{sec:g-and-k}. Finally, when ABC is implemented within an MCMC sampler, there is a further layer of practical issues that are usually of difficult management for the non-expert practitioner, such as coding an appropriate adaptive MCMC method for the generation of parameter proposals, also noting that the frequency of the adaptation affects results.
It is fair to say that calibration of ABC algorithms is often not trivial.

A more recent alternative to ABC is given by the synthetic likelihoods (SL) methodology introduced in \cite{wood2010statistical} and illustrated in section \ref{sec:synlik}. Same as ABC, synthetic likelihoods is a plug-and-play method for models with an intractable likelihood. SL requires the specification of data summaries, but no threshold parameter is introduced and the weighting of the summaries is implictly handled via simulation, thus the method is very easy to implement. However, while ABC sets no assumptions on said summaries, SL assumes a multivariate Gaussian distribution: hence, SL is less general than ABC and as discussed in \cite{price2016bayesian} significant departures from the assumed Gaussianity can have a negative impact on inference results.
In general, interest for likelihood-free methods is high, as the appeal of general methods only requiring forward simulations from a generative model allows the researchers to move beyond tractable (and often oversimplistic) mathematical representations of reality, to explore challenging models relevant for real world applications. Moreover, as shown in \cite{fasiolo2016comparison}, state-of-art (exact) methods based on sequential Monte Carlo can fail when the model is near-chaotic and when the algorithms are initialised at bad starting parameter values. In some cases where exact methods for SSM fail, approximate approaches can succeed, see \cite{fasiolo2016comparison} and our illustrative codes at \url{https://github.com/umbertopicchini/pomp-ricker}.

Given the above arguments, in the present work we consider the idea underlying the synthetic likelihood approach, and embed this into a stochastic EM algorithm for maximum likelihood estimation, the stochastic approximation (SAEM) algorithm of \cite{Delyon1999}, summarised in section \ref{sec:simulations}. The resulting SAEM-SL algorithm is a likelihood-free version of SAEM which is easy to code, requires minimal tuning and appeals a general class of incomplete-data models, either ``static'' (time-independent) and dynamic models. We compare our method with several other approaches, including iterated filtering, particle marginal methods, approximate Bayesian computation and SAEM incorporating a sequential Monte Carlo step. The comparison of our approach with alternatives is performed to highlight advantages and deficiencies of each approach. While SAEM-SL is primarily a likelihood-free version of SAEM, we also show that for those models where alternative methods are possible, these might produce suboptimal results. For example, in section \ref{sec:theophylline-results} we show that while SAEM-SL performs well, instead a pseudo-marginal Bayesian synthetic likelihood approach (via MCMC) does not allow the identification of two specific parameters. Although the latter is likely due to the specific priors we selected (uniforms, as in the methodology presented in \cite{wood2010statistical}), this is a way to highlight difficulties that cannot sometimes be handled without modifying the statistical model.

In the rest of this section, we illustrate the reasons for pursuing this approach. Since two of our simulation studies use state-space models, our notation introduces quantities that are time-indexed, however we emphasize that the methodology is suited for dynamic models that are not SSM and also for static models.

State-space models (SSM), sometimes denoted hidden Markov models (\citealp{cappe2005inference}), are used in many fields, such as biology, chemistry, ecology, signal processing etc.
We now introduce some notation. Consider a stochastic process $\{\mathbf{Y}_t\}_{t> t_0}$, $\mathbf{Y}_t\in\mathsf{Y}\subseteq \mathbb{R}^{d_y}$, which is observed at discrete sampling times $t\in\{t_1,...,t_n\}$ with $t_1>t_0\geq 0$, and we denote with $\mathbf{Y}_{1:n}=(\mathbf{Y}_1,...,\mathbf{Y}_n)$ the corresponding observations (data) from $\{\mathbf{Y}_t\}_{t> t_0}$ collected at said time points, where $\mathbf{Y}_{t_j}\equiv \mathbf{Y}_j$ for $j=1,...,n$.
Consider also a latent (unobservable) continuous-time stochastic process $\{\mathbf{X}_t\}_{t\geq t_0}$, $\mathbf{X}_t\in\mathsf{X}\subseteq \mathbb{R}^{d_x}$.
Process $\mathbf{X}_t\sim p(\mathbf{x}_t|\mathbf{x}_s,\boldsymbol{\theta}_x)$ is assumed Markov with transition densities $p(\mathbf{x}_t|\mathbf{x}_s,\cdot)$, $t_0\leq s<t$.
Denote with $\mathbf{X}_{1:n}=(\mathbf{X}_1,...,\mathbf{X}_n)$ the unobserved values for $\{\mathbf{X}_t\}_{t\geq t_0}$ at times $\{t_1,...,t_n\}$ and set $\mathbf{X}_{0:n}=(\mathbf{X}_0,\mathbf{X}_{1:n})$, where $\mathbf{X}_0$ is the (random or fixed) initial state for $\{\mathbf{X}_t\}$ at time $t_0$.
Both processes $\{\mathbf{X}_t\}$ and $\{\mathbf{Y}_t\}$ depend on their own (assumed unknown) vector-parameters $\boldsymbol{\theta}_x$ and $\boldsymbol{\theta}_y$ respectively.
We think at $\{\mathbf{Y}_t\}$ as a measurement-error-corrupted version of $\{\mathbf{X}_t\}$ and assume that observations for $\{\mathbf{Y}_t\}$ are conditionally independent given $\{\mathbf{X}_t\}$.
The SSM can be summarised as
\begin{equation}
\begin{cases}
\mathbf{Y}_j\sim p(\mathbf{y}_j|\mathbf{X}_j,\boldsymbol{\theta}_y),\quad j=1,...,n\\
\mathbf{X}_t\sim p(\mathbf{x}_t|\mathbf{x}_s,\boldsymbol{\theta}_x),\quad t_0\leq s<t.
\end{cases}
\label{eq:state-space-general}
\end{equation}
Typically $p(\mathbf{y}_j|\mathbf{X}_j,\cdot)$ is a known density (or probability mass) function, though this is not always the case (see for example $\alpha$-stable stochastic volatility, \citealp{calvet2014accurate}). Regarding the transition density $p(\mathbf{x}_t|\mathbf{x}_s,\cdot)$, this is typically unknown except for very simple toy models.

Goal of our work is to estimate the parameters $(\boldsymbol{\theta}_x,\boldsymbol{\theta}_y)$ using data $\mathbf{Y}_{1:n}=(\mathbf{Y}_1,...,\mathbf{Y}_n)$. For ease of notation we refer to the vector $\boldsymbol{\theta}:=(\boldsymbol{\theta}_x,\boldsymbol{\theta}_y)$ as the object of our inference. As previously remarked, the SAEM-SL methodology we introduce does not require data generated from a SSM, hence conditional independence of observations and Markovianity of $\{\mathbf{X}_t\}$ are not necessary for SAEM-SL to work.

Bayesian estimation for state-space models has been widely developed in the literature. There the goal is to derive analytically the posterior distribution $\pi(\boldsymbol{\theta}|\mathbf{Y}_{1:n})$ or, most frequently, implement an algorithm for sampling draws from the posterior. Sampling procedures are often carried out using Markov chain Monte Carlo (MCMC) or Sequential Monte Carlo (SMC) embedded in MCMC procedures (\citealp{andrieu2009pseudo}, \citealp{andrieu2010particle}).
In this work we instead aim at maximum likelihood estimation of $\boldsymbol{\theta}$. Several methods have been proposed in the literature, and some are based on the well-known EM algorithm \citep{Dempster1977}. The EM algorithm computes the conditional expectation of the complete-likelihood for the pair $(\mathbf{Y}_{1:n},\mathbf{X}_{0:n})$ and then produces a (local) maximizer for the data likelihood function based on observations $\mathbf{Y}_{1:n}$. One of the difficulties with EM is to compute the conditional expectation of the state $\{\mathbf{X}_t\}$ given the observations  $\mathbf{Y}_{1:n}$. This conditional expectation can be computed exactly with the Kalman filter when the state-space is linear and Gaussian \citep{cappe2005inference}, and otherwise it has to be approximated. In this work we focus on a stochastic approximation of the EM algorithm, namely the Stochastic Approximation EM (SAEM) \citep{Delyon1999}. The problem with implementing SAEM is at least two-fold: (i) it is necessary to generate  an appropriate ``proposal'' for the state $\{\mathbf{X}_t\}$, conditionally on the current value of $\boldsymbol{\theta}$. Sequential Monte Carlo (SMC) algorithms \citep{doucet2001sequential} have already been coupled to stochastic EM algorithms (see e.g. \cite{HuysPaninski2006, HuysPaninski2009, Lindsten2013, Ditlevsen2014} and references therein). However (ii) a second difficulty is that in order to use SAEM the complete likelihood of $\boldsymbol{\theta}$ based on the joint distribution of $(\mathbf{Y}_{1:n},\mathbf{X}_{0:n})$ must be tractable. With ``tractable'' we mean that the model at hand has a complete likelihood that it is possible to write in closed-form, and that additionally it is possible to derive analytically essential quantities, such as the corresponding sufficient statistics: this is because the convergence of SAEM to the maximizer of the data likelihood is ensured only for observations belonging to the exponential family. These requirements are usually  very difficult to satisfy, or result impossible for most realistic models. Even when these can be satisfied, the required analytic work is at best a tedious, difficult and error-prone task. This is particularly true when several candidate models are under screening. Also, such difficulties force the modeller to formulate oversimplified, tractable models so that SAEM can be implemented. However realistic models call for more complex formulations which are usually not amenable to closed form analytic computations.

\section{The complete likelihood and stochastic approximation EM}
\label{sec:complete-likelihood}

Recall that $\mathbf{Y}_{1:n}=(\mathbf{Y}_1,...,\mathbf{Y}_n)$ denotes the available data collected at times $(t_1,...,t_n)$ and denote with $\mathbf{X}_{1:n}=(\mathbf{X}_1,...,\mathbf{X}_n)$ the corresponding unobserved states. We additionally set $\mathbf{X}_{0:n}=(\mathbf{X}_0,\mathbf{X}_{1:n})$ for the vector including an initial (fixed or random) state $\mathbf{X}_0$, that is $\mathbf{X}_1$ is generated as $\mathbf{X}_1\sim p(\mathbf{x}_1|\mathbf{x}_0)$. When the transition densities between sampling times $p(\mathbf{x}_j|\mathbf{x}_{j-1})$ are available in closed form ($j=1,...,n$), the ``data likelihood'' function for $\boldsymbol{\theta}$ (sometimes denoted ``incomplete data likelihood'') can be written as
\begin{align}
p(\mathbf{Y}_{1:n};\boldsymbol{\theta}) &= \int p_{\mathbf{Y},\mathbf{X}}(\mathbf{Y}_{1:n},\mathbf{X}_{0:n}\, ; \boldsymbol{\theta})\,d\mathbf{X}_{0:n} = \int p_{\mathbf{Y}|\mathbf{X}}(\mathbf{Y}_{1:n}|\mathbf{X}_{0:n}\, ; \boldsymbol{\theta})p_{\mathbf{X}}(\mathbf{X}_{0:n}; \boldsymbol{\theta})\,d\mathbf{X}_{0:n}\nonumber\\
&=\int p(\mathbf{X}_0)\biggl\{\prod_{j=1}^n p(\mathbf{Y}_j|\mathbf{X}_j;\boldsymbol{\theta})p(\mathbf{X}_j|\mathbf{X}_{j-1};\boldsymbol{\theta})\biggr\}d\mathbf{X}_0\cdots d\mathbf{X}_n  \label{eq:likelihood}
\end{align}
where we have assumed a random initial state with density $p(\mathbf{X}_0)$.
Here $p_{\mathbf{Y},\mathbf{X}}$ is the ``complete data likelihood'', $p(\mathbf{Y}_j|\mathbf{X}_j)$ the conditional density of $\mathbf{Y}_j$ and $p_{\mathbf{X}}(\mathbf{X}_{0:n}; \boldsymbol{\theta})$ the joint density of $\mathbf{X}_{0:n}$. The last equality in \eqref{eq:likelihood} exploits the notion of conditional independence of observations given latent states and the Markovian property of $\{\mathbf{X}_t\}$.
In general the likelihood \eqref{eq:likelihood} is not explicitly known either because the integral is multidimensional or because  expressions for transition densities are typically not available.
In addition, when an exact simulator for the solution of the dynamical process associated with the Markov process $\{\mathbf{X}_t\}$ is unavailable, hence it is not possible to sample from $p(\mathbf{X}_j|\mathbf{X}_{j-1};\boldsymbol{\theta})$, numerical discretisation methods are required, see the example in section \ref{sec:theoph}. Without loss of generality, say that we have equispaced sampling times such that $t_j=t_{j-1}+\Delta$, with $\Delta>0$. Now introduce a discretisation for the interval $[t_1,t_n]$ given by $\{\tau_1,\tau_h,...,\tau_{Gh},...,\tau_{nGh}\}$ where $h=\Delta/G$ and $G\geq 1$. We take $\tau_1=t_1$, $\tau_{nGh}=t_n$ and therefore $\tau_{i}\in \{t_1,....,t_n\}$ for $i=1,Gh,2Gh,...,nGh$. We denote with $N$ the number of elements in the discretisation $\{\tau_1,\tau_h,...,\tau_{Gh},...,\tau_{nGh}\}$ and with $\mathbf{X}_{1:N}= (\mathbf{X}_{\tau_1}, \ldots, \mathbf{X}_{\tau_N})$ the corresponding values of $\{\mathbf{X}_t\}$ obtained when using a given numerical/approximated method of choice. Then the likelihood function becomes
\begin{align*}
p(\mathbf{Y}_{1:n};\boldsymbol{\theta}) &= \int p_{\mathbf{Y},\mathbf{X}}(\mathbf{Y}_{1:n},\mathbf{X}_{0:N}\, ; \boldsymbol{\theta})\,d\mathbf{X}_{0:N} = \int p_{\mathbf{Y}|\mathbf{X}}(\mathbf{Y}_{1:n}|\mathbf{X}_{0:N}\, ; \boldsymbol{\theta})p_{\mathbf{X}}(\mathbf{X}_{0:N}; \boldsymbol{\theta})\,d\mathbf{X}_{0:N}\\
&=\int \biggl\{\prod_{j=1}^n p(\mathbf{Y}_j|\mathbf{X}_j;\boldsymbol{\theta})\biggr\}p(\mathbf{X}_0)\prod_{i=1}^N p(\mathbf{X}_i|\mathbf{X}_{i-1};\boldsymbol{\theta})d\mathbf{X}_0\cdots d\mathbf{X}_N,
\end{align*}
where the product having index $j$ is over the $\mathbf{X}_{t_j}$'s and the product having index $i$ is over the $\mathbf{X}_{\tau_i}$'s.

\subsection{The standard SAEM algorithm}\label{sec:standard-SAEM}
The EM algorithm introduced by  \cite{Dempster1977} is a classical approach to estimate parameters of models with non-observed or incomplete data. Let us briefly cover the EM principle.
 The complete data of the model is  $(\mathbf{Y}_{1:n},\mathbf{X}_{0:N})$, where $\mathbf{X}_{0:N}\equiv \mathbf{X}_{0:n}$ if numerical discretisation is not required, and for ease of writing we denote this as $(\mathbf{Y},\mathbf{X})\equiv (\mathbf{Y}_{1:n},\mathbf{X}_{0:N})$ for the remaining of this section.
The EM algorithm maximizes the function
 $Q(\boldsymbol{\theta}|\boldsymbol{\theta}')=\mathbb{E}(L_c(\mathbf{Y},\mathbf{X};\boldsymbol{\theta})|\mathbf{Y};\boldsymbol{\theta}')$
in two steps, where $L_c(\mathbf{Y},\mathbf{X};\boldsymbol{\theta}):=\log p_{\mathbf{Y},\mathbf{X}}$ is the log-likelihood of the \textit{complete} data and $\mathbb{E}$ is the conditional expectation under the conditional distribution $p_{\mathbf{X}|\mathbf{Y}} (\cdot ; \boldsymbol{\theta}')$.
More explicitly, by denoting with $\hat{\boldsymbol{\theta}}^{(k-1)}$ the parameter estimate  obtained at iteration $k-1$ of EM, at $k$th iteration of EM the E-step  computes $Q(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}}^{(k-1)})=\int \log p_{\mathbf{Y},\mathbf{X}}(\mathbf{Y},\mathbf{X};\boldsymbol{\theta})p_{\mathbf{X}|\mathbf{Y}}(\mathbf{X}|\mathbf{Y};\hat{\boldsymbol{\theta}}^{(k-1)})d\mathbf{X}$. The M-step computes $\hat{\boldsymbol{\theta}}^{(k)}=\arg\max_{\boldsymbol{\theta}\in\Theta} Q(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}}^{(k-1)})$. The resulting sequence $\{\hat{\boldsymbol{\theta}}^{(k)}\}_k$ converges to a stationary point of the data likelihood $p(\mathbf{Y};\boldsymbol{\theta})$, under weak assumptions \citep{Dempster1977}. In most cases the E-step is difficult to perform, while the M-step can be considered relatively straightforward, meaning that standard optimization procedures for the M-step can be implemented, or closed form solutions are possible.

Important strategies for dealing with an intractable E-step are MCEM \citep{wei1990monte} and SAEM \citep{Delyon1999}, see also \cite{Lindsten2013} for a synthetic review.
In SAEM the integral in $Q(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}}^{(k-1)})$ is approximated using a stochastic procedure.
SAEM is proved to converge under general
conditions if $L_c(\mathbf{Y},\mathbf{X};\boldsymbol{\theta})$ belongs to the regular  exponential family
\begin{equation}
L_c(\mathbf{Y},\mathbf{X};\boldsymbol{\theta})= -\Lambda(\boldsymbol{\theta})+\langle \mathbf{S}_c(\mathbf{Y},\mathbf{X}),\Gamma(\boldsymbol{\theta})\rangle, \label{eq:complete-likelihood}
\end{equation}
where $\left\langle ...\right\rangle$ is the scalar product, $\Lambda$ and $\Gamma$ are two functions of $\boldsymbol{\theta}$ and  $\mathbf{S}_c(\mathbf{Y},\mathbf{X})$ is  the minimal sufficient statistic of the complete model.  The E-step is then divided into a
simulation step (S-step) of the missing data $\mathbf{X}^{(k)}$ under the
conditional distribution $p_{\mathbf{X}|\mathbf{Y}}(\cdot;\widehat{\boldsymbol{\theta}}^{(k-1)})$ and a stochastic
approximation step (SA-step) of the conditional expectation, using $(\gamma_k)_{k\geq 1}$ a sequence of real numbers in $[0,1]$, such that $\sum_{k=1}^\infty\gamma_k=\infty$ and $\sum_{k=1}^\infty\gamma_k^2<\infty$. This SA-step approximates  $\mathbb{E}\left\lbrack  \mathbf{S}_c(\mathbf{Y},\mathbf{X}) \vert  \widehat{\boldsymbol{\theta}}^{(k-1)} \right\rbrack$ at each iteration by the value $\mathbf{s}_k$ defined recursively as follows
$$
\mathbf{s}_{k}=\mathbf{s}_{k-1}+\gamma_k( \mathbf{S}_c(\mathbf{Y},\mathbf{X}^{(k)})-\mathbf{s}_{k-1}).
$$
The M-step is thus the update of the estimates $\widehat{\boldsymbol{\theta}}^{(k-1)}$
\begin{equation}
\widehat{\boldsymbol{\theta}}^{(k)}= \arg \max_{\boldsymbol{\theta} \in \Theta} \left(-\Lambda(\boldsymbol{\theta})+\langle \mathbf{s}_{k},\Gamma(\boldsymbol{\theta})\rangle \right).\label{eq:M-step}
\end{equation}
A schematic description of the SAEM procedure (coupled with a bootstrap filter) is in algorithm \ref{alg:saem-smc}, see also \cite{picchini2016coupling}. Also, notice that it is possible to implement a version of SAEM without the need to find the sufficient statistics $\mathbf{S}_c(\cdot)$, see \cite{Delyon1999} and \cite{Lindsten2013}, but then again convergence of the algorithm is only ensured when the complete loglikelihood is written as in \eqref{eq:complete-likelihood}. Moreover, when it is possible to parametrize the complete loglikelihood in terms of $\mathbf{S}_c(\cdot)$ as in \eqref{eq:complete-likelihood}, then it is sometimes possible to determine the $\widehat{\boldsymbol{\theta}}^{(k)}$ in \eqref{eq:M-step} explicitly, see sections \ref{sec:nonlingauss}--\ref{sec:theoph}, which has an obvious computational advantage.

Usually, the simulation step of the hidden trajectory $\mathbf{X}^{(k)}$ conditionally to the observations $\mathbf{Y}$ cannot be performed directly. A standard possibility is to use $M$ ``particles'' from sequential Monte Carlo filters, such as the bootstrap filter \citep{gordon1993novel}, see algorithm \ref{alg:bf}.
\begin{algorithm}
\small
\caption{SAEM with a bootstrap filter}\label{alg:saem-smc}
\begin{algorithmic}
\State Step 0. Set parameters starting values $\hat{\boldsymbol{\theta}}^{(0)}$, then set $M$, $\bar{M}$ and $k:=1$.
\State Step 1. For fixed $\hat{\boldsymbol{\theta}}^{(k-1)}$ apply the bootstrap filter in algorithm \ref{alg:bf} with $M$ particles and particles threshold $\bar{M}$.
\State 2 Sample an index $m'$ from the probability distribution $\{w_n^{(1)},...,w_n^{(m)}\}$ on $\{1,...,M\}$ and form the path $\mathbf{X}^{(k)}$ resulting from the genealogy of $m'$.
\State Step 3. {\bf Stochastic Approximation step :} update of the sufficient statistics
	$$\mathbf{s}_{k}  =  \mathbf{s}_{k-1} + \gamma_k\,\left(\mathbf{S}_c(\mathbf{Y}, \mathbf{X}^{(k)}) -\mathbf{s}_{k-1}\right)$$
\State Step 4. {\bf Maximisation step:} update $\boldsymbol{\theta}$
	$$\widehat{\boldsymbol{\theta}}^{(k)}= \arg \max_{\boldsymbol{\theta} \in \Theta} \left(-\Lambda(\boldsymbol{\theta})+\langle \mathbf{s}_{k},\Gamma(\boldsymbol{\theta})\rangle \right)$$
Set $k:=k+1$ and go to step 1.
\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\small
\caption{Bootstrap filter}\label{alg:bf}
\begin{algorithmic}
\State Step 0. Set $j=1$: for $m=1,...,M$ sample $\mathbf{X}_1^{(m)}\sim p(\mathbf{X}_0)$, compute weights $W_1^{(m)}=f(\mathbf{Y}_1|\mathbf{X}_1^{(m)})$ and normalize weights $w_1^{(m)}:=W_1^{(m)}/\sum_{m=1}^M W_1^{(m)}$.
\State Step 1.
   \If{$ESS(\{w_j^{(m)}\})<\bar{M}$}
\State resample $M$   particles $\{\mathbf{X}_j^{(m)},w_j^{(m)}\}$ and set $W_j^{(m)}=1/M$.
   \EndIf
\State Set $j:=j+1$ and if $j=n+1$, stop and return all constructed weights $\{W_j^{(m)}\}_{j=1:n}^{m=1:M}$ to sample a single path (see main text). Otherwise go to step 2.
\State Step 2. For $m=1,...,M$ sample $\mathbf{X}_j^{(m)}\sim p(\cdot|\mathbf{X}_{j-1}^{(m)})$. Compute
$$W_j^{(m)}:=w_{j-1}^{(m)}p(\mathbf{Y}_j|\mathbf{X}_j^{(m)})$$
normalize weights $w_j^{(m)}:=W_{j}^{(m)}/\sum_{m=1}^M W_j^{(m)}$ and go to step 1.
\end{algorithmic}
\end{algorithm}
The quantity ESS in algorithm \ref{alg:bf} is the effective sample size (e.g. \citealp{liu2008monte}) often estimated as $ESS(\{w_j^{(m)}\})=1/\sum_{m=1}^M(w_j^{(m)})^2$ and taking values between 1 and $M$, while $\bar{M}\leq M$ is a  threshold value that ``activates'' the resampling step, see \cite{Cappe2007} for an introduction to particle filters.
In addition to the procedure outlined in algorithm \ref{alg:bf}, once the set of normalised weights $\{w_n^{(1)},...,w_n^{(m)}\}$ is available at the end of the bootstrap filter, we sample a single index from the set $\{1,...,M\}$ having associated probabilities $\{w_n^{(1)},...,w_n^{(m)}\}$. Denote with $m'$ such index and with $a_j^m$ the ``ancestor'' of the generic $m$th particle sampled at time $t_{j+1}$, with $1\leq a_j^m \leq M$ ($m=1,...,M$, $j=1,...,n$). Then we have that particle $m'$ has ancestor $a_{n-1}^{m'}$ and in general particle $m''$ at time $t_{j+1}$ has ancestor $b_j^{m''}:=a_j^{b^{m''}_{j+1}}$, with $b_n^{m'}:=m'$. Hence, at the end of algorithm \ref{alg:bf} we can sample $m'$ and construct its genealogy (see also \citealp{andrieu2010particle}): the sequence of states $\{\mathbf{X}_t\}$ resulting from the genealogy of $m'$ is the chosen path that will be passed to SAEM in algorithm
\ref{alg:saem-smc}.

However, as explained in the Introduction and self-evident in the application in section \ref{sec:theoph}, constructing the SAEM machinery is a challenging task for most realistic models as typically the sufficient statistics $\mathbf{S}_c$ for the complete loglikelihood need to be available, for computational efficiency. Moreover, for state-space models it is necessary to know the expression of the transition densities, to construct the complete loglikelihood. For most stochastic nonlinear models, transition densities are typically unavailable in closed form.  In section \ref{sec:synlik} we propose a new, likelihood-free version of SAEM. But first, it is necessary to introduce the synthetic likelihoods methodology, due to \cite{wood2010statistical}.

\section{Synthetic likelihoods}\label{sec:synlik}

Approximate inference methods for complex models with an intractable likelihood function is an important topic that has generated enormous interest in the past fifteen years, especially thanks to developments in the approximate Bayesian computation (ABC) framework (see \cite{marin-et-al(2011)} for a review). The present work is not focussed on ABC methods, however we give a short account in section \ref{sec:g-and-k}, where we compare results from an  ABC-MCMC algorithm against our method. Briefly, in ABC parameters $\boldsymbol{\theta}^*$ are generated many times (either from their prior of from an MCMC kernel), and conditionally on these parameters artificial data $\mathbf{Y}^*$ are simulated from the data generating model. Then the corresponding (user defined) artificial summaries $\mathbf{S}^*:=\mathbf{S}(\mathbf{Y}^*)$ are created. By denoting with $\mathbf{S}(\mathbf{Y})$ the summaries for the observed data, when $\parallel \mathbf{S}^*-\mathbf{S}(\mathbf{Y}))\parallel <\delta$ for some distance $\parallel\cdot\parallel$ and a small positive constant $\delta,$ the corresponding  $\boldsymbol{\theta}^*$ is accepted (or receives a high ``weight'') and becomes a sample from the approximate posterior $\pi_\delta(\boldsymbol{\theta}|\mathbf{S}(\mathbf{Y}))$. For $\delta$ small enough $\pi_\delta(\boldsymbol{\theta}|\mathbf{S}(\mathbf{Y}))\approx \pi(\boldsymbol{\theta}|\mathbf{S}(\mathbf{Y}))$. In general however, this simple scheme is difficult to tune, as discussed in section \ref{sec:g-and-k}.
 An alternative, more recent and therefore less studied ``information reduction strategy'' is the synthetic likelihoods methodology \citep{wood2010statistical}. Briefly, synthetic likelihoods can be considered a ``parametric method'' as it assumes the data summaries to be jointly multivariate Gaussian as $\mathbf{S}(\mathbf{Y})\sim \mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{\theta}},\boldsymbol{\Sigma}_{\boldsymbol{\theta}})$, with unknown mean $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and unknown covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$. Instead, ABC does not make any parametric assumption on the summaries. Notation-wise we make explicit the dependence of the mean and covariance on $\boldsymbol{\theta}$, as later on it will be important to highlight this fact when estimating $\boldsymbol{\theta}$ (e.g. in equation \eqref{eq:synlik-maximization}).

Estimators for $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$ are found by simulating $R$ datasets independently from the assumed data-generating model, conditionally on some $\boldsymbol{\theta}$.
We denote the artificial datasets simulated from model \eqref{eq:state-space-general} with $\mathbf{y}_1^{*},...,\mathbf{y}_R^{*}$. These are such that $\dim(\mathbf{y}^*_r)=\dim(\mathbf{Y})$, $r=1,...,R$. For each dataset \cite{wood2010statistical} constructs the corresponding (vector valued) summary $\mathbf{s}_r^*$, such that $\dim(\mathbf{s}_r^{*})=\dim(\mathbf{S}(\mathbf{Y}))$. Then he computes the following estimators:
\[
\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}} = \frac{1}{R}\sum_{r=1}^{R}\mathbf{s}_r^{*},\qquad
\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}} = \frac{1}{R-1}\sum_{r=1}^{R}(\mathbf{s}_r^{*}-\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}})(\mathbf{s}_r^{*}-\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}})'.
\]
Thanks to the Gaussian assumption, a ``synthetic likelihood'' based on the summaries for the observed data is defined as $\hat{p}(\mathbf{S}(\mathbf{Y})|\boldsymbol{\theta}):=\mathcal{N}(\mathbf{S}(\mathbf{Y});\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$. It is then possible to numerically maximize $\hat{p}(\mathbf{S}(\mathbf{Y})|\boldsymbol{\theta})$ with respect to $\boldsymbol{\theta}$ or compute the MAP (maximum a posteriory) for the associated posterior distribution using MCMC, by considering uniform priors for the parameters. In order to construct synthetic likelihoods the only parameter that needs to be set is $R$, which is certainly easier to tune than the $\delta$ in ABC methods, also because in the latter case we want a small $\delta$ and this has an immediate impact on the algorithm acceptance rate (the smaller the $\delta$ the higher the number of rejections).

\section{SAEM with synthetic likelihoods }\label{sec:saem-sl}

We now use synthetic likelihoods (SL) to develop a likelihood-free version of SAEM. The two main consequences of our approach are (i) we allow the SAEM optimizer to be implemented for complex/intractable models and (ii) the algorithm does not require any ``advanced'' tuning, that is the user only needs to set two parameters ($L$ and $R$, the larger the better, given a limited computational budget). With specific reference to existing synthetic likelihoods approaches, with SAEM-SL the user does not need to set-up an MCMC implementation, as instead required in \cite{wood2010statistical} and \cite{price2016bayesian} and this usually comes with a need for expert tuning, as discussed in the introduction. Denote with $\mathbf{S}(\mathbf{Y})$ and $\mathbf{S}(\mathbf{X})$ \textit{user-defined} summary statistics for $\mathbf{Y}$ and $\mathbf{X}$ respectively. These are meant to encode information regarding $\boldsymbol{\theta}$, but should not be confused with the minimal sufficient statistics $\mathbf{S}_c(\cdot)$ of the complete likelihood.

Assume the complete likelihood for $\mathbf{s}=(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X})$ to be a multivariate Gaussian with mean $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and covariance $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$. That is for the corresponding complete synthetic log-likelihood evaluated at $\mathbf{s}$ we set
\begin{equation}
L_c(\mathbf{s};\boldsymbol{\theta}):=L_c(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X});\boldsymbol{\theta})=\log \mathcal{N}(\mathbf{s};\boldsymbol{\mu}_{\boldsymbol{\theta}},\boldsymbol{\Sigma}_{\boldsymbol{\theta}}).
\label{eq:gauss-complete}
\end{equation}
Of course $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$ are in general unknown. Also, here $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$ are not the same moments defined for the data likelihood in section \ref{sec:synlik}, as the latter is based solely on $\mathbf{S}(\mathbf{Y})$.

Here we illustrate an instance of SL for the current $\boldsymbol{\theta}$, this returning estimators $\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}}$ and $\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}}$. We call this procedure ``internal SAEM-SL'' to be distinguished from an ``external'' procedure described later. Crucially,
 thanks to the Gaussian assumption set on the summaries it is known that $(\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$ is jointly sufficient for $({\boldsymbol{\mu}}_{\boldsymbol{\theta}},{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$. Hence we can set ${\mathbf{S}}_c(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X})):=(\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$ without the need to perform analytic calculations. Then we plug the obtained moment estimates into the ``external SAEM-SL''. While below we describe the several steps of our approach, the complete procedure is illustrated in algorithm \ref{alg:saem-sl}.

\subsection*{Internal SAEM-SL}
Assume a value for $\boldsymbol{\theta}$ is given.
\begin{enumerate}
\item Simulate independently from the model $R$ realizations of processes $\{\mathbf{X}_t\}$ and $\{\mathbf{Y}_t\}$: $\mathbf{x}_r^*\sim p_\mathbf{X}(\mathbf{X}_{0:N};\boldsymbol{\theta})$ and $\mathbf{y}_r^*\sim p_{\mathbf{Y}|\mathbf{X}}(\mathbf{Y}_{1:n}|\mathbf{X}_{1:n};\boldsymbol{\theta})$, $r=1,...,R$.
\item compute user-defined summaries $\mathbf{s}_r^*=(\mathbf{S}(\mathbf{y}_r^*),\mathbf{S}(\mathbf{x}_r^*))$ for each $r$.
\item estimate moments (sufficient statistics for $(\boldsymbol{\mu}_{\boldsymbol{\theta}},\boldsymbol{\Sigma}_{\boldsymbol{\theta}})$)
\begin{equation}
\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}} = \frac{1}{R}\sum_{r=1}^{R}\mathbf{s}_r^{*},\qquad
\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}} = \frac{1}{R-1}\sum_{r=1}^{R}(\mathbf{s}_r^{*}-\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}})(\mathbf{s}_r^{*}-\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}})'\label{eq:synthetic-sample-moments}.
\end{equation}
\end{enumerate}

\subsection*{External SAEM-SL}

A generic iteration of SAEM is executed using the estimators $(\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$ obtained above. At iteration $k$ we update separately the moments for the complete loglikelihood as
\begin{align}
\hat{\boldsymbol{\mu}}^{(k)}_{\boldsymbol{\theta}}&=\hat{\boldsymbol{\mu}}^{(k-1)}_{\boldsymbol{\theta}}+\gamma^{(k)}( \hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}}-\hat{\boldsymbol{\mu}}^{(k-1)}_{\boldsymbol{\theta}})\label{eq:update-mean}\\
\hat{\boldsymbol{\Sigma}}^{(k)}_{\boldsymbol{\theta}}&=\hat{\boldsymbol{\Sigma}}^{(k-1)}_{\boldsymbol{\theta}}+\gamma^{(k)}( \hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}}-\hat{\boldsymbol{\Sigma}}^{(k-1)}_{\boldsymbol{\theta}})\label{eq:update-cov}.
\end{align}
From the quantities computed in \eqref{eq:update-mean}-\eqref{eq:update-cov} extract the corresponding mean and covariances for the two simulated processes, that is set $\hat{\boldsymbol{\mu}}^{(k)}\equiv \hat{\boldsymbol{\mu}}^{(k)}_{\boldsymbol{\theta}}=(\hat{\boldsymbol{\mu}}_{x},\hat{\boldsymbol{\mu}}_{y})$ and
\[
\hat{\boldsymbol{\Sigma}}^{(k)}\equiv \hat{\boldsymbol{\Sigma}}^{(k)}_{\boldsymbol{\theta}} = \left[\begin{array}{cc}
\hat{\boldsymbol{\Sigma}}_{x} & \hat{\boldsymbol{\Sigma}}_{xy}   \\
\hat{\boldsymbol{\Sigma}}_{yx} & \hat{\boldsymbol{\Sigma}}_{y} \\
\end{array}\right].
\]
We now sample $\mathbf{S}(\mathbf{X}^{(k)})$ conditionally on $\mathbf{S}(\mathbf{Y})$ by using well known properties of Gaussian distributions: we have $\mathbf{S}(\mathbf{X}^{(k)})|\mathbf{S}(\mathbf{Y})\sim \mathcal{N}(\hat{\boldsymbol{\mu}}_{x|y,{\boldsymbol{\theta}}}^{(k)},\hat{\boldsymbol{\Sigma}}_{x|y,{\boldsymbol{\theta}}}^{(k)})$ where (here we drop the index $k$ and subscript ${\boldsymbol{\theta}}$ for ease of reading)
\begin{align}
\hat{\boldsymbol{\mu}}_{x|y} & = \hat{\boldsymbol{\mu}}_{x}+\hat{\boldsymbol{\Sigma}}_{xy}\hat{\boldsymbol{\Sigma}}_{y}^{-1}(\mathbf{S}(\mathbf{Y})-\hat{\boldsymbol{\mu}}_{y})\label{eq:condmean}\\
\hat{\boldsymbol{\Sigma}}_{x|y} & = \hat{\boldsymbol{\Sigma}}_{x}-\hat{\boldsymbol{\Sigma}}_{xy}\hat{\boldsymbol{\Sigma}}_{y}^{-1}\hat{\boldsymbol{\Sigma}}_{yx}.\label{eq:condcov}
\end{align}
Some care should be used with the covariance matrix $\hat{\boldsymbol{\Sigma}}_{x|y}^{(k)}$ when sampling $\mathbf{S}(\mathbf{X}^{(k)})|\mathbf{S}(\mathbf{Y})$, as such covariance must be positive semi-definite. In fact $\hat{\boldsymbol{\Sigma}}_{x|y}^{(k)}$ is extracted from $\hat{\boldsymbol{\Sigma}}^{(k)}$, however while it is known that a linear combination (via \eqref{eq:update-cov}) of semi-positive definite matrices is a semi-positive definite matrix and while the sample covariance created in the Internal SAEM-SL is by definition semi-positive definite, in numerical calculations it can still happen that the resulting matrix has negative eigenvalues due to round-off errors in floating point approximations. Therefore, before using $\hat{\boldsymbol{\Sigma}}_{x|y}^{(k)}$ in our conditional sampling, we first check whether this is a positive definite matrix. If it turns out to be positive definite, by using the Cholesky decomposition of $\hat{\boldsymbol{\Sigma}}_{x|y}^{(k)}$, then we proceed with the sampling, that is we obtain the lower triangular matrix $\mathbf{M}$ such that $\mathbf{M}\bM'=\hat{\boldsymbol{\Sigma}}_{x|y}^{(k)}$ and then sample $\mathbf{S}(\mathbf{X}^{(k)})|\mathbf{S}(\mathbf{Y})$ using $\mathbf{S}(\mathbf{X}^{(k)}):=\hat{\boldsymbol{\mu}}_{x|y}^{(k)}+\mathbf{M}\mathbf{z}$, where $\mathbf{z}$ is a vector of independent draws from the standard normal distribution. For those rare instances where it is not positive definite (and not even semi-positive definite) it is possible to compute a ``nearest
semi-positive definite matrix'' (e.g. \citealp{higham1988computing}) and use this one for the sampling.

With the $\mathbf{S}(\mathbf{X}^{(k)})$ we have just sampled,  set $\mathbf{s}^{(k)}:=(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X}^{(k)}))$
and compute the M-step
\begin{equation}
\hat{\boldsymbol{\theta}}^{(k)}=\arg \max_{\boldsymbol{\theta} \in \Theta} \hat{L}_c(\mathbf{s}^{(k)};\boldsymbol{\theta})=\arg \max_{\boldsymbol{\theta} \in \Theta} \log \mathcal{N}(\mathbf{s}^{(k)};\boldsymbol{\mu}_{\boldsymbol{\theta}},\boldsymbol{\Sigma}_{\boldsymbol{\theta}})\label{eq:synlik-maximization}
\end{equation}
where maximization is obtained numerically, for example using $L$ iterations of a Nelder--Mead simplex. Each iteration of the maximizer used for \eqref{eq:synlik-maximization} tests a different value of $\boldsymbol{\theta}$ by invoking the Internal-SL procedure, hence each call evaluates the complete synthetic loglikelihood using a different set of simulated moments $(\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$ produced using the synthetic likelihoods approach.
At the end of the M-step \eqref{eq:synlik-maximization}, besides $\hat{\boldsymbol{\theta}}^{(k)}$ we also retrieve the corresponding ``optimal moments'' $(\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$. Optimal moments are passed to \eqref{eq:update-mean}-\eqref{eq:update-cov} for a further iteration of the External SAEM-SL procedure. Algorithm \ref{alg:saem-sl} details a single iteration of the SAEM-SL procedure, which should be executed for $k=1,...,K$ iterations, with quantities having $k=0$ denoting input/starting values. The generality of the algorithm implies that to implement our three case studies we did not really need to produce significant changes to our test code.

We initialize algorithm \ref{alg:saem-sl} by setting $\boldsymbol{\mu}^{(0)}$ and $\boldsymbol{\Sigma}^{(0)}$ to a vector of zeros and to a diagonal matrix with positive entries $\delta \mathbf{I}_d$ respectively, with $\delta=10^{-12}$ and $\mathbf{I}_d$ the $d$-dimensional identity matrix with $d$ the length of vector $(\mathbf{S}(\mathbf{X}),\mathbf{S}(\mathbf{Y}))$. Notice that each time a numeric maximizer evaluates \eqref{eq:synlik-maximization} for the current candidate parameters $\boldsymbol{\theta}^c$ the vector $\mathbf{s}^{(k)}$ does not vary within the Internal SL: $\mathbf{s}^{(k)}$ contains both the observed summaries and the summaries for the latent state $\mathbf{S}(\mathbf{X}^{(k)})$, which should not be altered when \eqref{eq:synlik-maximization} is executed. Also notice that while in step 2 of the Internal-SL procedure the quantity $\mathbf{S}(\mathbf{x}_r^*)$ is computed from the user defined set of summaries, the $\mathbf{S}(\mathbf{X}^{(k)})$ that is plugged into $\mathbf{s}^{(k)}$ is instead sampled from a multivariate Gaussian distribution.

For the sake of discussion, here we illustrate an ideal scenario which in practice cannot be attained for most realistic models, namely assuming that (a) the user defined sumaries $\mathbf{s}=(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X}))$ are jointly sufficient statistics for $\boldsymbol{\theta}$, and that (b) $\mathbf{s}$ is distributed according to a multivariate Gaussian, though (b) is much easier to obtain than (a). Then under (a)--(b) SAEM-SL does not result in any approximation and converges to a (local) maximizer of the data likelihood function under the same assumptions set for SAEM in \cite{Delyon1999}. In fact, if $\mathbf{s}$ is sufficient for $\boldsymbol{\theta}$ then $\mathbf{s}$ encodes the same amount of information regarding $\boldsymbol{\theta}$ as the couple $(\mathbf{Y},\mathbf{X})$, hence $\mathbf{S}_c(\mathbf{Y},\mathbf{X})\equiv \mathbf{S}_c(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X}))$. Then, under the additional Gaussian assumption, we have $\mathbf{S}_c(\mathbf{Y},\mathbf{X})\equiv \mathbf{S}_c(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X}))=(\hat{\boldsymbol{\mu}}_{\boldsymbol{\theta}},\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}})$. Therefore, since the synthetic complete loglikelihood \eqref{eq:gauss-complete} is  a member of the exponential family and can thus be written as \eqref{eq:complete-likelihood}, the two assumptions for the ``ideal'' scenario fit within the SAEM approach in section \ref{sec:standard-SAEM}. Even if the two assumptions (a)--(b) are met, deviations from what is expected from the theory is due to the non-availability of an explicit M-step, meaning that while in SAEM equation $\eqref{eq:M-step}$ can be sometimes solved analytically, with SAEM-SL the corresponding \eqref{eq:synlik-maximization} has to be solved numerically. Hence, for a finite computational budget and depending on the difficulty in exploring the surface in \eqref{eq:synlik-maximization}, we might not really obtain the exact maximizer from the M-step.

\begin{algorithm}
\small
\caption{A single iteration of SAEM-SL}\label{alg:saem-sl}
\begin{algorithmic}
\State \textbf{Input:} observed summaries $\mathbf{S}(\mathbf{Y})$, positive integers $L$ and $R$. Values for $\hat{\boldsymbol{\theta}}^{(k-1)}$, $\hat{\boldsymbol{\mu}}^{(k-1)}$ and $\hat{\boldsymbol{\Sigma}}^{(k-1)}$.
\State \textbf{Output:} $\hat{\boldsymbol{\theta}}^{(k)}$.\\

\State At iteration $k$ of External SAEM-SL:
\State 1. Extract $\hat{\boldsymbol{\mu}}_{x}$, $\hat{\boldsymbol{\mu}}_{y}$, $\hat{\boldsymbol{\Sigma}}_{x}$, $\hat{\boldsymbol{\Sigma}}_{y}$, $\hat{\boldsymbol{\Sigma}}_{xy}$ and $\hat{\boldsymbol{\Sigma}}_{yx}$ from $\hat{\boldsymbol{\mu}}^{(k-1)}$ and $\hat{\boldsymbol{\Sigma}}^{(k-1)}$. Compute conditional moments  $\hat{\boldsymbol{\mu}}_{x|y}$, $\hat{\boldsymbol{\Sigma}}_{x|y}$ using \eqref{eq:condmean}--\eqref{eq:condcov}.
\State 2. Sample $\mathbf{S}(\mathbf{X}^{(k-1)})|\mathbf{S}(\mathbf{Y})\sim \mathcal{N}(\hat{\boldsymbol{\mu}}_{x|y}^{(k-1)},\hat{\boldsymbol{\Sigma}}_{x|y}^{(k-1)})$ and form $\mathbf{s}^{(k-1)}:=(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X}^{(k-1)}))$.
\State 3. Obtain $(\boldsymbol{\theta}^{(k)},\boldsymbol{\mu}^{(k)},\boldsymbol{\Sigma}^{(k)})$ from \texttt{InternalSL($\mathbf{s}^{(k-1)},\hat{\boldsymbol{\theta}}^{(k-1)},R$)} starting at $\hat{\boldsymbol{\theta}}^{(k-1)}$.
\State 4. Increase $k:=k+1$ and go to step 1.
\\\hrulefill\\

\textbf{Function \texttt{InternalSL($\mathbf{s}^{(k-1)},\boldsymbol{\theta}^{\mathrm{start}},R$)}:}\\
\State \textbf{Input:} $\mathbf{s}^{(k-1)}$, starting parameters $\boldsymbol{\theta}^{\mathrm{start}}$, a positive integer $R$. Functions to compute simulated summaries $\mathbf{S}(\mathbf{y}^*)$ and $\mathbf{S}(\mathbf{x}^*)$ must be available.
\State \textbf{Output:} the best found $\boldsymbol{\theta}^{*}$ maximizing $\log \mathcal{N}(\mathbf{s}^{(k)};\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\Sigma}})$ and corresponding $(\boldsymbol{\mu}^*,\boldsymbol{\Sigma}^*)$.\\

\State Here $\boldsymbol{\theta}^c$ denotes a generic candidate value.
\State i. Simulate $\mathbf{x}_r^*\sim p_\mathbf{X}(\mathbf{X}_{0:N};\boldsymbol{\theta}^c)$, $\mathbf{y}_r^*\sim p_{\mathbf{Y}|\mathbf{X}}(\mathbf{Y}_{1:n}|\mathbf{X}_{1:n};\boldsymbol{\theta}^c)$ for $r=1,...,R$.
\State ii. Compute user-defined summaries $\mathbf{s}_r^*=(\mathbf{S}(\mathbf{y}_r^*),\mathbf{S}(\mathbf{x}_r^*))$ for $r=1,...,R$. Construct the corresponding $(\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\Sigma}})$.
\State iii. Evaluate $\log \mathcal{N}(\mathbf{s}^{(k)};\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\Sigma}})$. \\
Use a numerical procedure that performs (i)--(iii) $L$ times to find the best $\boldsymbol{\theta}^*$ maximizing $\log \mathcal{N}(\mathbf{s}^{(k)};\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\Sigma}})$ for varying $\boldsymbol{\theta}^c$.
Denote with $({\boldsymbol{\mu}}^*,\hat{\boldsymbol{\Sigma}}^*)$ the simulated moments corresponding to the best found $\boldsymbol{\theta}^*$. Set  $\boldsymbol{\theta}^{(k)}:=\boldsymbol{\theta}^*$.
\State iv. Update moments:
\begin{align*}
\hat{\boldsymbol{\mu}}^{(k)}&=\hat{\boldsymbol{\mu}}^{(k-1)}+\gamma^{(k)}( \hat{\boldsymbol{\mu}}^*-\hat{\boldsymbol{\mu}}^{(k-1)})\\
\hat{\boldsymbol{\Sigma}}^{(k)}&=\hat{\boldsymbol{\Sigma}}^{(k-1)}+\gamma^{(k)}( \hat{\boldsymbol{\Sigma}}^*-\hat{\boldsymbol{\Sigma}}^{(k-1)}).
\end{align*}
Return $(\boldsymbol{\theta}^{(k)},\hat{\boldsymbol{\mu}}^{(k)},\hat{\boldsymbol{\Sigma}}^{(k)})$.
\end{algorithmic}
\end{algorithm}

Notice that if we were to consider a maximization of \eqref{eq:synlik-maximization} with $\mathbf{s}^{(k)}$ containing exclusively the $\mathbf{S}(\mathbf{Y})$ component, then $\log \mathcal{N}(\mathbf{s}^{(k)};\boldsymbol{\mu}_{\boldsymbol{\theta}},\boldsymbol{\Sigma}_{\boldsymbol{\theta}})$ would be the synthetic log-likelihood function as defined in \cite{wood2010statistical}.

The advantage of the proposed method, which we call SAEM-SL (SAEM using synthetic likelihoods) is that (i) unlike the ``standard'' SAEM, SAEM-SL is completely plug-and-play, only the ability to simulate from the model is required.  The disadvantage is the requirement for the user to specify a set of summaries $\mathbf{s}=(\mathbf{S}(\mathbf{Y}),\mathbf{S}(\mathbf{X}))$ and that for each iteration of SAEM-SL the maximization of the loglikelihood \eqref{eq:synlik-maximization} consists of an iterative procedure. The choice of $\mathbf{s}$ is delicate as it is the only way to encode information about the model and the data, hence a poor choice for $\mathbf{s}$ will result in poor inference, even when the data generating model is correctly chosen. Unfortunately, outside simulation studies there is no way to assess whether  the statistics we set are appropriate. On the other hand the methodology expands considerably the set of problems that is possible to treat with SAEM. The standard SAEM itself is unable to deal with very complex models, unless it is possible to derive the necessary constructs (sufficient statistics for the complete likelihood and corresponding updating equations for the M-step), which is often a difficult task. If the model has an intractable complete likelihood, the task is actually impossible. However, even when it is tractable, deriving the necessary quantities to run SAEM is at best a tedious and prone to error task, when a range of models has to be screened. We believe that by merging the synthetic likelihoods idea with SAEM, the resulting SAEM-SL algorithm provides the modeller with a practical tool for statistical exploration.

\subsection{Direct maximization of the data synthetic likelihood and pseudo-marginal Bayesian SL}\label{sec:wood-synlik}
It is natural to wonder whether using SAEM-SL is worth at all, considering that it is possible to maximize the \textit{data synthetic likelihood} in \eqref{eq:wood-synloglik}, as opposed to the \textit{complete data synthetic loglikelihood} \eqref{eq:synlik-maximization}, as well as sampling via MCMC from the posterior $\pi(\boldsymbol{\theta}|\mathbf{S}(\mathbf{Y})$. Direct maximization of the data synthetic likelihood is given by
\begin{equation}
\hat{\boldsymbol{\theta}}=\arg \max_{\boldsymbol{\theta} \in \Theta} \log \mathcal{N}(\mathbf{S}(\mathbf{Y});\boldsymbol{\mu}(\boldsymbol{\theta}),\boldsymbol{\Sigma}(\boldsymbol{\theta})).\label{eq:wood-synloglik}
\end{equation}
However, as remarked in \cite{wood2010statistical}, the roughness of the surface implied by \eqref{eq:wood-synloglik} will make the optimization difficult. Therefore he recommends using a Markov chain Monte Carlo (MCMC) approach, setting independent uniform priors on each parameter, then returning the maximum a-posteriori (MAP) value. A MAP obtained this way is numerically equivalent to the maximum likelihood estimator. However, when considering an MCMC approach we prefer to follow \cite{price2016bayesian}, as the authors first construct an unbiased estimate $\hat{p}_U(\mathbf{S}(\mathbf{Y});\boldsymbol{\theta})$ of the (unknown) synthetic likelihood density $p(\mathbf{S}(\mathbf{Y});\boldsymbol{\theta})$ (the latter depends on the unknown $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$). This enables the use of a pseudo-marginal approach (in the spirit of \citealp{andrieu2009pseudo}), where parameters are sampled exactly from $\pi(\boldsymbol{\theta}|\mathbf{S}(\mathbf{Y}))\propto p(\mathbf{S}(\mathbf{Y});\boldsymbol{\theta})\pi(\boldsymbol{\theta})$, regardless the number of model simulations $R$ required to estimate $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$. In section \ref{sec:theophylline-results} we consider both the pseudo-marginal Bayesian SL approach and the optimization of \eqref{eq:wood-synloglik}.

\section{Simulation studies}\label{sec:simulations}
Simulations were coded in MATLAB (except for examples using the R \texttt{pomp} package) and executed on a Intel Core i7-4790 CPU 3.60 GhZ. In SAEM we always set $\gamma_k=1$ for the first $K_1$ iterations and $\gamma_k=(k-K_1)^{-1}$ for $k\geq K_1$ as in \cite{lavielle2014mixed}. However, we found that small modifications to this setup do not affect results significantly, that is using $\gamma_k=(k-K_1)^{-\beta}$ for $k\geq K_1$ and some $\beta \in (0.5,1]$ is also valid.
The numerical maximization of \eqref{eq:synlik-maximization} is performed using the Nelder-Mead simplex as implemented in the \textsc{Matlab} function \texttt{fminsearch}.
We compare our results with standard algorithms for Bayesian and ``classical'' inference: for Bayesian inference we considered a particle marginal methods (PMM) (\citealp{andrieu2009pseudo}, also including a version incorporating synthetic likelihoods, \citealp{price2016bayesian}). For the classical approach we considered the improved iterated filtering (denoted in literature as IF2) found in \cite{ionides2015inference} and SAEM incorporating a SMC step. To perform a fair comparison we make use of well tested code to fit models with PMM and IF2 via the R \texttt{pomp} package \citep{king2015statistical}. For a comparison between IF2, PMM, standard synthetic likelihoods and beyond the reader is referred to \cite{fasiolo2016comparison}. IF2 and PMM as implemented in \texttt{pomp} use the bootstrap filter sequential Monte Carlo (SMC) algorithm to approximate the likelihood function. The reader interested in using more advanced SMC algorithms suitable for a specific problem are referred to \cite{Cappe2007} and \cite{kantas2015particle}. For example in section \ref{sec:theoph} we use a filter specific for SDE models proposed in \cite{golightly2011bayesian}.

\subsection{Non-linear Gaussian state-space model}\label{sec:nonlingauss}

Here we study a simple non-linear model, using a setup similar to \cite{jasra2012filtering}. See also \cite{picchini2016coupling} for inference using algorithm \ref{alg:saem-smc} as well as SAEM coupled with an ABC filter.
We have
\begin{align}
\begin{cases}
Y_j = X_j + \sigma_y\nu_j,\qquad j\geq 1\\
X_j = 2\sin(e^{X_{j-1}})+\sigma_x\tau_j,
\end{cases}
\label{eq:nonlingauss-ssm}
\end{align}
with $\nu_j,\tau_j\sim N(0,1)$ i.i.d. and $X_0=0$. We assume  $\sigma_x,\sigma_y>0$ as the only unknowns and therefore conduct inference for $\boldsymbol{\theta} = (\sigma^2_x,\sigma^2_y)$.

We first consider the standard SAEM methodology outlined in section \ref{sec:standard-SAEM}, and therefore construct the set of sufficient statistics corresponding to the complete log-likelihood $L_c(\mathbf{Y},\mathbf{X})$. For this model the task is simple since $Y_j|X_j\sim N(X_j,\sigma^2_y)$ and $X_j|X_{j-1}\sim N(2\sin(e^{X_{j-1}}),\sigma^2_x)$ and it is easy to show that $S_{\sigma^2_x}=\sum_{j=1}^n (X_j-2\sin(e^{X_{j-1}}))^2$ and $S_{\sigma^2_y}=\sum_{j=1}^n (Y_j-X_{j})^2$ are sufficient for $\sigma^2_x$ and $\sigma^2_y$ respectively. By plugging these statistics into $L_c(\mathbf{Y},\mathbf{X})$ and equating to zero the gradient of $L_c$ with respect to $(\sigma^2_x,\sigma^2_y)$, we find that the M-step of SAEM results in updated values for $\sigma^2_x$ and $\sigma^2_y$ given by $S_{\sigma^2_x}/n$ and $S_{\sigma^2_y}/n$ respectively.
In the following, we write SAEM-SMC to refer to Algorithm \ref{alg:saem-smc}.

We generate $n=50$ observations for $\{Y_j\}$ using model \eqref{eq:nonlingauss-ssm}  with $\sigma^2_x=\sigma^2_y=5$. Our setup consists in running 30 independent experiments with SAEM-SMC: for each experiment we simulate parameter starting values for $(\log\sigma_x,\log\sigma_y)$ independently generated from a bivariate Gaussian distribution with mean the true value of the parameter, i.e. $(\log\sqrt{5},\log\sqrt{5})$, and diagonal covariance matrix having (2,2) on its diagonal. Hence the starting values are very spread.
We take $K_1=300$ as the number of warmup iterations (see beginning of section \ref{sec:simulations}) and use different numbers of particles $M$ in our simulation studies, see Table \ref{tab:nonlingauss}. We impose resampling when the effective sample size ESS gets smaller than $\bar{M}=200$, for any value of $M$. In summary, for all 30 simulations we use the same data and the same setup except that in each simulation we use different starting values for the parameters. Table \ref{tab:nonlingauss} reports the median of the 30 estimates and their $1^{st}-3^{rd}$ quartiles. Simulations for $\sigma_y$ converge to completely wrong values. We also experimented with $M=5,000$ using $\bar{M}=2,000$ but this does not solve the problem with SAEM-SMC, even if we let the algorithm start at the true parameter values. However, in \cite{picchini2016coupling} we learned that SAEM-SMC is very sensitive to the value of $\bar{M}$, and in particular that it is better to set a very low $\bar{M}$: in fact, when using $\bar{M}=20$ with $M=1,000$ results improve sensibly, see Table \ref{tab:nonlingauss}, though estimation of $\sigma_y$ is still unsatisfactory.

\begin{table}
\centering
\begin{tabular}{lcccc}
\hline

 $(M,\bar{M})$ & (500,200) & (1000,200) & (2000,200) & (1000,20)  \\
\hline
\multicolumn{4}{l}{$\sigma_x$ (true value 2.23)}\\
SAEM-SMC & 2.54 [2.53,2.54] & 2.55 [2.54,2.56] & 2.55 [2.54,2.56] & 1.99 [1.85,2.14] \\
IF2* & 1.26 [1.21,1.41] & 1.35 [1.28,1.41] & 1.33 [1.28,1.40] & --\\
\hline
\multicolumn{4}{l}{$\sigma_y$ (true value 2.23)}\\
SAEM-SMC & 0.11 [0.10,0.13] & 0.06 [0.06,0.07] & 0.04 [0.03,0.04] & 1.23 [1.00,1.39]\\
IF2* & 1.62 [1.56,1.75] & 1.64 [1.58,1.67] & 1.63 [1.59,1.67] & --\\
\hline
\end{tabular}
\caption{\footnotesize{Non-linear Gaussian model: medians and $1^{st}-3^{rd}$ quartiles for estimates obtained on 30 independent simulations, using different number of particles $M$ and different methods. (*)The IF2 method resamples at every time point, while SAEM-SMC resamples only when $ESS<\bar{M}$. Hence for IF2 it is always $\bar{M}\equiv
M$. }}
\label{tab:nonlingauss}
\end{table}

\begin{table}
\centering
\begin{tabular}{lccc}
\hline

 $R$ & 500 & 1000 & 2000   \\
\hline
\multicolumn{4}{l}{$\sigma_x$ (true value 2.23)}\\
SAEM-SL & 1.96 [1.27,2.52] & 1.90 [1.13,2.39] & 2.07 [1.57,2.18]\\
\hline
\multicolumn{4}{l}{$\sigma_y$ (true value 2.23)}\\
SAEM-SL & 2.35 [1.40,2.77] & 1.94 [1.30,2.44] & 1.70 [1.44,2.22]\\
\hline
\end{tabular}
\caption{\footnotesize{Non-linear Gaussian model: medians and $1^{st}-3^{rd}$ quartiles for estimates obtained on 30 independent simulations using SAEM-SL. }}
\label{tab:nonlingauss-SL}
\end{table}

We now compare the results above with the iterated filtering IF2 \citep{ionides2015inference} using the R package \texttt{pomp}. We do not provide a detailed description of IF2 here: it suffices to say that in IF2 particles are generated for both $\boldsymbol{\theta}$ (e.g. via perturbations using random walks) and for the systems state (using the bootstrap filter). Moreover a ``temperature'' parameter (to use an analogy with the simulated annealing optimization method) is let decrease until the algorithm ``freezes'' around an approximated MLE. This parameter that here we denote with $\rho$ is let decrease in $\rho\in\{0.9,0.7,0.4,0.3,0.2\}$ where the first value is used for the first 500 iterations of IF2, then   each of the remaining values is used for 100 iterations, for a total of 900 iterations. Notice that the tested version of \texttt{pomp} (v. 1.4.1.1) uses a bootstrap filter that resamples at each time point, and therefore results obtained with IF2 are not directly comparable with SAEM-SMC, hence the asterisk in Table \ref{tab:nonlingauss}. The output from one of the experiments obtained with $M=1,000$ is in Figure \ref{fig:nonlingauss-IF2}. From Figure \ref{fig:nonlingauss-IF2} we notice that the last major improvement for the loglikelihood maximization takes place  at iteration 600 when $\rho$ becomes $\rho=0.7$, and reducing $\rho$ further does not give any significant benefit (we have verified this in a number of experiments with this model), therefore we are confident about our setup. With IF2 the estimation of $\sigma_y$ is much improved compared to SAEM-SMC, however inference for $\sigma_x$ is more biased than with SAEM-SMC.

\begin{figure}
\centering
\includegraphics[scale=0.6]{nonlingauss-IF2}
\caption{\footnotesize{Non-linear Gaussian model: traces obtained for a single experiment with IF2 when using $M=1,000$ particles. (Top left) evolution of the loglikelihood function; (bottom left) evolution of $\sigma_x$; (top right) evolution of $\sigma_y$.}}
\label{fig:nonlingauss-IF2}
\end{figure}

We now consider a particle marginal method (PMM, \citealp{andrieu2009pseudo}) on a single simulation (instead of thirty), as PMM is a full Bayesian methodology and results are not directly comparable with SAEM nor IF2. Once more we make use of tools provided in \texttt{pomp}. We set wide uniform priors $U(0.1,15)$ for both $\sigma_x$ and $\sigma_y$ and use $M=2,000$ particles. Also, we set the algorithm in the most favourable way, by starting it at the true parameter values (here we are only interested in using PMM to obtain exact Bayesian inference, not as a competitor to the other frequentist approaches we have illustrated).
Parameters are proposed using an adaptive MCMC algorithm, and the algorithm is tuned to achieve the optimal 7\% acceptance rate \citep{sherlock2015efficiency}. We obtained the following posterior means and 95\% intervals: $\hat{ \sigma}_x=1.46$ [0.49,2.46], $\hat{ \sigma}_y=1.61$ [0.49,2.40]. Therefore, PMM seems to return values not very different from the ranges provided by IF2.

Finally, we consider inference with SAEM-SL. We performed simulations using $R=500$, 1,000 and 2,000 simulated summaries and $L=40$ iterations for the numerical maximization step. We used the same data as for SAEM-SMC and IF2, however we decide to make the estimation procedure more challenging, so we let the parameter start at random locations sampled from a Gaussian centred at $(\log \sqrt{16},\log \sqrt{16})$ and having diagonal covariance with variances $(2,2)$. Here we need to set a vector of summaries $(\mathbf{s}(\mathbf{x}_r^*),\mathbf{s}(\mathbf{y}_r^*))$. Vector $\mathbf{s}(\mathbf{x}_r^*)$ contains (i) the median value of $\mathbf{x}_r^*$; (ii) the median absolute deviation of $\mathbf{x}_r^*$ and (iii)  the 10th, 20th, 75th and 90th percentile of $\mathbf{x}_r^*$. Vector $\mathbf{s}(\mathbf{y}_r^*)$ contains the same summary functions, except that these are applied to $\mathbf{y}_r^*$. Of course summary functions for observed data $\mathbf{s}(\mathbf{y}):=\mathbf{S}(\mathbf{Y})$ are the same functions considered for $\mathbf{s}(\mathbf{y}^*_r)$ except that now they are evaluated at $\mathbf{y}$.
Same as before we consider thirty repetitions of our experiment: for each experiment we run a warmup of $K_1=10$ iterations and a total number of $K=20$ SAEM-SL iterations. Results are in Table \ref{tab:nonlingauss-SL} and trace plots for the case $R=1,000$ are in Figure \ref{fig:nonlingauss-traceplots_SAEM-SL}. As from Figure \ref{fig:nonlingauss-traceplots_SAEM-SL} we notice that those parameters initialized at much higher values than the true parameter values decay rapidly to approach the true values. As shown in Table \ref{tab:nonlingauss-SL}, the majority of them converges to reasonable values. SAEM-SL produces excellent inference for all tested values of $R$, and convergence is very rapid, well within 10 iterations, corresponding to about 10 seconds on a computer desktop when $R=1,000$.

For  one of the thirty repetitions, Figure \ref{fig:nonlingauss-qqplots} shows the normal qq-plots for the twelve chosen summary statistics (the six statistics in $\mathbf{s}(\mathbf{x}_r)$ and the six in $\mathbf{s}(\mathbf{y}_r)$) for the case $R=2,000$, generated at the optimum returned by SAEM-SL. Clearly there are no major departures from normality. Interestingly, we reach the same conclusion for the case $R=500$ (plots not reported).

\begin{figure}
\centering
\includegraphics[width=17cm,height=7cm]{nonlingauss-SAEM-SL_numsim1000-iter}
\caption{\footnotesize{Non-linear Gaussian model: trace plots for SAEM-SL ($\sigma_x$, left; $\sigma_y$, right) when $R=1,000$ for the thirty estimation procedures. Dashed lines denote the true parameter values.}}
\label{fig:nonlingauss-traceplots_SAEM-SL}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=15cm,height=8cm]{nonlingauss_qqplots_numsim2000}
\caption{\footnotesize{Non-linear Gaussian model: normal qq-plots for $\mathbf{s}(\mathbf{x}_r)$ and $\mathbf{s}(\mathbf{y}_r)$ corresponding to the last iteration of SAEM-SL when $M=2,000$.}}
\label{fig:nonlingauss-qqplots}
\end{figure}

\subsection{A pharmacokinetics model}\label{sec:theoph}

Here we consider a model for pharmacokinetics dynamics. For example, we may imagine to formulate a model to study the Theophylline drug pharmacokinetics. This example has often been described in literature devoted to longitudinal data modelling with random parameters (mixed--effects models), e.g. \cite{pinheiro1995approximations}. Same as in \cite{picchini2014inference} here we do not consider a mixed--effects model. It will be evident that in order to apply a standard SAEM it is required some preliminary analytic effort from the modeller, while running SAEM-SL is basically trivial and plug-and-play, provided ``good'' summary statistics are chosen.
We denote with $X_t$ the level of Theophylline drug concentration in blood at time $t$ (hrs).
Consider the following non-authonomous stochastic differential equation (SDE):
\begin{equation}
dX_t = \biggl(\frac{Dose\cdot K_a \cdot K_e}{Cl}e^{-K_a t}-K_eX_t\biggr)dt + \sigma \sqrt{X_t} dW_t, \qquad t\geq t_0
\label{eq:theophylline-sde}
\end{equation}
where $Dose$ is the known drug oral dose received by a subject, $K_e$ is the elimination rate constant, $K_a$ the absorption rate constant, $Cl$ the clearance of the drug and $\sigma$ the intensity of intrinsic stochastic noise.
We simulate data measured at $n=30$ equispaced sampling times $\{t_1,t_{\Delta},...,t_{30\Delta}\}=\{1,2,...,30\}$ where $\Delta=t_j-t_{j-1}=1$. The drug oral dose is chosen to be 4 mg.
After the drug is administered, we consider as $t_0=0$ the time when the concentration first reaches $X_{t_0}=X_0=8$. The error model is assumed to be linear, $Y_j=X_j+\varepsilon_j$ where the $\varepsilon_j\sim N(0,\sigma_{\varepsilon}^2)$ are i.i.d., $j=1,...,30$. Inference is based on data $\{Y_1,...,Y_{30}\}$ collected at corresponding sampling times. Parameter $K_a$ is assumed known as it is not possible to determine the sufficient statistic for $K_a$ analytically, hence parameters of interest are $\boldsymbol{\theta}=(K_e,Cl,\sigma^2,\sigma_{\varepsilon}^2)$ as $X_0$ is also assumed known.

Equation \eqref{eq:theophylline-sde} has no available closed-form solution, hence simulated data are created in the following way. We first simulate numerically a solution to \eqref{eq:theophylline-sde} using the Euler--Maruyama discretization with stepsize $h=0.05$ on the time interval $[t_0,30]$. The Euler-Maruyama scheme is defined as
\[
X_{t+h} = X_{t} + \biggl(\frac{Dose\cdot K_a \cdot K_e}{Cl}e^{-K_a t}-K_eX_t\biggr)h + \sigma \sqrt{X_t}Z_{t+h},
\]
where the $\{Z_t\}$ are i.i.d. $\mathcal{N}(0,h)$ distributed. The grid of generated values $\mathbf{X}_{0:N}$ is then linearly interpolated at sampling times $\{t_1,...,t_{30}\}$ to give $\mathbf{X}_{1:n}$, and finally residual error is added to $\mathbf{X}_{1:n}$ according to the error model $Y_j=X_j+\varepsilon_j$ as explained above. Data $\{Y_j\}$ are conditionally independent given the latent process $\{X_t\}$ and are generated with $(K_e,K_a,Cl,\sigma^2, \sigma_{\varepsilon}^2)=(0.05,1.492,0.04,0.01,0.102)$. The construction of the sufficient statistics to implement the standard SAEM approach is given in appendix, and this should make evident how applying SAEM can be laborious, even for a one-dimensional model. Considering that the model under study is not particularly complex (when compared with large, multidimensional stochastic models) the amount of analytic work typically required for the application of SAEM can be discouraging. In the results section below we show the simplicity of application of SAEM-SL for this specific example and compare SAEM-SL with a number of alternative approaches.

\subsubsection{Results}\label{sec:theophylline-results}

 Same as in section \ref{sec:nonlingauss}, for SAEM-SMC we run a number of independent repetitions of the estimation procedure: the dataset is shorter than in section \ref{sec:nonlingauss} and despite the need to resort to numerical integration of the SDE, we are able to run 100 estimation procedures in about 300 seconds overall. Each repetition generates a different dataset  using the true parameter values, then for each repetition SAEM-SMC is initialized at the same  parameter values $K_e=0.15$, $Cl=0.135$, $\sigma=0.135$ and $\sigma_{\varepsilon}=0.502$. We always use a warmup of $K_1=80$ iterations, $K=300$, $M=500$ particles and $\bar{M}=100$. We observed an $ESS>100$ at the last time point for each simulation. See Table \ref{tab:theophylline} and Figure \ref{fig:theophylline-SAEM-SMC-iter} for results: clearly $Cl$ and $\sigma$ are not identified. However these results can be improved, at least for $Cl$: for state-space models having additive Gaussian noise and an SDE model discretised using Euler-Maruyama, \cite{golightly2011bayesian} propose a ``fully adapted'' SMC filter where forward simulation of the particles is not blind to data (unlike the bootstrap filter). We refer the reader to \cite{golightly2011bayesian} for details and report results using their approach as SAEM-GW in Table \ref{tab:theophylline}. While $Cl$ is very well identified, the system noise $\sigma$ is still elusive.

With SAEM-SL we only need to set the vector of summaries $(\mathbf{s}(\mathbf{x}_r^*),\mathbf{s}(\mathbf{y}_r^*))$. The vector $\mathbf{s}(\mathbf{x}_r^*)$ contains (i) the median values of $\mathbf{X}_{0:N}^*$ ; (ii) the median absolute deviation of $\mathbf{X}_{0:N}^*$, (iii) a statistic for $\sigma$  computed from $\mathbf{X}_{0:N}^*$  (see below) and (iv) $(\sum_j(Y_{j}^*-X_{j}^*)^2/n)^{1/2}$ with $X_{j}$ the $j$th element of the $n$ interpolated values $\mathbf{X}_{1:n}$. Vector $\mathbf{s}(\mathbf{y}_r^*)$ contains: (i) the median value of $\mathbf{y}_r^*$; (ii) its median absolute deviation; (iii) the slope of the line connecting the first and last simulated observation $(Y_n^*-Y_1^*)/(t_n-t_1)$, since concentrations here show have a markedly decaying behaviour.

In \cite{miao} it is given that, for an SDE of the type $dX_t=\mu(X_t)dt+\sigma g(X_t)dW_t$ with $t\in [0,T]$, we have
\[
\frac{\sum_{\Gamma}|X_{i+1}-X_i|^2}{\sum_{\Gamma}g(X_i)(t_{i+1}-t_i)}\rightarrow \sigma^2 \qquad as\quad |\Gamma|\rightarrow 0
\]
where the convergence is in probability and $\Gamma$ a partition of $[0,T]$.
Therefore we deduce that using the discretization $\{X_0,X_1,...,X_N\}$  produced by the Euler-Maruyama scheme, we can take the square root of the left hand side in the limit above, which should be informative for $\sigma$. We use this as the third summary statistic in $\mathbf{s}(\mathbf{x}_r^*)$.

We used SAEM-SL on the same simulated data produced when implementing SAEM-SMC. We considered  $R=200$ simulated summaries and, since for this example SAEM-SL is computationally more intense than SAEM-SMC, we consider $K_1=50$ and $K=80$, with $L=30$ for the number of iterations in the maximization step. Notice for this example we found benefit in using robust methods for the computation of sample means and covariances, downweight summaries falling in the tails of the multivariate Gaussian synthetic likelihood (this is also recommended in \citealp{wood2010statistical}). Specifically, here we compute the moments \eqref{eq:synthetic-sample-moments} using the method in \cite{olive2010robust}.
See Table \ref{tab:theophylline} and Figure \ref{fig:theophylline-SAEM-SL-iter} for results. Notice that simulations for SAEM-SMC and SAEM-SL start at the same parameter values, even though from Figures \ref{fig:theophylline-SAEM-SMC-iter}--\ref{fig:theophylline-SAEM-SL-iter} it may seem otherwise (that is because SAEM-SMC reaches almost immediately the final values while SAEM-SL converges more slowly).
 SAEM-SL produces satisfactory results on all parameters.
For  one of the one-hundred repetitions, Figure \ref{fig:theophylline-qqplots} shows the normal qq-plots for the seven summary statistics (the four statistics in $\mathbf{s}(\mathbf{x}_r)$ and the three in $\mathbf{s}(\mathbf{y}_r)$), generated at the optimum returned by SAEM-SL. Also for this example, there are no major departures from normality.

We now compare SAEM-SL with the direct maximization of the data synthetic loglikelihood, as illustrated in section \ref{sec:wood-synlik}. In order to perform a fair comparison with SAEM-SL, we run the 100 optimization procedures with $R=200$ anew, using the same data as in SAEM-SMC and SAEM-SL and the same parameters starting value. Also, since with SAEM-SL we had $K\times L=2,400$, here to maximize $\log \mathcal{N}(\mathbf{S}(\mathbf{Y});\boldsymbol{\mu}_{\boldsymbol{\theta}},\boldsymbol{\Sigma}_{\boldsymbol{\theta}})$ we run 2,400 iterations of the Nelder-Mead (NM) simplex. Results in Table \ref{tab:theophylline} are clearly favourable to SAEM-SL. However, as discussed in section \ref{sec:wood-synlik}, this might likely be the consequence of using a deterministic optimizer for the direct optimization of the SL, which ends-up trapped in local optima.  Therefore we also run a single instance of the pseudo-marginal Bayesian SL algorithm of \cite{price2016bayesian}, see section \ref{sec:wood-synlik}. We impose independent uniform  priors $K_e\sim U(0.01,1)$, $Cl\sim U(0.01,20)$, $\sigma\sim U(0.01,0.2)$ and $\sigma_\varepsilon\sim U(0.05,1)$ and run 5000 MCMC iterations. Parameters were proposed using the adaptive Gaussian random walk of \cite{haario-saksman-tamminen}, obtaining an acceptance rate of about 25\%.  We first consider $R=200$, same as for SAEM-SL. Posterior means and 95\% posterior intervals for each parameter are: $\hat{K}_e=0.052$ [0.027,0.074], $\hat{Cl}=0.048$ [0.027,0.091], $\hat{\sigma}=0.105$ [0.024,0.195], $\hat{\sigma}_\varepsilon=0.541$ [0.087,0.9969].
We notice the first two parameters are correctly identified while the latter two parameters are essentially unidentified. This can also be noticed from their MCMC trace plots, spanning the support of the corresponding priors (plots not reported for brevity). Results can be partially improved using $R=2,000$, producing better identification for the first two parameters but not for the latter two, see Figure  \ref{fig:theophylline-bayesSL-hist}. Therefore, using uniform priors as suggested in \cite{wood2010statistical} when the MAP is the only object of interest is not appropriate here and strongly informative priors for $(\sigma,\sigma_\varepsilon)$ might be needed.

\begin{figure}
\centering
\includegraphics[width=14cm,height=6cm]{theophylline_qqplots_numsim500}
\caption{\footnotesize{Theophylline model: normal qq-plots for $\mathbf{s}(\mathbf{x}_r)$ and $\mathbf{s}(\mathbf{y}_r)$ corresponding to the last iteration of SAEM-SL.}}
\label{fig:theophylline-qqplots}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.6]{saemSMC_trajectories_numsim500_100repetitions}
\caption{\footnotesize{Theophylline model: $K=100$ iterations of 100 repetitions of SAEM-SMC using $M=500$. Top: $K_e$ (left) and $Cl$ (right). Bottom: $\sigma$ (left) and $\sigma_{\epsilon}$ (right). Horizontal lines are the true parameter values.}}
\label{fig:theophylline-SAEM-SMC-iter}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=16cm,height=7cm]{saemSL_trajectories_numsim500_100repetitions}
\caption{\footnotesize{Theophylline model: $K=80$ iterations of 100 repetitions of SAEM-SL using $R=200$. Top: $K_e$ (left) and $Cl$ (right). Bottom: $\sigma$ (left) and $\sigma_{\epsilon}$ (right). Dashed lines are the true parameter values.}}
\label{fig:theophylline-SAEM-SL-iter}
\end{figure}

\begin{table}
\centering
\begin{tabular}{lllll}
\hline
 & $K_e$  &  $Cl$ & $\sigma$&  $\sigma_{\varepsilon}$\\
\hline
true values & 0.050& 0.040&   0.100& 0.319\\
SAEM-SMC&  0.045  [0.042,0.049]& 0.085 [0.078,0.094] &   0.171 [0.158,0.184] & 0.395 [0.329,0.465] \\
SAEM-GW & 0.053 [0.049,0.058] & 0.039 [0.035,0.043] & 0.704 [0.549,0.963] & 0.175 [0.119,0.304]\\
SAEM-SL & 0.045 [0.037,0.049]  &  0.032 [0.027,0.038]  &  0.113 [0.088,0.144]  &  0.241 [0.200,0.294]\\
optim. SL & 0.062 [0.031,0.126] & 0.099 [0.045,0.209] & 0.305 [0.133,0.554] & 0.523 [0.396,0.771]\\
\hline
\end{tabular}
\caption{\footnotesize{Theophylline: medians and $1^{st}-3^{rd}$ quartiles for estimates obtained on 100 independent simulations using SAEM-SMC, SAEM-GW, SAEM-SL and the numerical optimization of the data synthetic loglikelihood \eqref{eq:wood-synloglik}.}}
\label{tab:theophylline}
\end{table}

\begin{figure}
\centering
\includegraphics[width=12cm,height=8cm]{theophylline_synlikMCMC_hist}
\caption{\footnotesize{Theophylline model: marginal posteriors from the peudo-marginal Bayesian SL using $R=2,000$. Vertical lines are true parameter values.}}
\label{fig:theophylline-bayesSL-hist}
\end{figure}

\subsection{Noisy $g$-and-$k$ distribution}\label{sec:g-and-k}

We now consider a ``static'' model, namely a $g$-and-$k$ distribution corrupted with noise. Noise-free versions of samples from  $g$-and-$k$ distributions have been considered numerous times in the ABC literature (e.g. \citealp{allingham2009bayesian,fearnhead-prangle(2011),picchini2016approximate}). This is a flexibly shaped distribution that is used to model non-standard data
through a small number of parameters. It is defined by its inverse distribution
function, but has no closed form density hence it is an example of model with an intractable likelihood. Therefore it cannot be dealt with using, say, standard SAEM methods, as the explicit computation of the complete likelihood (and its sufficient statistics) is impossible. However it is trivial to sample from a $g$-and-$k$ distribution and therefore ABC is an appealing methodology for this problem. The quantile function (inverse distribution function) is given by
\begin{equation}
F^{-1}(z;A,B,c,g,k)= A+B\biggl[1+c\frac{1-\exp(-g\cdot r(z))}{1+\exp(-g\cdot r(z))}\biggr](1+r^2(z))^kr(z)
\label{eq:g-k-inverse}
\end{equation}
where $r(z)$ is the $z$th standard normal quantile, $A$ and $B$ are location and scale parameters and $g$ and $k$ are related to skewness and kurtosis. Parameters restrictions are $B>0$ and $k>-0.5$.
An evaluation of \eqref{eq:g-k-inverse} returns a draw ($z$th quantile) from the $g$-and-$k$ distribution or, in other words, the $j$th sample $r_j:=r_j(z)\sim N(0,1)$ produces a draw $x_j:=F^{-1}(\cdot;A,B,c,g,k)$ from the $g$-and-$k$ distribution. However, unlike in previously mentioned references, we consider as data the vector $\mathbf{y}=(y_1,...,y_n)$, where $y_j=x_j+\varepsilon_i$, with i.i.d. noise $\varepsilon_j\sim N(0,\sigma^2_\varepsilon)$, where the $\varepsilon_j$'s are independent of the $x_j$'s, $j=1,...,n$. Also, denote $\mathbf{x}=(x_1,...,x_n)$. Notice that because SAEM-SL is an EM-type algorithm, and therefore it is suitable for ``incomplete data'', we would not be able to apply SAEM-SL to data observed directly as realizations from \eqref{eq:g-k-inverse}. That is while ABC methods can in principle accommodate inference based on either noisy data $\mathbf{y}$ and noise-free data $\mathbf{x}$, SAEM-SL can only deal with the former. We found the parameter $g$ to be of difficult identification and in the following we keep it fixed at its true value (see below): hence we assume $\boldsymbol{\theta}=(A,B,k,\sigma_\varepsilon)$ as parameter of interest, by noting that it is customary to keep $c$ fixed to $c=0.8$ (\citealp{drovandi2011likelihood,rayner2002numerical}).

 We initially consider the summaries $\tilde{\mathbf{S}}(\mathbf{x})=(S_{A,\mathbf{x}},S_{B,\mathbf{x}},S_{g,\mathbf{x}},S_{k,\mathbf{x}},P_{20,\mathbf{x}},P_{30,\mathbf{x}},P_{70,\mathbf{x}},P_{80,\mathbf{x}})$, where $P_{q,\mathbf{x}}$ is the $q$th empirical percentile of $\mathbf{x}$, whereas the remaining summaries are from \cite{drovandi2011likelihood}:
\begin{align*}
S_{A,\mathbf{x}}&=P_{50,\mathbf{x}} & S_{B,x}&=P_{75,\mathbf{x}}-P_{25,\mathbf{x}},\\
S_{g,\mathbf{x}}&=(P_{75,\mathbf{x}}+P_{25,\mathbf{x}}-2S_{A,\mathbf{x}})/S_{B,\mathbf{x}} & S_{k,\mathbf{x}}&=(P_{87.5,\mathbf{x}}-P_{62.5,\mathbf{x}}+P_{37.5,\mathbf{x}}-P_{12.5,\mathbf{x}})/S_{B,\mathbf{x}}.
\end{align*}
That is $S_{A,\mathbf{x}}$ and $S_{B,\mathbf{x}}$ are the median and the inter-quartile range of $\mathbf{x}$ respectively.
We define summaries $\tilde{\mathbf{S}}(\mathbf{y})$ for observed data in the analogous way as for $\mathbf{x}$, that is by plugging $\mathbf{y}$ in place of $\mathbf{x}$ in the summaries above. However we found that working with $\tilde{\mathbf{S}}(\mathbf{x})$ and $\tilde{\mathbf{S}}(\mathbf{y})$ produces unsatisfactory results, because the distributions of some of the simulated summaries are markedly asymmetric, i.e. far from being even approximately Gaussian. Therefore in practice we work with $\mathbf{S}(\mathbf{y}):=\log(\tilde{\mathbf{S}}(\mathbf{y})+\nu)$ and $\mathbf{S}(\mathbf{x}):=\log(\tilde{\mathbf{S}}(\mathbf{x})+\nu)$, where $\nu>0$ is a constant set so that the argument of the logarithms is strictly positive, and of course the same $\nu$ has to be used for $\mathbf{S}(\mathbf{x})$ and $\mathbf{S}(\mathbf{y})$ during the execution of SAEM-SL.
Therefore SAEM-SL is implemented with $\mathbf{s}=(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{x}))$.  For the specific data $\mathbf{y}$ simulated with the setting given below, $\nu=50$ was found to be appropriate.

Here we intend to compare SAEM-SL with an ABC algorithm. Therefore we produce a single dataset having length $n=500$, generated with $\boldsymbol{\theta}=(A,B,k,\sigma_\varepsilon)=(3,1,0.5,1)$ (we keep $c=0.8$ and $g=2$ fixed).
Starting values for SAEM-SL are $A=10$,          $B=10$,       $k= 4$ and      $\sigma_\varepsilon=0.3$. We run SAEM-SL with $R=3,000$, $K_1=10$ and $K=20$ and use $L=40$ iterations for the M-step. The result is in Figure \ref{fig:gk_saem_trajectories}. The simulation is relatively computer intensive, as computing the summaries (hence the percentiles) requires sorting procedures on each of the $R$ simulated data. Our SAEM-SL estimation required about 10 minutes of computation.

\begin{figure}
\centering
\includegraphics[width=17cm,height=8cm]{gk_saem_trajectories}
\caption{\footnotesize{g-and-k distribution: evolution of SAEM-SL. Dashed lines are the true parameter values.}}\label{fig:gk_saem_trajectories}
\end{figure}

\subsubsection{Bayesian estimation via ABC-MCMC}\label{sec:abc-mcmc}

Here we consider a comparison with the ``gold standard'' methodology for intractable likelihoods, that is approximate Bayesian computation (ABC). Several possible ABC methods could be considered: we choose an ABC-MCMC sampler, essentially a trivial modification of the one proposed in \cite{marjoram2003markov}, see for example \cite{sisson-fan(2011)}. As  shown in e.g. \cite{picchini2016approximate} it is possible to estimate parameters of  noise-free data from $g$-and-$k$ distributions using ABC-MCMC, and we now consider the case of noisy data. Briefly, with ABC the goal is to sample from an approximate posterior $\pi_{\delta}(\boldsymbol{\theta},\mathbf{z}|\mathbf{y})$ defined on the space of $\boldsymbol{\theta}$ augmented with the space of $\mathbf{z}$. Here $\mathbf{z}$ denotes synthetic observations defined on the same space as the actual observations $\mathbf{y}$, that is if $\mathbf{y}\in\mathcal{Y}$ are noisy observations then so are the $\mathbf{z}\in\mathcal{Y}$, and $\mathbf{z}$ should be simulated with the same generating model assumed for $\mathbf{y}$. However, typically in ABC studies a set of summary statistics is introduced to break the curse-of-dimensionality, and the resulting posterior is (by disregarding normalizing factors)
\begin{equation}
\pi_{\delta}(\boldsymbol{\theta},\mathbf{z}|\rho(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z}))\propto J_{\delta}(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z})){p(\mathbf{z}|\boldsymbol{\theta})\pi(\boldsymbol{\theta})},
\label{eq:abc-posterior}
\end{equation}
with $p(\mathbf{z}|\boldsymbol{\theta})$ the likelihood function based on $\mathbf{z}$ and $\pi(\theta)$ the prior for $\boldsymbol{\theta}$.
Here $\delta>0$ is a threshold value and $J_{\delta}(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z}))$ is a positive function assigning larger weights to values of $\boldsymbol{\theta}$ such that $\rho(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z}))<\delta$ for some appropriate distance $\rho(\cdot,\cdot)$. It can be shown that for a small enough $\delta$ the marginal ABC posterior $\pi_{\delta}(\boldsymbol{\theta}|\mathbf{S}(\mathbf{y}))=\int \pi_{\delta}(\boldsymbol{\theta},\mathbf{z}|\rho(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z}))) d\mathbf{z} $ is ``close'' to the true marginal $\pi(\boldsymbol{\theta}|\mathbf{S}(\mathbf{y}))$, if the summary statistics are informative for $\boldsymbol{\theta}$. Essentially, an ABC-MCMC algorithm produces a Markov chain for $\boldsymbol{\theta}$ having stationary distribution $\pi_{\delta}(\boldsymbol{\theta}|\mathbf{S}(\mathbf{y}))$. It should be remarked that in ABC $\mathbf{S}(\mathbf{y})$ and $\mathbf{S}(\mathbf{z})$ are the same set of summary functions, only applied to different arguments, as $\mathbf{z}$ and $\mathbf{y}$ are assumed to be defined on the same space and generated with the same underlying mechanism. For SAEM-SL the summaries we denoted with $\mathbf{S}(\mathbf{y})$ and $\mathbf{S}(\mathbf{x})$ in general \textit{do not} have to be the same functions, as $\mathbf{x}$ is a noise-free version of $\mathbf{y}$ hence these are defined on different spaces; however for this example we chose $\mathbf{S}(\mathbf{y})$ and $\mathbf{S}(\mathbf{x})$ to be the same set of functions.

To implement ABC-MCMC we choose a Gaussian kernel for $J_{\delta}(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z}))$, given by
\[
J_{\delta}(\mathbf{S}(\mathbf{y}),\mathbf{S}(\mathbf{z}))\propto \exp\{-(\mathbf{S}(\mathbf{z})-\mathbf{S}(\mathbf{y}))'\mathbf{\Omega}^{-1}(\mathbf{S}(\mathbf{z})-\mathbf{S}(\mathbf{y}))/2\delta^2\}
\]
where $'$ denotes transposition and $\mathbf{\Omega}$ is a positive definite matrix. For simplicity we assume a diagonal $\mathbf{\Omega}$ with elements $\mathbf{\Omega}=\mathrm{diag}\{\omega^{2}_1,...,\omega^{2}_{d_s}\}$, with $d_s=\dim \mathbf{S}(\mathbf{y})=\dim \mathbf{S}(\mathbf{z})$. When the elements in vector $\mathbf{S}(\mathbf{y})$ are  varying approximately on the same range of values it is possible to consider $(\omega^{2}_1,...,\omega^{2}_{d_s})=(1,...,1)$, however in general the variability of the statistics is unknown and, depending on the type of data and the underlying model, these can have very different magnitude.
The interested reader is referred to section 3.1 in \cite{picchini2016approximate} for further details (and disregarding the ``data cloning'' approach there exposed).

For $\mathbf{S}(\mathbf{y})$ and $\mathbf{S}(\mathbf{z})$ we consider the same set of summaries used with SAEM-SL and the same starting values for the parameters.
We run two attempts of an ABC-MCMC algorithm, with independent uniform priors $U(0,1)$ for $A$, $B$ and $k$ while we set $\sigma_\varepsilon\sim \Gamma(2,1)$, that is a Gamma distribution with mean 2. Parameters were proposed using an adaptive Metropolis algorithm with Gaussian innovations \citep{haario-saksman-tamminen}. At the first (pilot) attempt we use $(\omega^{2}_1,...,\omega^{2}_{d_s})=(1,...,1)$,
and let $\delta$ decrease every 20,000 iterations in $\delta\in\{0.03, 0.007, 0.003\}$, for a total of 60,000 iterations, where the $\delta$'s were chosen to target an acceptance rate of 1--3\% at the smallest $\delta$, usually considered a good compromise between accuracy and computational budget. Results were not encouraging, because the summaries vary on different scales but we assigned unit weight to each of them. However, we also collect the 20,000 summary statistics simulated at the smallest $\delta$, i.e. at $\delta=0.003$ and from these statistics we compute the median absolute deviation MAD for each coordinate of the accepted $\mathbf{S}(\mathbf{z})$ and define $(\omega_1,...,\omega_{d_s}):=(\mathrm{MAD}_1,...,\mathrm{MAD}_{d_s})$. We plug the new weights into $\mathbf{\Omega}$ for a further run of ABC-MCMC, this time using $\delta\in\{8, 3, 1, 0.3\}$, and the $\delta$'s had to be modified as a consequence of the different weights introduced, again targeting an acceptance rate of 1--3\% at the smallest $\delta$.
We use the parameter draws simulated in correspondence of $\delta=0.3$ to calculate the parameters posterior means and $95\%$ posterior intervals, and these result in: $\hat{A}=3.03$ [2.71,3.36], $\hat{B}=1.32$ [0.48,2.36], $k=0.51$ [0.05,1.45], $\sigma_\varepsilon=0.90$ [0.49,1.38].

Here the strength of ABC methods is on full display: ABC is not constrained by any parametric assumption regarding the distribution of the summaries, and when these are informative ABC is probably the go-to choice.
The essence of the comparison is that tuning ABC algorithms is not trivial. However, for each iteration of ABC-MCMC we only need to simulate a single realization of $\mathbf{z}$, while for each iteration of SAEM-SL we need at least $L\times R$ simulations from the model. However, a proper comparison between SAEM-SL and ABC is not problem independent. For example, in stochastic dynamical modelling an ABC-MCMC sampler will seldom produce accurate results and an ABC-SMC approach will usually be preferred (see e.g. \citealp{toni2009approximate}), this increasing the computational effort considerably.

\section{Summary}\label{sec:summary}

We have introduced a new method for approximate maximum likelihood estimation of the parameters of partially observed models, or equivalently for models having so-called incomplete data, following the terminology given in \cite{Dempster1977}. Under this framework, our method is able to deal with a large class of modelling scenarios, and both ``static'' (example in section \ref{sec:g-and-k}) and ``dynamic'' models (examples in sections \ref{sec:nonlingauss}--\ref{sec:theoph}) can be accommodated. We first illustrate the stochastic approximation EM algorithm (SAEM, \citealp{Delyon1999}) as one of the possible ways to implement an EM algorithm. To fully exploit the computational benefits of SAEM, namely convergence to a (local) maximizer of the data likelihood, it is required to analytically compute the complete likelihood of the model and derive the corresponding sufficient statistics. The latter step is far from being trivial (if at all possible) for most models of realistic complexity. Our SAEM-SL method makes use of the synthetic likelihoods (SL) approach proposed in \cite{wood2010statistical}: SL requires from the modeller the specification of ``appropriate'' (informative) summary statistics encoding the information about the parameter $\boldsymbol{\theta}$ that is contained in the available data. These summaries are assumed to follow a Gaussian distribution and we find that this assumption is convenient for exploitation in a SAEM context, as Gaussian likelihoods have trivial to compute sufficient statistics, which we obtain from SL simulations. Our approach constructs a version of SL for the ``complete synthetic loglikelihood'' and plugs it within SAEM. As a result, it bypasses the analytic calculation of the complete likelihood of the model by introducing a Gaussian approximation. SAEM-SL results in a plug-and-play, likelihood-free, approximated version of SAEM.
Under ideal scenarios, where the user-specified summaries are sufficient statistics for $\boldsymbol{\theta}$ \textit{and} are also Gaussian distributed, then SAEM-SL is equivalent to the standard SAEM and therefore should return a  stationary point of the true data-likelihood.

In three simulation studies we have shown the good performance of the method, which requires minimal tuning. From the computational complexity point of view the disadvantage of using SAEM-SL is that the M-step (maximization step) is performed using a numerical maximizer. This increases the computational time considerably, when compared to the standard SAEM, as with SAEM it is sometimes possible to derive analytically the maximizers, provided that the often challenging analytic work of deriving the sufficient statistics has been performed with success. Further research effort should be devoted to construct a more efficient optimizer for the M-step.
Moreover SAEM-SL requires from the modeller a set of summary statistics: this operation is clearly subjective and delicate. A possibility to automatize the process of selection of the statistics is to run a semi-automatic summaries selection algorithm as described in \cite{fearnhead-prangle(2011)}, from within an approximate Bayesian computation framework, then plug the constructed summaries into SAEM-SL. We have not considered the possibility to use the semi-automatic selection approach in the present work, and a study of the implications is left for future research. In conclusion SAEM-SL is an appealing application of SAEM for intractable models. We have performed several comparisons with well established methodologies (such as iterated filtering, particle marginal methods, approximate Bayesian computation and SAEM incorporating a sequential Monte Carlo step) and while SAEM-SL performs satisfactorily against available alternatives, in challenging settings the ``best'' approach for a specific problem is often a compromise between computational feasibility and statistical efficiency.

\section*{Acknowledgements}
Research was partially funded by the Swedish Research Council (VR grant 2013-5167). We thank Christopher Drovandi (Queensland University of Technology) for valuable discussion and for suggesting the Cholesky factorization check when sampling from a multivariate Gaussian.