Bayesian inference for Stable Levy driven Stochastic Differential Equations with high-frequency data
Ajay Jasra Kengo Kamatani Hiroki Masuda
math.ST stat.TH




1]Ajay Jasrastaja@nus.edu.sg
2]Kengo Kamatanikamatani@sigmath.es.osaka-u.ac.jp
3]Hiroki Masudahiroki@math.kyushu-u.ac.jp
[1]Department of Statistics  Applied Probability, National University of Singapore, SG
[2]Department of Engineering Science, Osaka University, JP
[3]Faculty of Mathematics, Kyushu University, JP




Bayesian inference for Stable Lévy driven Stochastic Differential Equations with high-frequency data
    [
    This version: December 30, 2023
====================================================================================================

 
In this article we consider parametric Bayesian inference for stochastic differential equations (SDE) driven by a pure-jump stable Lévy process, which is observed at high frequency. In most cases of practical interest, the likelihood function is not available, so we use a quasi-likelihood and place an associated prior on the unknown parameters. It is shown under regularity conditions that  there is a Bernstein-von Mises theorem associated to the posterior. We then develop a Markov chain Monte Carlo (MCMC) algorithm for Bayesian inference and assisted by our theoretical results, we show how to scale Metropolis-Hastings proposals when the frequency of the data grows, in order to prevent the acceptance ratio going to zero in the large data limit. Our algorithm is presented on numerical examples that help to verify our theoretical findings.

Keywords: Markov Chain; Monte Carlo; Lévy process; Bayesian inference; high-frequency data





§ INTRODUCTION


Stochastic differential equations (SDE) are found in a wide variety of real applications, such as finance and econometrics (see for instance <cit.> and the references therein), mathematical biology, movement ecology, turbulence, signal processing, to mention just a few. In this article we are concerned with Bayesian parameter estimation of discretely observed SDE driven by a pure-jump stable Lévy process with high-frequency data, which is often found in the afore-mentioned applications.
Among others, we refer to <cit.>, <cit.>, and <cit.> for a comprehensive account of stable distribution and stochastic processes driven by a stable .

Often, the main challenge with Bayesian or classical inference with Lévy driven SDE is the lack of tractability of the transition density of the process, hence one cannot evaluate a non-negative and unbiased estimate of the transition density. The latter is often required for inference methods.
In such scenarios, one often has to resort to time-discretization. We focus on the time discretization induced by the quasi-likelihood
approach, which, under suitable regularity conditions, enables us to deduce several desirable properties, such as consistency and asymptotic (mixed) normality; see <cit.> and the references therein.
As a result, one expects rather favorable properties of the associated posterior distribution in similar settings. In addition,
one must construct inference techniques, in this article MCMC, which will turn out to be robust in the high-frequency limit. That is, algorithms that will not collapse in some sense.

Bayesian inference for diffusions and jump diffusions have been considered in many articles.
This includes the fully observed (jump) diffusion case <cit.> and the partially observed jump diffusion 
case <cit.>. In particular the work <cit.> considered some related, but different classes of models to this article,
except using approximate Bayesian computation <cit.> methods with MCMC. 
We note that most of the previous works do not consider quasi-likelihood and/or the large data limit of the performance of MCMC.
Our contributions of this article are roughly summarized as follows:


  * To construct an approximate posterior distribution with favorable theoretical properties.
That is, under assumptions, that there is a Bernstein-von Mises theorem associated to the posterior.

  * To develop an MCMC algorithm and assisted by our theoretical results, show how to scale
Metropolis-Hastings proposals when the frequency of the data grows, in order to prevent the acceptance ratio going to zero in the large data limit.












This article is structured as follows. In Section <ref>, the model setup and an associated Bernstein-von Mises theorem are given.
Our MCMC algorithm is also described. In Section <ref>,
our theoretical results for MCMC are stated. Section <ref> provides some numerical simulations which confirm our theoretical findings.
We also analyze a real data set from the NYSE 
which features properties often captured well by the processes we study in this article. 
The appendix features a variety of proofs for our technical results.






§ STABLE-LÉVY SDES, MODEL AND ALGORITHM




 §.§ Setup


Let (Ω,ℱ,(_t)_t∈[0,T],ℙ) be a complete stochastic basis. 
Let {X_t;t∈ [0,T]} be a solution to the one-dimensional stochastic differential equation

    X_t=a(X_t,α) t+c(X_t-,γ) J_t,

where {J_t}_t is the stable Lévy-process independent of the initial variable X_0 and such that

    𝔼(e^iuJ_1)=e^-|u|^β,

where it is assumed throughout that _t=(J_s: s≤ t)∨(X_0) and that 1≤β <2.
Denote by ϕ_β the β-stable density of (J_1).
We will write ℙ_θ for the image measure of (X_t) associated with the parameter value

    θ:=(,)∈Θ_×Θ_=Θ⊂^p,

with Θ_⊂^p_ and Θ_⊂^p_ being bounded convex domains. The true value is denoted by =(_0,_0)∈Θ. For brevity, we will also write  for _θ_0 as well.
Let Θ denote the closure of Θ, and write a≲ b if a≤ Cb for some universal constant C>0. 

[Regularity of the coefficients]


  * The functions a(·,_0) and c(·,γ_0) are globally Lipschitz and of class ^2(), and c(x,γ) is positive for every (x,γ);

  * a(x,·)∈^3(Θ_) and c(x,·)∈^3(Θ_) for each x∈;

  * sup_θ∈Θ{max_0≤ k≤ 3max_0≤ l≤ 2(
|_^k_x^la(x,)| 
+ |_^k_x^lc(x,)|) + c^-1(x,)}≲ 1+ |x|^C.





[Identifiability conditions]
The random functions t↦(a(X_t,), c(X_t,)) and t↦(a(X_t,_0), c(X_t,_0)) 
on [0,T] a.s. coincide if and only if θ=.









The process X is observed only at discrete time t=nh (j=0,…, N) where h=T/N with the terminal sampling time T>0 being fixed.
Let X^N:={X_nh;n=0,…, N}, the available data set.
The Euler-Maruyama discretization under ℙ_θ is then given by

    Δ^N_nX ≈ a(X_(n-1)h,α)h+c(X_(n-1)h,γ)Δ^N_n J,

where Δ^N_n X:=X_nh-X_(n-1)h and Δ^N_nJ:=J_nh-J_(n-1)h.
This suggests to consider the following stable quasi-likelihood <cit.> (the multiplicative constant “h^-N/β” omitted):

    𝕃_N(θ|X^N)=∏_n=1^N1/c(X_(n-1)h,γ)ϕ_β(Δ^N_nX-a(X_(n-1)h,α)h/c(X_(n-1)h,γ)h^1/β).

We define the stable quasi-maximum likelihood estimator by any element

    =(,)∈ argmax_θ∈Θlog_N(θ|X^N).

Note that 𝕃_N(θ|X^N) is the true likelihood only when both a(x,) and c(x,) are constants.

It follows from <cit.> that under Assumptions <ref> and <ref> the estimator  is asymptotically mixed-normally distributed: on a suitably extended probability space (Ω̃,ℱ̃,ℙ̃) carrying a p-dimensional standard Gaussian random variable  η independent of , we have

    (√(n)h^1-1/β(-_0), √(n)(-_0)) _st→ I()^-1/2η,

where _st→ stands for the -stable convergence in law (see <cit.> for details)
and where I()=( I_(),  I_(_0)) denotes the a.s. positive definite quasi Fisher-information matrix specified by

    I_()=C_(β)Σ_T,(θ_0),      I_(_0)= C_(β)Σ_T,(_0)

with

    C_(β) := ∫(∂ϕ_β/ϕ_β(y))^2ϕ_β(y) y,      
    C_(β):=∫(1+ y∂ϕ_β/ϕ_β(y))^2ϕ_β(y) y, 
    
       Σ_T,(θ_0) :=1/T∫_0^T{∂_a(X_t,_0)}^⊗ 2/c^2(X_t,_0) t,     Σ_T,(_0) :=1/T∫_0^T{∂_c(X_t,_0)}^⊗ 2/c^2(X_t,_0) t.

It is known that the stable quasi-maximum likelihood estimator θ̂_N attains the Hajék-Le Cam minimax lower bound in some special cases, see <cit.> and <cit.>.
The convergence (<ref>) shows that we do have the conventional Studentization (asymptotic standard normality) result without a finite-variance property as well as any stability condition such as ergodicity.





 §.§ Posterior and Markov chain Monte Carlo methods




  §.§.§ Quasi-posterior distribution





















The quasi-posterior distribution is given by

    Π_N(θ|X^N)∝𝕃_N(θ|X^N)Π(θ),

where Π(θ)=π(θ)θ denotes the prior distribution of θ. The β-stable density ϕ_β does not have a closed form for β∈(1,2), hence cannot be computed pointwise and neither then can the posterior up-to a constant; this is required for stochastic simulation algorithms such as MCMC. In general, ϕ_β can be only computed via certain numerical integration, the iteration of which may be rather time-consuming <cit.>. However, we have access to a non-negative unbiased estimate of the un-normalized posterior, which does suffice to apply (e.g.) MCMC and it may be constructed as follows.
Let F_β( v)=f_β(v) v denote the positive β/2-stable distribution with Laplace transform 

    ∫_0^∞ e^-vtF_β( v)=e^-|2t|^β/2   (t≥ 0).

Then, the stable distribution (J_1) is the law of √(V)W where W∼ N(0,1) and V∼ F_β, 
and we have the normal variance-mixture representation:

    ϕ_β(x)=∫_0^∞1/√(2π v)exp(-x^2/2v)F_β( v).

If V_1,…, V_N ∼ F_β, 
and ξ_1,…, ξ_N∼ N(0,1), then ℒ(h^1/β√(V_n)ξ_n)=ℒ(Δ^N_nJ).
Invoking (<ref>) we may expect that the formal distributional approximation

    Δ^N_nX    ≈ a(X_(n-1)h,α)h+c(X_(n-1)h,γ)h^1/β√(V_n)ξ_n

under ℙ_θ would be meaningful.

Building on the above observation, we will make use of the “complete” quasi-likelihood 

    𝕃_N(θ|X^N,V^N)=∏_n=1^N1/c(X_(n-1)h,γ)√(V_n)ϕ(Δ^N_nX-a(X_(n-1)h,α)h/c(X_(n-1)h,γ)√(V_n)h^1/β)f_β(V_n),

where ϕ denotes the standard normal density. The corresponding posterior distribution is

    Π_N(θ,V^N|X^N)∝𝕃_N(θ|X^N,V^N)Π(θ).

Note that we can still not compute Π_N(θ,V^N|X^N) up to a constant. Nevertheless, we will still be able to devise a Metropolis-within-Gibbs MCMC algorithm which mitigates this issue.



  §.§.§ MCMC Algorithm


Let θ∈Θ. For n=1,…, N,  generate V_n from F_β( v|ϵ_n(θ))
where

    F_β( v|x)=
    	.1/√(2π v)exp(-x^2/2v)F_β( v)/ϕ_β(x)

and 

    ϵ_n(θ)= Δ^N_nX-a(X_(n-1)h,α)h/c(X_(n-1)h,γ)h^1/β.

Then update θ^*←θ+(noise)
and accept θ^* with probability 

    α(θ,θ^*)=min{1,𝕃_N(θ^*|X^N,V^N)π(θ^*)/𝕃_N(θ|X^N,V^N)π(θ)}.

This noise should be properly scaled by the rate matrix

    D_N := diag( √(N)h^1-1/βI_p_, √(N)I_p_),

 as specified in (<ref>). 
 The random variable V_n∼ F_β( v|ϵ_n(θ)) can be sampled via rejection sampling for each n∈{1,…,N}:
Generate v∼ F_β( v) and accept it with probability

    .{v^-1/2exp(-x^2/2v)}/ {|x|^-1/2exp(-|x|/2)}.

The algorithm is presented in Algorithm <ref>
(We write N_p(μ,Σ) for the p-dimensional Gaussian distribution with mean μ and covariance matrix Σ).





A minor modification is to consider a (correlated) pseudo-marginal algorithm <cit.> to avoid rejection sampling. 
In this case, update θ as in Algorithm in <ref>
together with

    V_n^*←ρ^2/β V_n + (1-ρ)^2/βξ_n,  ξ_n∼ F_β,

as in <cit.> and <cit.>. 
Then accept (θ^*, V^N*)  where V^N*={V_n^*}_n=1,… N with probability 

    α((θ,V^N),(θ^*,V^N*))=min{1,𝕃_N(θ^*|X^N,V^N*)π(θ^*)/𝕃_N(θ|X^N,V^N)π(θ)}.

Here, ρ is a tuning parameter and we recommend to take ρ close to 1. 


























 §.§ Bernstein-von Mises theorem



In this section we give a Bernstein-von Mises theorem associated with the stable quasi-likelihood (<ref>). Different from the classical version <cit.>, we will look at the convergence of the posterior distribution of the rescaled parameter centered not at , but at .


The prior distribution Π admits a bounded Lebesgue density π which is continuous and positive at .



Let 

    _N(θ):=log_N(θ|X^N)

be the logarithmic quasi-likelihood function. Following <cit.>, we introduce the quasi-likelihood ratio random field

    _N(u):=exp{_N(θ_0+D_N^-1u) - _N(θ_0) }

defined on the set _N=_N():=D_N(Θ-)⊂^p. For convenience, we extend the domain of _N into the whole ^p with _N≡ 0 on ^p∖_N. 
Then, the quasi-posterior distribution Π̅_N(du|X^N) of the rescaled parameter u=D_N(θ-) admits the density

    u ↦_N(u)π(+D_N^-1u)/∫_N(v)π(+D_N^-1v) v.

Let Δ_N():=D_N^-1_θ_N() and denote by  the convergence in probability.
For an -measurable random variable (μ,Σ)∈^p×^p× p, we keep denoting by N_p(μ,Σ) the -conditional Gaussian distribution associated with the characteristic function ξ↦[exp{iμ[ξ]-(1/2)Σ [ξ,ξ] }].
Let ν_TV=sup_|f|≤ 1|∫ f(x)ν( x)| denote the total variation norm of a signed measure ν. 


Under Assumptions <ref> to <ref>, we have

    Π̅_N(·|X^N)-N_p(I()^-1Δ_N(),  I()^-1)_TV 0.





The proof of Theorem <ref> is given in Appendix <ref>.






















§ LARGE SAMPLE PROPERTIES OF THE MARKOV CHAIN MONTE CARLO METHOD


In this section, we will first introduce a general framework for developing some stability properties of MCMC algorithms, and then apply it to the  MCMC algorithms proposed in Section <ref>.
In particular, we show how to scale the proposals when the frequency of the data grows.



 §.§ Local consistency property


For a moment we step away from the main context. Notations and terminologies in this section generally follow those of <cit.>.



  §.§.§ Some general notions


Let E be a normed space equipped with the Borel σ-algebra ℰ. 
For a probability measures μ and ν on (E,ℰ), let 

    μ-ν_BL=sup_f∈BL_1|∫ f(x)μ( x)-∫ f(x)ν( x)|

 be the bounded Lipschitz distance between μ and ν where 
BL_1 denotes the set of any ℰ-measurable real-valued function f such that 
sup_x∈ E|f(x)|≤ 1 and |f(x)-f(y)|≤x-y (x, y∈ E). 

Let 
ℰ_N={𝒳^(N),𝔄^(N),𝐏_θ^(N),θ∈Θ}  be a family of statistical experiments and x^N be the corresponding observation. Here, the parameter space 
Θ⊂ℝ^p be an open set equipped with the Borel σ-algebra Ξ. 
Let 𝕃_N(θ|x^N) be the 𝒳^(N)⊗Ξ-measurable quasi-likelihood function which approximates the true likelihood 𝐏_θ^(N)(x^N). 
Let Π(θ) be the prior distribution for the experiments, and Π_N(θ|x^N) be the quasi-posterior distribution defined by 

    Π_N(θ|x^N)=𝕃_N(θ|x^N)Π(θ)/∫_Θ𝕃_N(θ|x^N)Π(θ).



We shall refer to ℳ_N the Markov chain Monte Carlo method for the experiment ℰ_N
if the algorithm generates a Markov chain {θ^N_m;m=0,1,…} for each x^N.
We shall refer to F^N_M the empirical distribution of  {θ^N_m;0≤ m≤ M-1}, that is, 

    F^N_M(A)=1/M∑_m=0^M-11_A(θ_m^N)

for each Borel set A of ℝ^p. Let D_N be a  non-degenerate p× p-matrix for each N∈ℕ. 
In practice, we expect that Markov chain Monte Carlo methods are robust and does not collapse as N→∞. 
We formalize this favourable property as follows. 


	A family ℳ_N is called locally consistent if  as N, M→∞ 
	
    Π_N(θ_0 + D_N^-1 u|x^N)-F^N_M(θ_0 + D_N^-1 u)_BL→ 0

	in 𝐏_θ_0^(N)-probability. 
 

In the next subsection, we introduce a sufficient condition for this property under the general framework. 




  §.§.§ Sufficient conditions


Let Q_N(θ,ϑ|x^N) (n∈ℕ) be a probability transition kernel 
from 𝒳^(N)×Θ to Θ, 
and let Π_N(θ|x^N) be a probability transition kernel from 𝒳^(N) to Θ. 
Let A_N:𝒳^(N)×Θ^2→[0,1] be a measurable map. Let 

    P_N(θ,ϑ|x^N)   =Q_N(θ,ϑ|x^N)A_N(θ,ϑ|x^N)+R_N(θ|x^N)δ_θ(ϑ),

where 

    R_N(θ|x^N)=1-∫_Θ Q_N(θ,ϑ|x^N)A_N(θ,ϑ|x^N).

As probability measures on Θ^2, we assume that for each x^N∈𝒳^(N), 

    Π_N(θ|x^N)Q_N(θ,ϑ|x^N)A_N(θ,ϑ|x^N)
    	=
    	Π_N(ϑ|x^N)Q_N(ϑ,θ|x^N)A_N(ϑ,θ|x^N).

Then for each x^N, the probability transition kernel P_N(θ,ϑ|x^N) is Π_N(θ|x^N)-reversible, that is, 

    Π_N(θ|x^N)P_N(θ,ϑ|x^N)=Π_N(ϑ|x^N)P_N(ϑ,θ|x^N)

as probability measures on Θ^2. 
Let ℳ_N be the Markov chain Monte Carlo method associated with the transition kernel P_N, 
that is, for each x^N, the algorithm generates a Markov chain  {θ_m^N}_m  with respect to the probability transition kernel P_N(θ,ϑ|x^N). 
For simplicity, the initial point is generated from the quasi-posterior distribution, that is, 
θ_0^N∼Π_N(θ|x^N). 



Let q∈ℕ. Let Π( u|s)=π(u|s) u be a probability transition kernel from an open set 𝒮⊂ℝ^q to ℝ^p
and let Q(u,  v|s)=q(u,v|s) v be a probability transition kernel from 𝒮×ℝ^p to ℝ^p. 
Let A(u,v|s) be a [0,1]-valued function.
For each s∈𝒮 we define a probability transition kernel

    P(u, v|s)   =Q(u, v|s)A(u,v|s)+R(u|s)δ_u( v), 
    
    R(u|s)   =1-∫ Q(u, v|s)A(u,v|s).


[Regularity of the limit kernel]
For each u,v∈ℝ^p, π(u|s), q(u,v|s) and A(u,v|s) are continuous in s, and 

    sup_s∈ Ksup_u,v∈ℝ^pπ(u|s)+q(u,v|s)<∞

for any compact set K of 𝒮. 




Let 
s_N:𝒳^(N)→𝒮 be a sequence of random variables. The following proposition illustrates a sufficient condition for 
local consistency. See Appendix <ref> for the proof. 
In the following, we call that a transition probability kernel P is ergodic if it is irreducible and positive recurrent. 


Suppose that P(u, v|s) is ergodic with invariant probability measure Π( u|s) for each s∈𝒮. 
Suppose that s_N is 𝐏_θ_0^(N)-tight for each θ_0∈Θ. 
	Then, under Assumption <ref>, the family ℳ_N is locally consistent if A_N(θ_0 + D_N^-1 u,θ_0 + D_N^-1 v|x^N)- A(u,v|s_N)→ 0 in probability for each u, v∈ℝ^p and 

    Π_N(θ_0 + D_N^-1 u|x^N)Q_N(θ_0 + D_N^-1 u,θ_0 + D_N^-1 v|x^N)-Π( u|s_N)Q(u, v|s_N)_TV

converges in 𝐏_θ_0^(N)-probability to 0. 






 §.§ Main result


We now go back to the main context given in Section <ref>.
In Algorithm <ref>, the proposal distribution for observation X^N is 

    Q_N(θ_0 + D_N^-1 u,θ_0 + D_N^-1 v|X^N)=Q(u, v|s_N)=N_p(0,Σ)

for some non-degenerate matrix Σ∈ℝ^p× p. 
We will set s_N=(Δ_N(θ_0), I(θ_0), I^*(θ_0)) where I^*() is defined in Section <ref>.
By this setting of the proposal distribution, the following convergence is equivalent to  the condition (<ref>):  

    Π_N(θ_0 + D_N^-1 u|X^N)-Π( u|s_N)_TV=o_p(1).

The acceptance probability of  Algorithm <ref> is 

    A_N(θ,ϑ|X^N)=𝔼[.min{1,𝕃(ϑ|X^N,V^N)π(ϑ)/𝕃_N(θ|X^N,V^N)π(θ)}|X^N,θ].

Let s=(Δ, I, I^*)∈ℝ^p×ℝ^p× p×ℝ^p× p. 
We assume that I and I^* are p× p-symmetric matrices. 
In the following theorem, we will show that the scaled version of the acceptance ratio converges to

    A(u,u^*|Δ,I,I^*)=𝔼[min{1,η}]

where 

    η=Δ[u^*-u]+W[u^*-u]-I[(u^*)^⊗ 2-u^⊗ 2]-1/2I^*[(u^*-u)^⊗ 2]

with W∼ N(0, I^*). 
By showing this, we prove that Algorithm <ref> generates a locally consistent family of Markov chain Monte Carlo methods. 
In particular, the algorithm will not collapse in the large data limit. 
	



Under Assumptions <ref> to <ref>, 
	
    A_N(θ_0+D_N^-1u,θ_0+D_N^-1 v|X^N)-A(u,v|Δ_N(θ_0), I(θ_0), I^*(θ_0))→ 0

	in probability for each u, v∈ℝ^p. 
	In particular, the family of Markov chain Monte Carlo methods ℳ_N defined by Algorithm <ref> is locally consistent. 


The proof of the theorem is given in Section <ref>. 














































































































§ SIMULATIONS




 §.§ Simulated Example



We generated seven data sets from the time discretized models with
a(x,α)=α_1(x-α_2), c(x,γ)=exp(γcos(x)),
β=1.5, h_N=1/N and N∈{10,50,100,250,500,1000,2000}.
α_1, α_2, γ are independently standard normal variables in the prior.
Algorithm <ref> was run
for 10,000 iterations. All simulations are repeated 100 times.
The results are given in Figure <ref>.

Figure <ref> illustrates the average acceptance rate against the number of data.
As can be seen the average acceptance rate is stable from small data to large data.
In particular, the data increases the algorithm does not collapse and the acceptance rate is very reasonable.
The run time of the algorithm for N=2000 is only about five minutes and was coded in R (code is available upon request); the
code could be substantially improved to further reduce the computation time. Note that the update on V^N|⋯
can be parallelized to improve the running speed.





 §.§ Application to IBM Stock Data


We return to the data of NYSE. 
In the model, we set the same model with above, and 
β=1.411, h_N=T/N, where T is measured in minutes (T=390× 3) and
N=1156. The stable index β was estimated by applying `stableFit' function in R-package `fBasics' <cit.> to {Δ_n^N X}_n=1,…, N. The data contains some missing values which are ignored for the purposes of the analysis.
α_1, α_2,γ are independently normal variables in the prior with mean 0 and standard deviation 1.
Algorithm <ref> was run
for 100,000 iterations. 
The algorithm ran in approximately 1 hour
and the acceptance rate for the move on the parameters was 0.34.

In Figures <ref> and <ref> we can observe our results. 
Figure <ref> is a log-likelihood and 
time average drift and  jump coefficients: 

    1/N∑_n=1^Na(X_(n-1),h,α), 1/N∑_n=1^Nc(X_(n-1),h,γ).


For
this reasonable size data set, the algorithm performs well, with good mixing over the parameters.
The acceptance rate is very reasonable as is the run-time - recall one can improve the code or coding language.

Figure <ref> is a p-p plot of the posterior expected value of the standardized residual:

    Δ_n^N X-ha(X_(n-1)h,α)h/c(X_(n-1)h,γ)h^1/β

against β-stable distribution. This provides an idea of the ability of this model to fit the real NYSE data. We can see that the model, to some extent, can exhibit the behaviour
in the real data.









  §.§.§ Acknowledgements

We acknowledge JST CREST Grant Number JPMJCR14D7, Japan, for supporting the research.
AJ was additionally supported by Ministry of Education AcRF tier 2 grant,
R-155-000-161-112, and KK was additionally supported by JSPS KAKENHI Grant Number JP16K00046.







§ PROOFS FOR THE BERNSTEIN-VON MISES THEOREM



The purpose of this section is to prove Theorem <ref>. 
First, we introduce the quadratic random field

    ^0_N(u) =exp( Δ_N()[u] - 1/2I()[u,u] ).

By integrating the Gaussian density,

    ( ∫_N^0(u)π() u )^-1 = 1/π()exp( 1/2I()^-1[Δ_N()^⊗ 2] ) = O_p(1).

To complete the proof it suffices to show

    _N := ∫|_N(u)π(+D_N^-1u) - _N^0(u)π() | u  0,

since the left-hand side of (<ref>) can be bounded from above by the quantity

    ∫|_N(u)π(+D_N^-1u)/∫_N(u)π(+D_N^-1u) u - _N^0(u)π()/∫^0_N(u)π() u|  u
    ≤ 2 ( ∫_N^0(u)π()  u )^-1_N.

We only prove (<ref>) for β∈(1,2), since the remaining case of β=1, where we have the single localization rate √(N)(=√(N)h^1-1/β), is completely analogous and simpler.




First we will introduce a good event G_N∈ whose probability can get arbitrarily close to 1 for N→∞ (see (<ref>) below).
Let

    _1,N(u_1;)    := exp{_N(, _0+u_1/√(N)) - _N(, _0) }, 
    _2,N(u_2)    := exp{_N(_0 + u_2/√(N)h^1-1/β, _0) - _N(_0, _0) },

both being defined to be zero for u=(u_1,u_2)∈^p∖_N.
Also define

    _1,N(θ)    :=1/N{_N(,) - _N(,_0)}, 
    _1()    := 1/T∫_0^T∫[ 
    log{c(X_t,_0)/c(X_t,)ϕ_β(c(X_t,_0)/c(X_t,)z)}
    -logϕ_β(z)]ϕ_β(z) z t,
    
    _2,N()    := 1/Nh^2(1-1/β){_N(,_0) - _N(_0,_0)}, 
    _2()    := -1/2T∫_0^T( a(X_t,_0)-a(X_t,)/c(X_t,_0))^2 t ·∫(ϕ_β/ϕ_β(z))^2ϕ_β(z) z.

The random functions _1 and _2 represent quasi-Kullback-Leibler divergences for estimating  and , respectively;
we can estimate _0 more quickly than _0 in case of β∈(1,2).
For any matrix A we will write |A| for the Frobenius norm of A.
For later use we mention the following statements, which are given in <cit.> or can be directly deduced from the arguments therein.




  * We have (recall (<ref>)) under 

    (Δ_N(),  -D_N^-1_θ^2_N()D_N^-1) →( I()^1/2η,  I()),
    
    
       r̅_N:=sup_θ| D_N^-1_θ^3_N(θ)D_N^-1| = O_p(1).



  * There exists an a.s. positive random variable χ_0 such that for each κ>0,

    sup_;  |-_0|≥κ_1() ∨sup_;  |-_0|≥κ_2() ≤ -χ_0κ^2    a.s.



  * In case of β>1, we have

    sup_| 1/√(N)__N(,_0) | = O_p(1), 
       sup_| -1/N_^2_N(,_0) - C_(β)Σ_T,(_0) |=o_p(1).



In addition, we will also need the following uniform laws of large numbers with convergence rates:




  * There exists a constant q∈(0,1) for which

    (√(N))^qsup_θ|_1,N(θ)-_1()| ∨ (√(N)h^1-1/β)^qsup_|_2,N()-_2()|  0.


A sketch of derivation of (<ref>) will be given at the end of this section.



Set

    _N=N^-c    for     0<c<q/2(1/β-1/2)

with q∈(0,1) given in (<ref>). Then, since

    √(N)h^1-1/β∼ N^1/β-1/2

up to a multiplicative constant, we see that _N√(N)h^1-1/β↑∞ and _N^-2(√(N)h^1-1/β)^-q=O(1);
since β≥ 1, the latter implies that _N^-2N^-q/2=O(1). It is straightforward to deduce from (<ref>) that

    _N^-2sup_θ|_1,N(θ)-_1()| ∨_N^-2sup_|_2,N()-_2()|  0.

Let _min(A) denote the minimum eigenvalue of a square matrix A.
We now introduce the event G_N=G_N(M,) for positive constants M and :

    G_N   := {
    |Δ_N()| ∨sup_| 1/√(N)__N(,_0) | ≤ M,
    
           | -D_N^-1_θ^2_N()D_N^-1 - I() | ∨sup_| -1/N_^2_N(,_0) - C_(β)Σ_T,(_0) | <,
           _N^-2sup_θ|_1,N(θ)-_1()| ∨_N^-2sup_|_2,N()-_2()| < /2,
    
           _min(I_()) ∧_min(I_()) ≥ 4,   χ_0≥ 4,   r_N<3/_N}.

Given any '>0, we can find a triple (M_1,_1,_1) and an N_1∈ such that

    sup_N≥ N_1{G_N(M,)^c}<'

holds for every M≥ M_1 and ∈(0,_1].
Since the objective here is the convergence in probability, we may and do focus on the event G_N with M and  being sufficiently large and small, respectively.



We divide the domain of the integration in the definition of _N: ^p=A_N⊔ A_N^c where

    A_N:={u;  |D_N^-1u|≤_N},

and then denote the associated integrals by '_N and ”_N, respectively:
_N = '_N + ”_N with

    '_N   := ∫_A_N|_N(u)π(+D_N^-1u) - _N^0(u)π() | u, 
    ”_N   := ∫_A_N^c|_N(u)π(+D_N^-1u) - _N^0(u)π() | u.

We will deal with these terms separately.



First we show that '_N 0. We have '_N≤'_1,N+'_2,N, where

    '_1,N   := ∫_A_N_N(u)| π(+D_N^-1u)-π()| u,
    
    '_2,N   := π()∫_A_N| _N(u)-_N^0(u) | u.

By the third order Taylor expansion

    _N(u)    = exp{Δ_N()[u] - 1/2( I() + ( -D_N^-1_θ^2_N()D_N^-1 - I() )
    
            -1/3(D_N^-1_θ^3_N(θ̌_n(u))D_N^-1)[D_N^-1u]
    )[u,u]}

for some random point θ̌_n(u) on the segment joining  and ,
we see that _N(u) ≤exp(M|u| - |u|^2) for u∈ A_N on G_N.
Hence, under Assumption <ref> we obtain

    '_1,N≤sup_u ∈ A_N| π(+D_N^-1u)-π()|
    ∫ e^M|u| - |u|^2 u  0.

To handle '_2,N, we introduce the random function

    M_N(u) := I_G_N(ω) I_A_N(u) | _N(u)-_N^0(u) |.

To deduce the required convergence

    M_N:=∫ M_N(u) u  0,

we make use of the subsequence argument: fix any infinite sequence '⊂. In view of the estimate

    M_N(u)    ≤ I_G_N(ω) I_A_N(u)( _N(u) + _N^0(u)) 
       ≤ e^M|u| - |u|^2 + e^M|u| - 2|u|^2≲ e^M|u| - |u|^2    a.s.

and the dominated convergence theorem, it suffices to show that there exists a further subset ”={N”}⊂' along which M_N”(u) → 0 a.s. for each u.
We have


    | log_N(u) - log^0_N(u) |
    ≲ |u|^2(
    | -D_N^-1_θ^2_N()D_N^-1 - I() | + r̅_N_N) =: |u|^2R_N.

Since | -D_N^-1_θ^2_N()D_N^-1 - I() | ∨r̅_N_N→ 0, 
it is possible to pick a further subset ”={N”}⊂' along which R_N”→ 0 a.s.;
note that the random sequence (R_N) is free from the variable u.
With this ”, for each u we have

    M_N”(u)    ≤_N”^0(u)| exp{log_N(u) - log^0_N(u)}-1 | 
       ≤_N”^0(u) |u|^2R_N” e^|u|^2R_N”→ 0    a.s.,

followed by M_N”→ 0 a.s., hence by (<ref>). We conclude that '_N 0.



Now we turn to the proof of ”_N 0, still focusing on the event G_N of (<ref>).
Note that A_N^c↓∅ under (<ref>) since |D_N^-1u|>_N implies that |u| ≳_N√(N)h^1-1/β↑∞.
Under the condition (<ref>) we have _N/|D_N^-1| ∼ C_N√(N)h^1-1/β→∞, so that

    ∫_A_N^c_N^0(u) u ≤∫_|u| ≥_N/|D_N^-1|exp(M|u| - 2|u|^2)  u → 0

on G_N.
Since π is bounded, we are left to show that 

    ∫_A_N^c_N(u) u  0.

Write u=(u_2,u_1)∈^p_×^p_ and observe that (recall (<ref>) and (<ref>))

    _N(u) ≤(sup__1,N(u_1;)) _2,N(u_2) =: _1,N(u_1)_2,N(u_2).

Since A_N^c⊂{u; |u_2|≥ (_N/2)√(N)h^1-1/β}∪{u; |u_1|≥ (_N/2)√(N)}, we obtain

    ∫_A_N^c_N(u) u    ≤∫_|u_1|≥(_N/2)√(N)_1,N(u_1) u_1·∫_2,N(u_2) u_2
            + ∫_1,N(u_1) u_1·∫_|u_2|≥ (_N/2)√(N)h^1-1/β_2,N(u_2) u_2.

We will only show that

    ∫_|u_1|≥(_N/2)√(N)_1,N(u_1) u_1    0,
    
    ∫_1,N(u_1) u_1   = O_p(1).

Indeed, in view of the definition (<ref>) of the good event G_N, in order to deduce

    ∫_|u_2|≥(_N/2)√(N)h^1-1/β_2,N(u_2) u_2 0   and  ∫_2,N(u_2) u_2 = O_p(1),

we can follow exactly the same route as in the proofs of (<ref>) and (<ref>) below,
with _1,N(θ) and _1() replaced by _2,N() and _2(), respectively.

On G_N we have _N^-2sup_;  |-_0|≥_N/2_1()≤ -, so that

    sup_|u_1|≥ (_N/2)√(N)log_1,N(u_1)
       ≲ N ( sup_θ|_1,N(θ)-_1()| +sup_;  |-_0|≥_N/2_1()) 
       ≤ -/2(√(N)_N)^2≲ -N^1-2c.

Therefore, recalling that _1,N(u_1)=0 outside the set _1,N:=√(N)(Θ_-_0)⊂^p_ with Θ_ being bounded, we obtain the following estimate for some positive constants C_0 and C_1: on G_N,

    ∫_|u_1|≥(_N/2)√(N)_1,N(u_1) u_1   ≤ e^-C_0N^1-2c∫_{|u_1|≥ (_N/2)√(N)}∩_1,N u_1
       ≲ e^-C_0N^1-2c∫_|u_1|≤ C_1√(N) u_1≲ N^p_/2e^-C_0N^1-2c→ 0.

Hence (<ref>) is obtained.

It follows from (<ref>) that there exist N_1∈ and K>0 such that for every M_1>1,

    sup_N≥ N_1( G_N∩{∫_1,N(u_1) u_1 > M_1}) 
       ≤sup_N≥ N_1( G_N∩{∫_|u_1|≥ K_1,N(u_1) u_1 > M_1/2}) 
                + sup_N≥ N_1( G_N∩{sup_|u_1|≤ K_1,N(u_1)≳M_1/2K^p_}) 
       ≤sup_N≥ N_1( sup_|u_1|≤ K_1,N(u_1)≳M_1/2K^p_).

For every K>0 the sequence {sup_|u_1|≤ K_1,N(u_1)}_N is tight in , hence we can make the last upper bound arbitrarily small by taking a sufficiently large M_1.
This verifies (<ref>), and we are left to deduce (<ref>).



Proof of (<ref>).
Let

    a_n-1(α) = a(X _ (n - 1) h, α), 
    c_n-1(γ) = c(X _ (n - 1) h, γ),  Δ^N_n X = X _ nh - X _ (n-1)h,
       ϵ_n(θ) = ϵ_N,n(θ) := Δ^N_nX - h a_n-1(α)/h ^ 1 / βc_n-1(γ).

Let z_n:=h^-1/β(J_nh-J_(n-1)h)∼ J_1 and b_n-1(θ):=c^-1_n-1(){a_n-1(_0)-a_n-1()}.
Below we will repeatedly make use of several statements in <cit.>, hence at this stage it should be noted that, by the localization procedure <cit.>, without loss of generality we may and do suppose that

    𝔼(sup_t≤ T|X_t|^q)≲ 1,  sup_t∈ [s,s+h]∩ [0,T]𝔼(|X_t-X_s|^q|ℱ_s)≲ h(1+|X_s|^C)

for any q≥ 2 and s∈ [0,T].







This fact in particular implies that 

    sup_t∈ [s,s+h]∩ [0,T]𝔼(|X_t-X_s|)≲√(h).

For convenience, for a sequence of random functions {f_N(·)} on Θ and a positive sequence (b_N) we will write f_N(θ)=O^∗_p(b_N) if sup_θ|f_N(θ)|=O_p(b_N).



We can write _1,N(θ) as

    _1,N(θ)    = 1/N∑_n=1^Nlogc_n-1(_0)/c_n-1()
    +1/N∑_n=1^Nlogϕ_β(_n(θ))
    -1/N∑_n=1^Nlogϕ_β(_n(,_0)).

The function y↦logϕ_β(y) fulfills the conditions on η in Lemmas 6.2 and 6.3 of <cit.>.
Applying the two lemmas to the second and third terms in the right-hand side of (<ref>), we can deduce that

    _1,N(θ)    = 1/N∑_n=1^Nlogc_n-1(_0)/c_n-1()
    +1/N∑_n=1^N{logϕ_β(
    c_n-1(_0)/c_n-1()z_n+h^1-1/βb_n-1(θ))|_(n-1)h}
            -1/N∑_n=1^N{logϕ_β(
    z_n+h^1-1/βb_n-1(,_0))|_(n-1)h} 
    +O^∗_p(1/√(N)∨ h^2-1/β).

Proceeding as in the second equality in Eq.(6.19) of <cit.>, we obtain

    _1,N(θ)    = 1/N∑_n=1^N[
    logc_n-1(_0)/c_n-1()
    +∫{logϕ_β(c_n-1(_0)/c_n-1()z) - logϕ_β(z)}ϕ_β(z)dz]
    
            +O^∗_p(h^1-1/β∨1/√(N)∨ h^2-1/β).

Using the estimate (<ref>) combined with the inequality given in the proof of <cit.>, we can deduce from (<ref>) that

    sup_θ|_1,N(θ)-_1()|
       ≤ O_p(√(h)∨ h^1-1/β∨1/√(N)∨ h^2-1/β) 
       = O_p(h^1-1/β)=O_p(N^-(1-1/β)).

Hence (√(N))^qsup_θ|_1,N(θ)-_1()|  0 holds for q∈(0, 2(1-1/β)).

As for _2,N, following a similar line to the case of _1,N we can derive

    _2,N()    = 1/2N∑_n=1^Nb^2_n-1(,_0)
    (_y^2logϕ_β(y))|_y=_n()
    + O^∗_p(1/√(N)h^1-1/β∨ h^1-1/β).

Here, we also made use of <cit.> (the function y↦_ylogϕ_β(y) is odd) and the arguments in <cit.>.
Then it is not difficult to arrive at sup_|_2,N()-_2()|=O_p(N^-s) for any sufficiently small s>0. In view of (<ref>) we conclude that (√(N)h^1-1/β)^qsup_|_2,N()-_2()|  0 for q>0 small enough. The proof of (<ref>) is thus complete.






§ PROOF FOR THE LOCAL CONSISTENCY OF THE METROPOLIS-HASTINGS ALGORITHM



We prove Proposition <ref> in this section. 
Without loss of generality, we can assume s_N→ s in probability for some random variable s. 
For notational simplicity, we write F̅^N for the rescaled version of the function or measure F^N by u↦θ_0+D_N^-1u, and write 

    ν̅_N( u, v|x^N)=Π̅( u|x^N)Q̅_N(u, v|x^N),    μ̅_N( u, v|x^N)=Π̅( u|x^N)P̅_N(u, v|x^N)

and 

    ν( u, v|s)=Π( u|s)Q(u, v|s),    μ( u, v|s)=Π( u|s)P(u, v|s).

The equation (<ref>) becomes 

    ν̅_N(·|x^N)-ν(·|s_N)_TV→ 0.

By Lemmas 2 and 3 of <cit.> the following convergence is sufficient for local consistency:

    δ_N:=μ̅_N(·|x^N)-μ(·|s)_TV=o_p(1).

First, we prove the convergence of δ_N':=μ(·|s_N)-μ(·|s)_TV.
By triangular inequality, δ_N' is dominated above by the sum of 

    δ_N,1':=   ν(·|s)A(·|s)-ν(·|s_N)A(·|s_N)_TV,
    δ_N,2':=   Π( u|s)R(u|s)δ_u( v)-Π( u|s_N)R(u|s_N)δ_u( v)_TV.

For the former, by Assumption <ref>, 

    δ_N,1'=2∫(π(u|s)q(u,v|s)A(u,v|s)-π(u|s_N)q(u,v|s_N)A(u,v|s_N))^+ u v=o_p(1)

by the dominated convergence theorem, where x^+=max{0,x}. We can prove δ_N,3':=ν(·|s)-ν(·|s_N)_TV=o_p(1)
in the same way. 
On the other hand, by triangular inequality, 

    δ_N,2'   =Π( u|s)R(u|s)-Π( u|s_N)R(u|s_N)_TV
       ={ν(·×ℝ^p|s)-∫ν(·× v|s)A(·,v|s)}-{ν(·×ℝ^p|s_N)-∫ν(·× v|s_N)A(·,v|s_N)}_TV
       ≤δ_N,3'+δ_N,1'=o_p(1).

Hence, δ_N'≤δ_N,1'+δ_N,2'→ 0 in probability. 

Next, we prove the convergence of δ_N”:=μ̅_N(·|x^N)-μ(·|s_N)_TV.
By the same argument as above, it is sufficient to show the convergence of 
δ_N,1”:=ν̅_N(·|x^N)A̅_N(·|x^N)-ν(·|s_N)A(·|s_N)_TV
and δ_N,3”:=ν̅_N(·|x^N)-ν(·|s_N)_TV. 
We only proves the former. 
By triangular inequality, 

    δ_N,1”   ≤{ν̅_N(·|x^N)-ν(·|s_N)}A̅_N(·|x^N)_TV+ν(·|s_N)(A̅_N(·|x^N)-A(·|s_N))_TV
       ≤ν̅_N(·|x^N)-ν(·|s_N)_TV+∫ |A̅_N(u,v|x^N)-A(u,v|s_N)|ν( u, v|s_N).

The first term converges to 0 by assumption. 
Since ν(·|s_N) converges to ν(·|s), 
and s_N is tight, 
for any ϵ>0, we can choose a compact set K_p⊂ℝ^p×ℝ^p and K_q⊂𝒮 so that 

    lim sup_N→∞𝐄_θ_0^(N)[ν(K_p^c|s_N)]<ϵ/2,    lim sup_N→∞𝐏_θ_0^(N)(s_N∈ K_q^c)<ϵ/2.

Thus, 

    𝐄_θ_0^(N)[∫ |A̅_N(u,v|x^N)-A(u,v|s_N)|ν( u, v|s_N)]
       ≤ϵ + 𝐄_θ_0^(N)[∫_K_p |A̅_N(u,v|x^N)-A(u,v|s_N)|ν( u, v|s_N), s_N∈ K_q] + o(1)
       ≤ϵ + c∫_K_p𝐄_θ_0^(N)[|A̅_N(u,v|x^N)-A(u,v|s_N)|] u v+ o(1)=ϵ +o(1),

where c=sup_s∈ K_qsup_u,vπ(u|s)q(u,v|s) is the upper bound of the probability density function of ν( u, v|s_N) when s_N∈ K_q. 
This completes the proof of δ_N”→ 0 in probability, and hence 
δ_N≤δ_N'+δ_N”=o_p(1).  




§ CONVERGENCE OF THE ACCEPTANCE RATIO





 §.§ Setting and notation


We keep using some notation introduced in Section <ref>. 
We consider 
an extended probability space

    Ω̃ = Ω× A, ℱ̃ = ℱ⊗𝒜, ℙ̃ = ℙ(ω)ℚ_ω( a)

where (A,𝒜) is a measurable space, and ℚ_ω( a)
is a probability transition kernel from Ω to A. We consider a stochastic process

    X̃_t(ω, a) = X _ t(ω)=ω_t.

Fix u ∈ℝ ^ d. Set 
θ_N=(α_N,γ_N)=θ_0+D_N^-1u. 
We now consider random variables V_n, n=1,…, N, which correspond to the pseudo-data generated in Algorithm <ref> from parameter 
θ=θ_N.  Thus, for each ω∈Ω, the random variables V_n(ω,·):A→ℝ_+ (n = 1, …, N)
are independent, and 

    ℙ̃( V _ n ∈ v |ω) = ℚ_ω(V_n∈ v)=F_β( v|ϵ_n(θ_N)).








The log likelihood and the augmented-data (pseudo + observed data) log likelihood are

    ℍ_N(θ)=-∑_n=1^Nlog c_n-1(γ)+∑_n=1^Nlogϕ_β(ϵ_n(θ)), ℍ^†_N(θ)=
    -∑_n=1^Nlog c_n-1(γ)-1/2∑_n=1^Nϵ_n(θ)^2/V_n.

Here we omit terms that do not include θ. 
We define the Fisher information matrix for the augmented-data model by 

    I ^ †(θ _ 0)	:=
    diag(
    C _ α ^ † (β) Σ _ T, α(θ _ 0), 
    C _ γ ^ † (β) Σ _ T, γ(γ _ 0)
    ),

with

    C_α^†(β):=∫_ℝ_+1/vF_β( v),    C^†_γ(β):=2

and that for the pseudo-data model by 

    I ^ *(θ _ 0)	:=I ^ †(θ _ 0)-I(θ _ 0)=
    diag(
    C _ α ^ * (β) Σ _ T, α(θ _ 0), 
    C _ γ ^ * (β) Σ _ T, γ(γ _ 0)
    ),
 
with

    C_α^*(β):=C_α^†(β)-C_α(β),    C^*_γ(β):=C^†_γ(β)-C_γ(β).















In the next section, we will use the following law of large numbers. 


[Proposition 6.5 of <cit.>]
If η(x)∈𝒞^1(ℝ) and π(x,θ)∈𝒞^1(ℝ×Θ) satisfy sup_x∈ℝ|η(x)|+|η'(x)|≲ 1 and sup_θ|π(x,θ)|+|∂_θπ(x,θ)|≲ 1+|x|^C for some C>0, then

    N^-1∑_n=1^Nη(ϵ_n(θ_N))π(X_(n-1)h,θ_N)=O_p(1), 
        N^-1∑_n=1^Nη(ϵ_n(θ_N))π(X_(n-1)h,θ_N)= 1/T∫η(x)ϕ_β(x) x∫_0^Tπ(X_t,θ_0) t+o_p(1).












We also use the following fact: Let β∈ (1,2). If π(x,θ) satisfies the condition in Proposition <ref>, then

    1/Nh^1-1/β∑_n=1^N g_β(ϵ_n(θ_N))π(X_(n-1)h,θ_N)=O_p((√(N)h^1-1/β)^-1).

This convergence comes from Corollary 6.6. of <cit.> together with the estimate

    |g_β(ϵ_n(θ_N))-g_β(ϵ_n(α_0,γ_N))|≲ |ϵ_n(θ_N)-ϵ_n(α_0,γ_N))|≲1/√(N)h^1-1/β(1+|X_(n-1)h|^C).





 §.§ Some properties of F_β


