Weak approximations of Wright–Fisher equation

Gabriele Mongirdaite; Vigirdas Mackevičius

Articles

Esta obra está bajo una Licencia Creative Commons Atribución 4.0 Internacional.

Recepción: 10 Julio 2021

Publicación: 15 Diciembre 2021

DOI: https://doi.org/10.15388/LMR.2021.25217

Abstract: We construct weak approximations of the Wright-Fisher model and illustrate their accuracy by simulation examples.

Keywords: Wright–Fisher model, simulation, weak approximation .

Summary: Sukonstruota silpnoji pirmos eile˙s aproksimacija stochastinei Wright–Fisher lygčiai. Pavyzdžiais iliustruojamas jos tikslumas.

Keywords: Wright–Fisher modelis, modeliavimas, silpnoji aproksimacija.

Introduction

We consider Wright–Fisher process defined by the stochastic differential equation

d X_{t}^{x} = (a - b X_{t}^{x}) d t + σ \sqrt{X_{t}^{x} (1 - X_{t}^{x}) b B_{t},} X_{0}^{x} = x,

(1)

where $B$ is a standard Brownian motion, $0 ⩽ a ⩽ b, σ > 0,$ and $x \in [0, 1] .$

The Wright–Fisher model (Fisher 1930; Wright 1931) takes the values in the interval $[0, 1]$ and explicitly accounts for the effects of various evolutionary forces – random genetic drift, mutation, selection – on allele frequencies over time. This model can also accommodate the effect of demographic forces such as variation in population size through time and/or migration connecting populations [5].

In this note, we present a simple first-order weak approximation of the solution of Eq. (1) by discrete random variables that take two values at each approximation step.

Recall the definition of such an approximation. By a discretization scheme with time step $h > 0$ we mean any time-homogeneous Markov chain $\hat{X^{h}} = {\hat{X_{k h}^{h}}, k = 0, 1, . . .} .$

We say that a family of discretization schemes $\hat{X^{h}}, h > 0,$ is a first-order weak approximation of the solution $X^{x}$ of (1) in the interval $[0, T]$ if

| E f (\hat{X_{T}^{h}}) - E f (X_{T}^{x}) | ⩽ Ch, h = \frac{T}{N} ⩽ h o,

(2)

for a “sufficiently wide” class of functions $f : [0, 1] \to R$ and some constants $C$ and $h_{0} > 0$ (depending on the function $f$ ), where $N ∈ N .$ Note that because of the Markovity, the one-step approximation $\hat{X_{h}^{h}}$ completely defines (in distribution) a weak approximation $\hat{X_{k h}^{h}}, k = 0, 1, . . . .$ Thus, with some ambiguity, we also call it an approximation and denote it by $\hat{X_{h}^{x}},$ with $x$ indicating its starting point.

In our context, we introduce the following “sufficiently wide” function class of infinitely differentiable functions with “not too fast” growing derivatives:

C_{*}^{\infty} [0, 1] : = {f \in C^{\infty} [0, 1] : \underset{k \to ∞}{lim sup} \frac{1}{k!} \sup_{x \in [0, 1]} | f^{(k)} (x) | < \infty} .

We easily see that all functions from this class can be expanded by the Taylor series in the interval $[0, 1]$ around arbitrary $x_{0} \in [0, 1]$ (which, in fact, converges on the whole real line $R$ ) and contain, for example, all polynomials and exponential functions.

Approximation

Let us first construct an approximation for the “stochastic” part of Wright–Fisher equation, that is, the solution $S_{t}^{x}$ of Eq. (1) with $a = b = 0 .$ Similarly to [4] (see also [3]), we look for an approximation $\hat{X_{h}^{x}}$ as a two-valued discrete random variable taking values $x_{1, 2} ∈ [0, 1]$ with probabilities $p_{1, 2}$ such that

E (\hat{S_{h}^{x}} - x) = 0, x \in [0, 1],

(3)

E {(\hat{S_{h}^{x}} - x)}^{2} = σ^{2} x (1 - x) h + O (h^{2}), x \in [0, 1],

(4)

|E {(\hat{S_{h}^{x}} - x)}^{3} | = O (h^{2}), x \in [0, 1],

(5)

