This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

11affiliationtext: Universidad Adolfo Ibáñez. Faculty of Engineering & Sciences. Viña del Mar, Chile22affiliationtext: Data Observatory Foundation, ANID Technology Center No. DO210001, Eliodoro Yáñez 2990, 7510277, Providencia, Santiago, Chile33affiliationtext: Universidad de La República. IMERL. Uruguay44affiliationtext: McGill University, Canadathanks: gerardo.martinez@mail.mcgill.ca55affiliationtext: Escuela de Matemática UCV. Venezuela

Small-time approximation of the transition density for diffusions with singularities. Application to the Wright-Fisher model.

T.Roa tania.roa@uai.cl M. Fariello fariello@fing.edu.uy G. Martínez J.R. León rlramos@fing.edu.uy
Abstract

The Wright-Fisher (W-F) diffusion model serves as a foundational framework for interpreting population evolution through allele frequency dynamics over time. Despite the known transition probability between consecutive generations, an exact analytical expression for the transition density at arbitrary time intervals remains elusive. Commonly utilized distributions such as Gaussian or Beta inadequately address the fixation issue at extreme allele frequencies (0 or 1), particularly for short periods. In this study, we introduce two alternative parametric functions, namely the Asymptotic Expansion (AE) and the Gaussian approximation (GaussA), derived through probabilistic methodologies, aiming to better approximate this density.

The AE function provides a suitable density for allele frequency distributions, encompassing extreme values within the interval [0,1]. Additionally, we outline the range of validity for the GaussA approximation. While our primary focus is on W-F diffusion, we demonstrate how our findings extend to other diffusion models featuring singularities. Through simulations of allele frequencies under a W-F process and employing a recently developed adaptive density estimation method, we conduct a comparative analysis to assess the fit of the proposed densities against the Beta and Gaussian distributions.
Key words: diffusions with singularities, Wright-Fisher model, transition density.

Classification:

1 Introduction

The Wright-Fisher (W-F) model, initially proposed by Wright1931, offers a robust mathematical framework for studying the dynamics of allele frequencies within a population over successive generations.

To begin, we introduce the discrete W-F model alongside its diffusion approximation. The discrete W-F model, as delineated by Wright1931, operates within a population of constant size NN and two allelic types {a,A}\{a,A\}. It is represented as a Markov chain XnX_{n} with state space {0,1,2,,N}\{0,1,2,\ldots,N\}, where XnX_{n} denotes the count of allele aa at generation nn. Transition probabilities within this model are defined by a conditional Bernoulli law:

(Xn+1=j|Xn=i)=(Nj)(iN)j(1iN)Nj,j=0,,N.{\mathbb{P}}(X_{n+1}=j\,|X_{n}=i)=\binom{N}{j}\left(\frac{i}{N}\right)^{j}\left(1-\frac{i}{N}\right)^{N-j},\quad j=0,\ldots,N.

Scaling the number of generations relative to the population size permits approximate the model with a diffusion process. Let XnNX_{n}^{N} represent the Markov chain, emphasizing the population size, with 1NXN0x0\frac{1}{N}X^{N}_{0}\to x_{0} as NN\to\infty. Introducing the continuous-time process XN(t)=1NXN[Nt],t0,X_{N}(t)=\frac{1}{N}X^{N}_{[Nt]},t\geq 0, yields XN()X()X_{N}(\cdot)\to X(\cdot) as NN tends to infinity, being this convergence in distribution, as demonstrated in the book Ethier1986. The limiting process is a Markov diffusion with state space [0,1][0,1], governed by the following Stochastic Differential Equation (SDE):

dX(t)\displaystyle dX(t) =\displaystyle= X(t)(1X(t))dW(t)\displaystyle\sqrt{X(t)(1-X(t))}dW(t) (1.0.1)
X(0)\displaystyle X(0) =\displaystyle= x0\displaystyle x_{0}

Here, W(t)W(t) denotes a standard Wiener process.

The infinitesimal generator of this diffusion process is defined as:

(f)(x)=12(x(1x))d2fdx2(x),\mathcal{L}(f)(x)=\frac{1}{2}(x(1-x))\frac{d^{2}f}{dx^{2}}(x),

where ff is a function twice continuously differentiable and defined on (0,1)(0,1) with appropriate boundary conditions.

Let us introduce the main focus of our work. Since no simple analytical expression is available for the distribution of allele frequencies in a given generation, the application of the discrete W-F model in inference has relied on approximate results. An overview of these results can be found in Tataru2016, which compares several approximations of the W-F distribution used in statistical inference under different assumptions.

