Revisiting Feller Diffusion: Derivation and Simulation

TL;DR

This paper presents a simplified derivation of the Feller Diffusion's probability density function using Fourier Transform and Characteristics, along with simulation algorithms to verify hitting time properties.

Contribution

It introduces a more straightforward derivation method and discusses simulation techniques for Feller Diffusion, enhancing understanding and computational approaches.

Findings

Derived the Feller Diffusion PDF using Fourier Transform and Characteristics.

Validated hitting time properties through simulation.

Provided practical algorithms for simulating Feller Diffusion.

Abstract

We propose a simpler derivation of the probability density function of Feller Diffusion using the Fourier Transform and solving the resulting equation via the Method of Characteristics. We also discuss simulation algorithms and confirm key properties related to hitting time probabilities via the simulation.

Tables

Table 1. Table 1: Percentage of paths absorbed at origin over time for different values of σ 2 superscript 𝜎 2 \sigma^{2} with t 0 = 0 subscript 𝑡 0 0 t_{0}=0 , Δ ​ t i = 1 Δ subscript 𝑡 𝑖 1 \Delta t_{i}=1 , x 0 = 1000 subscript 𝑥 0 1000 x_{0}=1000 . (10,000 paths)

	Time (t)
${\sigma }^2$	100		1,000		5,000		10,000		20,000
	Theo.	Sim.	Theo.	Sim.	Theo.	Sim.	Theo.	Sim.	Theo.	Sim.
0.1	0.00	0.00	0.00	0.00	1.83	1.70	13.53	13.47	36.79	36.11
1	0.00	0.00	13.53	13.31	67.03	67.08	81.87	82.10	90.48	90.68
10	13.53	13.02	81.87	81.84	96.08	96.12	98.02	98.03	99.00	99.03
100	81.87	82.10	98.02	97.97	99.60	99.61	99.80	99.76	99.90	99.88

Table 2. Table 2: Percentage of paths absorbed at origin over time for different values of x 0 subscript 𝑥 0 x_{0} with t 0 = 0 subscript 𝑡 0 0 t_{0}=0 , σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 , Δ ​ t i = 1 Δ subscript 𝑡 𝑖 1 \Delta t_{i}=1 . (10,000 paths)

	Time (t)
$x_0$	100		1,000		5,000		10,000		20,000
	Theo.	Sim.	Theo.	Sim.	Theo.	Sim.	Theo.	Sim.	Theo.	Sim.
10	81.87	82.10	98.02	97.97	99.60	99.61	99.80	99.76	99.90	99.87
100	13.53	13.02	81.87	81.85	96.08	96.11	98.02	98.04	99.00	99.03
1000	0.00	0.00	13.53	13.31	67.03	67.08	81.87	82.10	90.48	90.68
10000	0.00	0.00	0.00	0.00	1.83	1.79	13.53	13.51	36.79	36.07

Table 3. Table 3: Average x 𝑥 x value for varying x 0 subscript 𝑥 0 x_{0} values with t 0 = 0 subscript 𝑡 0 0 t_{0}=0 , σ 2 = 10 superscript 𝜎 2 10 \sigma^{2}=10 , Δ ​ t i = 1 Δ subscript 𝑡 𝑖 1 \Delta t_{i}=1 . (10,000 paths)

	Time (t)
$x_0$	100	1000	5000	10000	20000
10	8.35	6.98	2.74	0.00	0.00
100	100.70	110.73	100.92	127.18	149.70
1000	1000.47	996.56	908.82	960.91	1011.71
10000	9956.68	9958.66	9752.99	9860.61	9775.20

Table 4. Table 4: Average x 𝑥 x value for varying x 0 subscript 𝑥 0 x_{0} values with t 0 = 0 subscript 𝑡 0 0 t_{0}=0 , σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 , Δ ​ t i = 1 Δ subscript 𝑡 𝑖 1 \Delta t_{i}=1 . (10,000 paths)

	Time (t)
$x_0$	100	1000	5000	10000	20000
10	10.07	11.07	10.29	12.68	14.92
100	100.05	99.72	90.96	96.19	99.01
1000	995.67	995.34	975.22	991.77	977.22
10000	9990.73	10000.84	10004.78	10049.72	9953.11

Table 5. Table 5: 95% Confidence Intervals of x ¯ ¯ 𝑥 \overline{x} values given in Table 3

$x_0$	Time (t)	$\overline{x}$	95% Confidence Interval
			Lower Bound	Upper Bound