Let 

    g_β(x):=∂/∂ xlogϕ _ β(x), 
    	h_β(x):=∂^2/∂ x^2logϕ_β(x)-1/x∂/∂ xlogϕ_β(x).

Recall that by the property of stable law (see pp.88–89 of <cit.>), we
have ϕ_β(x)∼ c|x|^-(β +1) as  |x|→∞ for some c>0.  
Moreover, the probability density function 
f_β(v) of F_β is bounded above, and as v→+∞, we have  f_β(v)∼ c|v|^-β/2-1 for some c>0, 
and as v→ 0 the density of the positive stable distribution f_β(v) converges to 0 exponentially fast. 
Thus ∫ v^-kF_β( v|x) is continuous at x=0 for k≥ 0. 
Moreover, we have the following identities and an estimate. 


Fix β∈ [1,2). 


  * 
    g_β(x)=-∫x/vF_β( v|x), 
    h_β(x)=∫(x/v+g_β(x))^2F_β( v|x).

	
  * 
    ∫ x g_β(x)ϕ_β(x) x=-1, ∫ h_β(x)ϕ_β(x) x=C_α^*(β), ∫ x^2h_β(x)ϕ_β(x) x=C_γ^*(β).

 
  * For any k≥ 0, 

    ∫_0^∞ v^-kF_β( v|x)≲ |x|^-2k.





