Online Smoothing for Diffusion Processes Observed with Noise
Shouto Yonekura Alexandros Beskos
stat.CO stat.ME


#1 1



0
  
  
  Online Smoothing for Diffusion Processes
  Observed with Noise


    Shouto Yonekura
    SY acknowledges funding from the Alan Turing Institute under grant number TU/C/000013.


    Department of Statistical Science, University College London, 

    Alan Turing Institute,

    and 

    Alexandros BeskosAB acknowledges funding via a Leverhulme Trust Prize. 

    Department of Statistical Science, University College London,

    Alan Turing Institute,
    December 30, 2023
=====================================================================================================================================================================================================================================================================================================================================================================================================

 

1
  

    Online Smoothing for Discretely Observed Jump Diffusions




   


We introduce a methodology for online estimation of smoothing expectations 
for a class of additive functionals, in the context of a rich family of diffusion processes 
(that may include jumps) – observed at discrete-time instances. We overcome the unavailability of the transition density of the underlying SDE by working on the augmented pathspace. 
The new method can be applied, for instance, to carry out online parameter inference for the designated class of models. 
Algorithms defined on the infinite-dimensional pathspace have been developed the last years mainly in the context of MCMC techniques. There, the main  benefit is the achievement of mesh-free mixing times for the practical time-discretised algorithm used on a PC.  Our own methodology sets up the framework for infinite-dimensional online filtering –  an important positive practical consequence is the construct of estimates with variance that does not increase with decreasing mesh-size.  
Besides regularity conditions, our method is, in principle, applicable under the weak assumption  – relatively to restrictive conditions often required in the MCMC or filtering literature of methods defined on  pathspace – that the SDE covariance matrix is invertible.



Keywords:  Jump Diffusion; Data Augmentation; Forward-Only Smoothing; Sequential Monte Carlo; Online Parameter Estimation.






















1.5 




§ INTRODUCTION


Research in Hidden Markov Models (HMMs) has – thus far – provided effective online algorithms for the estimation of expectations of the smoothing distribution for the case of a class of additive functionals of the underlying signal. Such methods necessitate knowledge of the transition density of the Markovian part of the model between observation times. We carry out a related exploration for the (common in applications) case when the signal corresponds to a diffusion process, thus we are faced with the challenge that such transition densities are typically unavailable. Standard data augmentation schemes that work with the multivariate density of a large enough number of imputed points of the continuous-time signal will lead to ineffective algorithms. The latter will have the abnormal characteristic that – for given Monte-Carlo iterates – the variability of the produced estimates will increase rapidly as the resolution of the imputation becomes finer. One of the ideas underpinning the work in this paper is that development of effective algorithms instead requires respecting the structural properties of the diffusion process, thus we build up imputation schemes on the infinite-dimensional diffusion pathspace itself. As a consequence, the time-discretised algorithm used in practice on a PC will be stable under mesh-refinement.    

We consider continuous-time jump-diffusion models observed at discrete-time instances. 
The d_x-dimensional process, X={X_t;t≥ 0}, d_x≥ 1, is defined via
the following time-homogeneous stochastic differential
equation (SDE), with X_t-:=lim_s↑ tX_t,


    dX_t =b(X_t-)dt+σ(X_t-)dW_t+dJ_t,   X_0=x_0∈ℝ^d_x,   t≥ 0.


Solution X is driven by the 
d_w-dimensional Brownian motion, W={W_t;t≥ 0}, d_w≥ 1, and the compound Poisson process,
J={J_t;t≥ 0}. The SDE involves a drift function b=b_θ:ℝ^d_x↦ℝ^d_x and coefficient matrix 
σ=σ_θ:ℝ^d_x↦ℝ^d_x× d_w, for parameter θ∈ℝ^p, p≥ 1.
Let {N_t;t≥ 0}
be a Poisson process with intensity function λ_θ(·), and 
 {ξ_k}_k≥ 1 i.i.d. sequence of  
 random variables with Lebesgue density h_θ(·); the càdlàg process J is determined as J_t= J_t,θ=∑_i=1^N_tξ_i.

We work under standard assumptions (e.g. linear growth, Lipschitz
continuity for b, σ) that guarantee a unique global solution of (<ref>), in a weak or strong sense,
see e.g. <cit.>.  
 

SDE (<ref>) is  observed with noise at 
 discrete-time instances
0 = t_0 < t_1<t_2<⋯ <t_n, n≥ 1. 
Without loss of generality, we assume equidistant observation times, 
with Δ:= t_1-t_0.
We consider data Y_t_0,…, Y_t_n, and for simplicity we set,


    x_i:=X_t_i,    y_i=Y_t_i,     0≤ i≤ n.


Let ℱ_i be a filtration generated by X_s for s∈[t_i-1,t_i], 0< i≤ n. We assume,


    [ Y_t_i | {Y_t_j;j< i}, {X_s;s∈[0,t_i]} ] ∼
    g_θ(dY_t_i | Y_t_i-1,ℱ_i),    0≤ i ≤ n,


for conditional distribution g_θ(·|Y_t_i-1,ℱ_i) on ℝ^d_y, d_y≥ 1, under convention Y_t_-1=y_-1=∅.


We write, 


    [ x_i |  x_i-1 ] ∼ f_θ(dx_i|x_i-1),

where f_θ(dx_i|x_i-1) is the transition distribution of the driving SDE process (<ref>).
We consider the density functions of g_θ(dy_i|y_i-1,ℱ_i) 
and f_θ(dx_i|x_i-1), and – with some abuse of notation – we  
write g_θ(dy_i|y_i-1,ℱ_i)= g_θ(y_i|y_i-1,ℱ_i)dy_i, f_θ(dx_i|x_i-1)= f_θ(x_i|x_i-1)dx_i, where dy_i, dx_i denote Lebesgue measures. Our work develops under the following regime.