10	100	8.35	6.39	10.31
	1000	6.98	0.78	13.18
	5000	2.74	-11.12	16.60
	10000	0.00	-19.60	19.60
	20000	0.00	-27.72	27.72
100	100	100.70	94.50	106.90
	1000	110.73	91.13	130.33
	5000	100.92	57.09	144.75
	10000	127.18	65.20	189.16
	20000	149.70	62.05	237.35
1000	100	1000.47	980.87	1020.07
	1000	996.56	934.58	1058.54
	5000	908.82	770.23	1047.41
	10000	960.91	764.91	1156.91
	20000	1011.71	734.53	1288.89
10000	100	9956.68	9894.70	10018.66
	1000	9958.66	9762.66	10154.66
	5000	9752.99	9314.73	10191.25
	10000	9860.61	9240.81	10480.41
	20000	9775.20	8898.68	10651.72

Table 6. Table 6: 95% Confidence Intervals of x ¯ ¯ 𝑥 \overline{x} values given in Table 4

$x_0$	Time (t)	$\overline{x}$	95% Confidence Interval
			Lower Bound	Upper Bound
10	100	10.07	8.11	12.03
	1000	11.07	4.87	17.27
	5000	10.29	-3.57	24.15
	10000	12.68	-6.92	32.28
	20000	14.92	-12.80	42.64
100	100	100.05	93.85	106.25
	1000	99.72	80.12	119.32
	5000	90.96	47.13	134.79
	10000	96.19	34.21	158.17
	20000	99.01	11.36	186.66
1000	100	995.67	976.07	1015.27
	1000	995.34	933.36	1057.32
	5000	975.22	836.63	1113.81
	10000	991.77	795.77	1187.77
	20000	977.22	700.04	1254.40
10000	100	9990.73	9928.75	10052.71
	1000	10000.84	9804.84	10196.84
	5000	10004.78	9566.52	10443.04
	10000	10049.72	9429.92	10669.52
	20000	9953.11	9076.59	10829.63

Equations

$$\dot{x}(t)$$

$$x(t_o)$$

$$⟨\xi (t){⟩}_{\xi }$$

$$⟨\xi (t)\xi (t^{\prime }){⟩}_{\xi }$$

$$\frac{\partial }{\partial t}P(x,t)$$

$$P(x,t_o)$$

$${\int }_0^{\infty }P(x,t)dx$$

$$P(x,t)\Delta x=\text{Probability walker}\in [x,x+\Delta x]\text{ at time }t$$

$$\hat{P}(k,t)={\int }_{-\infty }^{\infty }P(x,t)e^{-ikx}dx={\int }_0^{\infty }P(x,t)e^{-ikx}dx$$

$$\frac{\partial }{\partial t}\hat{P}(k,t)$$

$$\hat{P}(k,t_o)$$

$$\hat{P}(k=0,t)$$

$$\hat{P}(k,t)=\exp \left[\frac{-ik{x}_o}{1+\frac{ik{\sigma }^2}{2}(t-t_o)}\right]$$

$$\phi (k,t)=\hat{P}(-k,t)$$

$$\phi (k,t)=\exp \left[\frac{ik{x}_o}{1-\frac{ik{\sigma }^2}{2}(t-t_o)}\right]$$

$$F(t)=\exp \left[\frac{-2x_o}{{\sigma }^2(t-t_o)}\right]$$

$$c$$

$$\lambda$$

$$\overline{x}(t)=\frac{1}{N}\sum _{m=1}^N{x}_m(t)$$

$$S(t)$$

$$F(t)$$

$$T^{\ast }=\{\inf (t):x(t)=0\}$$

$$f(t)=-\frac{dS(t)}{dt},$$

$$\mathrm{Probability}{T}^{\ast }\in (t,t+dt)=f(t)dt$$

$$S(t)$$

$$F(t)$$

$$f(t)=\frac{x_o}{\sqrt{2}\sigma t^{3/2}}\exp \left[\frac{-x_o^2}{2{\sigma }^{2}t}\right],$$

$$t^{\ast }=\frac{x_o^2}{3{\sigma }^2}$$

$$f(t)=\frac{2x_o}{{\sigma }^2{t}^2}\exp \left[\frac{-2x_o}{{\sigma }^{2}t}\right]$$

$$⟨T^{\ast }⟩=⟨t⟩={\int }_0^{\infty }tf(t)dt=\frac{2x_o}{{\sigma }^2}{\int }_0^{\infty }\frac{e^{-v}}{v}dv$$