The expression of g_β(x) can be obtained via simple interchange of the derivative and the integral in the equation (<ref>). 
For the expression of h_β(x), we have

    h_β(x)   =ϕ_β”(x)/ϕ_β(x)-g_β(x)^2-1/xg_β(x).

By (<ref>), the second derivative of ϕ_β(x) is

    ∂^2/∂ x^2∫_ℝ_+1/√(2π v)exp(-x^2/2v)F_β( v)
       =
    ∫_ℝ_+(x^2/v^2-1/v)1/√(2π v)exp(-x^2/2v)F_β( v)
       =ϕ_β(x)∫_ℝ_+(x^2/v^2-1/v)F_β( v|x).

This equation yields the expression in (<ref>) by using the identity ∫ v^-1F_β( v|x)=-g_β(x)/x. 

Next, we prove identities (<ref>). Applying the change of variable u=x/√(v), we have

    ∫ x g_β(x)ϕ_β(x) x   =-∫ x{x/vF_β( v|x)}ϕ_β(x) x
       =-∫x^2/v1/√(2π v)exp(-x^2/2v) xF_β( v)=
    -∫ u^2ϕ(u) u=-1

where ϕ(u) is the probability density function of the standard normal distribution. 
In the same way, by the change of variable, we obtain

    ∫(x/v)^2F_β( v|x)ϕ_β(x) x=C_α^†(β), ∫ x^2(x/v)^2F_β( v|x)ϕ_β(x) x=3.

