###### Abstract

In this article we consider the development of unbiased estimators of the
Hessian, of the log-likelihood function with respect to parameters, for
partially observed diffusion processes. These processes arise in numerous
applications, where such diffusions require derivative information, either
through the Jacobian or Hessian matrix. As time-discretizations of diffusions
induce a bias, we provide an unbiased estimator of the Hessian. This is based
on using Girsanov’s Theorem and randomization schemes developed through
Mcleish [25] and Rhee & Glynn [27]. We demonstrate our developed estimator of
the Hessian is unbiased, and one of finite variance. We numerically test and
verify this by comparing the methodology here to that of a newly proposed
particle filtering methodology. We test this on a range of diffusion models,
which include different Ornstein–Uhlenbeck processes and the Fitzhugh–Nagumo
model, arising in neuroscience.

Key words: Partially Observed Diffusions, Randomization Methods, Hessian
Estimation,
Coupled Conditional Particle Filter
AMS subject classifications: 62C10, 60J60, 60J22, 65C40

Unbiased Estimation of the Hessian for Partially Observed Diffusions

BY NEIL K. CHADA, AJAY JASRA & FANGYUAN YU

Computer, Electrical and Mathematical Sciences and Engineering Division, King
Abdullah University of Science and Technology, Thuwal, 23955-6900, KSA.
E-Mail:  _neilchada123<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

In many scientific disciplines, diffusion processes [10] are used to model and
describe important phenomenon. Particular applications where such processes
arise include biological sciences, finance, signal processing and atmospheric
sciences [4, 24, 28, 30]. Mathematically, diffusion processes take the general
form

$dX_{t}=a_{\theta}(X_{t})dt+\sigma(X_{t})dW_{t},\quad
X_{0}=x_{\star}\in\mathbb{R}^{d},$ (1.1)

where $X_{t}\in\mathbb{R}^{d}$, $\theta\in\Theta$ is a parameter,
$X_{0}=x_{\star}$ is the initial condition with $x_{\star}$ given,
$a:\Theta\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ denotes the drift
term, $\sigma:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ denotes the
diffusion coefficient and $\\{W_{t}\\}_{t\geq 0}$ is a standard
$d-$dimensional Brownian motion. In practice it is often difficult to have
direct access to such continuous processes, where instead one has discrete-
time partial observations of the process $\\{X_{t}\\}_{t\geq 0}$, denoted as
$Y_{t_{1}},\ldots,Y_{t_{n}}$, where $0<t_{1}<\ldots<t_{n}=T$, such that
$Y_{t_{p}}\in\mathbb{R}^{d_{y}}$. Such processes are referred to as partially
observed diffusion processes (PODPs), where one is interested in doing
inference on the hidden process (1.1) given the observations. In order to do
such inference, one must time-discretize such a process which induces a
discretization bias. For (1.1) this can arise through common discretization
forms such as an Euler or Milstein scheme [23]. Therefore an important
question, related to inference, is how one can reduce, or remove the
discretization bias. Such a discussion motivates the development and
implementation of unbiased estimators.

The unbiased estimation of PODPs has been an important, yet challenging topic.
Some original seminal work on this has been the idea of exact simulation,
proposed in various works [5, 6, 13]. The underlying idea behind exact
simulation, is that through a particular transformation one can acquire an
unbiased estimator, subject to certain restrictions on the from of the
diffusion and its dimension. Since then there has been a number of extensions
aimed at going beyond this, w.r.t. to more general multidimensional diffusions
and continuous-time dynamics [8, 14]. However there has been recent attention
in unbiased estimation, for Bayesian computation, through the work of Rhee and
Glynn [17, 27], where they provide unbiased and finite variance estimators
through introducing randomization. In particular these methods allow to
unbiasedly estimate an expectation of a functional, by randomizing on the
level of the time-discretization in a type of multilevel Monte Carlo (MLMC)
approach [32], where there is a coupling between different levels. As a
result, this methodology has been considered in the context of both filtering
and Bayesian computation [11, 19, 21] and gradient estimation [18].

In this work we are interested in developing an unbiased estimator of the
Hessian for PODPs. This is of interest as current state of-the-art stochastic
gradient methodologies, exploit Hessian information for improved convergence,
such as Newton type methods [1, 9]. In order to develop an unbiased estimator,
our methodology will largely follow that described in [18], with the extension
of this from the score function to the Hessian. In particular we will exploit
the use of the conditional particle filter (CPF), first considered by Andrieu
et al. [2, 3]. We provide an expression for the Hessian of the likelihood,
while introducing an Euler time-discretization of the diffusion process in
order to implement our unbiased estimator. We then describe how one can attain
unbiased estimators, which is based on various couplings of the CPF. From this
we test this methodology to that of using the methods of [11, 22] for the
Hessian computation, as for a comparison, where we demonstrate the unbiased
estimator through both the variance and bias. This will be conducted on both a
single and multidimensional Ornstein–Uhlenbeck process, as well as a more
complicated model of the Fitzhugh–Nagumo model. We remark that our estimator
of the hessian is unbiased, but if the inverse hessian is required, it is
possible to adapt the forthcoming methodology to that context as well.

###  1.1 Outline

In Section 2 we present our setting for our diffusion process. We also present
a derived expression for the Hessian, with an appropriate time-discretization.
Then in Section 3 we describe our algorithm in detail for the unbiased
estimator of the Hessian. This will be primarily based on a coupling of a
coupled conditional particle filter. This will lead to Section 4 where we
present our numerical experiments, which provide variance and bias plots. We
compare the methodology of this work, with that of the Delta particle filter.
This comparison will be tested on a range of diffusion processes, which
include an Ornstein–Uhlenbeck process and the Fitzhugh–Nagumo model. We
summarize our findings in Section 5.

##  2 Model

In this section we introduce our setting and notation regarding our partially
observed diffusions. This will include a number of assumptions. We will then
provide an expression for the Hessian of the likelihood function, with a time-
discretization based on the Euler scheme. This will include a discussion on
the stochastic model where we define the marginal likelihood. Finally we
present a result indicating the approximation of the Hessian computation as we
take the limit of the discretization level.

###  2.1 Notation

Let $(\mathsf{X},\mathcal{X})$ be a measurable space. For
$\varphi:\mathsf{X}\rightarrow\mathbb{R}$ we write
$\mathcal{B}_{b}(\mathsf{X})$ as the collection of bounded measurable
functions, $\mathcal{C}^{j}(\mathsf{X})$ are the collection of $j-$times,
$j\in\mathbb{N}$ continuously differentiable functions and we omit the
subscript $j$ if the functions are simply continuous; if
$\varphi:\mathsf{X}\rightarrow\mathbb{R}^{d}$ we write
$\mathcal{C}_{d}^{j}(\mathsf{X})$ and $\mathcal{C}_{d}(\mathsf{X})$. Let
$\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}$,
$\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$ denotes the collection of real-
valued functions that are Lipschitz w.r.t. $\|\cdot\|_{2}$ ($\|\cdot\|_{p}$
denotes the $\mathbb{L}_{p}-$norm of a vector $x\in\mathbb{R}^{d}$). That is,
$\varphi\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$ if there exists a
$C<+\infty$ such that for any $(x,y)\in\mathbb{R}^{2d}$

$|\varphi(x)-\varphi(y)|\leq C\|x-y\|_{2}.$

We write $\|\varphi\|_{\textrm{Lip}}$ as the Lipschitz constant of a function
$\varphi\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$. For
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$, we write the supremum norm
$\|\varphi\|=\sup_{x\in\mathsf{X}}|\varphi(x)|$. $\mathcal{P}(\mathsf{X})$
denotes the collection of probability measures on $(\mathsf{X},\mathcal{X})$.
For a measure $\mu$ on $(\mathsf{X},\mathcal{X})$ and a
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$, the notation
$\mu(\varphi)=\int_{\mathsf{X}}\varphi(x)\mu(dx)$ is used. $B(\mathbb{R}^{d})$
denote the Borel sets on $\mathbb{R}^{d}$. $dx$ is used to denote the Lebesgue
measure. Let $K:\mathsf{X}\times\mathcal{X}\rightarrow[0,\infty)$ be a non-
negative operator and $\mu$ be a measure then we use the notations $\mu
K(dy)=\int_{\mathsf{X}}\mu(dx)K(x,dy)$ and for
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$,
$K(\varphi)(x)=\int_{\mathsf{X}}\varphi(y)K(x,dy).$ For $A\in\mathcal{X}$ the
indicator is written $\mathbb{I}_{A}(x)$. $\mathcal{U}_{A}$ denotes the
uniform distribution on the set $A$. $\mathcal{N}_{s}(\mu,\Sigma)$ (resp.
$\psi_{s}(x;\mu,\Sigma)$) denotes an $s-$dimensional Gaussian distribution
(density evaluated at $x\in\mathbb{R}^{s}$) of mean $\mu$ and covariance
$\Sigma$. If $s=1$ we omit the subscript $s$. For a vector/matrix $X$, $X^{*}$
is used to denote the transpose of $X$. For $A\in\mathcal{X}$,
$\delta_{A}(du)$ denotes the Dirac measure of $A$, and if $A=\\{x\\}$ with
$x\in\mathsf{X}$, we write $\delta_{x}(du)$. For a vector-valued function in
$d-$dimensions (resp. $d-$dimensional vector), $\varphi(x)$ (resp. $x$) say,
we write the $i^{\textrm{th}}-$component ($i\in\\{1,\dots,d\\}$) as
$\varphi(x)^{(i)}$ (resp. $x^{(i)}$). For a $d\times q$ matrix $x$ we write
the $(i,j)^{\textrm{th}}-$entry as $x^{(ij)}$. For
$\mu\in\mathcal{P}(\mathsf{X})$ and $X$ a random variable on $\mathsf{X}$ with
distribution associated to $\mu$ we use the notation $X\sim\mu(\cdot)$.

###  2.2 Diffusion Process

Let $\theta\in\Theta\subseteq\mathbb{R}^{d_{\theta}}$ be fixed and we consider
a diffusion process on the probability space
$(\Omega,\mathscr{F},\\{\mathscr{F}\\}_{t\geq 0},\mathbb{P}_{\theta})$, such
that

$dX_{t}=a_{\theta}(X_{t})dt+\sigma(X_{t})dW_{t},\quad
X_{0}=x_{\star}\in\mathbb{R}^{d},$ (2.1)

where $X_{t}\in\mathbb{R}^{d}$, $X_{0}=x_{\star}$ with $x_{\star}$ given,
$a:\Theta\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the drift term,
$\sigma:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ is the diffusion
coefficient and $\\{W_{t}\\}_{t\geq 0}$ is a standard $d-$dimensional Brownian
motion. We assume that for any fixed $\theta\in\Theta$,
$a_{\theta}^{(i)}\in\mathcal{C}^{2}(\mathbb{R}^{d})$ and
$\sigma^{(ij)}\in\mathcal{C}^{2}(\mathbb{R}^{d})$ for
$(i,j)\in\\{1,\ldots,d\\}^{2}$. For fixed $x\in\mathbb{R}^{d}$ we have
$a_{\theta}(x)^{(i)}\in\mathcal{C}(\Theta)$ for $i\in\\{1,\ldots,d\\}$.

Furthermore we make the following additional assumption, termed (D1).

>   1. 1.
>
> Uniform ellipticity: $\Sigma(x):=\sigma(x)\sigma(x)^{*}$ is uniformly
> positive definite over $x\in\mathbb{R}^{d}$.
>
>   2. 2.
>
> Globally Lipschitz: for any $\theta\in\Theta$, there exists a positive
> constant $C<\infty$ such that
>
>
> $|a_{\theta}(x)^{(i)}-a_{\theta}(x^{\prime})^{(i)}|+|\sigma(x)^{(ij)}-\sigma(x^{\prime})^{(ij)}|\leq
> C\|x-x^{\prime}\|_{2},$
>
> for all $(x,x^{\prime})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$,
> $(i,j)\in\\{1,\ldots,d\\}^{2}$.
>
>

Let $0<t_{1},<\ldots<t_{n}=T>0$ be a given collection of time points.
Following [18], by the use of Girsanov Theorem, for any
$\mathbb{P}_{\theta}$-integrable
$\varphi:\Theta\times\mathbb{R}^{nd}\rightarrow\mathbb{R}$,

$\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}},\ldots,X_{t_{n}})]=\mathbb{E}_{\mathbb{Q}}\bigg{[}\varphi_{\theta}(X_{t_{1}},\ldots,X_{t_{n}})\frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(\mathbf{X}_{T})\bigg{]},$
(2.2)

where $\mathbb{E}_{\theta}$ denotes the expectation w.r.t.
$\mathbb{P}_{\theta}$, set $\mathbf{X}_{T}=\\{X_{t}\\}_{t\in[0,T]}$, and the
change of measure is given by

$\frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(\mathbf{X}_{T})=\exp\bigg{\\{}-\frac{1}{2}\int^{T}_{0}\|b_{\theta}(X_{s})\|^{2}_{2}ds+\int^{T}_{0}b_{\theta}(X_{s})^{*}dW_{s}\bigg{\\}},$

with $b_{\theta}(x)=\Sigma(x)^{-1}\sigma(x)^{*}a_{\theta}(x)$ is a $d-$vector.
As it will be useful below, we can modify the above expression, by using that
$dX_{t}=\sigma(X_{t})dW_{t}$, where $X_{t}$ solves such a process,

$\frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(\mathbf{X}_{T})=\exp\bigg{\\{}-\frac{1}{2}\int^{T}_{0}\|b_{\theta}(X_{s})\|^{2}_{2}ds+\int^{T}_{0}b_{\theta}(X_{s})^{*}\Sigma(X_{s})^{-}1\sigma(X_{s})^{*}dX_{s}\bigg{\\}}.$

Set

$\rho_{\theta}(\mathbf{X}_{T})=\varphi_{\theta}(X_{t_{1}:t_{n}})\frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(\mathbf{X}_{T}).$