$$t^{\ast }=\frac{x_o}{{\sigma }^2}$$

$$\frac{\partial }{\partial t}P(x,t)=\frac{1}{2}{\sigma }^{2}x\frac{{\partial }^2}{\partial x^2}P(x,t)+{\sigma }^2\frac{\partial }{\partial x}P(x,t)$$

$$\frac{\partial }{\partial t}\hat{P}(k,t)=\frac{1}{2}{\sigma }^2\left[i\frac{\partial }{\partial k}(-k^2\hat{P}(k,t))\right]+ik{\sigma }^2\frac{\partial }{\partial k}\hat{P}(k,t)$$

$$\frac{d^{n}G(x)}{d{x}^n}\to (ik{)}^n\hat{G}(k)$$

$$x\cdot G(x)\to i\frac{d}{dk}\hat{G}(k)$$

Full text

Keywords— Feller Diffusion, Bessel Squared Process of Dimension Zero, Zero Degrees of Freedom, Square Root Diffusion, Simulation, Method of Characteristics, Absorbing Random Walks, Survival Probability, Hitting Time, First Passage Time.

1 Introduction

Feller Diffusion, also known as Bessel Squared Process of dimension zero [2, 3, 4, 5, 9], is given by the Îto stochastic differential equation (SDE):

[TABLE]

Here $x(t)$ denotes the position of a random walker following the trajectory given by equation (1) with initial condition $x_{o}$ at (initial) time $t_{o}$, and $\xi(t)$ is a Gaussian white noise independent of $x(t)$, with amplitude parameter $\sigma$, and the following correlation functions describing its statistics:

[TABLE]

The averaging $\langle\cdots\rangle_{\xi}$ denotes an average with respect to the realizations of the stochastic process $\xi(t)$. Due to presence of the square root term the dynamics are such that the random walker’s postion is always positive and if it reaches the origin it will stay there (absorbing barrier). Feller Diffusion is in the class of Îto SDEs known as square root diffusions [6] and is a special case in that it is driftless. In general, square root diffusions have a drift term in addition to the multiplicative square root noise term. Feller’s motivation in this area was genetics [3] with natural extensions to biological systems and population dynamics. Square root diffusions also have application in physics - in particular in the area of branching processes and percolation theory [1] and phase transitions with absorbing states [7]. In addition to the science applications listed above, there are applications of square root diffusions in quantitative finance - e.g. the CEV Model for asset pricing [6] and CIR Model for interest rate modeling [8, 14].

The Fokker-Planck equation corresponding to the SDE (1) and initial condition (2) is given by

[TABLE]

Here $P(x,t)=P(x,t|x_{o},t_{o})$ describes the probability density of the position of the random walker at a particular time via the following relation:

[TABLE]

2 Derivation of the Probability Density

The first step in deriving the probability density is to use the Fourier transform on $P(x,t)$:

[TABLE]

The Fokker-Planck equation (5) now becomes

[TABLE]

This partial differential equation (PDE) can be solved by the method of characteristics (please refer Appendix B) to yield

[TABLE]

Note that the initial condition given by equation (9) and normalization given by equation (10) are easily verified. The characteristic function $\varphi(k,t)$ can be easily derived by the transformation $k\rightarrow-k$ in the Fourier transform, i.e. using the relation:

[TABLE]

This gives the characteristic function for the probability density governing the SDE (1) with initial condition (2) as

[TABLE]

This is the characteristic function of a non-central chi-squared distribution with zero degrees of freedom (please refer to Appendix D) where

The random variable has been scaled such that $x\rightarrow\frac{4}{\sigma^{2}(t-t_{o})}x$ 2. 2.

The non-centrality parameter is $\lambda=\frac{4x_{o}}{\sigma^{2}(t-t_{o})}$ 3. 3.

Then $P(x,t|x_{o},t_{o})=g\left(\frac{4x}{\sigma^{2}(t-t_{o})};\frac{4x_{o}}{\sigma^{2}(t-t_{o})}\right)$, where $g(x;\lambda)$ is the density function of a non-central chi-squared distribution with zero degrees of freedom with non-centrality parameter $\lambda$.