Then we have

    ∫ h_β(x)ϕ_β(x) x=∫(x/v)^2F_β( v|x)ϕ_β(x) x-∫ g_β(x)^2ϕ_β(x) x= C_α^*(β)

and 

    ∫ x^2h_β(x)ϕ_β(x) x   =∫ x^2(x/v)^2F_β( v|x)ϕ_β(x) x-∫ x^2g_β(x)^2ϕ_β(x) x
       =3-∫ x^2g_β(x)^2ϕ_β(x) x=2-∫ (1+xg_β(x))^2ϕ_β(x) x=C_γ^*(β).

Finally we check (<ref>). By the property of stable law, 

    ∫_0^∞x^2k/v^kF_β( v|x)   =ϕ_β(x)^-1∫_0^∞x^2k/v^k1/√(2π v)exp(-x^2/2v)F_β( v)
       ≲∫_0^∞(x^2/v)^(2k+β+1)/2exp(-x^2/2v) v/v
       = ∫_0^∞(1/u)^(2k+β+1)/2exp(-1/2u) u/u<∞

by the change of variable u =v/x^2. 
Thus, the claim follows. 














 §.§ Estimates for the likelihood functions


Let Δ_N^†(θ_N)=D_N^-1∂_θℍ^†_N(θ_N)
and set Δ_N^*(θ_N)=Δ_N^†(θ_N)-Δ_N(θ_N). 


	
    ℒ(Δ_N^*(θ_N)|ω)- N_k(0, I^*(θ_0))_BL→ 0

	in ℙ-probability. 



	Observe that
	
    ∂_αϵ_n(θ)=-h^1-1/β∂_α a_n-1(α)/c_n-1(γ), ∂_γϵ_n(θ)=-ϵ_n(θ)∂_γ c_n-1(γ)/c_n-1(γ).

	By this fact, we have
	
    Δ_N^†(θ_N)=
    		N^-1/2∑_n=1^N
    		[             ϵ_n(θ_N)/V_n∂_α a_n-1(α_N)/c_n-1(γ_N); {ϵ_n(θ_N)^2/V_n - 1
    			}∂_γ c_n-1(γ_N)/c_n-1(γ_N) ],

		and hence 
	
    Δ_N^*(θ_N)=
    		N^-1/2∑_n=1^N
    		[            {ϵ_n(θ_N)/V_n	
    			+g_β(ϵ_n(θ_N))}∂_α a_n-1(α_N)/c_n-1(γ_N); {ϵ_n(θ_N)/V_n+g_β(ϵ_n(θ_N))
    			}ϵ_n(θ_N)
    			∂_γ c_n-1(γ_N)/c_n-1(γ_N) ]
    		=:N^-1/2∑_n=1^Nξ_n.

		Recall the definition of the extended probability space (<ref>). In this setting, the pseudo-data variables V_1,…, V_N are independent conditioned on ω∈Ω and generated from ℚ_ω. 
		Fix ω∈Ω. Since ∫ (x/v+g_β(x))F_β( v|x)=0,  the random variable Δ_N^*(θ_N)(ω,·):A→ℝ^p is a sum of independent variables  ξ_n(ω,·) in ℚ_ω. 
		By the expression of h_β, the covariance matrix of ξ_n conditioned on ω becomes 
		
    𝔼̃[ξ_n^⊗ 2|ω]=
    		[                                 h_β(ϵ_n(θ_N))
    			(∂_α a_n-1(α_N)/c_n-1(γ_N))^⊗ 2 h_β(ϵ_n(θ_N))ϵ_n(θ_N)
    			(∂_α a_n-1(α_N)/c_n-1(γ_N))⊗(∂_γ c_n-1(γ_N)/c_n-1(γ_N));                                                                             Sym.                      h_β(ϵ_n(θ_N))ϵ_n(θ_N)^2
    			 (∂_γ c_n-1(γ_N)/c_n-1(γ_N))^⊗ 2 ].

		This conditional covariance can be written as 𝔼̃[ξ_n^⊗ 2|ω]=η(ϵ_n(θ_N))[π(X_(n-1)h,θ_N)] where 
		
    η(x)=h_β(x)[   1   x;   x x^2 ], π(x,θ)=[∂_α a_n-1(α_N)/c_n-1(γ_N),∂_γ c_n-1(γ_N)/c_n-1(γ_N)]^⊗ 2.

		By (<ref>) and Assumption <ref>, η(x) and π(x,θ) satisfy the condition of Proposition <ref>. On the other hand,  for the forth moment on ξ_n, we have
		
    𝔼̃[|ξ_n|^4|ω]^1/4≲   {∫|ϵ_n(θ_N)/v|^4F_β( v|ϵ_n(θ_N))}^1/4|∂_α a_n-1(α_N)/c_n-1(γ_N)|
       +
    				{∫|ϵ_n(θ_N)^2/v|^4F_β( v|ϵ_n(θ_N))}^1/4|∂_γ c_n-1(α_N)/c_n-1(γ_N)|
    
    ≲    1+|X_(n-1)h|^C

		for some C>0 by (<ref>) and Assumption <ref> together with 
		the Minkowski and Jensen inequalities. 
		Therefore by Proposition <ref>, 
		
    𝔼̃[(Δ_N^*(θ_N))^⊗ 2|ω]=N^-1∑_n=1^N𝔼̃[ξ_n^⊗ 2|ω]→ I^*(θ_0), 
    			N^-2∑_n=1^N𝔼̃[|ξ_n|^4|ω]→ 0

		in ℙ̃-probability. For any infinite elements ℕ'⊂ℕ there exists a further subsequence ℕ”={N”}⊂ℕ'
		such that the above convergence satisfies in almost surely in ω∈Ω. Therefore, by the central limit theorem, 
		
    ℒ(Δ_N”^*(θ_N)|ω)→ N_k(0, I^*(θ_0))

		almost surely. Thus, the claim follows. 



	
    -D_N^-1∂^2_θℍ^†_N(θ_N)D_N^-1→ I^†(θ_0)

	in ℙ̃-probability, and 
	
    sup_Θmax_k=1,…, p|∂_θ_kD_N^-1∂^2_θℍ^†_N(θ)D_N^-1|=O_p(1).