The Gaussian distribution is commonly employed ( Bonhomme2010; Coop2010; Gautier2010; Pickrell2012; Lacerda2014; Terhorst2015), albeit with limitations due to the discrete nature of the model and the appearance of fixation or loss events at the extremes of the frequency spectrum. Most of the authors cite Nicholson2002 to justify the use of this distribution, but no proof is given in that article. This approximation is not always appropriate due to (i) the discrete nature of the model, (ii) the range of frequencies, which is restricted to the interval [0,1][0,1], and (iii) the appearance of atoms at 0 and 1 in the allele frequency, corresponding to fixation and loss of an allele, respectively. As the number of generations increases, given a constant population size, the atoms have an increasing probability of remaining fixed at 0 or 1. One possible solution, explored in the work of Tataru2015, is to use a Beta distribution, Beta with spikes at the extremes of the interval [0,1][0,1], or a truncated Gaussian distribution. The most promising proposal in terms of goodness of fit of the allele frequency distribution is to consider the diffusion model (1.0.1), but since its transition density does not have a closed analytical expression, it is not possible to perform statistical inference. These drawbacks make finding analytical expressions and/or valid ranges of the Gaussian distribution very important for statistical inference.

In this paper we propose a more general framework that allows obtaining a finite definite density function on the interval [0,1][0,1]. Our methodology, inspired by seminal works in the field, provides an asymptotic expansion of the probability density function solution for the W-F diffusion model, particularly effective when n2N\frac{n}{2N} is small. To find such a density, we generalize to a singular diffusion the technique developed in the work of DacunhaCastelle1986 and Voronka1975. Other works that complement the scope of the cited work of Voronka and Keller are those of antonelli1978, tier1978, and tier1981, where the Ohta–Kimura model with two-locus di-allelic systems with linkage are studied, for example.

As well as Voronka1975, our method provides an asymptotic expansion of the solution of the probability density function of the forward Kolmogorov equation for the W-F diffusion model that is valid when t=n2Nt=\frac{n}{2N} is small (where nn is the number of generations and NN is the constant population size). It is important to point out that our method is probabilistic in essence, instead of asymptotic methods in partial differential equations as is the case of the aforementioned work.

While our primary focus remains on the W-F diffusion, we present a generalized solution applicable to diffusions with singularities parameterized by two parameters, aa and bb. When these parameters are equal to 12\frac{1}{2}, we obtain W-F diffusion. Thus, we show how to develop a general procedure to find the transition density of a diffusion that is state dependent and singular.

To evaluate the efficacy of our proposed density and compare it with existing methods, we conduct simulations of discrete W-F trajectories. For a fixed time t,t[0,1]t,\,t\in[0,1], we use an adaptive estimation method to obtain a continuous empirical density from the simulated data (see, e.g., Bertin2011). Having a continuous density from the simulated data allows us to better evaluate the fit of different continuous densities to the data using the Hellinger distance and the 𝕃2\mathbb{L}^{2} norm.

Our work is organized into various sections. In Section 2, we introduce our notation, the underlying model that we are investigating, the process of transforming state dependent diffusion, significant findings in discovering transition densities of diffusion processes, and the lemma that provides the transition density. We provide a comprehensive proof of this lemma in Section 3. For Section 4, we present significant outcomes associated with the W-F model, as well as corresponding proofs and special scenarios, such as neutral, with mutation and selection. Proofs for these scenarios can be found in Section LABEL:apend. Additionally, in Section LABEL:section:_model_evaluation, we assess the suitability of the defined models and two other models commonly used: the Beta and Gaussian distribution.

2 Preliminaries

In this section we begin by defining a stochastic differential equation (SDE) that has a drift component, μ()\mu(\cdot), and a local variance, σ()\sigma(\cdot), depending on X(t)X(t). We then define the adaptation of Lemma 2, presented by DacunhaCastelle1986 paper, so that the transition density can be studied. In this lemma, both the drift term and the local variance term require certain regularity of these functions. Using the Lamperti transform (see, for example, moller2010), we study a transformed process and the conditions for the drift term of this process to fulfill the conditions of the proposed lemma.

We will find in subsection 3.3 (proposition 6) an exact formula for the transition density of the following diffusion model

dX(t)=X(t)a(1X(t))bdW(t);a,b,X(t)(0,1),dX(t)=X(t)^{a}(1-X(t))^{b}dW(t);\qquad a,b,X(t)\in(0,1), (2.0.1)

where WW is a standard Wiener process. The W-F diffusion is a particular case of this model, when a=b=12a=b=\frac{1}{2}.

As local variance term σ(x)=xa(1x)b\sigma(x)=x^{a}(1-x)^{b} of this SDE is not constant, so it is necessary to transform it.

2.1 Local variance transformation for a SDE

Through the Lamperti transformation we can convert a SDE, with drift and local variance term dependent on tt, into a SDE with a constant local variance term.