Now if we assume that $\varphi_{\theta}$ is differentiable w.r.t. $\theta$,
then one has for $i\in\\{1,\dots,d_{\theta}\\}$

$\mathfrak{G}_{\theta}^{(i)}:=\frac{\partial}{\partial\theta^{(i)}}\Big{(}\log\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}},\ldots,X_{t_{n}})]\right\\}\Big{)}=\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\bigg{[}\frac{\partial}{\partial\theta^{(i)}}\Big{(}\log\left\\{\rho_{\theta}(\mathbf{X}_{T})\right\\}\Big{)}\bigg{]},$
(2.3)

where
$\overline{\mathbb{P}}_{\theta}=\varphi_{\theta}\mathbb{P}_{\theta}/\mathbb{P}_{\theta}(\varphi_{\theta})$
and $\mathbb{P}_{\theta}[\varphi_{\theta}(X_{t_{1}},\ldots,X_{t_{n}})]$. From
herein we will use the short-hand notation
$\varphi_{\theta}(X_{t_{1}},\ldots,X_{t_{n}})=\varphi_{\theta}(X_{t_{1}:t_{n}})$
and also set, for $i\in\\{1,\dots,d_{\theta}\\}$,

$G_{\theta}(\mathbf{X}_{T})^{(i)}=\frac{\partial}{\partial\theta^{(i)}}\Big{(}\log\left\\{\rho_{\theta}(\mathbf{X}_{T})\right\\}\Big{)}.$

###  2.3 Hessian Expression

Given the expression (2.3) our objective is now to write the matrix of second
derivatives, for $(i,j)\in\\{1,\dots,d_{\theta}\\}^{2}$

$\mathfrak{H}_{\theta}^{(ij)}:=-\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]\right\\}\Big{)},$

in terms of expectations w.r.t. $\overline{\mathbb{P}}_{\theta}$.

We have the following simple calculation

$\displaystyle\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]\right\\}\Big{)}$
$\displaystyle=\frac{\partial}{\partial\theta^{(i)}}\Big{(}\frac{\frac{\partial}{\partial\theta^{(j)}}\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]\right\\}}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]}\Big{)}$
$\displaystyle=\frac{\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]\right\\}}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]}-\frac{\frac{\partial}{\partial\theta^{(i)}}\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]\right\\}\frac{\partial}{\partial\theta^{(j)}}\left\\{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]\right\\}}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{t_{1}:t_{n}})]^{2}}$
$\displaystyle=:T_{1}-T_{2}.$

Under relatively weak conditions, once can express $T_{1}$ and $T_{2}$ as

$\displaystyle T_{1}$
$\displaystyle=\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\Big{[}\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(i)}}\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(j)}}\Big{]}+\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}\Big{)}\right],$
$\displaystyle T_{2}$
$\displaystyle=\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(i)}}\right]\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(j)}}\right].$

Therefore we have the following expression

$\displaystyle\mathfrak{H}_{\theta}^{(ij)}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(i)}}\right]\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(j)}}\right]-\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\Big{[}\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(i)}}\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(j)}}\Big{]}-$
(2.4)
$\displaystyle\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}\Big{)}\right].$

Defining, for $(i,j)\in\\{1,\dots,d_{\theta}\\}^{2}$

$H_{\theta}(\mathbf{X}_{T})^{(ij)}:=\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}\Big{)},$

one can write more succinctly

$\mathfrak{H}_{\theta}^{(ij)}=\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[G_{\theta}(\mathbf{X}_{T})^{(i)}\right]\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[G_{\theta}(\mathbf{X}_{T})^{(j)}\right]-\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[G_{\theta}(\mathbf{X}_{T})^{(i)}G_{\theta}(\mathbf{X}_{T})^{(j)}\right]-\mathbb{E}_{\overline{\mathbb{P}}_{\theta}}\left[H_{\theta}(\mathbf{X}_{T})^{(ij)}\right].$

####  2.3.1 Stochastic Model

Consider a sequence of random variables $(Y_{t_{1}},\ldots,Y_{t_{n}})$ where
$0<t_{1}<\ldots<t_{n}=T$, where $Y_{t_{p}}\in\mathbb{R}^{d_{y}}$, which are
assumed to have the following joint Lebesgue density

$p_{\theta}(y_{t_{1}},\ldots,y_{t_{n}}|\\{x_{s}\\}_{0\leq s\leq
T})=\prod^{n}_{p=1}g_{\theta}(y_{t_{p}}|x_{t_{p}}),$

where
$g:\Theta\times\mathbb{R}^{d}\times\mathbb{R}^{d_{y}}\rightarrow\mathbb{R}^{+}$
for any
$(\theta,x)\in\Theta\times\mathbb{R}^{d},\int_{\mathbb{R}^{d_{y}}}g_{\theta}(y|x)dy=1$
such that $dy$ is the Lebesgue measure. Now if one considers instead
realizations of the random variables $(Y_{t_{1}},\ldots,Y_{t_{n}})$, then we
have a state-space model with marginal likelihood

$p_{\theta}(y_{t_{1}},\ldots,y_{t_{n}}):=\mathbb{E}_{\theta}\bigg{[}\prod^{n}_{p=1}g_{\theta}(y_{t_{p}}|X_{t_{p}})\bigg{]}.$

Note that the framework to be investigated in this article is not restricted
to this special case, but, we shall focus on it for the rest of the paper. So
to clarify
$\varphi_{\theta}(x_{t_{1}},\ldots,x_{t_{n}})=\prod^{n}_{p=1}g_{\theta}(y_{t_{p}}|x_{t_{p}})$
from herein.

In reference to (2.3) and (2.4), we have that

$\displaystyle\frac{\partial\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}}{\partial\theta^{(i)}}$
$\displaystyle=$
$\displaystyle\sum^{n}_{p=1}\frac{\partial}{\partial\theta^{(i)}}\Big{(}\log\\{g_{\theta}(y_{t_{p}}|x_{t_{p}})\\}\Big{)}-\frac{1}{2}\int^{T}_{0}\frac{\partial}{\partial\theta^{(i)}}\Big{(}\|b_{\theta}(X_{s})\|^{2}_{2}\Big{)}ds+$
$\displaystyle\frac{\partial}{\partial\theta^{(i)}}\Big{(}\int^{T}_{0}b_{\theta}(X_{s})^{*}\Sigma(X_{s})^{-1}\sigma(X_{s})^{*}dX_{s}\Big{)}$
$\displaystyle\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\\{\rho_{\theta}(\mathbf{X}_{T})\\}\Big{)}$
$\displaystyle=$
$\displaystyle\sum^{n}_{p=1}\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\\{g_{\theta}(y_{t_{p}}|x_{t_{p}})\\}\Big{)}-\frac{1}{2}\int^{T}_{0}\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\|b_{\theta}(X_{s})\|^{2}_{2}\Big{)}ds+$
$\displaystyle\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\int^{T}_{0}b_{\theta}(X_{s})^{*}\Sigma(X_{s})^{-1}\sigma(X_{s})^{*}dX_{s}\Big{)}.$

###  2.4 Time-Discretization

From herein, we take the simplification that
$t_{p}=p,p\in\\{1,\ldots,n\\},T=n$. Let $l\in\mathbb{N}_{0}$ be given and
consider the Euler discretization of step size
$\Delta_{l}=2^{-l},k\in\\{1,2,\ldots,\Delta^{-1}_{l}T\\}$ with
$\tilde{X}_{0}=x_{\star}$:

$\tilde{X}_{k\Delta_{l}}=\tilde{X}_{(k-1)\Delta_{l}}+a_{\theta}(\tilde{X}_{(k-1)\Delta_{l}})\Delta_{l}+\sigma(\tilde{X}_{(k-1)\Delta_{l}})[W_{k\Delta
l}-W_{(k-1)\Delta_{l}}].$ (2.5)

Set
$\mathbf{x}_{T}^{l}=(x_{\star},\tilde{x}_{\Delta_{l}},\ldots,\tilde{x}_{T})$.
We then consider the vector-valued function
$G^{l}:\Theta\times(\mathbb{R}^{d})^{\Delta_{l}^{-1}T+1}\rightarrow\mathbb{R}^{d_{\theta}}$
and the matrix-valued function
$H^{l}:\Theta\times(\mathbb{R}^{d})^{\Delta_{l}^{-1}T+1}\rightarrow\mathbb{R}^{d_{\theta}\times
d_{\theta}}$ defined as, for $(i,j)\in\\{1,\dots,d_{\theta}\\}^{2}$

$\displaystyle G_{\theta}^{l}(\mathbf{x}_{T}^{l})^{(i)}$ $\displaystyle=$
$\displaystyle\sum^{n}_{p=1}\frac{\partial}{\partial\theta^{(i)}}\Big{(}\log\\{g_{\theta}(y_{p}|\tilde{x}_{p})\\}\Big{)}-\frac{\Delta_{l}}{2}\sum^{\Delta_{l}^{-1}T-1}_{k=0}\frac{\partial}{\partial\theta^{(i)}}\Big{(}\|b_{\theta}(\tilde{x}_{k\Delta_{l}})\|^{2}_{2}\Big{)}+$
$\displaystyle\sum^{\Delta_{l}^{-1}T-1}_{k=0}\frac{\partial}{\partial\theta^{(i)}}\Big{(}b_{\theta}(\tilde{x}_{k\Delta_{l}})^{*}\Sigma(\tilde{x}_{k\Delta_{l}})^{-1}\sigma(\tilde{x}_{k\Delta_{l}})^{*}[\tilde{x}_{(k+1)\Delta_{l}}-\tilde{x}_{k\Delta_{l}}]\Big{)},$
$\displaystyle H_{\theta}^{l}(\mathbf{x}_{T}^{l})^{(ij)}$ $\displaystyle=$
$\displaystyle\sum^{n}_{p=1}\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\log\\{g_{\theta}(y_{p}|\tilde{x}_{p})\\}\Big{)}-\frac{\Delta_{l}}{2}\sum^{\Delta_{l}^{-1}T-1}_{k=0}\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}\|b_{\theta}(\tilde{x}_{k\Delta_{l}})\|^{2}_{2}\Big{)}+$
$\displaystyle\sum^{\Delta_{l}^{-1}T-1}_{k=0}\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(j)}}\Big{(}b_{\theta}(\tilde{x}_{k\Delta_{l}})^{*}\Sigma(\tilde{x}_{k\Delta_{l}})^{-1}\sigma(\tilde{x}_{k\Delta_{l}})^{*}[\tilde{x}_{(k+1)\Delta_{l}}-\tilde{x}_{k\Delta_{l}}]\Big{)}.$

Then, noting (2.4), we have an Euler approximation of the Hessian

$\displaystyle\mathfrak{H}_{\theta}^{l,(ij)}$ $\displaystyle:=$
$\displaystyle\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(i)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(i)}G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-$
$\displaystyle\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})H_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(ij)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}.$

In the context of the model in Section 2.3.1, if one sets

$\pi_{\theta}^{l}(d\mathbf{x}_{T}^{l})\propto\left\\{\prod^{n}_{p=1}g_{\theta}(y_{p}|\tilde{x}_{p})p_{\theta}^{l}(\tilde{x}_{p-1},\tilde{x}_{p})\right\\}d\mathbf{x}_{T}^{l},$

where $p_{\theta}^{l}$ is the transition density induced by (2.5) (over unit
time) and we use the abuse of notation that $d\mathbf{x}_{T}^{l}$ is the
Lebesgue measure on the co-ordinates
$(\tilde{x}_{\Delta_{l}},\ldots,\tilde{x}_{T})$, then one has that

$\mathfrak{H}_{\theta}^{l,(ij)}=\pi_{\theta}^{l}(G_{\theta}^{l,(i)})\pi_{\theta}^{l}(G_{\theta}^{l,(j)})-\pi_{\theta}^{l}(G_{\theta}^{l,(i)}G_{\theta}^{l,(j)})-\pi_{\theta}^{l}(H_{\theta}^{l,(ij)}),$
(2.6)

where we are using the short-hand
$G_{\theta}^{l,(i)}=G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(i)}$ and
$H_{\theta}^{l,(ij)}=H_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(ij)}$ etc.

We have the following result whose proof and assumption (D2) is in Appendix A.

######  Proposition 2.1.

Assume (D1-D2). Then for any $(i,j)\in\\{1,\dots,d_{\theta}\\}^{2}$ we have

$\lim_{l\rightarrow\infty}\mathfrak{H}_{\theta}^{l,(ij)}=\mathfrak{H}_{\theta}^{(ij)}.$

The main strategy of the proof is by strong convergence, which means that one
can characterize an upper-bound on
$|\mathfrak{H}_{\theta}^{l,(ij)}-\mathfrak{H}_{\theta}^{(ij)}|$ of
$\mathcal{O}(\Delta_{l}^{1/2})$ but that rate is most likely not sharp, as one
expects $\mathcal{O}(\Delta_{l})$.

##  3 Algorithm

The objective of this Section is, using only approximations of (2.6), is to
obtain an unbiased estimate of $\mathfrak{H}_{\theta}^{(ij)}$ for any fixed
$\theta\in\Theta$ and $(i,j)\in\\{1,\dots,d_{\theta}\\}^{2}$. Our approach is
essentially an application of the methodology in [18] and so we provide a
review of that approach in the sequel.

###  3.1 Strategy

To focus our description, we shall suppose that we are interested in computing
an unbiased estimate of $\mathfrak{G}_{\theta}^{(i)}$ for some fixed $i$; we
remark that this specialization is not needed and is only used for notational
convenience. An Euler approximation of $\mathfrak{G}_{\theta}^{(i)}$ is
$\pi_{\theta}^{l}(G_{\theta}^{l,(i)})=:\mathfrak{G}_{\theta}^{l,(i)}$. To
further simplify the notation we will simply write $G_{\theta}^{l}$ instead of
$G_{\theta}^{l,(i)}$.