By calculation, 

    (Nh^2(1-β))^-1∂_α^2 ℍ_N^† (θ_N)	
    =1/Nh^1-1/β∑_n=1^Nϵ_n(θ_N)/V_n(∂_α^2 a_n-1(α_N)/c_n-1(γ_N))-1/N∑_n=1^N1/V_n(∂_α a_n-1(α_N)/c_n-1(γ_N))^⊗ 2.

The first term of the right-hand side of (<ref>) is 

    -1/Nh^1-1/β∑_n=1^Ng_β(ϵ_n(θ_N))(∂_α^2 a_n-1(α_N)/c_n-1(γ_N))
    +1/Nh^1-1/β∑_n=1^N(ϵ_n(θ_N)/V_n + g_β(ϵ_n(θ_N)))(∂_α^2 a_n-1(α_N)/c_n-1(γ_N)).

Then, the first term is negligible by Proposition <ref> for β=1 and by (<ref>) for β∈ (1,2), and the second term is also negligible since it is a sum of independent variables conditioned on ω, and its variance is o_p(1). 
Also, the second term of (<ref>) is 

    1/N∑_n=1^N(g_β(ϵ_n(θ_N))/ϵ_n(θ_N))(∂_α a_n-1(α_N)/c_n-1(γ_N))^⊗ 2
    -
    1/N∑_n=1^N(1/V_n+g_β(ϵ_n(θ_N))/ϵ_n(θ_N))(∂_α a_n-1(α_N)/c_n-1(γ_N))^⊗ 2,