The non-centrality parameter [13] enables us to estimate the probability that the random variable takes the value 0 as $e^{-\lambda/2}$. In terms of the stochastic process (1) this is the probability the random walker gets absorbed [9] at the origin at time $t$ given the starting point $x_{o}$. We express these in terms of the survival function444$S(t)$ is the probability the random walker has survived up to time $t$., $S(t)$ and the corresponding absorption probability555$F(t)=1-S(t)$ and is the probability the random walker has been absorbed by time $t$. $F(t)$:

[TABLE]

Note that in the limit $t\rightarrow\infty$ this probability will tend to 1 which means that all random walkers following process (1) will eventually be absorbed666Alternatively $\lim_{t\rightarrow\infty}S(t)\rightarrow 0$. at the origin.

3 Sample Path Simulation

Using the fact that a non-central chi-squared distribution can be expressed as a Poisson compound mixture of central chi-square distribution (refer Appendix D) with even degrees of freedom [13] leads to a very simple formulation for the simulation of a random walk [6] corresponding to (1). We make use of items 1 and 2 along with the Markov property [9, 11] of the random walkers to describe the algorithm for generating a sample path:

For $i=1,2,\ldots,n$:

[TABLE] 2. 2.

Generate $N\sim$ Poisson ($\lambda/2$) 3. 3.

If $N=0$ set $X=0$ otherwise generate $X\sim\chi^{2}_{2N}$ 4. 4.

Set $x_{i}=cX$

The algorithm illustrated above was implemented in Rcpp (an implementation of C++ in R) and R was used in generating plots and data manipulation. The Rcpp implementation is given in Appendix E. A simple simulation of 100 such generated paths is shown in Figure 1.

3.1 Examining Survival/Absorption Probabilities

Using our simulations we can explore several properties of Feller Diffusion paths777Henceforth we set $t_{o}=0$ in both our simulation and computation.. We begin by looking at the absorption probability (13) and estimate this via the evolution of the (number of) simulated surviving paths. Our simulations match the intuition - i.e. more volatile walkers (i.e. larger values of $\sigma$) have lower survival rate and for large times the survival drops. Further the values are in close agreement with the theoretical values (Table 1 - using 10,000 simulations).

We also observe that the further away from the origin the starting point (larger values of $x_{o}$) the better the chance of survival - as per the intuition from equation (13). These results are displayed in (Table 2 - using 10,000 simulations) and show close agreement with the theoretical values:

3.2 Martingale Property and Average Position

As equation (1) is a martingale888An SDE with no drift term., the expected position999The averaging is with respect to $P(x,t)$. of a walker at time $t$ to be $x_{o}$ - calculation can be found in the Appendix C. We can examine this via our simulation by evaluating the average position of several walkers at a particular time via:

[TABLE]

The fact that $\langle x(t)\rangle=x_{o}$ (given $x(0)=x_{o}$) might seem paradoxical given that the Feller process is one where eventually all walkers are absorbed at the origin. This can be understood through the simulation where the walkers that survive (even though this number will be small) will most likely be very far away from the origin and thus driving up the average given by equation (14). Results are displayed in Tables 3 and 4. Note the error bounds (confidence intervals) for the Monte Carlo estimates are in agreement with the theory [6] - details can be found in Appendix C.

3.3 Hitting Time Distribution

From equation (13) we defined the absorption probability $F(t)$ and the corresponding survival probability $S(t)$, with the relation $S(t)=1-F(t)$. Given the nature of the Feller random walker, we can define these probabilities as:

[TABLE]

We define a new random variable, the hitting time or first-passage time, denoted $T^{*}$:

[TABLE]

We can then define the hitting (first passage) time density by [12]:

[TABLE]

where the relation between $f(t)$ and $T^{*}$ is given by

[TABLE]

This in turn allows us to express the relations (15) and (16) in terns of relationships between $T^{*}$ and $S(t)$, $F(t)$:

[TABLE]

Let us first turn our attention to a Brownian motion originating at $x_{o}$ described by the SDE $\dot{x}(t)=\sigma\xi(t)$ with the origin as an absorbing boundary. In this case, the hitting time density is given by [12]:

[TABLE]

whose integral over $t\in[0,\infty)$ is 1, indicating that absorption is certain. Paradoxically however the expected hitting time101010The averaging is with respect to the hitting time density $f(t)$. $\langle t\rangle=\langle T^{*}\rangle$ is infinite due to the large time power law behavior $f(t)\sim t^{-3/2}$ [10]. Setting $\frac{df}{dt}=0$ for equation (20) yields the most likely (typical) hitting time

[TABLE]