Suppose that one can construct a sequence of random variables
$(\widehat{\pi}_{\theta}^{l}(G_{\theta}^{l}))_{l\in\mathbb{N}_{0}}$ on a
potentially extended probability space with expectation operator
$\overline{\mathbb{E}}_{\theta}$, such that for each $l\in\mathbb{N}_{0}$,
$\overline{\mathbb{E}}_{\theta}[\widehat{\pi}_{\theta}^{l}(G_{\theta}^{l})]=\pi_{\theta}^{l}(G_{\theta}^{l})$.
Moreover, consider the independent sequence of random variables,
$(\Xi_{\theta}^{l})_{l\in\mathbb{N}_{0}}$ which are constructed so that for
$l\in\mathbb{N}_{0}$

$\overline{\mathbb{E}}_{\theta}[\Xi_{\theta}^{l}]:=\overline{\mathbb{E}}_{\theta}[\widehat{\pi}_{\theta}^{l}(G_{\theta}^{l})]-\overline{\mathbb{E}}_{\theta}[\widehat{\pi}_{\theta}^{l-1}(G_{\theta}^{l-1})]=\pi_{\theta}^{l}(G_{\theta}^{l})-\pi_{\theta}^{l-1}(G_{\theta}^{l-1}),$
(3.1)

with
$\overline{\mathbb{E}}_{\theta}[\widehat{\pi}_{\theta}^{-1}(G_{\theta}^{-1})]:=\pi_{\theta}^{-1}(G_{\theta}^{-1}):=0$.
Now let $\mathbb{P}_{L}$ be a positive probability mass function on
$\mathbb{N}_{0}$ and set
$\overline{\mathbb{P}}_{L}(l)=\sum_{p=l}^{\infty}\mathbb{P}_{L}(p)$. Now if,

$\sum_{l\in\mathbb{N}_{0}}\frac{1}{\overline{\mathbb{P}}_{L}(l)}\Big{\\{}\overline{\mathbb{V}\textrm{ar}}_{\theta}[\Xi_{\theta}^{l}]+(\mathfrak{G}_{\theta}^{l,(i)}-\mathfrak{G}_{\theta}^{(i)})^{2}\Big{\\}}<+\infty,$
(3.2)

then if one samples $L$ from $\mathbb{P}_{L}$ independently of the sequence
$(\Xi_{l})_{l\in\mathbb{N}_{0}}$ then by e.g. [32, Theorem 5] the estimate

$\widehat{\mathfrak{G}}_{\theta}^{(i)}:=\sum_{l=0}^{L}\frac{\Xi_{\theta}^{l}}{\overline{\mathbb{P}}_{L}(l)},$
(3.3)

is an unbiased and finite variance estimator of $\mathfrak{G}_{\theta}^{(i)}$.
The main issue is to construct the sequence of independent random variables
$(\Xi_{\theta}^{l})_{l\in\mathbb{N}_{0}}$ such that (3.1) and (3.2) hold and
that the expected computational cost for doing so is not unreasonable as a
functional of $\Delta_{l}$: a method for doing this is in [18] as we will now
describe.

###  3.2 Computing $\Xi_{\theta}^{0}$

The computation of $\Xi_{\theta}^{0}$ is performed by using exactly the
coupled conditional particle filter (CCPF) that has been introduced in [20].
This is an algorithm which allows one to construct a random variable
$\widehat{\pi}_{\theta}^{0}(G_{\theta}^{0})$ such that
$\overline{\mathbb{E}}_{\theta}[\widehat{\pi}_{\theta}^{0}(G_{\theta}^{0})]=\pi_{\theta}^{0}(G_{\theta}^{0})$
and we will set $\Xi_{\theta}^{0}=\widehat{\pi}_{\theta}^{0}(G_{\theta}^{0})$.

  1. 1.

Input $x_{\Delta_{l}:n}^{\prime}\in\mathsf{X}^{l}$. Set $k=1$,
$x_{0}^{i}=x_{\star}$, $a_{0}^{i}=i$ for $i\in\\{1,\dots,N-1\\}$.

  2. 2.

Sampling: for $i\in\\{1,\dots,N-1\\}$ sample
$x_{k-1+\Delta_{l}:k}^{i}|x_{k-1}^{a_{k-1}^{i}}$ using the Markov kernel
$p_{\theta}^{l}$. Set $x_{k-1+\Delta_{l}:k}^{N}=x_{k}^{\prime}$ and for
$i\in\\{1,\dots,N-1\\}$,
$x_{\Delta_{l}:k}^{i}=(x_{\Delta_{l}:k-1}^{a_{k-1}^{i}},x_{k-1+\Delta_{l}:k}^{i})$.
If $k=n$ go to 4..

  3. 3.

Resampling: Construct the probability mass function on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{k}|x_{k}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|x_{k}^{j})}.$

For $i\in\\{1,\dots,N-1\\}$ sample $a_{k}^{i}$ from $r_{1}^{i}$. Set $k=k+1$
and return to the start of 2..

  4. 4.

Construct the probability mass function on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{n}|x_{n}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|x_{n}^{j})}.$

Sample $i\in\\{1,\dots,N\\}$ using this mass function and return
$x_{\Delta_{l}:n}^{i}$.

Algorithm 1 Conditional Particle Filter at level $l\in\mathbb{N}_{0}$.

We begin by introducing the Markov kernel
$C^{l}:\Theta\times\mathsf{X}^{l}\rightarrow\mathcal{P}(\mathsf{X}^{l})$ in
Algorithm 1. To that end we will use the notation
$x^{i,l}_{\Delta_{l}:k}\in(\mathbb{R}^{d})^{k\Delta_{l}^{-1}}$, where
$l\in\mathbb{N}_{0}$ is the level of discretization, $i\in\\{1,\dots,N\\}$ is
a particle (sample) indicator, $k\in\\{1,\dots,n\\}$ is a time parameter and
$x^{i,l}_{\Delta_{l}:k}=(x_{\Delta_{l}}^{i,l},x_{2\Delta_{l}}^{i,l},\dots,x_{k}^{i,l})$.
The kernel described in Algorithm 1 is called the called the conditional
particle filter, as developed in [2] and allows one to generate, under minor
conditions, an ergodic Markov chain of invariant measure $\pi_{\theta}^{l}$.
By itself, it does not provide unbiased estimates of expectations w.r.t.
$\pi_{\theta}^{l}$, unless $\pi_{\theta}^{l}$ is the initial distribution of
the Markov chain. However, the kernel will be of use in our subsequent
discussion.

Our approach generates a Markov chain $\\{Z_{m}\\}_{m\in\mathbb{N}_{0}}$ on
the space $\mathsf{Z}^{0}:=\mathsf{X}^{0}\times\mathsf{X}^{0}$,
$Z_{m}\in\mathsf{Z}^{0}$. In order to describe how one can simulate this
Markov chain, we introduce several objects which will be needed. The first of
which is the kernel
$\check{p}^{l}:\Theta\times\mathbb{R}^{2d}\rightarrow\mathcal{P}(\mathbb{R}^{2d})$,
which we need in the case $l=0$ and its simulation is described in Algorithm
2. We will also need to simulate the maximal coupling of two probability mass
functions on $\\{1,\dots,N\\}$, for some $N\in\mathbb{N}$, and this is
described in Algorithm 3.

######  Remark 3.1.

Step 4. of Algorithm 3 can be modified to the case where one generates the
pair $(i,j)\in\\{1,\dots,N\\}^{2}$ from any coupling of the two probability
mass functions $(r_{4},r_{5})$. In our simulations in Section 4 we will do
this by sampling by inversion from $(r_{4},r_{5})$, using the same uniform
random variable. However, to simplify the mathematical analysis that we will
give in the Appendix, we consider exactly Algorithm 3 in our calculations.

To describe the CCPF kernel, we must first introduce a driving coupled
conditional particle filter, which is presented in Algorithm 4. The driving
coupled conditional particle filter is nothing more than an ordinary coupled
particle filter, except the final pair of trajectories is ‘frozen’ as is given
to the algorithm (that is $(x_{1:n}^{\prime},\bar{x}_{1:n}^{\prime})$ as in
step 1\. of Algorithm 4) and allowed to interact with the rest of the particle
system. Given the ingredients in Algorithms 2-4 we are now in a position to
describe the CCPF kernel, which is a Markov kernel
$K^{0}:\Theta\times\mathsf{Z}^{0}\rightarrow\mathcal{P}(\mathsf{Z}^{0})$,
whose simulation is presented in Algorithm 5. We will consider the Markov
chain $\\{Z_{m}\\}_{m\in\mathbb{N}_{0}}$,
$Z_{m}=(X_{1:n}(m),\bar{X}_{1:n}(m))$, generated by the CCPF kernel in
Algorithm 5 and with initial distribution

$\nu_{\theta}^{0}\Big{(}d(x_{1:n},\bar{x}_{1:n})\Big{)}=\int_{\mathsf{Z}^{0}}\Big{(}\prod_{k=1}^{n}p_{\theta}^{0}(x_{k-1}^{\prime},x_{k}^{\prime})\Big{)}\Big{(}\prod_{k=1}^{n}p_{\theta}^{0}(\bar{x}_{k-1}^{\prime},\bar{x}_{k}^{\prime})\Big{)}C_{\theta}^{0}(x_{1:n}^{\prime},dx_{1:n})\delta_{\\{\bar{x}_{1:n}^{\prime}\\}}(d\bar{x}_{1:n}),d(x_{1:n},\bar{x}_{1:n})$
(3.4)

where $x_{0}^{\prime}=\bar{x}_{0}^{\prime}=x_{\star}$.

We remark that in Algorithm 5, marginally, $x_{1:n}^{i}$ (resp.
$\bar{x}_{1:n}^{j}$) has been generated according to
$C_{\theta}^{0}(x_{1:n},\cdot)$ (resp. $C_{\theta}^{0}(\bar{x}_{1:n},\cdot)$).
A rather important point is that if the two input trajectories in step 1. of
Algorithm 5 are equal, i.e. $x_{1:n}=\bar{x}_{1:n}$, then the output
trajectories will also be equal. To that end, define the stopping time
associated to the given Markov chain

$\tau^{0}=\inf\\{m\geq 1:x_{1:n}(m)=\bar{x}_{1:n}(m)\\}.$

Then, setting $m^{*}\in\\{2,3,\dots\\}$ one has the following estimator

$\widehat{\pi}_{\theta}^{0}(G_{\theta}^{0}):=G_{\theta}^{0}(X_{1:n}(m^{*}))+\sum_{m=m^{*}+1}^{\tau^{0}-1}\\{G_{\theta}^{0}(X_{1:n}(m))-G_{\theta}^{0}(\bar{X}_{1:n}(m))\\},$
(3.5)

and one sets $\Xi_{\theta}^{0}=\widehat{\pi}_{\theta}^{0}(G_{\theta}^{0})$.
The procedure for computing $\Xi_{\theta}^{0}$ is summarized in Algorithm 6.

  1. 1.

Input $(x_{0},\bar{x}_{0})\in\mathbb{R}^{2d}$ and the level
$l\in\mathbb{N}_{0}$.

  2. 2.

Generate
$V_{k\Delta_{l}}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d}(0,\Delta_{l}I_{d})$,
for $k\in\\{1,\dots,\Delta_{l}^{-1}\\}$.

  3. 3.

Run the two recursions, for $k\in\\{1,\dots,\Delta_{l}^{-1}\\}$:

$\displaystyle X_{k\Delta_{l}}$ $\displaystyle=$ $\displaystyle
X_{(k-1)\Delta_{l}}+a_{\theta}(X_{(k-1)\Delta_{l}})\Delta_{l}+\sigma(X_{(k-1)\Delta_{l}})V_{k\Delta_{l}}$
$\displaystyle\bar{X}_{k\Delta_{l}}$ $\displaystyle=$
$\displaystyle\bar{X}_{(k-1)\Delta_{l}}+a_{\theta}(\bar{X}_{(k-1)\Delta_{l}})\Delta_{l}+\sigma(\bar{X}_{(k-1)\Delta_{l}})V_{k\Delta_{l}}.$

  4. 4.

Return $(x_{1},\bar{x}_{1})\in\mathbb{R}^{2d}$.

Algorithm 2 Simulating the Kernel $\check{p}_{\theta}^{l}$.

  1. 1.

Input: Two probability mass functions (PMFs) $(r_{1}^{1},\dots,r_{1}^{N})$ and
$(r_{2}^{1},\dots,r_{2}^{N})$ on $\\{1,\dots,N\\}$.

  2. 2.

Generate $U\sim\mathcal{U}_{[0,1]}$.

  3. 3.

If $U<\sum_{i=1}^{N}\min\\{r_{1}^{i},r_{2}^{i}\\}=:\bar{r}$ then generate
$i\in\\{1,\dots,N\\}$ according to the probability mass function:

$r_{3}^{i}=\frac{1}{\bar{r}}\min\\{r_{1}^{i},r_{2}^{i}\\},$

and set $j=i$.

  4. 4.

Otherwise generate $i\in\\{1,\dots,N\\}$ and $j\in\\{1,\dots,N\\}$
independently according to the probability mass functions

$r_{4}^{i}=\frac{1}{1-\bar{r}}(r_{1}^{i}-\min\\{r_{1}^{i},r_{2}^{i}\\}),$

and

$r_{5}^{j}=\frac{1}{1-\bar{r}}(r_{2}^{j}-\min\\{r_{1}^{j},r_{2}^{j}\\}),$

respectively.

  5. 5.

Output: $(i,j)\in\\{1,\dots,N\\}^{2}$. $i$, marginally has PMF $r_{1}^{i}$ and
$j$, marginally has PMF $r_{2}^{j}$.

Algorithm 3 Simulating a Maximal Coupling of Two Probability Mass Functions on
$\\{1,\dots,N\\}$.

  1. 1.