and again, the second term is negligible since it is a sum of conditionally independent random variables. Thus, we obtain

    -(Nh^2(1-β))^-1∂_α^2 ℍ_N^† (θ_N)	
       = -1/N∑_n=1^N(g_β(ϵ_n(θ_N))/ϵ_n(θ_N))(∂_α a_n-1(α_N)/c_n-1(γ_N))^⊗ 2+o_p(1)
    
       = C_α^†(β)Σ_T,α(θ_0)+o_p(1)

by Proposition <ref>. 

Similarly, we have 

    (Nh^(1-β))^-1∂_α∂_γℍ_N^† (θ_N)	
    =   -2/N∑_n=1^Nϵ_n(θ_N)/V_n(∂_α a_n-1(α_N)/c_n-1(γ_N))⊗(∂_γ c_n-1(γ_N)/c_n-1(γ_N))
    
    =   2/N∑_n=1^Ng_β(ϵ_N(θ_N))(∂_α a_n-1(α_N)/c_n-1(γ_N))⊗(∂_γ c_n-1(γ_N)/c_n-1(γ_N))
       -2/N∑_n=1^N(g_β(ϵ_N(θ_N))+ϵ_n(θ_N)/V_n)(∂_α a_n-1(α_N)/c_n-1(γ_N))⊗(∂_γ c_n-1(γ_N)/c_n-1(γ_N)).