Let us consider the following stochastic differential equation,

dX(t)=b(X(t))dt+σ(X(t))dW(t),dX(t)=b(X(t))dt+\sigma(X(t))dW(t), (2.1.1)

where its drift b()b(\cdot) component is s smooth function, and the local varainace σ()>0\sigma(\cdot)>0 is twice differentiable in [0,1][0,1], the continuity of this function implies that it is lower bounded by a constant c>0c>0, and WW ia a standard Wiener process.

Let Y(t)=F(X(t))Y(t)=F(X(t)) be a transformation of X(t)X(t), where F()F(\cdot) is a twice differentiable function; Itô lemma yields

dY(t)=[F(X(t))b(X(t))+12F(X(t))σ2(X(t))]dt+F(X(t))σ(X(t))dW(t).dY(t)=\left[F^{\prime}(X(t))b(X(t))+\frac{1}{2}F^{\prime\prime}(X(t))\sigma^{2}(X(t))\right]dt+F^{\prime}(X(t))\sigma(X(t))\,dW(t).

Defining F()F(\cdot) as

F(x)=0x1σ(u)du, we have,\displaystyle F(x)=\int_{0}^{x}\frac{1}{\sigma(u)}du,\quad\text{ we have,}
F(x)=1σ(x),and,\displaystyle F^{\prime}(x)=\frac{1}{\sigma(x)},\text{and,} (2.1.2)
F(x)=σ(x)σ2(x).\displaystyle F^{\prime\prime}(x)=-\frac{\sigma^{\prime}(x)}{\sigma^{2}(x)}.

Then the result for the transformed process is

dY(t)=[F(X(t))b(X(t))+12F(X(t))σ2(X(t))]dt+dW(t).dY(t)=\left[F^{\prime}(X(t))b(X(t))+\frac{1}{2}F^{\prime\prime}(X(t))\sigma^{2}(X(t))\right]dt+dW(t). (2.1.3)

The process Y(t)Y(t) is now the solution of a SDE with a constant local variance term equal to one. With the additional hypothesis that F()F(\cdot) has an inverse, i.e x=F1(y)x=F^{-1}(y), we can re-write the equation (2.1.3) as

dY(t)=[F(F1(Y(t)))b(F1(Y(t)))+12F(F1(Y(t)))σ2(F1(Y(t)))]dt+dW(t).dY(t)=\left[F^{\prime}(F^{-1}(Y(t)))b(F^{-1}(Y(t)))\right.\\ \left.+\frac{1}{2}F^{\prime\prime}(F^{-1}(Y(t)))\sigma^{2}(F^{-1}(Y(t)))\right]dt+dW(t). (2.1.4)

Now, we have a SDE, resulting from the Lamperti transformation, with constant local variance equal to one and drift term equal to

μ(y)=F(F1(y))b(F1(y))+12F(F1(y))σ2(F1(y)).\mu(y)=F^{\prime}(F^{-1}(y))b(F^{-1}(y))+\frac{1}{2}F^{\prime\prime}(F^{-1}(y))\sigma^{2}(F^{-1}(y)).

In the particular case when b=0b=0 we obtain

dY(t)=12σ(F1(Y(t)))dt+dW(t).dY(t)=-\dfrac{1}{2}\sigma^{\prime}(F^{-1}(Y(t)))dt+dW(t). (2.1.5)

2.2 Transition density for a diffusion with constant local variance term

The above transformation allows us to work, without loss of generality, with a diffusion that has a constant term of local variance. Now based on Lemma 2, of DacunhaCastelle1986, we present the following lemma that allows obtaining a exact formula for the transition density.

Lemma 1.

[Transition density] Let us consider the following stochastic differential equation

dX(t)=μ(X(t))dt+dW(t)dX(t)=\mu(X(t))dt+dW(t) (2.2.1)

where W(t)W(t) is the standard Wiener process and μ()\mu(\cdot) is a smooth drift function and verifies

μ2(x)+μ(x)=O(|x|2) when |x|.\mu^{2}(x)+\mu^{\prime}(x)=O(|x|^{2})\mbox{ when }|x|\to\infty. (2.2.2)

Let ptp_{t} be the transition density of the process defined in (2.2.1), and {B(t)}t0\{B(t)\}_{t\geq 0} be a standard Brownian bridge and let us denote M(t)=0tμ(s)dsM(t)=\int_{0}^{t}\mu(s)\,ds. Then, ptp_{t} can be written as

pt(x,y)=12πtexp[(xy)22t+(M(y)M(x))]×𝔼exp[t201ν(tB(s)+(1s)x+sy)ds],p_{t}(x,y)=\frac{1}{\sqrt{2\pi t}}\exp\left[-\frac{(x-y)^{2}}{2t}+(M(y)-M(x))\right]\\ \times\,{\mathbb{E}}\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\right], (2.2.3)