Input $(x_{1:n}^{\prime},\bar{x}_{1:n}^{\prime})\in\mathsf{Z}^{0}$. Set $k=1$,
$(x_{0}^{i},\bar{x}_{0}^{i})=(x_{\star},x_{\star})$,
$a_{0}^{i}=\bar{a}_{0}^{i}=i$ for $i\in\\{1,\dots,N-1\\}$.

  2. 2.

Sampling: for $i\in\\{1,\dots,N-1\\}$ sample
$(x_{k}^{i},\bar{x}_{k}^{i})|(x_{k-1}^{a_{k-1}^{i}},\bar{x}_{k-1}^{\bar{a}_{k-1}^{i}})$
using the Markov kernel $\check{p}_{\theta}^{0}$ in Algorithm 2. Set
$(x_{k}^{N},\bar{x}_{k}^{N})=(x_{k}^{\prime},\bar{x}_{k}^{\prime})$ and for
$i\in\\{1,\dots,N-1\\}$,
$(x_{1:k}^{i},\bar{x}_{1:k}^{i})=((x_{1:k-1}^{a_{k-1}^{i}},x_{k}^{i}),(\bar{x}_{1:k-1}^{\bar{a}_{k-1}^{i}},\bar{x}_{k}^{i}))$.
If $k=n$ stop.

  3. 3.

Resampling: Construct the two probability mass functions on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{k}|x_{k}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|x_{k}^{j})}\quad
r_{2}^{i}=\frac{g_{\theta}(y_{k}|\bar{x}_{k}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|\bar{x}_{k}^{j})}\quad
i\in\\{1,\dots,N\\}.$

For $i\in\\{1,\dots,N-1\\}$ sample $(a_{k}^{i},\bar{a}_{k}^{i})$ from the
maximum coupling of the two given probability mass functions, using Algorithm
3. Set $k=k+1$ and return to the start of 2..

Algorithm 4 Driving Coupled Conditional Particle Filter at level $0$.

  1. 1.

Input $(x_{1:n},\bar{x}_{1:n})\in\mathsf{Z}^{0}$.

  2. 2.

Run Algorithm 4.

  3. 3.

Construct the two probability mass functions on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{n}|x_{n}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|x_{n}^{j})}\quad
r_{2}^{i}=\frac{g_{\theta}(y_{n}|\bar{x}_{n}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|\bar{x}_{n}^{j})}\quad
i\in\\{1,\dots,N\\}.$

Sample $(i,j)\in\\{1,\dots,N\\}^{2}$ from the maximum coupling of the two
given probability mass functions, using Algorithm 3. Return
$(x_{1:n}^{i},\bar{x}_{1:n}^{j})$ which are the path of samples at indices $i$
and $j$ in step 2. of Algorithm 4 when $k=n$.

Algorithm 5 The Coupled Conditional Particle Filter at level $0$.

  1. 1.

Initialize the Markov chain by generating $Z_{0}$ using (3.4). Set $m=1$

  2. 2.

Generate $Z_{m}|Z_{m-1}$ using the Markov kernel described in Algorithm 5. If
$x_{1:n}(m)=\bar{x}_{1:n}(m)$ stop and return
$\Xi_{0}=\widehat{\pi}_{\theta}^{0}(G_{\theta}^{0})$ as in (3.5). Otherwise
set $m=m+1$ and return to the start of 2..

Algorithm 6 Computing $\Xi_{\theta}^{0}$.

###  3.3 Computing $(\Xi_{\theta}^{l})_{l\in\mathbb{N}}$

We are now concerned with the task of computing
$(\Xi_{\theta}^{l})_{l\in\mathbb{N}}$ such that (3.1)-(3.2) are satisfied.
Throughout the section $l\in\mathbb{N}$ is fixed. We will generate a Markov
chain $\\{\check{Z}_{m}^{l}\\}_{m\in\mathbb{N}_{0}}$ on the space
$\mathsf{Z}^{l}\times\mathsf{Z}^{l-1}$, where
$\mathsf{Z}^{l}=\mathsf{X}^{l}\times\mathsf{X}^{l}$ and
$\check{Z}_{m}^{l}\in\mathsf{Z}^{l}\times\mathsf{Z}^{l-1}$. In order to
construct our Markov chain kernel, as in the previous section, we will need to
provide some algorithms. We begin with the Markov kernel
$\check{q}^{l}:\Theta\times\mathbb{R}^{4d}\rightarrow\mathcal{P}(\mathbb{R}^{\Delta_{l}^{-1}2d}\times\mathbb{R}^{\Delta_{l-1}^{-1}2d})$
which will be needed and whose simulation is described in Algorithm 7. We will
also need to sample a coupling for four probability mass functions on
$\\{1,\dots,N\\}$ and this is presented in Algorithm 8.

  1. 1.

Input
$(x_{0}^{l},\bar{x}_{0}^{l},x_{0}^{l-1},\bar{x}_{0}^{l-1})\in\mathbb{R}^{4d}$
and the level $l\in\mathbb{N}_{0}$.

  2. 2.

Generate
$V_{k\Delta_{l}}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d}(0,\Delta_{l}I_{d})$,
for $k\in\\{1,\dots,\Delta_{l}^{-1}\\}$.

  3. 3.

Run the two recursions, for $k\in\\{1,\dots,\Delta_{l}^{-1}\\}$:

$\displaystyle X_{k\Delta_{l}}^{l}$ $\displaystyle=$ $\displaystyle
X_{(k-1)\Delta_{l}}^{l}+a_{\theta}(X_{(k-1)\Delta_{l}}^{l})\Delta_{l}+\sigma(X_{(k-1)\Delta_{l}}^{l})V_{k\Delta_{l}}$
$\displaystyle\bar{X}_{k\Delta_{l}}^{l}$ $\displaystyle=$
$\displaystyle\bar{X}_{(k-1)\Delta_{l}}^{l}+a_{\theta}(\bar{X}_{(k-1)\Delta_{l}}^{l})\Delta_{l}+\sigma(\bar{X}_{(k-1)\Delta_{l}}^{l})V_{k\Delta_{l}}.$

  4. 4.

Run the two recursions, for $k\in\\{1,\dots,\Delta_{l-1}^{-1}\\}$:

$\displaystyle X_{k\Delta_{l}}^{l-1}$ $\displaystyle=$ $\displaystyle
X_{(k-1)\Delta_{l}}^{l-1}+a_{\theta}(X_{(k-1)\Delta_{l}}^{l-1})\Delta_{l}+\sigma(X_{(k-1)\Delta_{l}}^{l-1})[V_{(2k-1)\Delta_{l}}+V_{2k\Delta_{l}}]$
$\displaystyle\bar{X}_{k\Delta_{l}}^{l-1}$ $\displaystyle=$
$\displaystyle\bar{X}_{(k-1)\Delta_{l}}^{l-1}+a_{\theta}(\bar{X}_{(k-1)\Delta_{l}}^{l-1})\Delta_{l}+\sigma(\bar{X}_{(k-1)\Delta_{l}}^{l-1})[V_{(2k-1)\Delta_{l}}+V_{2k\Delta_{l}}].$

  5. 5.

Return
$(x_{\Delta_{l}:1}^{l},\bar{x}_{\Delta_{l}:1}^{l},x_{\Delta_{l-1}:1}^{l-1},\bar{x}_{\Delta_{l-1}:1}^{l-1})\in\mathbb{R}^{\Delta_{l}^{-1}2d}\times\mathbb{R}^{\Delta_{l-1}^{-1}2d}$.

Algorithm 7 Simulating the Kernel $\check{q}_{\theta}^{l}$.

  1. 1.

Input: Four PMFs
$(r_{1}^{1},\dots,r_{1}^{N}),\dots,(r_{4}^{1},\dots,r_{4}^{N})$ on
$\\{1,\dots,N\\}$.

  2. 2.

If $r_{1}^{i}=r_{3}^{i}$ for every $i\in\\{1,\dots,N\\}$ and there exists at
least one $i\in\\{1,\dots,N\\}$ such that $r_{2}^{i}\neq r_{4}^{i}$ then
sample $(i_{1},i_{2})$ according to the maximal coupling of
$(r_{1}^{i_{1}},r_{2}^{i_{2}})$ in Algorithm 3. Implement 5. with
$r_{5}^{i}=r_{1}^{i}$ and $r_{6}^{i}=r_{4}^{i}$, $i\in\\{1,\dots,N\\}$ and
$i_{5}=i_{1}$. Set $(i_{3},i_{4})=(i_{5},i_{6})$ where $(i_{5},i_{6})$ have
been computed from step 5. Go to 7..

  3. 3.

If $r_{2}^{i}=r_{4}^{i}$ for every $i\in\\{1,\dots,N\\}$ and there exists at
least one $i\in\\{1,\dots,N\\}$ such that $r_{1}^{i}\neq r_{3}^{i}$ then
sample $(i_{1},i_{2})$ according to the maximal coupling of
$(r_{1}^{i_{1}},r_{2}^{i_{2}})$ in Algorithm 3. Implement 5. with
$r_{5}^{i}=r_{2}^{i}$ and $r_{6}^{i}=r_{3}^{i}$, $i\in\\{1,\dots,N\\}$ and
$i_{5}=i_{2}$. Set $(i_{3},i_{4})=(i_{6},i_{5})$ where $(i_{5},i_{6})$ have
been computed from step 5. Go to 7..

  4. 4.

Otherwise implement 6. with $r_{j+6}^{i}=r_{j}^{i}$,
$(i,j)\in\\{1,\dots,N\\}\times\\{1,\dots,4\\}$. Set
$(i_{1},\dots,i_{4})=(i_{7},\dots,i_{10})$ where $(i_{7},\dots,i_{10})$ have
been computed from step 6. Go to 7..

  5. 5.

Conditional Algorithm based on [31].

    1. (a)

Input two PMFs $(r_{5}^{1},\dots,r_{5}^{N}),(r_{6}^{1},\dots,r_{6}^{N})$ on
$\\{1,\dots,N\\}$ and $i_{5}\in\\{1,\dots,N\\}$ drawn according to $r_{5}$.

    2. (b)

Sample $U\sim\mathcal{U}_{[0,r_{5}^{i_{5}}]}$. If $U<r_{6}^{i_{5}}$ set
$i_{6}=i_{5}$ and go to (c). Otherwise go to (b).

    3. (c)

Sample $i_{6}^{\prime}$ from $r_{6}^{i_{6}^{\prime}}$. Sample
$U^{\prime}\sim\mathcal{U}_{[0,r_{6}^{i_{6}^{\prime}}]}$. If
$U^{\prime}>r_{5}^{i_{6}^{\prime}}$ set $i_{6}=i_{6}^{\prime}$ and go to (c).
Otherwise start (b) again.

    4. (d)

Output: $(i_{5},i_{6})$.

  6. 6.

Sampling Maximal Couplings of Maximal Couplings

    1. (a)

Input four PMFs
$(r_{7}^{1},\dots,r_{7}^{N}),\dots,(r_{10}^{1},\dots,r_{10}^{N})$ on
$\\{1,\dots,N\\}$. For $j\in\\{7,9\\}$ define the PMFs

$\check{r}_{j}(i_{j},i_{j+1})=r_{j}^{i_{j}}\wedge
r_{j+1}^{i_{j+1}}+\frac{r_{j}^{i_{j}}-r_{j}^{i_{j}}\wedge
r_{j+1}^{i_{j}}}{1-\sum_{i=1}^{N}r_{j}^{i}\wedge
r_{j+1}^{i}}\Big{\\{}r_{j+1}^{i_{j+1}}-r_{j}^{i_{j+1}}\wedge
r_{j+1}^{i_{j+1}}\Big{\\}}.$

    2. (b)

Sample $(i_{7},i_{8})$ according to the maximal coupling of
$(r_{7}^{i_{7}},r_{8}^{i_{8}})$ in Algorithm 3. Generate
$U\sim\mathcal{U}_{[0,\check{r}_{7}(i_{7},i_{8})]}$. If
$U<\check{r}_{9}(i_{7},i_{8})$ set $(i_{9},i_{10})=(i_{7},i_{8})$ and go to
(d). Otherwise go to (c).

    3. (c)

Sample $(i_{9}^{\prime},i_{10}^{\prime})$ according to the maximal coupling of
$(r_{9}^{i_{9}^{\prime}},r_{10}^{i_{10}^{\prime}})$ in Algorithm 3. Sample
$U^{\prime}\sim\mathcal{U}_{[0,\check{r}_{9}(i_{9}^{\prime},i_{10}^{\prime})]}$.
If $U^{\prime}>\check{r}_{7}(i_{9}^{\prime},i_{10}^{\prime})$ set
$(i_{9}^{\prime},i_{10}^{\prime})=(i_{9},i_{10})$ and go to (d). Otherwise
start (c) again.

    4. (d)

Output: $(i_{7},i_{8},i_{9},i_{10})$.

  7. 7.

Output: $(i_{1},i_{2},i_{3},i_{4})\in\\{1,\dots,N\\}^{4}$. $i_{j}$, marginally
has PMF $r_{j}^{i}$, $j\in\\{1,\dots,4\\}$.

Algorithm 8 Simulating a Maximal Coupling of Maximal Couplings associated to
Four Probability Mass Functions on $\\{1,\dots,N\\}$.

To continue onwards, we will consider a generalization of that in Algorithm 4.
The driving coupled conditional particle filter at level $l$ is described in
Algorithm 10. Now given Algorithms 7-10 we are in a position to give our
Markov kernel,
$\check{K}^{l}:\Theta\times\mathsf{Z}^{l}\times\mathsf{Z}^{l-1}\rightarrow\mathcal{P}(\mathsf{Z}^{l}\times\mathsf{Z}^{l-1})$,
which we shall call the coupled-CCPF (C-CCPF) and it is given in Algorithm 11.
To assist the subsequent discussion, we will introduce the marginal Markov
kernel:

$\check{q}_{\theta}^{(l)}\Big{(}[x_{0}^{l},x_{0}^{l-1}],d[x_{\Delta_{l}:1}^{l},x_{\Delta_{l-1}:1}^{l-1}]\Big{)}:=$
$\int_{(\mathbb{R}^{d})^{\Delta_{l}^{-1}}\times(\mathbb{R}^{d})^{\Delta_{l-1}^{-1}}}\check{q}_{\theta}^{l}\Big{(}[(x_{0}^{l},\bar{x}_{0}^{l}),(x_{0}^{l-1},\bar{x}_{0}^{l-1})],d[(x_{\Delta_{l}:1}^{l},\bar{x}_{\Delta_{l}:1}^{l}),(x_{\Delta_{l-1}:1}^{l-1},\bar{x}_{\Delta_{l-1}:1}^{l-1})]\Big{)}.$
(3.6)

Given this kernel, one can describe the CCPF at two different levels $l,l-1$
in Algorithm 9. Algorithm 9 details a Markov kernel
$\check{C}^{l}:\Theta\times\mathsf{X}^{l}\times\mathsf{X}^{l-1}\rightarrow\mathcal{P}(\mathsf{X}^{l}\times\mathsf{X}^{l-1})$
which we will use in the initialization of our Markov chain to be described
below.

  1. 1.

Input
$(x_{\Delta_{l}:n}^{\prime},\bar{x}_{\Delta_{l-1}:n}^{\prime})\in\mathsf{X}^{l}\times\mathsf{X}^{l-1}$.
Set $k=1$, $(x_{0}^{i},\bar{x}_{0}^{i})=(x_{\star},x_{\star})$,
$a_{0}^{i}=\bar{a}_{0}^{i}=i$ for $i\in\\{1,\dots,N-1\\}$.

  2. 2.

Sampling: for $i\in\\{1,\dots,N-1\\}$ sample
$(x_{k-1+\Delta_{l}:k}^{i},\bar{x}_{k-1+\Delta_{l-1}:k}^{i})|(x_{k-1}^{a_{k-1}^{i}},\bar{x}_{k-1}^{\bar{a}_{k-1}^{i}})$
using the Markov kernel $\check{q}_{\theta}^{(l)}$ in (3.6). Set
$(x_{k-1+\Delta_{l}:k}^{N},\bar{x}_{k-1+\Delta_{l-1}:k}^{N})=(x_{k-1+\Delta_{l}:k}^{\prime},\bar{x}_{k-1+\Delta_{l-1}:k}^{\prime})$
and for $i\in\\{1,\dots,N-1\\}$,
$(x_{\Delta_{l}:k}^{i},\bar{x}_{\Delta_{l-1}:k}^{i})=((x_{\Delta_{l}:k-1}^{a_{k-1}^{i}},x_{k-1+\Delta_{l}:k}^{i}),(\bar{x}_{\Delta_{l-1}:k-1}^{\bar{a}_{k-1}^{i}},\bar{x}_{k-1+\Delta_{l-1}:k}^{i}))$.
If $k=n$ go to 4..

  3. 3.

Resampling: Construct the two probability mass functions on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{k}|x_{k}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|x_{k}^{j})}\quad
r_{2}^{i}=\frac{g_{\theta}(y_{k}|\bar{x}_{k}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|\bar{x}_{k}^{j})}\quad
i\in\\{1,\dots,N\\}.$

For $i\in\\{1,\dots,N-1\\}$ sample $(a_{k}^{i},\bar{a}_{k}^{i})$ from the
maximum coupling of the two given probability mass functions, using Algorithm
3. Set $k=k+1$ and return to the start of 2..

  4. 4.

Construct the two probability mass functions on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{n}|x_{n}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|x_{n}^{j})}\quad
r_{2}^{i}=\frac{g_{\theta}(y_{n}|\bar{x}_{n}^{i})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|\bar{x}_{n}^{j})}\quad
i\in\\{1,\dots,N\\}.$

Sample $(i,j)\in\\{1,\dots,N\\}^{2}$ from the maximum coupling of the two
given probability mass functions, using Algorithm 3. Return
$(x_{\Delta_{l}:n}^{i},\bar{x}_{\Delta_{l-1}:n}^{j})$ which are the path of
samples at indices $i$ and $j$ in step 2. when $k=n$.

Algorithm 9 Coupled Conditional Particle Filter at level $l,l-1$,
$l\in\mathbb{N}$.

We will consider the Markov chain
$\\{\check{Z}_{m}^{l}\\}_{m\in\mathbb{N}_{0}}$, with

$\check{Z}_{m}^{l}=\Big{(}(X_{\Delta_{l}:n}^{l}(m),\bar{X}_{\Delta_{l}:n}^{l}(m)),(X_{\Delta_{l-1}:n}^{l-1}(m),\bar{X}_{\Delta_{l-1}:n}^{l-1}(m))\Big{)},$

generated by the C-CCPF kernel in Algorithm 11 and with initial distribution

$\check{\nu}_{\theta}^{l}\Big{(}d[(x_{\Delta_{l}:n}^{l},\bar{x}_{\Delta_{l}:n}^{l}),(x_{\Delta_{l-1}:n}^{l-1},\bar{x}_{\Delta_{l-1}:n}^{l-1})]\Big{)}=\int_{\mathsf{Z}^{l}\times\mathsf{Z}^{l-1}}\Big{(}\prod_{k=1}^{n}\check{q}_{\theta}^{(l)}\Big{(}[x_{k-1}^{l,\prime},x_{k-1}^{l-1,\prime}],d[x_{k-1+\Delta_{l}:k}^{l,\prime},x_{k-1+\Delta_{l-1}:k}^{l-1,\prime}]\Big{)}\Big{)}\times$
$\Big{(}\prod_{k=1}^{n}\check{q}_{\theta}^{(l)}\Big{(}[\bar{x}_{k-1}^{l,\prime},\bar{x}_{k-1}^{l-1,\prime}],d[\bar{x}_{k-1+\Delta_{l}:k}^{l,\prime},\bar{x}_{k-1+\Delta_{l-1}:k}^{l-1,\prime}]\Big{)}\check{C}_{\theta}^{l}\Big{(}[x_{\Delta_{l}:n}^{l,\prime},x_{\Delta_{l-1}:n}^{l-1,\prime}],d[x_{\Delta_{l}:n}^{l},x_{\Delta_{l-1}:n}^{l-1}]\Big{)}\times$
$\delta_{\\{\bar{x}_{\Delta_{l}:n}^{l,\prime},\bar{x}_{\Delta_{l-1}:n}^{l-1,\prime}\\}}(d[\bar{x}_{\Delta_{l}:n}^{l},\bar{x}_{\Delta_{l-1}:n}^{l-1}]),$
(3.7)

where
$x_{0}^{l,\prime}=x_{0}^{l-1,\prime}=\bar{x}_{0}^{l,\prime}=\bar{x}_{0}^{l-1,\prime}=x_{\star}$.
An important point, as in the case of Algorithm 5, is that if the two input
trajectories in step 1. of Algorithm 11 are equal, i.e.
$x_{\Delta_{l}:n}^{l}=\bar{x}_{\Delta_{l}:n}^{l}$, or
$x_{\Delta_{l-1}:n}^{l-1}=\bar{x}_{\Delta_{l-1}:n}^{l-1}$, then the associated
output trajectories will also be equal. As before, we define the stopping
times associated to the given Markov chain
$(\check{Z}_{m}^{l})_{m\in\mathbb{N}_{0}}$, $s\in\\{l,l-1\\}$

$\tau^{s}=\inf\\{m\geq
1:X_{\Delta_{s}:n}^{s}(m)=\bar{X}_{\Delta_{s}:n}^{s}(m)\\}.$

Then, setting $m^{*}\in\\{2,3,\dots\\}$ one has the following estimator

$\widehat{\pi}_{\theta}^{s}(G_{\theta}^{s}):=G_{\theta}^{s}(X_{\Delta_{s}:n}^{s}(m^{*}))+\sum_{m=m^{*}+1}^{\tau^{s}-1}\\{G_{\theta}^{s}(X_{\Delta_{s}:n}^{s}(m))-G_{\theta}^{s}(\bar{X}_{\Delta_{s}:n}^{s}(m))\\},$
(3.8)

and one sets
$\Xi_{\theta}^{l}=\widehat{\pi}_{\theta}^{l}(G_{\theta}^{l})-\widehat{\pi}_{\theta}^{l-1}(G_{\theta}^{l-1})$.
The procedure for computing $\Xi_{\theta}^{l}$ is summarized in Algorithm 12.

  1. 1.

Input
$((x_{\Delta_{l}:n}^{l},\bar{x}_{\Delta_{l}:n}^{l}),(x_{\Delta_{l-1}:n}^{l-1},\bar{x}_{\Delta_{l-1}:n}^{l-1}))\in\mathsf{Z}^{l}\times\mathsf{Z}^{l-1}$.
Set $k=1$,
$(x_{0}^{i,l},\bar{x}_{0}^{i,l})=(x_{\star},x_{\star})=(x_{0}^{i,l-1},\bar{x}_{0}^{i,l-1})$,
$a_{0}^{i,l}=\bar{a}_{0}^{i,l}=a_{0}^{i,l-1}=\bar{a}_{0}^{i,l-1}=i$ for
$i\in\\{1,\dots,N-1\\}$.

  2. 2.

Sampling: for $i\in\\{1,\dots,N-1\\}$ sample

$\Big{(}(x_{k-1+\Delta_{l}:k}^{i,l},\bar{x}_{k-1+\Delta_{l}:k}^{i,l}),(x_{k-1+\Delta_{l-1}:k}^{i,l-1},\bar{x}_{k-1+\Delta_{l-1}:k}^{i,l-1})\Big{)}\Big{|}\Big{(}(x_{k-1}^{a_{k-1}^{i,l},l},\bar{x}_{k-1}^{\bar{a}_{k-1}^{i,l},l}),(x_{k-1}^{a_{k-1}^{i,l-1},l-1},\bar{x}_{k-1}^{\bar{a}_{k-1}^{i,l-1},l-1})\Big{)}$

using the Markov kernel $\check{q}_{\theta}^{l}$ in Algorithm 7. Set

$\Big{(}(x_{k-1+\Delta_{l}:k}^{N,l},\bar{x}_{k-1+\Delta_{l}:k}^{N,l}),(x_{k-1+\Delta_{l-1}:k}^{N,l-1},\bar{x}_{k-1+\Delta_{l-1}:k}^{N,l-1})\Big{)}=$
$\Big{(}(x_{k-1+\Delta_{l}:k}^{l},\bar{x}_{k-1+\Delta_{l}:k}^{l}),(x_{k-1+\Delta_{l-1}:k}^{l-1},\bar{x}_{k-1+\Delta_{l-1}:k}^{l-1})\Big{)}$

and for $i\in\\{1,\dots,N-1\\}$

$\displaystyle(x_{\Delta_{l}:k}^{i,l},\bar{x}_{\Delta_{l}:k}^{i,l})$
$\displaystyle=$
$\displaystyle\Big{(}(x_{\Delta_{l}:k-1}^{a_{k-1}^{i,l},l},x_{k-1+\Delta_{l}:k}^{i,l}),(\bar{x}_{\Delta_{l}:k-1}^{\bar{a}_{k-1}^{i,l},l-1},\bar{x}_{k-1+\Delta_{l}:k}^{i,l})\Big{)},$
$\displaystyle(x_{\Delta_{l-1}:k}^{i,l-1},\bar{x}_{\Delta_{l-1}:k}^{i,l-1})$
$\displaystyle=$
$\displaystyle\Big{(}(x_{\Delta_{l-1}:k-1}^{a_{k-1}^{i,l-1},l},x_{k-1+\Delta_{l-1}:k}^{i,l-1}),(\bar{x}_{\Delta_{l-1}:k-1}^{\bar{a}_{k-1}^{i,l-1},l-1},\bar{x}_{k-1+\Delta_{l}:k}^{i,l-1})\Big{)}.$

If $k=n$ stop.

  3. 3.

Resampling: Construct the four probability mass functions on
$\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{k}|x_{k}^{i,l})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|x_{k}^{j,l})}\quad
r_{3}^{i}=\frac{g_{\theta}(y_{k}|\bar{x}_{k}^{i,l})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|\bar{x}_{k}^{j,l})}\quad
i\in\\{1,\dots,N\\},$

and

$r_{2}^{i}=\frac{g_{\theta}(y_{k}|x_{k}^{i,l-1})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|x_{k}^{j,l-1})}\quad
r_{4}^{i}=\frac{g_{\theta}(y_{k}|\bar{x}_{k}^{i,l-1})}{\sum_{j=1}^{N}g_{\theta}(y_{k}|\bar{x}_{k}^{j,l-1})}\quad
i\in\\{1,\dots,N\\}.$

For $i\in\\{1,\dots,N-1\\}$ sample
$(a_{k}^{i,l},a_{k}^{i,l-1},\bar{a}_{k}^{i,l},\bar{a}_{k}^{i,l-1})$ using
Algorithm 8. Set $k=k+1$ and return to the start of 2..

Algorithm 10 Driving Coupled Conditional Particle Filter at level
$l\in\mathbb{N}$.

  1. 1.

Input
$((x_{\Delta_{l}:n}^{l},\bar{x}_{\Delta_{l}:n}^{l}),(x_{\Delta_{l-1}:n}^{l-1},\bar{x}_{\Delta_{l-1}:n}^{l-1}))\in\mathsf{Z}^{l}\times\mathsf{Z}^{l-1}$.

  2. 2.

Run Algorithm 10.

  3. 3.

Construct the four probability mass functions on $\\{1,\dots,N\\}$:

$r_{1}^{i}=\frac{g_{\theta}(y_{n}|x_{n}^{i,l})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|x_{n}^{j,l})}\quad
r_{3}^{i}=\frac{g_{\theta}(y_{n}|\bar{x}_{n}^{i,l})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|\bar{x}_{n}^{j,l})}\quad
i\in\\{1,\dots,N\\},$