The first term is o_p(1) by Proposition <ref> together with the fact ∫ g_β(x)ϕ_β(x) x=0, 
and the second term converges to 0 in L^2 since it is a sum of conditionally independent random variables. Then 
(Nh^(1-β))^-1∂_α∂_γℍ_N^† (θ_N)	→ 0 in probability. 

In the same way, decomposing the sum into the main term and the sum of independent variables, we have

    N^-1∂_γ^2 ℍ_N^† (θ_N)	
    =   -2/N∑_n=1^Nϵ_n(θ_N)^2/V_n(∂_γ c_n-1(γ_N)/c_n-1(γ_N))^⊗ 2
       +
    1/N∑_n=1^N(ϵ_n(θ_N)^2/V_n-1)
    {(∂_γ^2 c_n-1(γ_N)/c_n-1(γ_N))
    -(∂_γ c_n-1(γ_N)/c_n-1(γ_N))^⊗ 2}
    
    =   2/N∑_n=1^Nϵ_n(θ_N)g_β(ϵ_n(θ_N))
    (∂_γ c_n-1(γ_N)/c_n-1(γ_N))^⊗ 2
       +
    1/N∑_n=1^N(-ϵ_n(θ_N)g_β(ϵ_n(θ_N))-1)
    {(∂_γ^2 c_n-1(γ_N)/c_n-1(γ_N))
    -(∂_γ c_n-1(γ_N)/c_n-1(γ_N))^⊗ 2}
       +o_p(1).