E [{(\hat{S_{h}^{x}} - x)}^{4} = O (h^{2}), x \in [0, 1] .

(6)

By solving the equation system (3)–(4) with respect to $x_{1}, x_{2}, p_{1}, p_{2},$ we get the solution

x_{1} = x + (1 - x) σ^{2} h - \sqrt{(x + (1 - x) σ^{2} h) (1 - x) σ^{2} h,} x \in [0, 1],

(7)

x_{2} = x + (1 - x) σ^{2} h + \sqrt{(x + (1 - x) σ^{2} h) (1 - x) σ^{2} h,} x \in [0, 1]

(8)

with $p_{1, 2} = \frac{x}{2_{x_{1, 2}}} .$ It also satisfies conditions (5)–(6). However, for the values of $x$ near 1, the values of $x_{2}$ a slightly greater than 1, which is unacceptable. We overcome this problem by using the symmetry of the solution of the stochastic part with respect to the point $\frac{1}{2};$ to be precise, $S_{t}^{x} \overset{d}{=} 1 - S_{t}^{1 - x} .$ Therefore, in the interval $[0, 1 / 2],$ we can use the values $x_{1, 2}$ defined by (7)-(8), whereas in the interval $(1 / 2, 1),$ we use the values corresponding to the process $1 - \hat{S_{t}^{1 - x}}$ , that is,

\hat{x_{1, 2}} = \hat{x_{1, 2}} (x, h) : = 1 - x_{1, 2} (1 - x, h) = x - x σ^{2} h \pm \sqrt{(1 - x + x σ^{2} h) x σ^{2} h}

(9)

with probabilities

with probabilities $\hat{p_{1, 2}} = \frac{1 - x}{2_{x_{1, 2}} (1 - x, h)} .$ Thus we obtain a correct (i.e., with values in $[0, 1]$ ) approximation $\hat{S_{h}^{x}}$ taking the values

\tilde{x_{1, 2}} : = {\binom{x_{1, 2} (x, h) With probabilities p_{1, 2} = \frac{x}{2_{x_{1, 2}} (x, h)}, x \in [0, 1 / 2],}{1 - x_{1, 2} (1 - x, h) with probabilities p_{1, 2} = \frac{1 - x}{2_{x_{1, 2}} (1 - x, h)}, x \in (1 / 2, 1] .}

Now for the initial equation (1), we obtain an approximation $\hat{X_{h}^{x}}$ by a simple "splitstep" procedure (again, see, e.g., [4] or [3]):

\hat{X_{h}^{x}} : = \hat{S_{h}^{x}} e^{- b h} + \frac{a}{b} (1 - e^{- b h}) .

(10)

Now we can state the following:

Theorem 1.Let $\hat{X_{t}^{x}}$ be the discretization scheme defined by one-step approximation (10). Then $\hat{X_{t}^{x}}$ is a first-order weak approximation of equation (1) for functions $f \in C_{*}^{\infty} [0, 1] .$

Backward Kolmogorov equation

The constructed approximation is in fact a so-called potential first-order weak approximation of Eq. (1) (for a definition, see, e.g., Alfonsi [1], Section 2.3.1). The proof that, indeed, it is a first-order weak approximation, is based on the following:

Theorem 2. Let $f \in C_{*}^{∞} [0, 1] . The u (t, x) : = E f (X_{ι}^{x}) i s a C^{∞}$ function on $[0, 1] x R$ that solves the backward Kolmogorov equation

\partial_{t u} (t, x) = A u (t, x), x \in [0, 1], t ⩾ 0 .

In particular,

∀ T > 0, ∀ l, m ∈ N, \exists C_{l, m} : | \partial_{l} \partial_{m} u (t, x) | ⩽ C_{l, m,} t \in [0, T], x \in [0, 1] .

Such theorem is stated for $f \in C^{∞} [0, 1]$ in $[1, Thm . 6.1.12],$ based on the results of [2]. Our class of functions $f$ is slightly narrower, but our proof of the theorem is significantly simpler and is based on the estimates of the moments of $X_{t}^{x},$ which show that they grow slower than factorials. The recurrent relations of the moments $E [{(X_{t}^{x})}^{k}]$ show that they are infinitely differentiable with respect to $t$ and $x,$ which allows us infinitely differentiate the series

u (t, x) = E f (X_{t}^{x}) = \sum_{k = 0}^{\infty} c_{k} E [{(X_{t}^{x})}^{k}]

termwise with respect to $t$ and $x,$ where $f (x) = \sum_{k = 0}^{\infty} c_{{k x}^{k}}$ is the Taylor expansion of $f .$

Simulation examples

We illustrate our approximation for $f (x) = x^{4}$ and $f (x) = \exp {- x} .$ Since we do not explicitly know the moments $E \exp {- X_{t}^{x}},$ we use the approximate equality $\exp {- x} \approx 1 - x + \frac{x^{2}}{2} - \frac{x^{3}}{6} + \frac{x^{4}}{24} .$ In Figs. 1 and 2, we compare the moments $E f (\hat{X_{t}^{x}})$ and $E f (x_{t}^{x})$ as functions of $t$ $(left plots, h = 0.001)$ and as functions of discretization step $h$ $(right plots, t = 1) .$ As expected, the approximations agree with exact values pretty well.

References

1. A. Alfonsi. Affine Diffusions and Related Processes: Simulation, Theory and Applications. Springer, 2015.

2. C. Epstein, R. Mazzeo. Wright–Fisher diffusion in one dimension. SIAM J. Math. Anal., 42:568–608, 2010.

3. G. Lileika, V. Mackevičius. Weak approximation of CKLS and CEV processes by discrete random variables. Lith. Math. J., 20(2):208–224, 2020.

4. V. Mackevičius. Weak approximation of CIR equation by discrete random variables. Lith. Math. J., 51(3):385–401, 2011.

5. P. Tataru, M. Simonsen, T. Bataillon, A. Hobolth. Statistical inference in the Wright– Fisher model using allele frequency data. Syst. Biol., 66(1):e30–e46, 2016.