and

$r_{2}^{i}=\frac{g_{\theta}(y_{n}|x_{n}^{i,l-1})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|x_{n}^{j,l-1})}\quad
r_{4}^{i}=\frac{g_{\theta}(y_{n}|\bar{x}_{n}^{i,l-1})}{\sum_{j=1}^{N}g_{\theta}(y_{n}|\bar{x}_{n}^{j,l-1})}\quad
i\in\\{1,\dots,N\\}.$

Sample $(i_{1},i_{2},i_{3},i_{4})\in\\{1,\dots,N\\}^{4}$ using Algorithm 8.
Return
$((x_{\Delta_{l}:n}^{i_{1},l},\bar{x}_{\Delta_{l}:n}^{i_{3},l}),(x_{\Delta_{l-1}:n}^{i_{2},l-1},\bar{x}_{\Delta_{l-1}:n}^{i_{4},l-1}))$
which are the path of samples at indices $i_{1:4}$ in step 2. of Algorithm 10
when $k=n$.

Algorithm 11 The Coupled-CCPF at level $l\in\mathbb{N}$.

  1. 1.

Initialize the Markov chain by generating $\check{Z}_{0}^{l}$ using (3.7). Set
$m=1$

  2. 2.

Generate $\check{Z}_{m}^{l}|\check{Z}^{l}_{m-1}$ using the Markov kernel
described in Algorithm 11. If
$x_{\Delta_{l}:n}^{l}(m)=\bar{x}_{\Delta_{l}:n}^{l}(m)$ and
$x_{\Delta_{l-1}:n}^{l-1}(m)=\bar{x}_{\Delta_{l-1}:n}^{l-1}(m)$ stop and
return
$\Xi_{\theta}^{l}=\widehat{\pi}_{\theta}^{l}(G_{\theta}^{l})-\widehat{\pi}_{\theta}^{l-1}(G_{\theta}^{l-1})$
where $\widehat{\pi}_{\theta}^{s}(G_{\theta}^{s})$, $s\in\\{l,l-1\\}$, is as
(3.8). Otherwise set $m=m+1$ and return to the start of 2..

Algorithm 12 Computing $\Xi_{\theta}^{l}$.

###  3.4 Estimate and Remarks

Given the commentary above we are ready to present the procedure for our
unbiased estimate of $\mathfrak{H}_{\theta}^{(ij)}$ for each $(i,j)\in
D_{\theta}:=\\{(i,j)\in\\{1,\dots,d_{\theta}\\}^{2}:i\leq j\\}$;
$\mathfrak{H}_{\theta}$ is a symmetric matrix. The two main Algorithms we will
use (Algorithms 6 and 12) are stated in terms of providing $\Xi_{\theta}^{l}$
in terms of one specified function $G_{\theta}^{l,(i)}$ (recall that $i$ was
suppressed from the notation). However, the algorithms can be run once and
provide an unbiased estimate of
$\pi_{\theta}^{l}(G_{\theta}^{l,(i)})-\pi_{\theta}^{l-1}(G_{\theta}^{l-1,(i)})$
for every $i\in\\{1,\dots,d_{\theta}\\}$, of
$\pi_{\theta}^{l}(G_{\theta}^{l,(i)}G_{\theta}^{l,(j)})-\pi_{\theta}^{l-1}(G_{\theta}^{l-1,(i)}G_{\theta}^{l-1,(j)})$
and
$\pi_{\theta}^{l}(H_{\theta}^{l,(ij)})-\pi_{\theta}^{l-1}(H_{\theta}^{l-1,(ij)})$
for every $(i,j)\in D_{\theta}$. To that end we will write
$\Xi_{\theta}^{l}(G_{\theta}^{l,(i)})$,
$\Xi_{\theta}^{l}(G_{\theta}^{l,(i)}G_{\theta}^{l,(j)})$ and
$\Xi_{\theta}^{l}(H_{\theta}^{l,(ij)})$ to denote the appropriate estimators.

Our approach consists of the following steps, repeated for
$k\in\\{1,\dots,M\\}$:

  1. 1.

Generate $(L_{k},\tilde{L}_{k})$ according to
$\mathbb{P}_{L}\otimes\mathbb{P}_{L}$.

  2. 2.

Compute $\Xi_{\theta}^{k,0}(G_{\theta}^{0,(i)})$ for every
$i\in\\{1,\dots,d_{\theta}\\}$ and
$\Xi_{\theta}^{k,0}(G_{\theta}^{0,(i)}G_{\theta}^{0,(j)})$,
$\Xi_{\theta}^{k,0}(H_{\theta}^{0,(ij)})$ for every $(i,j)\in D_{\theta}$
using Algorithm 6. Independently, compute
$\tilde{\Xi}_{\theta}^{k,0}(G_{\theta}^{0,(i)})$ for every
$i\in\\{1,\dots,d_{\theta}\\}$ using Algorithm 6.

  3. 3.

If $L_{k}>0$ then independently for each $l\in\\{1,\dots,L_{k}\\}$ and
independently of step 2. calculate $\Xi_{\theta}^{k,l}(G_{\theta}^{l,(i)})$
for every $i\in\\{1,\dots,d_{\theta}\\}$ and
$\Xi_{\theta}^{k,l}(G_{\theta}^{l,(i)}G_{\theta}^{l,(j)})$,
$\Xi_{\theta}^{k,l}(H_{\theta}^{l,(ij)})$ for every $(i,j)\in D_{\theta}$
using Algorithm 12.

  4. 4.

If $\tilde{L}_{k}>0$ then independently for each
$l\in\\{1,\dots,\tilde{L}_{k}\\}$ and independently of steps 2. and 3.
calculate $\tilde{\Xi}_{\theta}^{k,l}(G_{\theta}^{l,(i)})$ for every
$i\in\\{1,\dots,d_{\theta}\\}$ using Algorithm 12.

  5. 5.

Compute for every $(i,j)\in D_{\theta}$

$\widehat{\mathfrak{H}}_{\theta}^{k,(ij)}=\Big{(}\sum_{l=0}^{L_{k}}\frac{\Xi_{\theta}^{k,l}(G_{\theta}^{l,(i)})}{\overline{\mathbb{P}}_{L}(l)}\Big{)}\Big{(}\sum_{l=0}^{\tilde{L}_{k}}\frac{\tilde{\Xi}_{\theta}^{k,l}(G_{\theta}^{l,(i)})}{\overline{\mathbb{P}}_{L}(l)}\Big{)}-\sum_{l=0}^{L_{k}}\frac{\Xi_{\theta}^{k,l}(G_{\theta}^{l,(i)}G_{\theta}^{l,(j)})}{\overline{\mathbb{P}}_{L}(l)}-\sum_{l=0}^{L_{k}}\frac{\Xi_{\theta}^{k,l}(H_{\theta}^{l,(ij)})}{\overline{\mathbb{P}}_{L}(l)}.$

Then our estimator is for each $(i,j)\in D_{\theta}$

$\widehat{\mathfrak{H}}_{\theta}^{(ij)}=\frac{1}{M}\sum_{k=1}^{M}\widehat{\mathfrak{H}}_{\theta}^{k,(ij)}.$
(3.9)

The algorithm and the various settings are described and investigated in
details in [18] as well as enhanced estimators. We do not discuss the
methodology further in this work.

######  Proposition 3.1.

Assume (D1-2). Then there exists choices of $\mathbb{P}_{L}$ so that (3.9) is
an unbiased and finite variance estimator of $\mathfrak{H}_{\theta}^{(ij)}$
for each $(i,j)\in D_{\theta}$.

###### Proof.

This is the same as [18, Theorem 2], except one must repeat the arguments of
that paper given Lemma A.3 in the appendix and given the rate in the proof of
Proposition 2.1. Since the arguments and calculations are almost identical,
they are omitted in their entirety. ∎

The main point is that the choice of $\mathbb{P}_{L}$ is as in [18], which is:
In the case that $\sigma$ is constant
$\mathbb{P}_{L}(l)\propto\Delta_{l}(l+1)\log_{2}(2+l)^{2}$ and in the non-
constant case
$\mathbb{P}_{L}(l)\propto\Delta_{l}^{1/2}(l+1)\log_{2}(2+l)^{2}$; both choices
achieve finite variance and costs to achieve an error of
$\mathcal{O}(\epsilon)$ with high-probability as in [27, Propositions 4 and
5].

As we will see in the succeeding section, we will compare our methodology
which is based on the C-CCPF to that of another methodology, which is the
$\Delta$PF, within particle Markov chain Monte Carlo. Specifically it will be
a particle marginal Metropolis Hastings algorithm. We omit such a description
of the later, as we only use it as a comparison, but we refer the reader to
[11] for a more concrete description. However we emphasis with it, that it is
only asymptotically unbiased, in relation to the Hessian identity (2.4).

######  Remark 3.2.

It is important to emphasis that with inverse Hessian, which is required for
Newton methodologies, we can debias both the C-CCPF and the $\Delta$PF. This
can be achieved by using the same techniques which are presented in the work
of Jasra et al. [21].

##  4 Numerical experiments

In this section we demonstrate that our estimate of the Hessian is unbiased
through various different experiments. We consider testing this through the
study of the variance and bias of the mean square error, while also provided
plots related to the Newton-type learning. Our results will be demonstrated on
three underlying diffusion process. Firstly that of a univariate
Ornstein–Uhlenbeck process, a multivariate OU process and the Fitzhugh–Nagumo
model. We compare our methodology to that of using the $\Delta$PF instead of
the coupled-CCPF within our unbiased estimator.

A python package which allows implementation of Hessian estimate (3.9), as
well as score estimate in [18] for general partially observed diffusion model
can be found in: https://github.com/fangyuan-ksgk/Hessian_Estimate.

###  4.1 Ornstein–Uhlenbeck process

Our first set of numerical experiments will be conducted on a univariate
Ornstein–Uhlenbeck (OU) process, which takes the form

$\displaystyle dX_{t}$ $\displaystyle=-\theta_{1}X_{t}dt+\sigma dW_{t},$
$\displaystyle X_{0}$ $\displaystyle=x_{0},$

where $x_{0}\in\mathbb{R}^{+}$ is our initial condition,
$\theta_{1}\in\mathbb{R}$ is a parameter of interest and $\sigma\in\mathbb{R}$
is the diffusion coefficient. For our discrete observations, we assume that we
have Gaussian measurement errors, $Y_{t}|X_{t}\sim
g_{\theta}(\cdot|X_{t})=\mathcal{N}(X_{t},\theta_{2})$ for
$t\in\\{1,\ldots,T\\}$ and for some $\theta_{2}\in\mathbb{R}^{+}$. Our
observations will be generated with parameter choices defined as
$\theta=(\theta_{1},\theta_{2})=(0.46,0.38)$, $x_{0}=1$ and $T=500$.
Throughout the simulation, one observation sequence
$\\{Y_{1},Y_{2},...,Y_{T}\\}$ is used. The true distribution of observations
can be computed analytically, therefore the Hessian is known. In Figure 1, we
present the surface plots comparing the true Hessian with the estimated
Hessian, obtained by the Rhee & Glynn estimator (3.9) truncated at
discretization level $L=8$, this is done by letting
$\mathbb{P}_{L}(l)\propto\Delta_{l}\mathbb{I}(l\leq L)$. We use $M=10^{4}$ to
obtain the estimate Hessian surface plot. Both surface plots are evaluated at
$\theta_{1},\theta_{2}\in\\{0.2,0.3,0.4,\ldots,1.0\\}$. In Figure 2, we test
out the convergence of bias of the Hessian estimate (3.9) with respect to its
truncated discretization level. This essentially tests the result in Lemma
A.2. We uses $L=\\{2,3,4,5,6,7\\}$ and plot the bias against $\Delta_{L}$.

The bias is obtained by using $M=10^{4}$ i.i.d. samples, and taking its entry-
wise difference with the true Hessian entry-wise value. Note that the Hessian
estimate here is evaluated with true parameter choice. As the parameter
$\theta$ is two-dimensional, we present four $\log$ $\log$ plots where the
rate represents the fitted slope of $\log$-scaled bias against $\log$-scaled
$\Delta_{L}$. We observe that the Hessian estimate bias is of order
$\Delta_{L}^{\alpha}$ where $\alpha\in\\{0.9629,0.7536,0.7361,0.8949\\}$
respectively for the four entries, this verifies our result in Lemma A.2. We
also compare the wall-clock time cost of obtaining one realization of Hessian
estimate (3.9) with the cost of obtaining one realization of score estimate
(see [18]), both truncated at same discretization levels
$L=\\{2,3,4,5,6,7\\}$, here $M=100$. The comparison result is provided on top
of Figure 3.

We observe that the cost of obtaining the Hessian estimate. is on average 3
times more expensive than obtaining a Score estimate. The reason for this is
that we needs to simulate three CCPF paths in order to obtain one summand in
the Hessian estimate, while to estimate the score function, we need only one
path. We also record the fitted slope of $\log$-scaled Cost against
$\log$-scaled $\Delta_{L}$ for both estimates, the cost for Hessian estimates
is roughly proportional to $\Delta_{L}^{-0.9}$. To verify the rate obtained in
Lemma A.3, we compare the variance of the Hessian incremental estimate with
respect to discretization level $L\in\\{1,2,3,4,5,6\\}$. The incremental
variance is approximated with the sample variance over $10^{3}$ repetitions,
and we sum over all $2\times 2$ entries and present the $\log$ $\log$ plot of
the summed variance against $\Delta_{L}$ on the bottom of Figure 3. We observe
that the incremental variance is of order $\Delta_{L}^{1.15}$ for the OU
process model. This verifies the result obtained in Lemma A.3. It is known
that when truncated, the Rhee $\&$ Glynn method essentially serves as a
variance reduction method. As a result, compared to the discrete Hessian
estimate (2.6), the truncated Hessian estimate (3.9) will require less cost to
achieve the same MSE target.