The first term in the right-hand side converges to -C_γ^†(β)Σ_T,γ(θ_0)
and the second term is o_p(1) by Lemma <ref> together with the fact that ∫ xg(x)ϕ_β(x) x=-1. 
Then, we obtain

    -N^-1∂_γ^2 ℍ_N^† (θ_N)		→ C_γ^†(β)Σ_T,γ(θ_0)

in probability. Almost the same arguments give convergence of the thrice derivatives. We omit the detail. 






 §.§ Proof of Theorem <ref>


The proof of Theorem <ref> is an application of Proposition <ref>. Since Π(u|s)=N_p(I(θ_0)^-1s, I(θ_0)^-1) and Q(u, v|s)=N_p(u,Σ),  Assumption <ref> is easy to check. Ergodicity is also obvious by irreducibility. 
The Bernstein-von Mises theorem is already proved in Theorem <ref> which is a sufficient condition for (<ref>) in this case. 
We will complete the proof  by showing the convergence of the acceptance ratio A_N. 


Let θ_N=θ_0+D_N^-1u and θ_N^*=θ_0+D_N^-1v. 
By (<ref>) and (<ref>) together with Taylor's expansion, we have the difference of the log quasi-posterior density satisfies 

    ξ_N    :={ℍ_N^†(θ^*_N)+logπ(θ_N^*)}-{ℍ_N^†(θ_N)+logπ(θ_N)}
       =∂_θℍ^†_N(θ_N^*-θ_N)
    	+1/2∂_θ^2ℍ_N^†[(θ_N^*-θ)^⊗ 2]+o_p(1)
       =Δ^†_N(θ_N)[v-u]-1/2I^†(θ_0)[(v-u)^⊗ 2]+o_p(1).

Now we rewrite the two terms in the right-hand side of the above equation. 
First, observe that difference of the score function at θ_N and θ_0 is 

    Δ_N(θ_N)-Δ_N(θ_0)=D_N^-1∂_θ^2ℍ_N(θ_0)D_n^-1[u]+o_p(1)=-I(θ_0)[u]+o_p(1)

by Taylor's expansion. Also, by I^†(θ_0)=I^*(θ_0)+I(θ_0), we have the following identity among Fisher information matrices: 

    I(θ_0)[u, v-u] +1/2I^†(θ_0)[(v-u)^⊗ 2]=I(θ_0)[v^⊗ 2-u^⊗ 2]+1/2I^*(θ_0)[(v-u)^⊗ 2].

By these facts, we can rewrite ξ_N by

    ξ_N=   Δ_N(θ_N)[v-u]+Δ(θ_N)[v-u]-1/2I^†(θ_0)[(v-u)^⊗ 2]+o_p(1)
    
    =   Δ_N(θ_N)[v-u]+Δ(θ_0)[v-u]-I(θ_0)[u,v-u]-1/2I^†(θ_0)[(v-u)^⊗ 2]+o_p(1)
    
    =   η_N(Δ_N^*(θ_N))+o_p(1)

where 

    η_N(w)=Δ_N(θ_0)[v-u]+w[v-u]-I(θ_0)[v^⊗ 2-u^⊗ 2]-1/2I^*(θ_0)[(v-u)^⊗ 2]

for w∈ℝ^p. By the expression of the difference of the log quasi-posterior density together with the function ψ(x)=min{1,exp(x)}, we have the expression of the acceptance probabilities: 

    A_N(θ_N,θ^*_N|X^N)=𝔼̃[ψ(ξ_N)|ω], A(u, v|Δ_N(θ_0), I(θ_0), I^*(θ_0))=𝔼[ψ(η_N(W))|w]

where W∼ N(0, I^*(θ_0)). 
On the other hand,  we have

    |𝔼̃[ψ(ξ_N)|ω]-𝔼[ψ(η_N(W))|w]|
    	≤   |𝔼̃[ψ(ξ_N)|ω]-𝔼[ψ(η_N(Δ_N^*(θ_N))|w]|
       + |𝔼̃[ψ(η_N(Δ_N^*(θ_N))|ω]-𝔼[ψ(η_N(W))|w]|
    ≤    o_p(1)+ℒ(Δ^*_N,T(θ_N)|w)-N(0, I^*(θ_0))_BL=o_p(1).

Hence the claim follows. 



99    

bns
Barndorff-Nielsen, O. E. & Shephard, N. (2001). Non-Gaussian OU-based models and some of their
uses in financial economics (with discussion). J. Roy. Statist. Soc. B, 63, 167–241.

BKP71
Borwanker, J., Kallianpur, G., & Rao, B. P. (1971). The Bernstein-von Mises theorem for Markov processes. The Annals of Mathematical Statistics, 43, 1241–1253.

CleGlo15
Clément, E. and Gloter, A. (2015). Local asymptotic mixed normality property for discretely observed stochastic differential equations driven by stable Lévy processes. Stochastic Process. Appl. 125, no. 6, 2316–2352.

cont
Cont, R. & Tankov, P. (2004). Financial modelling with Jump Processes. Chapman & Hall: London.

2015arXiv151104992D
Deligiannidis, G., Doucet, A., and Pitt, M. K. (2015). 
The Correlated Pseudo-Marginal Method.
arXiv:1511.04992v2 (v2). 

gander
Gander, M. P. S. & Stephens, D. A. (2007). Simulation and inference for stochastic volatility models
driven by Lévy processes. Biometrika,  94, 627–646.

golight
Golightly, A., & Wilkinson, D. J. (2008). Bayesian inference for nonlinear multivariate diffusion models observed with error.
Comp. Stat. Data Anal., 52, 1674–1693.

griffin
Griffin, J. & Steel, M. (2006). Inference with non-Gaussian Ornstein-Uhlenbeck processes for stochastic
volatility. J. Econom., 134, 605–644.

IH
Ibragimov, I. A. & Has'minskii, R. Z. (1981).
Statistical estimation, asymptotic theory.
New York : Springer-Verlag.

JacPro11
Jacod, J., and Protter, P. E. (2011). Discretization of processes. Springer.

JW94
Janicki, A. & Weron, A. (1994).
Simulation and chaotic behavior of -stable stochastic processes.
Monographs and Textbooks in Pure and Applied Mathematics, 178. Marcel Dekker, Inc., New York.

jasra
Jasra, A.,  Stephens, D. A., Doucet, A. & Tsagaris, T. (2011).
Inference for Lévy driven stochastic volatility models via adaptive
sequential Monte Carlo. Scand. J. Statist., 38, 1–22.

Kamatani10
Kamatani, K. (2014).
Local consistency of Markov chain Monte Carlo methods. Ann. Inst. Statist. Math., 66, 63–74.

2016arXiv160202889K
Kamatani, K. (2017).
Ergodicity of Markov chain Monte Carlo with reversible proposal.
Journal of Applied Probability 54.

MatTak06
Matsui, M. & Takemura, A. (2006). 
Some improvements in numerical evaluation of symmetric stable density
  and its derivatives.
Comm. Statist. Theory Methods, 35(1-3):149–172.

marin
Marin, J.-M., Pudlo, P., Robert, C.P. & Ryder, R. (2012). Approximate Bayesian computational methods. 
Statist. Comp., 22, 1167–1180.

Mas09
Masuda, H. (2009). Joint estimation of discretely observed stable Lévy processes with symmetric Lévy density. J. Japan Statist. Soc. 39, no. 1, 49–75.

masuda13
Masuda, H. (2013).
Convergence of Gaussian quasi-likelihood random fields for ergodic Lévy driven SDEs observed at high-frequency. Ann. Statist.,  41, 1593–1641.

masuda16.sqmle
Masuda, H. (2017).
Non-Gaussian quasi-likelihood estimation of SDE driven by locally stable Lévy process. arXiv:1608.06758 (v3)

MR1723510
Neal, R. M. (1999).
Regression and classification using Gaussian process priors.
In Bayesian statistics, 6 (Alcoceber, 1998).
Oxford Univ. Press, New York pp. 475–501.

peters
Peters G.W., Sisson S.A. & Fan Y. (2012). Likelihood-free Bayesian inference for α-stable models. 
Comp. Stat. Data Anal., 56, 3743–3756.




fbasics
Rmetrics Core Team, Wuertz, D., Setz, T., and Chalabi, Y.
fBasics: Rmetrics - Markets and Basic Statistics, 2014.
R package version 3011.87.

roberts
Roberts, G. O. & Stramer, O. (2001). On inference for partially observed nonlinear diffusion models using the Metropolis–Hastings algorithm.
Biometrika, 88, 603–621.

SamTaq
Samorodnitsky, G. & Taqqu, M. S. (1994). 
Stable non-Gaussian random processes. Stochastic models with infinite variance. Chapman & Hall, New York.

MR1739520
Sato, K. (1999).
Lévy processes and infinitely divisible distributions,
 volume 68 of Cambridge Studies in Advanced Mathematics.
Cambridge University Press, Cambridge.





Yos11
Yoshida, N. (2011). Polynomial type large deviation inequalities and quasi-likelihood analysis for stochastic differential equations. Ann. Inst. Statist. Math., 63, 431–479.