where ν(x)=(μ2(x)+μ(x))\nu(x)=(\mu^{2}(x)+\mu^{\prime}(x)).

A detailed proof of Lemma 1 can be found in the following section.

Furthermore, using the above lemma we can obtain an useful approximation for the transition density when t0t\to 0.

Lemma 2.

For t0t\to 0 and yxy\neq x, we have

limt01t𝔼exp[t201ν(tB(s)+(1s)x+sy)ds]=12(yx)xyν(u)du\lim_{t\to 0}\frac{1}{t}{\mathbb{E}}\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\right]=-\frac{1}{2(y-x)}\int_{x}^{y}\nu(u)\,du
Proof.

Using the approximation exp(f(t))1f(t)\exp(f(t))-1\approx f(t) when f(t)0f(t)\to 0, we get that

limt01t[exp(t201ν(tB(s)+(1s)x+sy)ds)1]limt01201ν(tB(s)+(1s)x+sy)ds=1201ν((1s)x+sy)ds=12(yx)xyν(u)du.\lim_{t\to 0}\frac{1}{t}\left[\exp\left(-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\right)-1\right]\\ \approx\lim_{t\to 0}-\frac{1}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\\ =-\frac{1}{2}\int_{0}^{1}\nu((1-s)x+sy)\,ds=-\frac{1}{2(y-x)}\int_{x}^{y}\nu(u)\,du.

Then

1t𝔼{exp[t201ν(tB(s)+(1s)x+sy)ds]1}t012(yx)xyν(u)du\frac{1}{t}{\mathbb{E}}\left\{\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\right]-1\right\}\\ \underset{t\to 0}{\longrightarrow}-\frac{1}{2(y-x)}\int_{x}^{y}\nu(u)\,du (2.2.4)

Corollary 3.

[Asymptotic expansion] The transition density, ptp_{t}, can be written as follows

pt(x,y)=12πtexp[(xy)22t+(M(y)M(x))]×(1t2(yx)xyν(u)du+o(t)).p_{t}(x,y)=\frac{1}{\sqrt{2\pi t}}\exp\left[-\frac{(x-y)^{2}}{2t}+(M(y)-M(x))\right]\\ \times\big{(}1-\frac{t}{2(y-x)}\int_{x}^{y}\nu(u)\,du+o(t)\big{)}. (2.2.5)

Above, we have used the Landau’s symbols o()o(\cdot) and O()O(\cdot), we say the d(t)=o(t)d(t)=o(t) or O(t))O(t)), if limt0d(t)t=0\lim_{t\to 0}\frac{d(t)}{t}=0 or a constant, respectively.

The equation defined in (2.2.5) gives us an approximation to the transition density that is valid when tt approaches 0. For this reason, we will refer to this approximation as the Asymptotic Expansion (AE).

Given that we are going to work with a transformed process Y(t)Y(t), it is necessary to establish a link between the results for the original process, X(t)X(t) and the transformed one. The following lemma is presented in order to state the relation between their two transition densities.

Lemma 4.

Let X(t)X(t) be a real-valued Markov process with transition density ptXp_{t}^{X} and FF an increasing and differentiable function. Let us define the process Y(t)=F(X(t))Y(t)=F(X(t)) with transition density ptYp_{t}^{Y}. Then, for any real numbers x0<xx_{0}<x, the transition density of X(t)X(t) can be written as

ptX(x0,x)=ptY(F(x0),F(x))F(x).p_{t}^{X}(x_{0},x)=p_{t}^{Y}(F(x_{0}),F(x))F^{\prime}(x). (2.2.6)
Proof.

Let c1,c2c_{1},c_{2} be real numbers such that c1<c2c_{1}<c_{2}. Then, {IEEEeqnarray*}rCl P^x_0(X(t)∈[c_1, c_2])=P^y_0(Y(t)∈[F(c_1),F(c_2)])&=∫_F(c_1)^F(c_2)p^Y_t(y_0,y)dy
=∫_c_1^c_2 p^Y_t(F(x_0),F(x))F’(x)dx, where y0=F(x0)y_{0}=F(x_{0}) and y=F(x)y=F(x). By the definition of transition density, we conclude that (2.2.6) holds. ∎

This lemma implies that if we have an exact or approximate analytical expression for the transition density of Y(t)Y(t), we also have one for X(t)X(t).

3 Proof for transition density lemma 1

To study the transition density defined in (2.2.3), we further analyze the terms involved in it. Starting with the associated semigroup and how it is related to the drift function, continuing with the associated drift function and finally and analysis of the function σ()\sigma(\cdot) considered.