We present on top of Figure 4 the $\log$-$\log$ plot of cost against MSE for
discrete Hessian estimate (2.6) and the Rhee $\&$ Glynn (R $\&$ G) Hessian
estimate (3.9). We observe that (2.6) requires much less cost for a MSE target
compared to (3.9). For (2.6), the cost is proportional to
$\mathcal{O}(\epsilon^{-2.974})$ for a MSE target of order
$\mathcal{O}(\epsilon^{2})$. While for (3.9), the cost is proportional to
$\mathcal{O}(\epsilon^{-2.428})$. The average cost ratio between (3.9) and
(2.6) under same MSE target is $5.605$. On the bottom of Figure 4, we present
the $\log$-$\log$ plot of cost against MSE for (3.9) and the hessian estimate
obtained by the $\Delta$PF method. We observe that under similar MSE target,
the latter method on average has costs $5.054$ times less than (3.9). In
Figure 5, we present the convergence plots for Stochastic Gradient Descent
(SGD) method with Score estimate and Newton method with Score $\&$ Hessian
estimate. For both method, the parameter is initialized at $(0.1,0.1)$,
learning rate for the SGD method is set to $0.002$.

Figure 1: Experiments for the OU model. Left: true values of the Hessian.
Right: estimated values of the Hessian. Figure 2: Experiments for the OU
model: bias values of Hessian estimate.

Figure 3: Experiments for the OU model. Top: cost of Hessian & score estimate.
Bottom: incremental Hessian estimate variance summed over all entries.

Figure 4: Hessian estimate cost against summed MSE for the OU model

Figure 5: Parameter estimate for the OU model. Top: SGD with score estimate.
Bottom: Newton method with score & Hessian estimate.

###  4.2 Multivariate Ornstein-Uhlenbeck process model

Our second model of interest is a two-dimensional OU process defined as

$\begin{bmatrix}dX^{(1)}_{t}\\\
dX^{(2)}_{t}\end{bmatrix}=\begin{bmatrix}\theta_{1}-\theta_{2}X_{t}^{(1)}\\\
-\theta_{3}X_{t}^{(2)}\end{bmatrix}dt+\begin{bmatrix}\sigma_{1}\\\
\sigma_{2}\\\ \end{bmatrix}dW_{t},\quad\quad X_{0}=x_{0}.$

where $x_{0}\in\mathbb{R}^{2}$ is the initial condition and
$(\sigma_{1},\sigma_{2})\in\mathbb{R}^{+}\times\mathbb{R}^{+}$ are the
diffusion coefficients. We assume Gaussian measurement errors,
$Y_{t}|X_{t}\sim
g_{\theta}(\cdot|X_{t})=\mathcal{N}_{2}(X_{t},\theta_{4}I_{2})$ where $I_{2}$
is a two-dimensional Identity matrix. We generate one sequence of observations
up to time $T=500$ with parameter choice
$\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})=(0.48,0.78,0.37,0.32)$,
$\sigma_{1}=0.8$, $\sigma_{2}=0.6$, $x_{0}=(1,1)^{T}$. As before, we study
various properties of (3.9) with the true parameter choice. In Figure 6, we
present the $\log$-$\log$ plot of bias against $\Delta_{L}$ for (3.9), where
the five points are evaluated with $L\in\\{2,3,4,5,6\\}$. The bias is
approximated by the difference between (3.9) and the true Hessian with
$M=10^{4}$, we sum over all entry-wise bias and present it on the plot. We
observe that the summed bias is of order $\Delta_{L}^{0.972}$. This verifies
result in Lemma A.2. On top of Figure 7, we present a $\log$-$\log$ plot of
cost against $\Delta_{L}$ for (3.9) and the R$\&$G score estimate both with
$M=10$. The experiments is done over $L\in\\{2,3,4,5,6\\}$. We observe that
the cost of (3.9) is proportional to $\Delta_{L}^{-0.866}$. This rate is
similar to that of the score estimate, on average the cost ratio between (3.9)
and the score estimate is $3.495$.

In Figure 7, we present on the bottom the $\log$-$\log$ plot of summed
incremental variance of Hessian estimate against $\Delta_{L}$ for
$L\in\\{2,3,4,5,6\\}$. We compute the entry-wise sample variance of the
incremental Hessian estimate for $10^{3}$ times, and plot the summed variance
against $\Delta_{L}$. We observe that the Hessian incremental variance is
proportional to $\Delta_{L}^{1.068}$. This verifies the result in Lemma A.3.
On top of Figure 8, we present the $\log$-$\log$ plot of the cost against MSE
for (3.9) and (2.6), where the MSE is approximated through averaging over
$10^{3}$ i.i.d. repetitions of both estimators. We observe that under a summed
MSE target of $\mathcal{O}(\epsilon^{2})$, the cost for (3.9) is of order
$\mathcal{O}(\epsilon^{-2.362})$, while the cost for (2.6) is of order
$\mathcal{O}(\epsilon^{-2.958})$. On average, the cost ratio between (2.6) and
(3.9) is $3.575$. This verifies the variance reduction effect of truncated
R$\&$G scheme. On bottom of Figure 8, we present the $\log$-$\log$ plot of the
cost against MSE for (3.9) and hessian estimate using $\Delta$PF. We observe
that under a similar MSE target, the latter method on average costs $3.165$
times less than that of (3.9). In Figure 9, we present the convergence plots
for SGD and Newton method. Both the score estimate and the Hessian estimate
(3.9) is obtained with $M=2\times 10^{3}$, truncated at level $L=8$. The
learning rate for the SGD is set to $0.005$. The training reach convergence
when the relative Euclidean distance between trained and true $\theta$ is no
bigger than $0.02$. We initialize the training parameter at
$(0.1,0.1,0.1,0.1)$, we observe that the SGD method reaches convergence with
$122$ iterations, while the Modified Newton method reaches convergence with
$4$ iterations. The actual training time until convergence for the Newton
method is roughly $7.6$ times faster than the SGD method.

Figure 6: Hessian estimate bias summed over all entries for the multivariate
OU diffusion model.

Figure 7: Experiments for the Multivariate OU diffusion model. Top: cost of
Hessian & score estimate. Bottom: incremental Hessian estimate variance summed
over all entries.

Figure 8: Hessian estimate cost against summed MSE for the multivariate OU
diffusion model.

Figure 9: Parameter estimate for the multivariate OU model. Top: SGD with the
score estimate. Bottom: Newton method with score & Hessian estimate.

###  4.3 FitzHugh–Nagumo model

Our next model will be a two-dimensional ordinary differential equation, which
arises in neuroscience, known as the FitzHugh–Nagumo (FHN) model [7, 12]. It
is concerned with the membrane potential of a neuron and a (latent) recovery
variable modeling the ion channel kinetics. We consider a stochastic perturbed
extended version, given as

$\begin{bmatrix}dX^{(1)}_{t}\\\
dX^{(2)}_{t}\end{bmatrix}=\begin{bmatrix}\theta_{1}(X^{(1)}_{t}-(X_{t}^{(1)})^{3}-X^{(2)}_{t})\\\
\theta_{2}X^{(1)}_{t}-X^{(2)}_{t}+\theta_{3}\end{bmatrix}dt+\begin{bmatrix}\sigma_{1}\\\
\sigma_{2}\\\ \end{bmatrix}dW_{t},\quad\quad X_{0}=u_{0}.$

for the discrete observations, we assume Gaussian measurement errors,
$Y_{t}|X_{t}\sim
g_{\theta}(\cdot|X_{t})=\mathcal{N}_{2}(X_{t},\theta_{4}I_{2})$, where
$(\theta_{1},\theta_{2},\theta_{3},\theta_{4})\in\mathbb{R}^{+}\times\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{+}$,
( $\sigma_{1},\sigma_{2})\in\mathbb{R}^{+}\times\mathbb{R}^{+}$ are the
diffusion coefficients and, as before, $\\{W_{t}\\}_{t\geq 0}$ is a Brownian
motion. We generate one observation sequence with parameter choices
$\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})=(0.89,0.98,0.5,0.79)$,
$\sigma=(0.2,0.4)$. As the true distribution of the observation is not
available analytically, we uses $L=10$ to simulate out
$\\{Y_{1},Y_{2},...,Y_{T}\\}$ where $T=500$. In Figure 10, we compared the
bias of (3.9), truncated at discretization level $L\in\\{2,3,4,5,6,7\\}$ and
plot it against $\Delta_{L}$ ($\log$-$\log$ plot). The summed bias is obtained
by taking element-wise difference between average of $10^{3}$ i.i.d.
realizations of the Hessian estimate and the true Hessian, then summed over
all the element-wise difference. The true Hessian is approximated by (3.9)
with $M=10^{4}$ and $L=10$. We observe that the summed bias is of order
$\mathcal{O}(\Delta_{L}^{1.402})$. This verifies the result in Lemma A.2. At
the top of Figure (11), we present the $\log$-$\log$ plot of cost against
$\Delta_{L}$ for (3.9) and the R$\&$G score estimate both with $M=10$. The
experiments is done over $L\in\\{2,3,4,5,6,7\\}$. We observe that the cost of
(3.9) is of order $\mathcal{O}(\Delta_{L}^{-0.866})$, while the cost for score
estimate is of order $\mathcal{O}(\Delta_{L}^{-0.867})$. The average cost
ratio between (3.9) and the score estimate is $3.495$. In the bottom of Figure
11, we present the $\log$-$\log$ plot of the summed incremental variance with
$\Delta_{L}$, where $L\in\\{2,3,4,5,6,7\\}$. We observe that the summed
incremental variance is of order $\mathcal{O}(\Delta_{L}^{1.118})$. This
verifies the result in Lemma A.3.

At the top of Figure 12 we present a $\log$-$\log$ plot of the cost against
the summed MSE of 3.9 over all entries for both (3.9) and (2.6). We observe
that under a MSE target of $\epsilon^{2}$, (3.9) requires cost of order
$\mathcal{O}(\epsilon^{-2.482})$, while (2.6) requires cost of order
$\mathcal{O}(\epsilon^{-2.97})$. The average cost ratio between (3.9) and
(2.6) under same MSE target is $3.987$. This verifies the variance reduction
effect of truncated R$\&$G scheme. On the bottom of Figure 12 we present a
$\log$-$\log$ plot of cost against summed MSE for (3.9) and Hessian estimate
using the $\Delta$PF. We observe that under similar MSE target, the latter
method on average costs $4.627$ times less than that of (3.9). In Figure 13,
we present the convergence plots of SGD and the modified Newton method. For
the modified Newton method, we set all the off-diagonal entries to zero for
the Hessian estimate, and add $0.0001$ to the diagonal entries to avoid
singularity. When the $L_{2}$ norm of the score is smaller than $0.1$, we
scale the searching step by a learning rate of $0.002$. Both the score
estimate and the Hessian estimate (3.9) is obtained with $M=2\times 10^{3}$,
truncated at level $L=8$. The learning rate for the SGD is set to $0.001$. The
training reach convergence when the relative Euclidean distance between
trained and true $\theta$ is no bigger than $0.02$. We initialize the training
parameter at $(0.8,0.8,0.8,0.8)$, we observe that the SGD method reaches
convergence with $247$ iterations, while the modified Newton method reaches
convergence with $5$ iterations.

Figure 10: Hessian estimate bias summed over all entries for the FHN model.

Figure 11: Experiments for the FHN model. Top: cost of Hessian & score
estimate. Bottom: incremental Hessian estimate variance summed over all
entries.

Figure 12: Hessian Estimate Cost against summed MSE for FHN model.

Figure 13: Parameter estimate for the FHN model. Top: SGD with score estimate.
Bottom: modified Newton method with score & Hessian estimate.

##  5 Summary

In this work we were interested in developing an unbiased estimator of the
Hessian, related to partially observed diffusion processes. This task is of
interest, as computing the Hessian is primarily biased, due to its
computational cost, but also its has improved convergence over the score
function. We presented a general expression for the Hessian and proved, in the
limit of discretization level, that it is consistent with the continuous-form.
We demonstrated that we were able to reduce the bias, arising from the
discretization. This was shown through various numerical experiments that were
tested on a range of diffusion processes. This not only highlighted the
reduction in bias, but that convergence is better compared to computing and
using the score function. In terms of research directions beyond what we have
done, it would be nice firstly to extend this to more complicated diffusion
models, such as ones arising in mathematical finance [15, 16]. Such diffusion
models would be rough volatility models. Another potential direction would be
to consider diffusion bridges, and analyze how one can could use the tools
here and adapt them. This has been of interest, with recent works such as [26,
29]. Finally one could aim to

## Acknowledgments

This work was supported by KAUST baseline funding.

##  Appendix A Proofs for Proposition 2.1

In this Section we will consider a diffusion process $\\{X_{t}^{x}\\}_{t\geq
0}=\mathbf{X}_{T}^{x}$ which follows (2.1) and has an initial condition
$X_{0}=x\in\mathbb{R}^{d}$ and we will also consider Euler discretizations
(2.5), at some given level $l$, which are driven by the same Brownian motion
as $\\{X_{t}^{x}\\}_{t\geq 0}$ and the same initial condition, written
$(\widetilde{X}_{\Delta_{l}}^{x},\widetilde{X}_{2\Delta_{l}}^{x},\dots)$. We
also consider another diffusion process $\\{X_{t}^{x_{\star}}\\}_{t\geq 0}$
which also follows (2.1), initial condition $X_{0}=x_{\star}\in\mathbb{R}^{d}$
with the same Brownian motion as $\\{X_{t}^{x}\\}_{t\geq 0}$ and associated
Euler discretizations, at level $l$, which are driven by the same Brownian
motion as $\\{X_{t}^{x_{\star}}\\}_{t\geq 0}$ and the same initial condition,
written
$(\widetilde{X}_{\Delta_{l}}^{x_{\star}},\widetilde{X}_{2\Delta_{l}}^{x_{\star}},\dots)$.
The use of signifying the initial condition will be made apparent later on in
the appendix. The expectation operator for the described process is written
$\mathbb{E}_{\theta}$.