Now returning to our Feller walker we see an analogous situation to that of the Brownian walker. In this case the hitting time density is:

[TABLE]

We saw earlier from equation (13) that absorption is certain, in terms of the density this is equivalent to the Brownian case with the integral of the hitting density over $t\in[0,\infty)$ equal to 1. As with the Brownian case, due to the large time power law (in this case $f(t)\sim t^{-2}$) the expected hitting time is infinite:

[TABLE]

Similarly setting $\frac{df}{dt}=0$ for equation (21) yields the typical hitting time for the Feller walker as

[TABLE]

We use simulations to demonstrate the hitting time distributions for a Feller walker with $t_{0}=0$, $t_{n}=20000$, $\Delta t_{i}=1$, $x_{0}=100$ and $so^{2}=1,10$ (Figure 3). A reference line (shown in red) has been marked on the plots to indicate the theoretical hitting time $t^{*}$. It can be seen that this appears at the peak of the histograms confirming this result. The theoretical distributions indicated in equation 21 are overlaid (in blue) as well to illustrate how closely the simulation agrees with the theoretical distribution.

4 Summary

In this article, we re-examine the derivation of the probability density of the Feller Diffusion equation. As opposed to the prior approaches [2, 3, 4, 5], we use the Fourier transform on the Fokker-Planck equation corresponding to the Feller Diffusion and solve the resulting advection-type equation in Fourier space via the method of characteristics. The solution is easily mapped to the characteristic function of a scaled version of a non-central chi-squared distribution with zero degrees of freedom. As the theory of these distributions is well-known [13], it is straight forward to interpret the results in terms of survival probabilities and further develop the algorithm for developing sample paths. We carry out such a simulation and i) find the simulated survival probability is very close to the theoretical value and ii) prove the martingale property by estimating the average position of the random walkers using Monte Carlo. We also find the hitting time distribution, including the typical hitting time shows good agreement between the theoretical results and simulation. For future work we plan to apply the Feller diffusion in the setting of quantitative finance where it has the potential to be used in the following modelling contexts: i) price dynamics of a risky asset or interest rate, ii) derivative pricing on risky underlying assets, iii) bond pricing with risky interest rates, iv) default probability estimation, v) survival analysis including hazard rate modelling and regression methods.

5 Acknowledgements

To the best of our knowledge this work and our approach has not been carried out previously. Our only aim is to add to the body of knowledge in this area. This work has been funded and supported by MIND Analytics & Management. Contribution to this work as follows:

• •

Motivation for Investigation: RM/LK

• •

Fourier Transform + Method of Characteristics Approach: RM

• •

Simulation Algorithm Development: RM/LK

• •

Simulation, Coding & Results: PK

• •

General Discussion & Future Direction: RM/LK/PK.

• •

Writing: RM/PK/LK

The authors declare that there is no conflict of interest regarding the publication of this article.

Appendix A Fourier Transform of Fokker Planck Equation

In this section we look at deriving equation (8). We begin by re-writing equation (5) as

[TABLE]

Using the Fourier transform we obtain

[TABLE]

Note the use of following properties of Fourier transforms in deriving equation (24):

[TABLE]

Expanding the right hand side of equation (24) and simplifying yields the result in equation (8):

[TABLE]

Appendix B Solution via Method of Characteristics

Let us solve the following PDE for $u(z,t)$ using the method of characteristics:

[TABLE]

The method of characteristics postulates a solution of the form $u(z(s),t(s))$ satisfying

[TABLE]

Now equation (27) implies that along the characteristic curve $(z(s),t(s))$, which is parametrized by $s$, the solution is constant, i.e.

[TABLE]

To be explicit the $(z(s),t(s))$ and $(z(0),t(0))$ lie on the same characteristic curve. Matching the terms in equation (25) with equation (28) leads to the following system of ordinary differential equations (ODEs) for the characteristics:

[TABLE]

Equation (30) is solved by integrating as follows:

[TABLE]

Equation (31) is solved by integrating as follows:

[TABLE]

Using equation (29) along with equations (32) and (33) enables us to obtain the solution

[TABLE]

We see that (34) satisfies the initial condition $u(z,t_{o})=\psi(z)$. For the particular choice of initial condition (as in our problem):

[TABLE]

We make the identifications to map back to our original problem:

[TABLE]

It is then straight forward to use (35)to arrive at the final answer given by equation (11):

[TABLE]

Appendix C Mean and Variance of Feller Diffusion Process