For ease of notation let us write the transition density of Y(t)Y(t) as

ptY(y0,y)=qt(y0,y)exp(M(y0)M(y))×𝔼exp[t201ν(tB(s)+(1s)x+sy)ds],p_{t}^{Y}(y_{0},y)=q_{t}(y_{0},y)\exp\left(M(y_{0})-M(y)\right)\\ \times\,{\mathbb{E}}\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\right], (3.0.1)

where

qt(y0,y)=12πtexp[(y0y)22t] and M(y)=0y12σ(F1(u))dudy.q_{t}(y_{0},y)=\frac{1}{\sqrt{2\pi t}}\exp\left[-\frac{(y_{0}-y)^{2}}{2t}\right]\text{ and }M(y)=\int_{0}^{y}-\dfrac{1}{2}\sigma^{\prime}(F^{-1}(u))du\,dy.

By equation (2.2.6), we obtain

ptX(x0,x)=qt(F(x0)),F(x))exp[M(F(x0))M(F(x))]F(x)×𝔼F(x0)exp[t201ν(tB(s)+(1s)x+sy)ds].p_{t}^{X}(x_{0},x)=q_{t}(F(x_{0})),F(x))\exp\left[M(F(x_{0}))-M(F(x))\right]F^{\prime}(x)\\ \times\,{\mathbb{E}}^{F(x_{0})}\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)x+sy)ds\right]. (3.0.2)

3.1 Semigroup associated

The semigroup associated to the Markov process XX defined in (2.2.1) is given by

Ptf(x)=𝔼x[f(X(t))]=Ipt(x,y)f(y)dy,P_{t}f(x)={\mathbb{E}}^{x}[f(X(t))]=\int_{I}p_{t}(x,y)f(y)dy, (3.1.1)

where f()f(\cdot) is a smooth function with compact support I[0,1]I\subset[0,1] and 𝔼x{\mathbb{E}}^{x} is the expectation conditioning on X=xX=x. The semigroup in (3.1.1) can be written using the Cameron-Martin-Girsanov formula (see Oksendal2000 pp. 123) (this formula holds under the condition μ2(x)=O(|x|2)\mu^{2}(x)=O(|x|^{2}) when |x||x|\to\infty) as {IEEEeqnarray}rCl P_tf(x) & = E^x[e^∫_0^t μ(W(s))dW(s)-12∫_0^t μ^2(W(s))dsf(W(t))]

Which, when it is applied to (2.1.5), gives