We require the following additional assumption called (D2) and all derivatives
are assumed to be well-defined.

  * •

$[\Sigma^{-1}]^{j,k}\in\mathcal{B}_{b}(\mathbb{R}^{d})\cap\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$
$(j,k)\in\\{1,\dots,d\\}^{2}$.

  * •

For any $\theta\in\Theta$, $a_{\theta}^{j}\in\mathcal{B}_{b}(\mathbb{R}^{d})$,
$\sigma^{j,k}\in\mathcal{B}_{b}(\mathbb{R}^{d})$,
$(j,k)\in\\{1,\dots,d\\}^{2}$.

  * •

For any $\theta\in\Theta$, there exists $0<\underline{C}<\overline{C}<+\infty$
such that for any $(x,y)\in\mathbb{R}^{d}\times\mathbb{R}^{d_{y}}$,
$\underline{C}\leq g_{\theta}(y|x)\leq\overline{C}$. In addition for any
$(\theta,y)\in\Theta\times\mathbb{R}^{d_{y}}$,
$g_{\theta}(y|\cdot)\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$.

  * •

For any $(\theta,y)\in\Theta\times\mathbb{R}^{d_{y}}$,
$\frac{\partial}{\partial\theta^{(i)}}(\log\\{g_{\theta}(y|\cdot)\\})\in\mathcal{B}_{b}(\mathbb{R}^{d})\cap\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$,
$i\in\\{1,\dots,d_{\theta}\\}$.

  * •

For any $\theta\in\Theta$,

$\frac{\partial}{\partial\theta^{(i)}}(b_{\theta}^{(j)}),\frac{\partial}{\partial\theta^{(i)}}(b_{\theta}^{(j)})^{2},\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(k)}}(b_{\theta}^{(j)}),\frac{\partial^{2}}{\partial\theta^{(i)}\partial\theta^{(k)}}(b_{\theta}^{(j)})^{2}\in\mathcal{B}_{b}(\mathbb{R}^{d})\cap\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$

$(i,k,j)\in\\{1,\dots,d_{\theta}\\}^{2}\times\\{1,\dots,d\\}$.

The following result is proved in [18] and is Lemma 1 of that article.

######  Lemma A.1.

Assume (D1-2). Then for any
$(n,r,\theta,i)\in\mathbb{N}\times[1,\infty)\times\Theta\times\\{1,\dots,d_{\theta}\\}$,
there exists a $C<\infty$ such that for any
$(l,x)\in\mathbb{N}_{0}\times\mathbb{R}^{d}$

$\mathbb{E}_{\theta}[|G_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(i)}-G_{\theta}(\mathbf{X}_{T}^{x})^{(i)}|^{r}]^{1/r}\leq
C\Delta_{l}^{1/2}.$

We have the following result which can be proved using very similar arguments
to Lemma A.1.

######  Lemma A.2.

Assume (D1-2). Then for any
$(n,r,\theta,i,j)\in\mathbb{N}\times[1,\infty)\times\Theta\times\\{1,\dots,d_{\theta}\\}^{2}$,
there exists a $C<\infty$ such that for any
$(l,x)\in\mathbb{N}_{0}\times\mathbb{R}^{d}$:

$\displaystyle\mathbb{E}_{\theta}[|\varphi_{\theta}(\widetilde{X}_{1:n}^{x})G_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(i)}-\varphi_{\theta}(X_{1:n}^{x})G_{\theta}(\mathbf{X}_{T}^{x})^{(i)}|^{r}]^{1/r}$
$\displaystyle\leq$ $\displaystyle C\Delta_{l}^{1/2},$
$\displaystyle\mathbb{E}_{\theta}[|\varphi_{\theta}(\widetilde{X}_{1:n}^{x})H_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(ij)}-\varphi_{\theta}(X_{1:n}^{x})H_{\theta}(\mathbf{X}_{T}^{x})^{(ij)}|^{r}]^{1/r}$
$\displaystyle\leq$ $\displaystyle C\Delta_{l}^{1/2},$
$\displaystyle\mathbb{E}_{\theta}[|\varphi_{\theta}(\widetilde{X}_{1:n}^{x})G_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(i)}G_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(j)}-\varphi_{\theta}(X_{1:n}^{x})G_{\theta}(\mathbf{X}_{T}^{x})^{(i)}G_{\theta}(\mathbf{X}_{T}^{x})^{(j)}|^{r}]^{1/r}$
$\displaystyle\leq$ $\displaystyle C\Delta_{l}^{1/2}.$

###### Proof of Proposition 2.1.

In the following proof we will suppress the initial condition from the
notation. We have

$\mathfrak{H}_{\theta}^{l,(ij)}-\mathfrak{H}_{\theta}^{(ij)}=\sum_{j=1}^{3}T_{j},$

where

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(i)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})G_{\theta}(\mathbf{X}_{T})^{(i)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]}\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})G_{\theta}(\mathbf{X}_{T})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]},$
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(i)}G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})G_{\theta}(\mathbf{X}_{T})^{(i)}G_{\theta}(\mathbf{X}_{T})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]},$
$\displaystyle T_{3}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})H_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(ij)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})H_{\theta}(\mathbf{X}_{T})^{(ij)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]}.$

We remark that

$|\mathbb{E}_{\theta}[\varphi_{\theta}(\widetilde{X}_{1:n})]-\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]|\leq
C\Delta_{l}^{1/2},$ (A.1)

by using (D2) and convergence of Euler approximations of diffusions. For
$T_{1}$ we have

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\Bigg{(}\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(i)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})G_{\theta}(\mathbf{X}_{T})^{(i)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]}\Bigg{)}\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}+$
$\displaystyle\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})G_{\theta}(\mathbf{X}_{T})^{(i)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]}\Bigg{(}\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})G_{\theta}^{l}(\mathbf{X}_{T}^{l})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(\tilde{X}_{1:n})]}-\frac{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})G_{\theta}(\mathbf{X}_{T})^{(j)}]}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_{1:n})]}\Bigg{)}.$

So one can easily deduce by [18, Theorem 1] and (D2) that for some $C$ that
does not depend upon $l$

$T_{1}\leq C\Delta_{l}^{1/2}.$

Now for any real numbers, $a,b,c,d$ with $b,d$ non-zero we have the simple
identity

$\frac{a}{b}-\frac{c}{d}=\frac{a}{bd}[d-b]+\frac{1}{d}[a-c].$

So for $T_{2},T_{3}$ combining this identity with Lemma A.2, (A.1) and (D2)
one can easily conclude that for some $C$ that does not depend upon $l$

$\max\\{T_{2},T_{3}\\}\leq C\Delta_{l}^{1/2}.$

From here the proof is easily concluded. ∎

We end the section with a couple of results which are more-or-less direct
Corollaries of [18, Remarks 1 & 2]. We do not prove them.

######  Lemma A.3.

Assume (D1-2). Then for any
$(i,j,r,\theta)\in\\{1,\dots,d_{\theta}\\}^{2}\times[1,\infty)\times\Theta$,
there exists a $C<\infty$ such that for any
$(l,x,x_{\star})\in\mathbb{N}\times\mathbb{R}^{2d}$

$\displaystyle\mathbb{E}_{\theta}[|H_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(ij)}-H_{\theta}^{l-1}(\mathbf{X}_{T}^{l-1,x_{\star}})^{(ij)}|^{r}]^{1/r}$
$\displaystyle\leq$ $\displaystyle
C\Big{(}\Delta_{l}^{1/2}+\|x-x_{\star}\|_{2}\Big{)}$
$\displaystyle\mathbb{E}_{\theta}[|G_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(i)}G_{\theta}^{l}(\mathbf{X}_{T}^{l,x})^{(j)}-G_{\theta}^{l-1}(\mathbf{X}_{T}^{l-1,x_{\star}})^{(i)}G_{\theta}^{l-1}(\mathbf{X}_{T}^{l-1,x_{\star}})^{(j)}|^{r}]^{1/r}$
$\displaystyle\leq$ $\displaystyle
C\Big{(}\Delta_{l}^{1/2}+\|x-x_{\star}\|_{2}\Big{)}.$

## References

  * [1] N. Agarwal, B. Bullins and E. Hazan. Second-order stochastic optimization for machine learning in linear time. Journal of Machine Learning Research, 18, 1–40, 2017.
  * [2] C. Andrieu, A. Doucet, and R. Holenstein. Particle Markov chain Monte Carlo methods. J. R. Stat. Soc. Ser. B Stat. Methodol., 72(3):269–342, 2010.
  * [3] C. Andrieu, A. Lee, and M. Vihola. Uniform ergodicity of the iterated conditional SMC and geometric ergodicity of particle Gibbs samplers. Bernoulli, 24(2):842–872, 2018.
  * [4] A. Bain and D. Crisan.  Fundamentals of Stochastic Filtering.  Springer, New York, 2009.
  * [5] A. Beskos and G. O. Roberts. Exact simulation of diffusions. The Annals of Applied Probability, 15(4): 2422–2444, 2005.
  * [6] A. Beskos, O. Papaspiliopoulos, and G. O. Roberts. Retrospective exact simulation of diffusion sample paths with applications.  Bernoulli, 12(6):1077–1098, 2006.
  * [7] J. Bierkens, F. van der Meulen and M. Schauer. Simulation of elliptic and hypo-elliptic conditional diffusions. Adv. Appl. Probab., 52, 173–212, 2020.
  * [8] J. Blanchet and F. Zhang. Exact simulation for multivariate Ito diffusions. arXiv preprint arXiv:1706.05124, 2017
  * [9] R. H. Byrd, G. M. Chin, W. Neveitt, and J. Nocedal. On the Use of stochastic Hessian information in optimization methods for machine learning. SIAM J. Optim., 21(3), 977–995, 2011.
  * [10] O. Cappé, E. Moulines, and T. Ryden.  Inference in Hidden Markov Models. 580 Springer, New York, 2005.
  * [11] N. K. Chada, J. Franks, A Jasra, K. J. H. Law and M. Vihola. Unbiased inference for discretely observed hidden Markov model diffusions. _SIAM/ASA J. Unc. Quant._ , 9 (2), 763–787, 2021.
  * [12] S. Ditlevsen and A. Samson. Hypoelliptic diffusions: filtering and inference from complete and partial observations. J. R. Stat. Soc. Ser. B Stat. Methodol.,81(2): 361–384, 2019.
  * [13] P. Fearnhead, O. Papaspiliopoulos, and G. O. Roberts. Particle filters for partially observed diffusions. J. R. Stat. Soc. Ser. B Stat. Methodol., 70(4):755–777, 2008
  * [14] P. Fearnhead, O. Papaspiliopoulos, G. O. Roberts, and A. Stuart. Random-weight particle filtering of continuous time processes. J. R. Stat. Soc. Ser. B Stat. Methodol., 72(4):497–512, 2010.
  * [15] J. Gatheral.  The volatility surface: a practitioner’s guide. John Wiley & Sons, 2006.
  * [16] J. Gatheral, T. Jaisson, and M. Rosenbaum. Volatility is Rough. _Quant. Finance_ , 933–949, 2018.
  * [17] P. W. Glynn and C.H. Rhee. Exact estimation for Markov chain equilibrium expectations. Journal of Applied Probability, 51(A):377–389, 2014
  * [18] J. Heng, J. Houssineau and A. Jasra. On unbiased score estimation for partially observed diffusions. arXiv: 2105:04912, 2021.
  * [19] J. Heng, A. Jasra, K. J. Law, and A. Tarakanov. On unbiased estimation for discretized models. arXiv preprint arXiv:2102.12230, 2021.
  * [20] P. Jacob, F. Lindstein and T. Schön. Smoothing with couplings of conditional particle filters. _J. Amer. Statist. Assoc._ 115, 721–729, 2021.
  * [21] A. Jasra, K. J. H Law and F. Yu. Unbiased filtering for a class of partially observed diffusion processes. Adv. Appl. Probab., (to appear), 2021.
  * [22] A. Jasra, K. Kamatani, K. J. H. Law and Y. Zhou. Bayesian static parameter estimation for partially observed diffusions via multilevel Monte Carlo. _SIAM J. Sci. Comp._ , 40, A887-A902, 2018.
  * [23] P. E. Kloeden and E. Platen. Numerical solution of stochastic differential equations. volume 23. Springer Science & Business Media, 2013
  * [24] A. Majda and X. Wang.  Non-linear Dynamics and Statistical Theories for Basic Geophysical Flows, Cambridge University Press,  2006\.
  * [25] D. Mcleish. A general method for debiasing a Monte Carlo. Monte Carlo Methods and Applications, 17: 301–315, 2011.
  * [26] F. van der Meulen and M. Schauer. Bayesian estimation of incompletely observed diffusions.  Stochastic, 90(5), 641–662, 2018.
  * [27] C. H. Rhee and P. Glynn. Unbiased estimation with square root convergence for SDE models. _Op. Res._ 63, 1026–1043, 2016.
  * [28] L. M. Ricciardi.  Diffusion Processes and Related Topics in Biology. Springer, Lecture Notes in Biomathematics, 1977.
  * [29] M. Schauer, F. Van Der Meulen, and H. Van Zanten. Guided proposals for simulating multi-dimensional diffusion bridges.  Bernoulli, 23(4A):2917–2950, 2017.
  * [30] S. E. Shreve.  Stochastic calculus for finance II: Continuous-time model. Springer Science & Business Media, 2004.
  * [31] H. Thorisson Coupling, stationarity, and regeneration. Springer: New York, 2002.
  * [32] M. Vihola. Unbiased estimators and multilevel Monte Carlo. _Op. Res._ , 66, 448–462, 2018.