We may use the cumulant generating function $H(k,t)$ to derive the cumulants of the Feller Diffusion, where

[TABLE]

In particular the mean is given by

[TABLE]

This matches our intuition as the random walk described by equation (1) is drift-less and thus a martingale. The variance of $x(t)$ is given by

[TABLE]

We use this value to obtain an error-bound on the Monte Carlo estimate of the average position, $\overline{x}(t)$ of a Feller random walker at time $t$ given by equation(14). By Monte Carlo theory [6] $\overline{x}(t)$ is distributed (approximately) as a normal with a mean $x_{o}$ (i.e. true value) and variance $\sigma^{2}x_{o}t/N$. This enables us to obtain a confidence interval for $x_{o}$ [6]:

[TABLE]

Here $\Phi^{-1}$ is the inverse cumulative normal distribution and $\alpha$ is the confidence level111111For the 95 percent confidence interval $\alpha=0.05$ and $\Phi^{-1}(0.975)=1.96$. For the 99 percent confidence interval $\alpha=0.01$ and $\Phi^{-1}(0.995)=2.58$.

To demonstrate this, we include here the results of the Monte Carlo simulations given in Section 3.2 with their 95% confidence intervals in Tables 5 and 6. It is noted that all the simulated results indicate that $\langle x(t)|x(0)=x_{o}\rangle=x_{o}$ within the error bounds.

Appendix D Zero Degrees of Freedom Non-central Chi-Squared Distribution

A non-central chi-squared distribution with zero degrees of freedom and non-centrality parameter $\lambda$ [13] is given by the density function:

[TABLE]

Here $q(x;\nu)$ is the density function of a central chi-squared distribution with $\nu$ degrees of freedom.

[TABLE]

Using the convention $q(x;0)=\delta(x)$ we may re-write equation (38) as

[TABLE]

Equations (38) and (40) describe a compound Poisson mixture of central chi-square distributions with even degrees of freedom [13]. The corresponding characteristic function $\varphi_{f}(k;\lambda)$ can be derived from (40) very easily:

[TABLE]

We have used the fact that the characteristic function of $q(x;2j)$ is $(1-2ik)^{-j}$ in the derivation. In addition, we note that we may write the density function in the alternative form [13, 5]

[TABLE]

Here $I_{1}(z)$ denotes the modified Bessel function of order 1 which is given by

[TABLE]

Substituting equations (39) and (43) in equation (40) yields

[TABLE]

Re-arranging the last line leads to equation (42). The distribution has a straight forward interpretation where the the probability that random variable takes the value 0 is $e^{-\lambda/2}$ and the probability that the random variable lies between $a(\neq 0)$ and $b$ is given by the integral of the second term in (42) over the corresponding limits. If $a=0$ we would include the $\delta$ function term in the integration. Finally, note for the Feller Diffusion process the density function (with $t_{o}$ set to 0) is given by [5]

[TABLE]

Appendix E Rcpp Implementation of Sample Path Simulation

The function CFellerPaths implemented below outputs a list containing time and the set of paths generated. The function allows to specify an initial seed, the number of paths required, a starting and ending time, $\Delta t_{i}$, the variance and the starting point.

⬇

#include<Rcpp.h>

#include<math.h>

#include<random>

//[[Rcpp::plugins(cpp11)]]

using namespace Rcpp;

//[[Rcpp::export]]

List CFellerPaths(int seed, int num, int t0, int tn, double dt, double var,

            **double** x0){

    **int** n = ceil(((**double**(tn - t0))/dt) + 1.0);

    NumericMatrix Paths(n, num);

    **static** std::mt19937 gen;

    NumericVector t(n);

    t[0] = t0;

    **for**(**int** j=0; j<num; j++){

            NumericVector x(n);

            x[0] = x0;

            gen.seed(seed);

            seed += 1;

            **for**(**int** i=1; i<n; i++){

                    **double** c = var/(4.0*dt);

                    **double** lambda = **float**(x[i-1])/c;

                    std::poisson_distribution<> pd((lambda/2.0));

                    **double** N = pd(gen);

                    **if**(N==0){

                            x[i] = 0.0;

                    }**else**{

                            std::chi_squared_distribution<> cd((2*N));

                            x[i] = c*cd(gen);

                    }

            }

    Paths(_,j) = x;

    }

    **for**(**int** i=1; i<n; i++){

            t[i] = t[i-1] + dt;

    }

    **return**(List::create(Named("t")=t,

    Named("paths")=Paths));

}