𝔼y0[f(Y(t))]=𝔼y0[e120tσ(F1(W(s)))dW(s)140t(σ(F1(W(s)))2dsf(W(t))].{\mathbb{E}}^{y_{0}}[f(Y(t))]={\mathbb{E}}^{y_{0}}[e^{-\frac{1}{2}\int_{0}^{t}\sigma^{\prime}(F^{-1}(W(s)))dW(s)-\frac{1}{4}\int_{0}^{t}(\sigma^{\prime}(F^{-1}(W(s)))^{2}ds}f(W(t))]. (3.1.2)

Now, we analyze the drift function involved.

3.2 Drift function

In our case, the drift function is defined by μ(z)=12σ(F1(z))\mu(z)=-\frac{1}{2}\sigma^{\prime}(F^{-1}(z)). Considering that M(z)=0zμ(s)dsM(z)=\int_{0}^{z}\mu(s)ds and F()F(\cdot) is an increasing function, we have

M(z)\displaystyle M(z) =12F1(0)F1(z)σ(x)σ(x)dx=12(logσ(F1(z))logσ(F1(0))),\displaystyle=-\dfrac{1}{2}\int_{F^{-1}(0)}^{{F^{-1}(z)}}\dfrac{\sigma^{\prime}(x)}{\sigma(x)}dx=-\dfrac{1}{2}\left(\log\sigma(F^{-1}(z))-\log\sigma(F^{-1}(0))\right),
M(z)\displaystyle M^{\prime}(z) =12σ(F1(z)),\displaystyle=-\frac{1}{2}\sigma^{\prime}(F^{-1}(z)),
M(z)\displaystyle M^{\prime\prime}(z) =12σ(F1(z))σ(F1(z)),\displaystyle=-\frac{1}{2}\sigma^{\prime\prime}(F^{-1}(z))\sigma(F^{-1}(z)),

Considering these results and by means of Itô lemma,

M(W(t))M(W(0))=120tσ(F1(W(s)))dW(s)120tσ(F1(W(s)))σ(F1(W(s)))ds\displaystyle M(W(t))-M(W(0))=-\frac{1}{2}\int_{0}^{t}\sigma^{\prime}(F^{-1}(W(s)))dW(s)-\frac{1}{2}\int_{0}^{t}\sigma^{\prime\prime}(F^{-1}(W(s)))\sigma(F^{-1}(W(s)))ds

Implying that,

M(W(t))M(W(0))+120tσ(F1(W(s)))σ(F1(W(s)))ds=120tσ(F1(W(s)))dW(s)\displaystyle M(W(t))-M(W(0))+\frac{1}{2}\int_{0}^{t}\sigma^{\prime\prime}(F^{-1}(W(s)))\sigma(F^{-1}(W(s)))ds=-\frac{1}{2}\int_{0}^{t}\sigma^{\prime}(F^{-1}(W(s)))dW(s)

Replacing this result in (3.1.2), it results a formula without terms involving stochastic integrals

𝔼y0[f(Y(t))]=eM(W(0))𝔼y0[e120tσ(F1(W(s)))σ(F1(W(s)))ds140t(σ(F1(W(t)))2dseM(W(t))f(W(t))].{\mathbb{E}}^{y_{0}}[f(Y(t))]=e^{-M(W(0))}{\mathbb{E}}^{y_{0}}[e^{\frac{1}{2}\int_{0}^{t}\sigma^{\prime\prime}(F^{-1}(W(s)))\sigma(F^{-1}(W(s)))ds-\frac{1}{4}\int_{0}^{t}(\sigma^{\prime}(F^{-1}(W(t)))^{2}ds}e^{M(W(t))}f(W(t))]. (3.2.1)
Remark 5.

It is important to point out that all the results obtained above can be obtained without difficulty from the paper of DacunhaCastelle1986.

3.3 Analysis of σ()\sigma(\cdot)

This subsection aims to extend the precedent results to more general local variance terms. First let us assume that σ:[0,1]+\sigma:[0,1]\to\mathbb{R}^{+} is concave and let us define

V(x)=12|σ(F1(x))|σ(F1(x)))+14(σ(F1(x))20.V(x)=\frac{1}{2}|\sigma^{\prime\prime}(F^{-1}(x))|\sigma(F^{-1}(x)))+\frac{1}{4}(\sigma^{\prime}(F^{-1}(x))^{2}\geq 0.

The above result enables us to write,

𝔼y0[f(Y(t))]=F(0)F(1)pYt(y0,y)f(y)dy=F(0)F(1)𝔼[e120tV(y0+W(s))ds|W(t)=yy0]qt(y0y)eM(y)M(y0)f(y)dy=F(0)F(1)𝔼[e120tV(y0+(W(s)stW(t))+stW(t))ds|W(t)=yy0]qt(y0y)eM(y)M(y0)f(y)dy=F(0)F(1)𝔼[e120tV((1st)y0+sty+(W(s)stW(t))ds]qt(y0y)eM(y)M(y0)f(y)dy=F(0)F(1)𝔼[et201V((1u)y0+uy+t(W(u)uW(1))du]qt(y0y)eM(y)M(y0)f(y)dy.{\mathbb{E}}^{y_{0}}[f(Y(t))]=\int_{F(0)}^{F(1)}p^{Y}_{t}(y_{0},y)f(y)dy\\ =\int_{F(0)}^{F(1)}{\mathbb{E}}[e^{-\frac{1}{2}\int_{0}^{t}V(y_{0}+W(s))ds}\,|W(t)=y-y_{0}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})}f(y)dy\\ =\int_{F(0)}^{F(1)}{\mathbb{E}}[e^{-\frac{1}{2}\int_{0}^{t}V(y_{0}+(W(s)-\frac{s}{t}W(t))+\frac{s}{t}W(t))ds}\,|W(t)=y-y_{0}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})}f(y)dy\\ =\int_{F(0)}^{F(1)}{\mathbb{E}}[e^{-\frac{1}{2}\int_{0}^{t}V((1-\frac{s}{t})y_{0}+\frac{s}{t}y+(W(s)-\frac{s}{t}W(t))ds}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})}f(y)dy\\ =\int_{F(0)}^{F(1)}{\mathbb{E}}[e^{-\frac{t}{2}\int_{0}^{1}V((1-u)y_{0}+uy+\sqrt{t}(W(u)-uW(1))du}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})}f(y)dy.

By duality, it results that for F(0)y0,yF(1)F(0)\leq y_{0},y\leq F(1), yields,

pYt(y0,y)=𝔼[et201V((1u)y0+uy+t(W(u)uW(1))du]qt(y0y)eM(y)M(y0).p^{Y}_{t}(y_{0},y)={\mathbb{E}}[e^{-\frac{t}{2}\int_{0}^{1}V((1-u)y_{0}+uy+\sqrt{t}(W(u)-uW(1))du}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})}.

In addition, for y0=F(x0)y_{0}=F(x_{0})