The transition density f_θ(x'|x) is intractable;
the density g_θ(y'|y,ℱ) is analytically available –
for appropriate x', x, y', y, ℱ consistent with  preceding definitions.


The intractability of the transition density f_θ(·|·) will pose challenges for the main  problems this paper aims to address.

Models defined via (<ref>)-(<ref>) are extensively used, e.g.,  in finance and
econometrics, for instance for capturing the market microstructure noise,
see <cit.>.
The above setting belongs to the general class of HMMs, with a signal defined in continuous-time. See <cit.>
for a general treatment of HMMs fully specified in discrete-time. A number of methods have been suggested in the literature 
for approximating the unavailable transition density – mainly in the case of processes without jumps – including:
asymptotic expansion techniques <cit.>;
martingale estimating functions <cit.>;
generalized method of moments <cit.>; Monte-Carlo approaches <cit.>.
See, e.g., <cit.> for a detailed review.


For a given sequence {a_m}_m, we use the notation a_i:j:=(a_i,…,a_j), 
for integers i≤ j. 
Let p_θ(y_0:n) denote the joint density of y_0:n. 
Throughout the paper, p_θ(·) is used generically to represent probability distributions or densities of random variables appearing 
as its arguments.

Consider the maximum likelihood estimator (MLE),


    θ̂_n    :=max_θ∈Θlog p_θ(y_0:n).


Except for limited cases, one cannot obtain the MLE analytically
for HMMs (even for discrete-time signal) due to the intractability of 
p_θ(y_0:n).



We have set up the modelling context for this work. 
The main contributions of the paper 
in this setting – several of which relate with overcoming the intractability of the transition density of the SDE, and developing a methodology that is well-posed on the infinite-dimensional pathspace – will be as follows:



  (i) We present an online algorithm that delivers Monte-Carlo estimators of smoothing expectations,


    𝒮_θ,n=∫ S_θ(𝐱_0:n) p_θ(d𝐱_0:n|y_0:n),    n≥ 1,


 for the class of additive functionals S_θ(·) of the structure, 
 

    S_θ(𝐱_0:n)=∑_k=0^ns_θ,k(x_k-1,𝐱_k),
 

under the conventions x_-1=∅, 𝐱_0 = x_0.
The bold type notation 𝐱_k, k≥ 0, is reserved for carefully defined variables involving elements of the infinite-dimensional pathspace. The specific construction of 𝐱_k will depend on the model at hand, and will be explained in the main part of the paper. We note that the class of additive functionals considered supposes additivity over time, and therefore includes the general case of evaluating the score function.
The online solution of the smoothing problem is often used as the means to solve some concrete inferential problems - see, e.g., (ii) below.


  (ii) We take advantage of the new approach to show numerical applications, with emphasis on
carrying out online parameter inference for the designated class of models via a gradient-ascent  approach (in a Robbins-Monro stochastic gradient framework). A critical aspect of this particular online algorithm (partly likelihood based, when concerned with parameter estimation; partly Bayesian, with regards to identification of filtering/smoothing expectations) is that it delivers estimates of the evolving score function,  of the model parameters, together with particle representations of the filtering distributions, through a single passage of the data. 
This is a unique favourable algorithmic characteristic, when constrasted with alternative algorithms with similar objectives, such as, e.g.,  Particle MCMC <cit.>, or SMC^2 <cit.>.   


  (iii) In this work, we will not characterise analytically the size of the
time-discretisation bias relevant to the SDE models at hand, and are content that: (I) the bias can be decreased by increasing the resolution of 
the numerical scheme (typically an Eyler-Maruyama one,
or some other Taylor scheme, see e.g. <cit.>);
(II) critically, the Monte-Carlo algorithms are developed in a manner that the variance does not increase (in the limit, up to infinity) when increasing the resolution of the time-discretisation method; to achieve such an effect, the  
algorithms are (purposely) defined on the infinite-dimensional pathspace, and SDE paths are only discretised when implementing the algorithm on a PC (to allow, necessarily, for finite computations).


  (iv) Our method draws inspiration from earlier works, in the context of online filtering for discrete-time HMMs and infinite-dimensional pathspace MCMC methods. The complete construct is novel; one consequence of this is that it is applicable, in principle, for a wide class of SDEs, under the following weak assumption (relatively to restrictive conditions often imposed in the literature of infinite-dimensional MCMC methods).



The diffusion covariance matrix function,


    Σ_θ(v):=σ_θ(v)σ_θ(v)^⊤∈ℝ^d_x× d_x,

 is invertible, for all relevant v, θ.



Thus, the methodology does not apply as defined here only for the class of hypoelliptic SDEs. 

An elegant solution to the online smoothing problem posed above in (i), for the case of a standard HMM with discrete-time 
signal of known transition density f_θ(x'|x), is given in <cit.>. Our own work overcomes 
 the unavailability of the transition density in the continuous-time scenario by following the above literature but augmenting the hidden state with the complete continuous-time SDE path. Related augmentation approaches in this setting – though for different inferential objectives – have appeared in 
 <cit.>, where the auxiliary variables are derived via the Poisson estimator of transition densities for SDEs  (under strict conditions on the class of SDEs; no jumps), introduced in <cit.>,
 and in <cit.> where the augmentation involves indeed the continuous-time path (the objective therein is to solve the filtering problem and the method is applicable for SDEs with additive Wiener noise; no jump processes are considered).   

A Motivating Example: 
Fig. <ref> (Top-Panel) shows
estimates of the score function, evaluated at the true parameter value θ=θ^†, for parameter θ_3 
of the
Ornstein–Uhlenbeck (O–U) process, 
dX_t = θ_1(θ_2-X_t)dt + θ_3 dW_t, X_0=0.0, for n=10 observations y_i = x_i + ϵ_i, ϵ_ii.i.d∼𝒩(0,0.1^2), 1≤ i ≤ n. Data were simulated from 
θ^†=(0.5,0.0,0.4) with an Euler-Maruyama scheme of M^†=1,000 grid points per unit of time.
Fig. <ref> (Top-Panel) illustrates the `abnormal' effect of a standard data-augmentation scheme, where for N=100 particles, the Monte-Carlo method (see later sections for details) produces estimates of increasing variability as algorithmic resolution increases with 
M=10, 50, 100, 200 – i.e., as it approaches the `true' resolution used for the data generation. In contrast, Fig. <ref> (Bottom Panel) shows results for the same score function as estimated by our proposed method in this paper. 
Our method is well-defined on the infinite-dimensional pathspace, and, consequently, it manages to robustly estimate  the score function under mesh-refinement.





The rest of the paper is organised as follows.
Section <ref> reviews the Forward-Only algorithm for the online approximation of expectations of a class of additive functionals described in <cit.>.
Section <ref> sets up the framework for the treatment of pathspace-valued SDEs, first for the conceptually simpler case of SDEs without jumps, and then proceeding with incorporating jumps.
Section <ref> provides the complete online approximation algorithm, constructed on the infinite-dimensional pathspace.
Section <ref> discusses the adaptation of the developed methodology for the purposes of online inference
for unknown parameters of the given SDE model. 
Section <ref> shows numerical applications of the developed methodology.
Section <ref> contains conclusions and directions for future research. 











§ FORWARD-ONLY SMOOTHING











A bootstrap filter <cit.> is applicable
in the continuous-time setting, as it only requires forward sampling of the 
underlying signal {X_t}_t≥ 0; this is trivially possible – under numerous approaches – and is typically associated with the introduction of some time-discretisation bias <cit.>. However, the transition density is still required for the smoothing problem we have posed in the Introduction. In this section, we assume 
a standard discrete-time HMM, with initial distribution p_θ(dx_0), transition density f_θ(x'|x), 
and likelihood g_θ(y|x), for appropriate x, x'∈ℝ^d_x, 
y∈ℝ^d_y,
and review the online algorithm developed in 
<cit.> for this setting. Implementation of the bootstrap filter provides 
an approximation of the smoothing distribution p_θ(dx_0:n|y_0:n)
by following the geneology of the particles. This method is studied, e.g., in <cit.>.

Let {x_0:n^(i),W_n^(i)}_i=1^N, N≥ 1,
be a particle approximation of the smoothing distribution p_θ(dx_0:n|y_0:n), in
the sense that we have the estimate,




    p_θ(dx_0:n| y_0:n)=∑_i=1^NW_n^(i)δ_x_0:n^(i)(dx_0:n),  ∑_i=1^NW_n^(i)=1,




with δ_x_0:n^(i)(dx_0:n) 
 the Dirac
measure with an atom at x_0:n^(i). 


Then, replacing  p_θ(x_0:n|y_0:n) with 
its estimate in (<ref>) provides consistent estimators of expectations of the HMM smoothing distributions. 

Though the
method is online and the computational cost per time step is 𝒪(N), it typically suffers from
the well-documented path-degeneracy problem –
as illustrated via theoretical results or numerically <cit.>. 
That is, as n increases, the particles 
representing p_θ(dx_0:n| y_0:n) obtained by the
above method will eventually all share the same ancestral particle due
to the resampling steps, and the approximation
collapses for big enough n. This is well-understood not to be a solution to the approximation of the smoothing distribution for practical applications. 

An approach which overcomes path-degeneracy is the  
Forward Filtering Backward Smoothing (FFBS) algorithm of <cit.>. 
We briefly review the method here, following closely the notation and development in 
<cit.>. In the forward direction, assume that a filtering algorithm (e.g. bootstrap) has 
provided a particle approximation of the filtering distribution p_θ(dx_k-1|y_0:k-1) – assuming a relevant k –,


    p_θ(dx_k-1|y_0:k-1)= ∑_i=1^NW_k-1^(i)δ_x_k-1^(i)(dx_k-1),


for weighted particles {x_k-1^(i),W_k-1^(i)}_i=1^N.
In the backward direction, assume that one is given the particle approximation 
of the marginal smoothing distribution p_θ(dx_k|y_0:n), 


    p_θ(dx_k|y_0:n) = ∑_i=1^NW_k|n^(i)δ_x_k^(i)(dx_k).


One has that <cit.>,


    p_θ(dx_k-1:k| y_0:n)    = p_θ(dx_k|y_0:n)⊗ p_θ(dx_k-1|x_k,y_0:k-1) 
       = p_θ(dx_k|y_0:n)⊗ p_θ(dx_k-1| y_0:k-1)f_θ(x_k| x_k-1)/∫ f_θ(x_k| x_k-1)p_θ(x_k-1|y_0:k-1)dx_k-1.


Using (<ref>)-(<ref>), and based on equation (<ref>),
we obtain the 
approximation, 


    p_θ(dx_k-1:k| y_0:n) = ∑_j=1^NW_k|n^(j)∑_i=1^Nf_θ(x_k^(j)|x_k-1^(i))W_k-1^(i)/∑_l=1^Nf_θ(x_k^(j)| x_k-1^(l))W_k-1^(l)δ_(x_k-1^(i),x_k^(j))(dx_k-1:k).


Recalling the expectation of additive functionals in (<ref>)-(<ref>) – where now, in the discrete-time setting, we can ignore the bold 𝐱_k elements, and 
simply use x_k instead – 
the above calculations give rise to the following estimator of the target quantity 𝒮_θ,n in (<ref>),


    𝒮_θ,n = ∑_k=0^n∫ s_θ,k(x_k-1,x_k)  p_θ(dx_k-1:k| y_0:n).


To be able to apply the above method, the marginal smoothing approximation in (<ref>) is 
obtained via a backward recursive approach. In particular, 
starting from k=n (where the approximation is provided by the standard forward particle filter), one proceeds as follows. 
Given k, the quantity for k-1 is directly obtained 
by integrating out x_k in (<ref>), 
thus we have, 


    p_θ(dx_k-1|y_0:n) = ∑_i=1^N W_k-1|n^(i)δ_x_k-1^(i)(dx_k-1),


for the normalised weights, 

    W_k-1|n^(i)∝∑_j=1^NW_k|n^(j) f_θ(x_k^(j)|x_k-1^(i))W_k-1^(i)/∑_l=1^Nf_θ(x_k^(j)| x_k-1^(l))W_k-1^(l).




Notice that – in this version of FFBS – 
the same particles {x_k^(i)}_i=1^N are used in both directions (the ones before resampling at the forward filter), but with different weights.

An important development made in <cit.> is transforming the  above offline algorithm 
 into an online one. This is achieved by consideration of the sequence of  instrumental functionals,   


    T_θ,0(x_0)= s_θ,0(x_0);       
    T_θ,n(x_n):=  ∫ S_θ,n(x_0:n)p_θ(dx_0:n-1|y_0:n-1,x_n), 
       n≥ 1.


Notice that, first,


    𝒮_θ,n
    
      =∫ T_θ,n(x_n)p_θ(dx_n| y_0:n).



We also have that, 


    T_θ,n(x_n)    =∫[T_θ,n-1(x_n-1)+s_θ,n(x_n-1,x_n)]p_θ(dx_n-1|y_0:n-1,x_n) 
       ≡∫[T_θ,n-1(x_n-1)+s_θ,n(x_n-1,x_n)]f_θ(x_n|x_n-1)p_θ(
    dx_n-1|y_0:n-1)/∫ f_θ(x_n|x_n-1)p_θ(dx_n-1|y_0:n-1),


– see Proposition 2.1 of <cit.> for a (simple) proof. 
This recursion provides an online – forward-only – 
advancement of FFBS for estimating the smoothing 
expectation of additive functionals. The complete method is summarised in  Algorithm <ref>: one key ingredient is that, during the recursion, 
values of the functional T_θ,n(x_n) 
are only required at the discrete positions x_n^(i) determined by the forward particle filter.





In the SDE context, under Assumption <ref>, the transition
density f_θ(·|·) is considered intractable, 
thus Algorithm <ref> – apart from serving as a review of the method in <cit.> – does not appear to be practical in 
the continuous-time case. 






§ DATA AUGMENTATION ON DIFFUSION PATHSPACE




To overcome the intractability of the transition
density f_θ(·|·) of the SDE, we will work with an algorithm that is defined in 
continuous-time and makes use of the complete SDE path-particles in its development. 
The new method has connections with earlier works in the literature. <cit.> 
focus on the filtering problem for a class of models related to (<ref>)-(<ref>), and 
come up with an approach that requires the complete SDE path, for a limited class of diffusions with additive noise and no jumps. <cit.> also deal with the filtering problem, and – equipped with an unbiased estimator of the unknown transition density – recast the 
problem as one of filtering over an augmented space that incorporates the randomness for the unbiased estimate. The method is accompanied by strict conditions on the drift and diffusion coefficient (one should be able to transform the SDE – no jumps – into one of unit diffusion coefficient; the drift of the SDE must have a gradient form). Our contribution requires, in principle, solely the diffusion coefficient invertibility Assumption <ref>; arguably, the weakened condition we require stems from the fact that our approach appears as the relatively 
most `natural' extension (compared to alternatives) of the standard discrete-time algorithm of <cit.>. 

The latter discrete-time 
method requires the density  f_θ(x'|x)=f_θ(dx'|x)/dx'.
In continuous-time, we obtain
an analytically available Radon-Nikodym derivative of p_θ(d𝐱'|x), 
for a properly defined variate 𝐱' that involves information about the 
continuous-time path for moving from x to x' within time Δ. 
We will give the complete algorithm in Section <ref>. 
In this section, 
we prepare the ground via carefully determining 𝐱' given x, and calculating the 
relevant densities to be later plugged in into our method.



 §.§ SDEs with Continuous Paths



We work first with the process with continuous sample-paths, i.e. of dynamics,

    dX_t =b_θ(X_t)dt+σ_θ(X_t)dW_t.


We adopt an approach motivated by techniques used for MCMC
algorithms <cit.>.
Assume we are given starting point x∈ℝ^d_x, ending point x'∈ℝ^d_x,
and the complete continuous-time path  for the signal process in (<ref>) on 
[0,T], for some T≥ 0. That is, we now work with the path process,


    [ {X_t;t∈[0,T]} |  X_0=x,X_T=x' ].


Let ℙ_θ,x,x' denote the probability distribution of the pathspace-valued variable in (<ref>).

We consider
the auxiliary bridge process X̃={X̃_t;t∈[0,T]} defined as, 


    dX̃_t    ={ b_θ(X̃_t)+x'-X̃_t/T-t} dt+σ_θ(X̃_t)dW_t,  X̃_0=x,   t∈[0,T],


with corresponding probability distribution ℚ_θ,x,x'. Critically,
a path of X̃ starts at point x and finishes at x', w.p. 1.

Under regularity conditions, <cit.> prove that probability measures ℙ_θ,x,x', ℚ_θ,x,y are absolutely continuous with respect to each other. We treat the auxiliary SDE (<ref>) as a transform
from the driving noise to the solution, whence
a sample path, X, of the process 
X̃={X̃_̃t̃; t∈[0,T]}, 
is produced by a mapping – determined by (<ref>) – of a corresponding sample path, say Z, of the 
Wiener process. That is, we have set up a map, and – under Assumption <ref> – its inverse,   


    Z ↦ X=: F_θ(Z;x,x'),    Z =F^-1_θ(X;x,x').


More analytically, F_θ^-1 is given via the transform,


    dZ_t    =σ_θ(X_t)^-1{dX_t-b_θ(X_t)dt-x'-X_t/T-tdt}.


In this case we define,


    𝐱' := (x',Z),


and the probability measure of interest is, 

    p_θ(d𝐱'|x) := f_θ(dx'|x)⊗ p_θ(dZ|x',x).
































Let 𝕎 be the standard
Wiener probability measure on [0,T]. 
Due to the 1–1 transform, we have that, 


    p_θ(dZ|x',x)/𝕎(dZ)≡dℙ_θ,x,x'/dℚ_θ,x,x'(F_θ(Z;x,x')),


so it remains to obtain the density dℙ_θ,x,x'/
dℚ_θ,x,x'. 



Such a Radom-Nikodym derivative has been object of interest in many works. <cit.> provided detailed conditions and a proof, but (seemingly) omit an expression for the normalising constant which is important in our case, as it involves the parameter θ – in our applications later in the paper, we aim to infer about θ.
<cit.> provide an expression based on a conditioning argument for the projection of the probability measures on [0,t), t<T, and passage to the limit t↑ T. 
The derivations in <cit.> are extremely rigorous,


so we will make use the expressions in that paper. Following carefully the proofs of some of their main results (Theorem 5, together with Lemmas 7, 8) one can indeed retrieve the constant in the deduced density. In particular,
<cit.> impose the following conditions.




  (i) SDE (<ref>) admits a strong solution;

  (ii) v↦σ_θ(v) is in  𝒞^2_b (i.e., twice continuously differentiable and bounded, with bounded first and second derivatives), and it is invertible, with bounded inverse;

  (iii) v↦ b_θ(v) is locally Lipschitz, locally bounded.




Under Assumption <ref>,
<cit.> prove that,  






    dℙ_θ,x,x'/dℚ_θ,x,x'(X) =  |Σ_θ(x')|^1/2/|Σ_θ(x)|^1/2×𝒩(x';x, TΣ_θ(x))/f_θ(x'|x)×φ_θ(X;x,x'),


where 𝒩(v;μ,V) is the density function of the Gaussian law on ℝ^d_x 
with mean μ, variance V, evaluated at v, 
|·| is matrix determinant, 
and φ_θ(X;x,x') is such that, 


    logφ_θ(X;x,x') = ∫_0^T   ⟨ b_θ(X_t),Σ^-1_θ(X_t)dX_t ⟩ 
     - 12∫_0^T⟨  b_θ(X_t),  Σ_θ^-1(X_t)b_θ(X_t)dt ⟩
       - 12∫_0^T⟨  (x'-X_t), dΣ^-1_θ(X_t)(x'-X_t) ⟩T-t 
    -12∫_0^T∑_i,j=1^d_xd [ Σ_θ,ij^-1, (x'_i-X_t,i)(x'_j-X_t,j) ]T-t.


Here, [·,·] denotes the quadratic variation process for semi-martingales; also, ⟨·, ·⟩ is the standard inner-product on ℝ^d_x. We note that 
transforms different from (<ref>) have been proposed in the literature  <cit.>
to achieve the same effect of obtaining an 1–1 mapping of the latent path that has a density with respect to a measure that does not depend on x, x' or θ.
However, such methods are mostly applicable for scalar diffusions <cit.>. Auxiliary variables involving a random, finite selection of points of the latent path, based on 
the (generalised) Poisson estimator of <cit.> are similarly
restrictive. In contrast to other attempts, our
methodology may be applied for a much more general class of SDEs, as determined by Assumption <ref> – and further regularity conditions, e.g. as in Assumption <ref>. Thus, continuing from (<ref>) we have obtained that, 


    p_θ(d𝐱'|x)/(Leb^⊗ d_x⊗𝕎)(d𝐱')   = φ_θ(F_θ(Z;x,x');x,x')×𝒩(x';x,TΣ_θ(x))×|Σ_θ(x')|^1/2/|Σ_θ(x)|^1/2
       =: p_θ(𝐱'|x) ≡ p_θ(𝐱'|x; T).


We have added the extra argument involving the length of path, T, in the last expression, as it will be of use in the next section. 


A critical point here is that the above density is analytically tractable, thus by working on pathspace we have overcome the 
unavailability of the transition density f_θ(x'|x).






 §.§ SDEs with Jumps



We extend the above developments to
the more general case of the  d_x-dimensional jump diffusion model given in (<ref>), which we re-write here for convenience,


    dX_t    =b_θ(X_t-)dt+σ_θ(X_t-)dW_t+dJ_t,   X_0=x∈ℝ^d_x,   t∈[0,T].


Recall that J_θ= J={ J_t}
denotes a compound Poisson process with jump intensity λ_θ(·) and
jump-size density h_θ(·).

Let 𝔽_θ,x(·) denote the law of the unconditional process
(<ref>) and 𝕃_θ the law of the involved compound Poisson process.






We write J=((τ_1,b_1), …, (τ_κ,b_κ))
to denote the jump process, 
where {τ_i} are the times of events, {b_i} the jump sizes
and κ≥ 0 the total number of events. 
In addition, we consider the reference 
measure 𝕃, corresponding to unit rate Poisson process measure on [0,T] multiplied with ⊗_i=1^κ+1Leb^⊗ d_x.







  §.§.§ Construct One

We consider the random variate,


    𝐱' = ( J, {x_τ_i-}_i=1^κ+1, {Z(i)}_i=1^κ+1),


under the conventions x_τ_0-≡ x, x_τ_κ+1-≡ x', 
where we have defined,


    Z(i) = F^-1_θ(X(i) ; x_τ_i-1,x_τ_i-),    
    X(i) := {X_t;t∈ [x_τ_i-1, x_τ_i)},   
    1≤ i ≤κ+1.


We have that, 

    p_θ(d𝐱'|x) := 
    𝕃_θ(dJ)⊗[ ⊗_i=1^κ+1{ f_θ(dx_τ_i-|x_τ_i-1)⊗ p_θ(dZ(i)|x_τ_i-1,x_τ_i-) }].


Using the results about SDEs without jumps in Section <ref>,
and in particular expression (<ref>), 
upon defining, 


    𝐱'(i):=(x_τ_i-, Z(i)),


we have that – under an apparent adaptation of Assumption <ref>,


    f_θ(dx_τ_i-|x_τ_i-1)⊗ p_θ(dZ(i)|x_τ_i-1,x_τ_i-)/Leb^⊗ d_x(dx_τ_i-)⊗𝕎(dZ(i)))
     = 
    p_θ(𝐱'(i)|x_τ_i-; τ_i-τ_i-1),


with the latter density p_θ(𝐱'(i)|x_τ_i-) 
determined as in (<ref>), given clear adjustments.
Thus, the density of p_θ(d𝐱'|x) with respect to the reference 
measure,

    μ(d𝐱'):=𝕃(dJ)⊗[  ⊗_i=1^κ+1{Leb^⊗ d_x(dx_τ_i-)⊗𝕎(dZ(i))}],


is equal to, 

    p_θ(d𝐱'|x)/μ(d𝐱')=
    e^-∫_0^Tλ_θ(t)dt/e^-T·∏_i=1^κ{λ_θ(τ_i))h_θ(b_i)}×∏_i=1^κ+1
    
    p_θ(𝐱'(i)|x_τ_i-;τ_i-τ_i-1).





































  §.§.§ Construct Two


We adopt an idea used – for a different problem – in 
<cit.>. 






Given x,x'∈ℝ^d_x, we define an auxiliary process X̃_t as follows,


    dX̃_t     ={ b_θ(X̃_t)+x'-J_T-X̃_t+J_t/T-t} dt+σ_θ(X̃_t)dW_t+dJ_t,   X_0=x,


so that X̃_T=x', w.p. 1.
As with (<ref>), we view (<ref>)
as a transform, projecting a path, Z of the Wiener process and the compound process, J,  onto a path, X, of the jump process. That is, we consider the 1–1 maps,


    (J,Z) ↦ X=: G_θ(J,Z;x,x'),    (J,Z) =G^-1_θ(X;x,x').


Notice that for the inverse transform, the J-part is obtained immediately, whereas for the Z-part one uses 
the expression
– well-defined due to Assumption <ref> –,



    dZ_t    =σ_θ(X_t)^-1{dX_t- dJ_t - b_θ(X_t)dt-x'-J_T-X_t+J_t/T-tdt}.



We denote by ℙ_θ,x,x'
the law of
the original process in (<ref>) conditionally 
on hitting x' at time T. 
Also, we denote the distribution on pathspace 
induced by (<ref>) as ℚ_θ,x,x'.  
Consider the variate, 


    𝐱' = ( x', J, Z ),


so that, 


    p_θ(d𝐱'|x) := f_θ(dx'|x)⊗ p_θ(d(J,Z)|x',x).


Due to the employed 1–1 transforms, we have that, 


    p_θ(d𝐱'|x)/(Leb^⊗ d_x⊗𝕃_θ⊗𝕎)(d𝐱')   = f_θ(x'|x)×dℙ_θ,x,x'/dℚ_θ,x,x'(G_θ(J,Z;x,x')).


Thus, using the parameter-free reference measure, 


    μ(d𝐱'):= Leb^⊗ d_x⊗𝕃⊗𝕎,


one obtains that, 

    p_θ(d𝐱'|x)/μ(d𝐱')    = f_θ(x'|x)×e^-∫_0^Tλ_θ(t)dt/e^-T·∏_i=1^κ{λ_θ(τ_i))h_θ(b_i)}
                                   ×dℙ_θ,x,x'/dℚ_θ,x,x'(F_θ(J,Z;x,x')).














































































<cit.> obtained the Radon-Nikodym derivative expression in  (<ref>) after a great amount of rigorous analysis.
A similar development for the case of conditioned jump diffusions does not follow from their work, and can only be subject of dedicated research at the scale of a separate paper. This is beyond the scope of our work. In practice, one can proceed as follows. 
For grid size M≥ 1, and δ=T/M, let  ℙ^M_θ,x,x'(X_δ,…, X_(M-1)δ | X_0=x, X_Mδ=x')
denote the time-discretised Lebesgue density of the (M-1)-positions of the 
conditioned diffusion with law ℙ_θ,x,x'.
Once (<ref>) is obtained, a time-discretisation approach will give,

    ℙ^M_θ,x,x'(X_δ,…,    X_(M-1)δ | X_0=x, X_Mδ=x')  
       =ℙ_θ,x,x'^M(X_δ,…, X_(M-1)δ, X_Mδ=x' | X_0=x)/ℙ^M_θ,x,x'
    (X_Mδ=x' | X_0=x)
    .

In this time-discretised setting, f_θ(x'|x) in (<ref>) will be
replaced by ℙ^M_θ,x,x'
(X_Mδ=x' | X_0=x). Thus, the intractable transition density over the complete time period will  cancel out, and one is left with an explicit expression to use on a PC. 
Compared to the method in Section <ref>, and the Construct One
in the current section, we do not have explicit theoretical evidence of a density on the pathspace. Yet, all numerical experiments we tried showed that the deduced algorithm was stable under mesh-refinement.  
We thus adopt the approach (or, conjecture) that the density in (<ref>) exists – under assumptions –, and can be obtained pending future research.

















§ FORWARD-ONLY SMOOTHING FOR SDES





 §.§ Pathspace Algorithm


We are ready to develop a forward-only particle smoothing 
method, under the scenario in (<ref>)-(<ref>) on
the pathspace setting.
We will work with the pairs of random elements

    (x_k-1,𝐱_k),    1≤ k≤ n,


with 𝐱_k as defined in Section <ref>, i.e. containing pathspace elements,
and given by an 1–1 transform of {X_s;s∈[t_k-1,t_k]}, such that we can obtain 
a density for p_θ(𝐱_k|x_k-1) with respect to 
a reference measure that does not involve θ. Recall that  p_θ(d𝐱_k|x_k-1) denotes the
probability law for the augmented variable 𝐱_k given 
x_k-1. 

We also write the corresponding density as, 


    p_θ(𝐱_k|x_k-1) := p_θ(d𝐱_k|x_k-1)/μ(d𝐱_k).



The quantity of interest is now,

    𝒮_θ,n=∫ S_θ(𝐱_0:n) p_θ(d𝐱_0:n|y_0:n),    n≥ 1,


for the class of additive functionals S(·) of the structure,


    S_θ(𝐱_0:n)=∑_k=0^ns_θ,k(x_k-1,𝐱_k),
 

under the convention that x_-1=∅.


Notice that we now allow s_k(·,·) to be a function of 
x_k-1 and 𝐱_𝐤; thus, s_k(·,·) can potentially correspond to integrals, or other pathspace functionals. 

We will work with a transition density on the enlarged space of 𝐱_k.


Similarly to the discrete-time case in Section <ref>, we define the functional,


    T_θ,n(𝐱_n)     :=∫ S_θ(𝐱_0:n) p_θ(d𝐱_0:n-1| y_0:n-1,𝐱_n).




We have that,

    𝒮_θ,n    =∫ T_θ,n(𝐱_n) p_θ(d𝐱_n| y_0:n).




See Appendix <ref>


Critically, as in (<ref>), we obtain the
following recursion.


 
For any n≥1, we have that, 

    T_θ,n(𝐱_n)     =∫[T_θ,n-1(𝐱_n-1)+s_θ,n(x_n-1,𝐱_n)] p_θ(d𝐱_n-1| y_0:n-1,𝐱_n)
    
       ≡∫[T_θ,n-1(𝐱_n-1)+s_θ,n(x_n-1,𝐱_n)] p_θ(𝐱_n|x_n-1) p_θ(
    d𝐱_n-1|y_0:n-1)/∫ p_θ(𝐱_n|x_n-1) p_θ(dx_n-1|y_0:n-1).





See Appendix <ref>.


Proposition <ref> gives rise to a Monte-Carlo methodology for a forward-only, online approximation of the smoothing expectation of interest. This is given in Algorithm <ref>.





Algorithm <ref> uses a simple bootstrap filter with multinomial resampling applied at each step. 
The variability of the Monte-Carlo estimates can be further reduced by incorporating: more effective resampling (e.g., systematic resampling <cit.>, stratified resampling <cit.>);
dynamic resampling via use of Effective Sample Size; non-blind proposals in the propagation of the particles.






 §.§ Pathspace versus Finite-Dimensional Construct

One can attempt to define an algorithm without reference to the underlying pathspace. That is, in the case of no jumps (for simplicity)
an alternative approach can involve working with  a regular grid
on the period [0,T], say {s_j=jδ}_j=0^M, with δ=T/M for chosen size M≥ 1. Then, defining 
𝐱'=(x_δ,x_2δ,…, x_Mδ),
and using, e.g., an Euler-Maruyama time-discretisation scheme to obtain the joint density of such an 𝐱' given x=x_0,
a Radon-Nikodym derivative, p_θ^M(𝐱'|x) on ℝ^M× d_x, can be obtained with respect to the Lebesgue reference measure Leb^⊗ (d_x× M), as a product of M conditionally Gaussian densities. 
As shown e.g. in the motivating example in the Introduction, such an approach would lead to estimates with variability that increases rapidly with M, for fixed Monte-Carlo particles N.





A central argument in this work is that one should develop the algorithm in a manner that respect the probabilistic properties of the SDE pathspace, before applying (necessarily) a time-discretisation for implementation on a PC. This procedure is not followed for purposes of mathematical rigour, but it has practical effects on algorithmic performance.




 §.§ Consistency


For completeness, we provide an asymptotic result for Algorithm <ref> following standard results from the literature. We consider the following assumptions.



Let 𝐗 and 𝖷
denote the state spaces of 𝐱 and x respectively.



  (i) For any relevant y', y, ℱ, x', we have that g_θ(y'| y,ℱ)≡ g_θ(y'|x'), 
where the latter 
is a positive function such that, for any y, sup_x∈𝖷 g_θ(y| x)<∞.

  (ii) sup_(x,𝐱')∈ X×𝐗 p_θ(𝐱'| x)<∞.










  (i) Under Assumption <ref>, for any n≥0, there exist constants
b_n,c_n>0, such that for any ϵ>0,


    Prob[ |𝒮_θ,n-𝒮̂_θ,n|>ϵ ] ≤ b_ne^-c_nNϵ^2.



  (ii) For any n≥0, 𝒮̂_θ,n→𝒮_θ,n
w.p.1, as N→∞.




See Appendix <ref>.




§ ONLINE PARAMETER/STATE ESTIMATION FOR SDES



In this section, we  derive an online gradient-ascent 
for partially observed SDEs. 









<cit.> use the score function estimation methodology to
propose an online gradient-ascent algorithm for obtaining an
MLE-type parameter estimate, following ideas in <cit.>.
In more detail, the method is based on the  Robbins-Monro-type of recursion,


    θ_n+1    =θ_n+γ_n+1∇log p_θ_0:n(y_n| y_0:n-1)
    
       =θ_n+γ_n+1{∇log p_θ_0:n(y_0:n)-∇log p_θ_0:n-1(y_0:n-1)}


where {γ_k}_k is a positive decreasing sequence with,

    ∑_k=1^∞γ_k=∞,   ∑_k=1^∞γ_k^2<∞.


The meaning of quantity ∇log p_θ_0:n(y_0:n) is that – given a recursive method (in n) for 
the estimation of θ↦∇log p_θ(y_0:n)
as we describe below and based on the methodology of Algorithm <ref> – one uses θ_n-1 
when incorporating y_n-1, then θ_n for y_n, 
and similarly for k>n. See <cit.> for analytical
studies of the convergence properties of the deduced algorithm, where 
under strong conditions the recursion
is shown to converge to the `true' parameter value, say θ^†, as n→∞. 


Observe that, from Fisher's
 identity (see e.g. <cit.>) we have that,
 

    ∇log p_θ(y_0:n)     =∫∇log p_θ(𝐱_0:n,y_0:n)p_θ(d𝐱_0:n|y_0:n).


Thus, in the context of Algorithm <ref>,  estimation of the score function corresponds to
the choice, 


    S_θ(𝐱_0:n) = ∇log p_θ(𝐱_0:n,y_0:n)≡∑_k=0^n∇log p_θ(𝐱_k,y_k|x_k-1).


and, 


    s_θ,k(x_k-1,𝐱_k) ≡∇log p_θ(𝐱_k,y_k|x_k-1).


Combination of the Robins-Morno recursion (<ref>) with the one in Algorithm <ref>, 
delivers Algorithm <ref>, which we have presented  here 
in some detail for clarity.





When the joint density of (𝐱_0:n,y_0:n) is in the exponential family, an online EM algorithm can also be developed;
see <cit.> for the discrete-time case.





§ NUMERICAL APPLICATIONS



All numerical experiments shown in this section were obtained on a standard PC, with code written in . Unless otherwise stated, executed algorithms discretised the pathspace using the Euler–Maruyama scheme with M=10 imputed points between consecutive observations, the latter assumed to be separated by 1 time unit. Algorithmic costs scale linearly with M.

To specify the schedule for the scaling parameters {γ_k}
in (<ref>), we use the standard adaptive method termed , by <cit.>. Assume that after n steps, we have 


c_n:=-(∇log p_θ̂_0:n(y_0:n)-∇log p_θ̂_0:n-1(y_0:n-1)).
 applies the following iterative procedure,


    m_n=m_n-1β_1+(1-β_1)c_n,     v_n=v_n-1β_2+(1-β_2)c_n^2,
    m̂_n=m_n/(1-β_1^n),  v̂_n=v_n/(1-β_2^n),
    θ̂_n+1=θ̂_n-αm̂_n/(√(v̂_n)+ϵ),


where (β_1,β_2,α,ϵ) are tuning parameters.
Convergence properties of  have been
widely studied <cit.>. Following <cit.>, in all cases  

below we set   (β_1,β_2,α,ϵ)=(0.9,0.999,0.001,10^-8).
 is nowadays a standard and very effective addition to the type of recursive inference algorithms we are considering here, even more so as for increasing dimension of unknown parameters. See the above references for more motivation
and details.

 






 §.§ Ornstein-Uhlenbeck SDE

We consider the model,


    dX_t  =θ_1(θ_2-X_t)+θ_3dW_t+dJ_t,   X_0=0,
    
    y_i=x_i+ϵ_i,    i≥ 1,   ϵ_ii.i.d.∼𝒩(0,0.1^2),   Δ=1,


where J_t=∑_i=1^N_tξ_i, 
N_t is a Poisson process of intensity λ≥ 0,
and ξ_ii.i.d.∼𝒰(-ζ,ζ), 
ζ> 0. 

O–U: Experiment 1. 
Here, we choose λ=0, so the score function is
analytically available; we also fix θ_2=0. 
We simulated datasets of size n= 2,500, 5,000, 7,500, 10,000 under true remaining parameter values (θ_1^†,θ_3^†) =(0.4,0.5).
We executed Algorithm <ref>
to approximate the score function with  N=50, 100, 150 particles,
for all above datasets.

Fig. <ref> summarises results of R=50 replicates 
of estimates of the bivariate score function, with the dashed lines showing the true values.












O–U: Experiment 2. 
We now present two scenarios: 
one still without jumps (i.e., λ=0); and one with jumps, such that we 
fix λ=0.5, ζ = 0.5.
In both scenarios, the parameters to be estimated are θ=(θ_1,θ_2,θ_3).
Beyond the mentioned fixed values for λ, ζ, 
in both scenarios we generated n=20,000 data points using the true parameters values
θ^†=(θ_1^†,θ_2^†,θ_3^†)=(0.2,0.0,0.2).
We applied Algorithm <ref> with N=100 particles, and initial parameter values θ̂_0=(1.0,1.0,1.0) for both scenarios. We used Construct 1 for the scenario with jumps.
Fig. <ref> summarises the results. 


O–U: Experiment 3. 
To compare the efficiency of Constructs 1 and 2 for the case of jump processes, as developed in Section <ref>,
we applied Algorithm <ref> – for each case – 
to approximate the score function (at the true values for (θ_1,θ_2,θ_3)), with N=50,100,150,200 particles,
with n=10 data generated from O–U process (<ref>) 
with fixed λ=0.5, ζ=0.5, and true remaining parameters
(θ_1^†,θ_2^†,θ_3^†)=(0.3,0.0,0.2).
Fig. <ref> shows boxplots summarising the results 
from R=50 replications of the obtained estimates.
We note that performance of Construct 2 is, in this 
case, much better than that of Construct 1. 












 §.§ Periodic Drift SDE






We consider the (highly) nonlinear
model,


    dX_t =sin(X_t-θ_1)dt+θ_2dW_t,   X_0=0,
    
    
    y_i=x_i+ϵ_i,    i≥ 1,   ϵ_ii.i.d.∼𝒩(0,0.1^2),   Δ=1,


for θ_1∈[0,2π), θ_2>0.
We used true values (θ_1^†,θ_2^†)=(π/4,0.9),  
generated n=10^4 data, and applied Algorithm <ref> with N=100 particles and initial
values (θ̂_1,0, θ̂_2,0) =(0.1,2). To illustrate the mesh-free performance of Algorithm <ref>, we choose two values, M=10 and M=100 and compare the obtained numerical results. Figure <ref> shows the data (Left Panel) and the trajectory of the estimated parameter values (Right Panel). Note that the SDE in (<ref>) is making infrequent jumps between a number of modes.  Indeed, in this case, the results are stable when  increasing M from 10 to 100.























 §.§ Heston Model


Consider the Heston model <cit.>,
for price S_t and volatility X_t, modelled as,  


    dS_t   =θ_4S_tdt+S_tX_t^1/2dB_t, 
    
    dX_t   =θ_1(θ_2-X_t)dt+θ_3X_t^1/2dW_t,   
    X_0= θ_2,


where B_t, W_t are independent standard Brownian motions. 
Define the log-price
U_t=log(S_t), so that application of Itô's lemma gives, 


    dU_t   =(θ_4-X_t2)dt+X_t^1/2dB_t,
    
    dX_t   =θ_1(θ_2-X_t)dt+θ_3X_t^1/2dW_t,   X_0=θ_2.

We make the standard assumption 2θ_1θ_2>θ_3^2, 
so that the CIR process X={X_t} will not hit 0; 
also, we have θ_1, θ_2, θ_3, θ_4>0. Process U_t is observed at discrete times,
so that y_i=U_t_i, i≥ 1. Thus, we have, 


    [ y_i | y_i-1, {X_s;s∈[t_i-1,t_i]} ]    ∼𝒩(y_i;μ_i,Σ_i),


where we have set, 


    μ_i=y_i-1+∫_t_i-1^t_i( θ_4-X_s2) ds,  
    Σ_i   =∫_t_i-1^t_iX_sds.


We chose true parameter value θ^†=(0.1,1.0,0.2,0.45),
and generated n=10^4 observations, with Δ=1.
We applied Algorithm <ref> with θ̂_0=(0.005,0.1,0.4,0.3) and N=100 particles.
The trajectories of the estimated parameters are given in Figure <ref>. 






 §.§ Sequential Model Selection – Real Data

We use our methodology to carry out sequential model selection, motivated by applications in <cit.>.
Recall that BIC <cit.> for model ℳ, with parameter vector θ, and data y_1:n
is given by,


    BIC(ℳ)     :=-2ℓ_θ̂_MLE(y_1:n)+
    dim(θ)log n,


where θ̂_MLE denotes the MLE, 
and ℓ(·) the log-likelihood.
BIC and the (closely related) Bayes Factor are known to have good asymptotic properties, e.g. they are consistent Model Selection Criteria, under the assumption of model-correctness, in the context of nested models, for specific classes of models (see e.g. <cit.>, and <cit.> for the case of discrete-time nested 
HMMs). One can plausibly conjecture such criteria will also perform well
for the type of continuous-time HMMs we consider here. See also <cit.> for
a rigorous analysis of BIC for diffusion-type models.
Given models ℳ_a, ℳ_b, with ℳ_a⊂ℳ_b, 
we use the difference, 


    BIC(ℳ_ab):=-2(ℓ_θ̂_MLE,a(y_1:n)+
    ℓ_θ̂_MLE,b(y_1:n))+
    (dim(θ_a)-dim(θ_b))log n,

with involved terms defined as obviously, to choose between ℳ_a and ℳ_b.


























 
We note that 
<cit.> use Bayes Factor for model selection 
in the context of SDE processes; that work 
uses MCMC to estimate parameters, and is not sequential (or online). 
We stress that our methodology allows for carrying out inference for parameters and the signal part of the model, and simultaneously allowing for model selection, all in an online fashion. Also, it is worth noting that <cit.> use
the sequential likelihood ratio for model comparison, but such quantity
can overshoot, i.e. the likelihood ratio will tend to
choose a large model. 
Also, that work uses fixed calibrated
parameters, so it does not relate to online parameter inference. 
















We will use the running estimate of the parameter vector as a proxy for the
MLE given the data already taken under consideration. 
Similarly, we use the weights of the particle filter to obtain a running 
proxy of the log-likelihood evaluated at the MLE, thus overall 
an online approximation of BIC. 
Such an approach, even not fully justified theoretically,  
can provide reasonable practical insights when performing model comparison 
in an online manner – particularly so, given when there is no alternative, to the best of our knowledge, for such an objective.










We consider the following family of nested SDE models,


    dX_t    =b_θ^(i)(X_t)+θ_4√(X_t)dW_t,


where,


    ℳ_1:        b_θ^(1)(X_t)=θ_0+θ_1X_t,
    ℳ_2:        b_θ^(2)(X_t)=θ_0+θ_1X_t+θ_2^2X_t^2,
    ℳ_3:        b_θ^(3)(X_t)=θ_0+θ_1X_t+θ_2^2X_t^2+θ_3/X_t.


Such models have been used for short-term interest rates; see <cit.>
and references therein for more details. Motivated by
<cit.>, we applied
our methodology to daily 3-month Treasury Bill rates,
from 2nd of January 1970, to 29th of December 2000, providing n=7,739 observations.
The data can be 
obtained from the Federal Reserve Bank of St. Louis, at webpage <https://fred.stlouisfed.org/series/TB3MS>. The dataset is shown
at Fig. <ref>.















Fig. <ref> shows the BIC differences, obtained online, for each of the three pairs of models. 


The results suggest that ℳ_2 does not
provide a good fit to the data – relatively to ℳ_1, ℳ_3 – for almost all period under consideration. 
In terms of ℳ_1 against ℳ_3, there is evidence 
that once the data from around 1993 onwards are taken under consideration, ℳ_3 is preferred to ℳ_1.
In general, one can claim that models with non-linear drift should be taken under consideration for fitting daily 3-month Treasury Bill rates,
without strong evidence in favour or against linearity. 
This non-definite conclusion is in some agreement with the empirical
studies in <cit.>.





§ CONCLUSIONS AND FUTURE DIRECTIONS



We have introduced an online particle smoothing methodology for discretely observed (jump) diffusions with intractable transition densities. Our approach overcomes 
such intractability by formulating the problem on
pathspace, thus delivering an algorithm that – besides regulatory conditions – requires only the 
weak invertibility Assumption <ref>. Thus, we have covered 
a rich family of SDE models, when related literature  
imposes quite strong restrictions. Combining our online
smoothing algorithm with a Robbins-Monro-type approach of Recursive Maximum-Likelihood, we set up an online stochastic gradient-ascent for the likelihood function of the SDEs under consideration. 
The algorithm provides a wealth of interesting output, that can provide a lot of useful insights in statistical applications. The numerical examples show a lot of promise for the performance of the methodology.
Our framework opens up a number of routes for insights and future research, 
including the ones described below.



  (i) In the case of SDEs of jumps, we have focused on jump dynamics driven by compound Poisson processes. There is great scope for generalisation here, and one can extend the algorithm to different cases of jump processes, also characterised by more complex dependencies between the jumps and the paths of the solution of the SDE, X={X_t}.  
Extensions to time-inhomogeneous cases with b(v)=b(t,v), σ(v) = σ(t,v),
are immediate; we have chosen the time-homogeneous models only for purposes of  presentation simplicity. 




  (ii) Since the seminal work of <cit.>, more `tuned' auxiliary bridge processes have appeared in the literature, see e.g. the works of <cit.>. Indicatively, the work in <cit.> considers bridges of the form (in one of the many options they consider),


    dX̃_t = {b(X̃_t) + Σ(X̃_t)Σ^-1(x')x'-X̃_t/T-t}dt + 
    σ(X̃_t)dW_t.


Auxiliary bridge processes that are closer in dynamics to the diffusion bridges of the given signal are expected to reduce the variability of Monte-Carlo algorithm, thus progress along the above direction can be immediately incorporated in our methodology and improve its performance. For instance, as noted in <cit.>, use of (<ref>) will give a Radon-Nikodym derivative where stochastic integrals cancel out. Such a setting is known to considerably reduce the variability of Monte-Carlo methods, see e.g. the numerical examples in <cit.> and the discussion in <cit.>.


  (iii) The exact specification of the recursion used for the online estimation of unknown parameters is in itself an problem of intensive research in the field of stochastic optimisation. One would ideally aim for the recursion procedure to provide parameter estimates which are as close to the unknown parameter as the data (considered thus far) permit. In our case, we have used a fairly `vanilla' recursion, maybe with the exception of the  variation. E.g., recent works in the Machine Learning community have pointed at the use of `velocity' components in the recursion to speed up convergence,
see, e.g., <cit.>.

  (iv) 
We mentioned in the main text several modifications that can
improve algorithmic performance: dynamic resampling, stratified resampling, non-blind proposals in the filtering steps, choice of auxiliary processes. Parallelisation and use of  
HPC are obvious additions in this list.


  (v) 
We stress that the algorithm involves a filtering step, and a step 
that approximates the values of the instrumental function. 
These two procedures should be thought of separately.  
One reason for including two approaches in the case of jump diffusions (Constructs 1 and 2) is indeed to highlight this point. The two Constructs are identical in terms of the filtering part.
Construct 1 incorporates in 𝐱' the location of the path at all times of jumps; thus, when the algorithm `mixes' all pairs of {x_k-1^(j)},  {𝐱_k^(i)},  at the update of the instrumental function (see Step (iv) of Algorithm <ref>), many of such pairs can be incompatible (simply consider any such pair with i≠ j, when {𝐱_k^(i)} has in fact been generated conditionally on {x_k-1^(i)}). Such an effect is even stronger in the case of the standard algorithm applied in the motivating example in the Introduction, and partially explains the inefficiency of that algorithm. In contrast, in Construct 2, 𝐱' contains less information about the underlying paths, 
thereby improving the compatibility of pairs selected from  particle populations {x_k-1^(j)},  {𝐱_k^(i)},
thus – not surprisingly – Construct 2 seems more effective 
than Construct 1. 

  (vi)
The experiments included in Section 6 of the paper aim to highlight the qualitative aspects of the introduced approach. 
The key property of robustness to mesh-refinement underlies all numerical applications, and is explicitly illustrated  in the experiment in Section 6.2, and within Figure 6.4.
In all cases, the code was written in , was not parallelised or optimised and was run on a PC, with executions times that were quite small due to the online nature of the algorithm (the most expensive experiment took about 30 minutes).
We have avoided a comparison with potential alternative approaches as we believe this would obscure the main messages we have tried to convey within the numerical results section. Such a detailed numerical study can be the subject of future work.




§ ACKNOWLEDGMENTS

We are grateful to the AE and the referees for several useful comments that helped us to improve the content of the paper.



















apalike







§ PROOF OF PROPOSITION 1
 

We have the integral, 


    ∫ T_θ,n(𝐱_n) p_θ(d𝐱_n| y_0:n) =  
    ∫ S_θ(𝐱_0:n) p_θ(d𝐱_0:n-1| y_0:n-1,𝐱_n)⊗ 
    p_θ(d𝐱_n| y_0:n).


Also, simple calculations give, 


    p_θ(d𝐱_0:n-1| y_0:n-1,   𝐱_n)
    ⊗ p_θ(𝐱_n| y_0:n) 
       ≡  p_θ(d𝐱_0:n-1| y_0:n,𝐱_n)⊗ p_θ(𝐱_n| y_0:n)
    = p_θ(d𝐱_0:n| y_0:n).


Using this expression in the integral on the right side of (<ref>) completes the proof.




§ PROOF OF PROPOSITION 2
 

Simply note that,


    ∫ T_θ,n-1(𝐱_n-1)
        p_θ(d𝐱_n-1| y_0:n-1,𝐱_n) =
       =
    ∫ S_θ,n-1(𝐱_0:n-1)  p_θ(d𝐱_0:n-2| y_0:n-2,𝐱_n-1) p_θ(d𝐱_n-1| y_0:n-1,𝐱_n).

Then, we have that,


    p_θ   (d𝐱_0:n-2| y_0:n-2,𝐱_n-1) p_θ(d𝐱_n-1| y_0:n-1,𝐱_n) 
           ≡ 
    p_θ(d𝐱_0:n-2| 𝐱_n-1, y_0:n-1,𝐱_n) p_θ(d𝐱_n-1| y_0:n-1,𝐱_n)
           ≡   p_θ(d𝐱_0:n-1| y_0:n-1,𝐱_n).


Replacing the probability measure on the left side of the above equality with its equal on the right side, and using the latter in the integral above completes the proof for the first equation in the statement of the proposition.
The second equation follows from trivial use of Bayes rule.




§ PROOF OF PROPOSITION 3
 

Part (i) follows via the same arguments as in <cit.>, based on <cit.> and <cit.>.
For part (ii), one can proceed as follows.
Given n≥0, we define the event  A_N(1/j):={|𝒮_θ,n-𝒮̂_θ,n^N|>1/j},
where we added superscript N to stress the dependency of the estimate
on the number of particles. Then,


    Prob [ lim_N→∞𝒮̂_θ,n^N=𝒮_θ,n ]=
    1-Prob    [ ∪_j=1^∞lim sup_N→∞A_N(1/j) ]
       ≥ 1-∑_j=1^∞Prob [ lim sup_N→∞A_N(1/j) ].


From (i), we get   Prob [ A_N(1/j) ]  ≤ b_ne^-c_nN(1/j)^2, so
∑_N=1^∞Prob [ A_N(1/j) ]<∞.
The Borel–Cantelli lemma gives Prob [ lim sup_N→∞A_N(1/j) ]=0,
therefore the result follows from (<ref>).