x0(X(t)[a,b])=y0(F(X(t))[F(a),F(b)])=F(a)F(b)pYt(y0,y)dy=abpYt(F(x0),F(x))1σ(x)dx.{\mathbb{P}}^{x_{0}}(X(t)\in[a,b])={\mathbb{P}}^{y_{0}}(F(X(t))\in[F(a),F(b)])=\int_{F(a)}^{F(b)}p^{Y}_{t}(y_{0},y)dy=\int_{a}^{b}p^{Y}_{t}(F(x_{0}),F(x))\frac{1}{\sigma(x)}dx.

This relationship means that

pXt(x0,x)=pYt(F(x0),F(x))1σ(x).p^{X}_{t}(x_{0},x)=p^{Y}_{t}(F(x_{0}),F(x))\frac{1}{\sigma(x)}. (3.3.1)

Below we extend the class of functions σ()\sigma(\cdot) for which the above result holds, based on the following assumptions

  • (A1)

    Let us consider that σ:[0,1]+\sigma:[0,1]\to{\mathbb{R}}^{+}, is concave, and it tends to zero at 0 and 11, such that 011σ(u)du<,\int_{0}^{1}\frac{1}{\sigma(u)}du<\infty,.

  • (A2)

    We can approximate σ\sigma and its two derivatives by a sequence of functions with two continuous derivatives such that σn0\sigma^{\prime\prime}_{n}\leq 0, σnσ\sigma_{n}\to\sigma, as well as its derivatives

These assumptions allow us to approximate the expression (3.2.1) by a similar expression written for the sequence instead of σ\sigma; then, by the dominated convergence theorem,we can exchange limit with expectation.

To clarify this point, let VnV_{n} be defined as VV, but replacing σn()\sigma_{n}(\cdot) instead of σ()\sigma(\cdot), and let us denote pYntp^{Y_{n}}_{t} the transition density. The Kolmogorov differential equation, implies that pYnt(y0,y)pYt(y0,y)p^{Y_{n}}_{t}(y_{0},y)\to p^{Y}_{t}(y_{0},y).

Additionally, we have

𝔼[et201Vn((1u)y0+uy+tBn(u)du]qt(y0y)eMn(y)Mn(y0)=𝔼[et201Vn((1u)y0+uy+tB(u)du]qt(y0y)eMn(y)Mn(y0).{\mathbb{E}}[e^{-\frac{t}{2}\int_{0}^{1}V_{n}((1-u)y_{0}+uy+\sqrt{t}B_{n}(u)du}]q_{t}(y_{0}-y)e^{M_{n}(y)-M_{n}(y_{0})}\\ ={\mathbb{E}}[e^{-\frac{t}{2}\int_{0}^{1}V_{n}((1-u)y_{0}+uy+\sqrt{t}B(u)du}]q_{t}(y_{0}-y)e^{M_{n}(y)-M_{n}(y_{0})}. (3.3.2)

where, by equality in distribution, we have substituted, the sequence of Brownian bridges, Bn(u)B_{n}(u), by a fixed Brownian bridge B(u)B(u).

Considering that VnV_{n} is a positive function, we use dominated convergence theorem to show that the term (3.3.2) converges to

𝔼[et201V((1u)y0+uy+tB(u)du]qt(y0y)eM(y)M(y0).{\mathbb{E}}[e^{-\frac{t}{2}\int_{0}^{1}V((1-u)y_{0}+uy+\sqrt{t}B(u)du}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})}.

Thus, we obtain

pYt(y0,y)=E[et201V((1u)y0+uy+tB(u)du]qt(y0y)eM(y)M(y0),p^{Y}_{t}(y_{0},y)=E[e^{-\frac{t}{2}\int_{0}^{1}V((1-u)y_{0}+uy+\sqrt{t}B(u)du}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})},

and also the relationship (3.3.1), holds.

These computations allow us to give the following result

Proposition 6.

Based on assumptions (A1) and (A2), we get

pYt(y0,y)=E[et201V((1u)y0+uy+tB(u)du]qt(y0y)eM(y)M(y0),p^{Y}_{t}(y_{0},y)=E[e^{-\frac{t}{2}\int_{0}^{1}V((1-u)y_{0}+uy+\sqrt{t}B(u)du}]q_{t}(y_{0}-y)e^{M(y)-M(y_{0})},
Remark 7.

We think this is a new result; we have not found a similar one in the literature.. The class of functions defined by σa,b(x)=xa(1x)b\sigma_{a,b}(x)=x^{a}(1-x)^{b}, for 0<a,b<10<a,b<1 satisfies the conditions (A1) and (A2). When a=b=12a=b=\frac{1}{2} it corresponds to the W-F diffusion.

4 Wright - Fisher diffusion: Asymptotic Expansion and Gaussian approximation

In the following we concentrate in the main subject of this paper: the W-F diffusion ( a=b=1/2a=b=1/2).

A parallel development can be carried out for other values of the parameters, although we will not consider this problem in the present work.

4.1 Analytical expression for the transition density

The goal of this section will be to prove an exact analytical expression for the transition density and its asymptotic expansion for small tt, of the W-F diffusion. When this diffusion is considered, it produces the following

dX(t)=X(t)(1X(t))dW(t)\displaystyle dX(t)=\sqrt{X(t)(1-X(t))}dW(t) (4.1.1)
dY(t)=12cot(Y(t))dt+dW(t)\displaystyle dY(t)=-\dfrac{1}{2}\cot(Y(t))dt+dW(t) (4.1.2)

Let us recall that in this case the following function writes ν(y)=μ2(F1(y))+μ(F1(y))=14cot2(y)+12csc2(y)=12+34cot2(y)\nu(y)=\mu^{2}(F^{-1}(y))+\mu^{\prime}(F^{-1}(y))=\frac{1}{4}\cot^{2}(y)+\frac{1}{2}\csc^{2}(y)=\frac{1}{2}+\frac{3}{4}\cot^{2}(y). The main result is

Theorem 8.

The transition density of the W-F diffusion can be written as

(i)pXt(x0,x)=(x0(1x0))14(x(1x))34×𝔼[et201ν((1u)F(x0)+uF(x)+tB(u))du]qt(F(x),F(x0)),(i)\qquad p^{X}_{t}(x_{0},x)=\frac{(x_{0}(1-x_{0}))^{\frac{1}{4}}}{(x(1-x))^{\frac{3}{4}}}\\ \times{\mathbb{E}}\left[e^{-\frac{t}{2}\int_{0}^{1}\nu((1-u)F(x_{0})+uF(x)+\sqrt{t}B(u))du}\right]q_{t}(F(x),F(x_{0})), (4.1.3)

and its asymptotical expansion as

(ii)pXt(x0,x)=(x0(xx0))14(x(1x))34exp[12t(F(x)F(x0))2]×[1t21F(x)F(x0)F(x0)F(x)ν(u)du+o(t)](ii)\qquad p^{X}_{t}(x_{0},x)=\frac{(x_{0}(x-x_{0}))^{\frac{1}{4}}}{(x(1-x))^{\frac{3}{4}}}\exp\left[-\frac{1}{2t}(F(x)-F(x_{0}))^{2}\right]\\ \times\left[1-\frac{t}{2}\,\frac{1}{F(x)-F(x_{0})}\int_{F(x_{0})}^{F(x)}\nu(u)du+o(t)\right] (4.1.4)

for all x0,x(0,1)x_{0},x\in(0,1).

Proof.

Let us start by proving (i)(i). For ease of notation let us write the transition density of Y(t)Y(t), based on Proposition 6, as

ptY(y0,y)=qt(y0,y)exp(M(y)M(y0))×𝔼exp[t201ν(tB(s)+(1s)y0+sy)ds],p_{t}^{Y}(y_{0},y)=q_{t}(y_{0},y)\exp\left(M(y)-M(y_{0})\right)\\ \times\,{\mathbb{E}}\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)y_{0}+sy)ds\right], (4.1.5)

where

qt(y0,y)=12πtexp[(y0y)22t] and M(y)=0y12cot(y)dy.q_{t}(y_{0},y)=\frac{1}{\sqrt{2\pi t}}\exp\left[-\frac{(y_{0}-y)^{2}}{2t}\right]\text{ and }M(y)=\int_{0}^{y}-\frac{1}{2}\cot(y)\,dy.

By equation (2.2.6), we obtain

ptX(x0,x)=qt(F(x0)),F(x))exp[M(F(x0))M(F(x))]1x(1x)×𝔼exp[t201ν(tB(s)+(1s)F(x0)+sF(x))ds].p_{t}^{X}(x_{0},x)=q_{t}(F(x_{0})),F(x))\exp\left[M(F(x_{0}))-M(F(x))\right]\frac{1}{\sqrt{x(1-x)}}\\ \times\,{\mathbb{E}}\exp\left[-\frac{t}{2}\int_{0}^{1}\nu(\sqrt{t}B(s)+(1-s)F(x_{0})+sF(x))ds\right]. (4.1.6)

Let us further develop the terms in equation (4.1.5). For every pair of positive real numbers yy and y0y_{0} such that y>y0y>y_{0}, we have that {IEEEeqnarray*}rCl M(y)-M(y_0) &= -12∫_y_0^ycot(u) du
= 12 (log—siny—- log—siny_0— )
= 14 [logsin^2(y_0)-logsin^2(y) ]
= 14[log