A Spectral Method Approach to Quadratic Normal Volatility Diffusions

Abstract

This paper is concerned with the quadratic volatility family of driftless stochastic differential equations (SDEs), also known in the literature as Quadratic Normal Volatility models (QNV), which have found applications primarily in mathematical finance, but can also model dynamics of stochastic processes in other fields such as mathematical biology and physics. These SDE models are characterized by a quadratic volatility term and can be reduced to one of four distinct possibilities depending on the roots of the quadratic volatility term and the position of initial value. We find explicit weak solutions for each case by a combination of Itô calculus and Fourier analysis, which can be described as a ’spectral method’. Furthermore, for all cases we also express the weak solutions as fairly simple functions of Brownian motion, which allows for efficient one-step Monte Carlo evaluation of functionals of the solutions $X_t$ of the form $\mathbb{E}\left(f\left(X_t\right)\right)$ and also more general functionals of the solution. The method used to compute the solutions may also be of interest itself more generally in other fields where SDEs play a fundamental role.

1. Introduction

QNV models have been studied, in financial mathematics in particular, for more than three decades. Several authors have suggested that a number of financial variables (such as foreign exchange rates, equity prices, or forward interest rates, see for instance [1,2,3]) can be modelled by SDEs of the type

$$d{X}_t=\sigma \left(X_t\right)d{W}_t$$

(1)

where the volatility term $\sigma \left(x\right)$ is a quadratic function of x and $W_t$ is standard Brownian motion. The condition in financial mathematics that there should exist a ’risk-neutral’ probability measure with regard to which asset process $X_t$ is a (local) martingale is trivially met by such equations since any solution to (1) is a local martingale by definition of the Itô integral, hence the interest in driftless SDEs as models for underlying variables of financial instruments. Of course, in the financial literature the focus is on the pricing of various derivatives (such as European call and put options), which are basically more or less complex functionals of the underlying process $X_t$, and not necessarily the dynamics of $X_t$ itself.

We suggest that dynamics and probabilistic properties of SDEs with quadratic volatility terms, driftless or not, can be interesting in themselves and can serve as viable models for stochastic processes not only in financial mathematics.

As a simple example in evolutionary dynamics, consider for instance a gene that has lost its evolutionary function, but is also not detrimental, meaning that no positive or negative evolutionary pressure is applied to the frequency of that gene in a population. What could be a plausible large population and continuous-time model for the frequency of such a gene over time? Reasonable assumptions for an equation that models the frequency $X_t\in [0,1]$ of such a gene would be that changes in the value of $X_t$ are (1) purely random and that (2) during a short period of time $X_t$ is just as likely to increase by a given amount as it is to decrease by that same amount. Put another way, the expected net change in $X_t$ after a while should be zero (which agrees with the notion of a martingale). Also, the amount of change in frequency during a short time interval should be (3) proportional both to the number of individuals that possess the gene as well as the number that do not possess the gene (although for this assumption to be reasonable a more specific description of the setup would be required). Then, the simplest possible SDE that fulfils these assumptions would be

$$\begin{array}{ccc}\hfill d{X}_t& =& a{X}_t(1-X_t)d{W}_t\hfill \\ \hfill X_0& =& x_0\in [0,1]\hfill \end{array}$$

(2)

for some constant a, which is an example of Equation (1).

In mathematical physics, there are also situations where processes of the type (1) or related SDEs can arise naturally. We feel that giving a sufficiently detailed example here somewhat exceeds the intended scope of this paper, but we are currently working on a separate paper where our results here can be applied to just such a case.

From a purely theoretical viewpoint, QNV models are also interesting as a perfect example of the difference between a true martingale and a strict local martingale, and how much more care is needed when dealing the latter.

The paper is organized as follows: in Section 2, we give a short summary of the existing literature on QNV models and the additional insight we hope to present here. The list of references here is not meant to be exhaustive. The intention is to provide some context for our paper, highlight the original results and aspects of our contribution, and contrast them to the results in the existing literature. In Section 3, we define the four classes of QNV models in terms of four ‘dimensionless’ SDEs, and summarize some qualitative probabilistic properties of our diffusions, which follow from general SDE theory and are used to justify the subsequent analysis. In Section 4, we describe the general procedure used to obtain our results and list the eigenfunctions of the diffusion ‘generators’ for these equations along with the transformations of the equations we will be using. In Section 5, we obtain expressions for the weak solutions for each case, along with alternative expressions relating the solutions to (possibly reflected) Brownian motion. Concise forms of the weak solutions to all four cases are then listed in Section 6, which also contains some comments and interpretations of the results. In Section 7, we present the results of numerical simulations, which demonstrate the accuracy of the obtained results and the efficiency of the proposed Monte Carlo schemes for evaluation of expectations. Finally, in Section 8, we present our conclusions and outline a few possible avenues for further research on this or related topics.

2. Our Results in Context with Existing QNV Model Literature

Initially, QNV models were approached mainly by Markov process and partial differential equation means (since SDEs can in principle be associated to Markov process diffusion generators) leading to heat equation problems for specific prices of financial derivatives of QNV modelled assets (i.e., functionals of solutions to (1)). This approach, modelled on the original Black–Scholes formulae derivation, is outlined in [4], and carried out for specific root configurations and option prices in for instance in [5,6].

However, as subsequently explained in [7], this approach has a fundamental flaw in the fact that for most root configurations solutions to (1) are local martingales that are not true martingales. This in practice means that the naive correspondence between such an SDE and its Markov process diffusion generator breaks down which in turn leads to incorrect pricing formulas derived by these methods. A more careful probabilistic treatment was therefore needed, which accounts for the strict local martingale nature of (1) (for most root configurations). This was carried out in [7], leading to correct pricing formulas for several choices of option prices and ’non-trivial’ root configurations.

Even more was achieved in [8] regarding a purely probabilistic (or more specifically, martingale theory) treatment of QNV models, where even a characterization of QNV models was obtained as the only local martingales that can be expressed as certain transformations of stopped Brownian motion, that also leads to pricing formulas for very general functionals of the underlying QNV model (see Theorem 1 and Corollary 1 in [8]). The comprehensive nature of the results in [8] probably explains the apparent lack of more recent significant theoretical results regarding the topic.

In this paper, we offer a more analytically minded alternative approach for dealing with the strict local martingale issue, which leads to a somewhat different but similarly complete probabilistic description of QNV models as to that in [8].

Roughly speaking, in both [7,8], the problem of taking expectations of local martingales was handled by localisation, i.e., stopping times and limiting arguments, a standard and natural technique was given the definition of a local martingale, and in [8] the QNV models were associated to stopped Brownian motion via a set of martingale measure changes using (generalisations of) the Girsanov theorem.

Here, the basic idea is to identify sufficiently large families of bounded functions of the solution $X_t$ to (1) that also happen to be local martingales, relying on the fundamental fact that a bounded local martingale is a true martingale. This leads to eigenfunction problems via the Itô formula, and it turns out that taking expectations of these families of true martingales actually produces variants of Fourier transforms of the transition probability measures from which weak solutions to (1) can be obtained.

The general idea and, to a degree, also the analysis is similar for all root configurations, but the properties of the solutions as well as the details of the analysis vary significantly. For instance, for the case of two real roots with initial values located between the roots, the transition measure of a transformed process can be expressed as a weighted average of two moving Gaussian measures, while for the case of two complex roots, the transition kernel can be expressed in terms of classical Jacobi theta functions and related to reflected Brownian motion.

We believe that the treatment presented here, which is more oriented towards classical analysis, has some merits. For one, our approach seems conceptually comparatively straightforward and yields fairly simple explicit expressions for the transition densities in the Markov process picture of QNV models. These transition densities are all expressed using elementary functions (with the exception of theta functions in one case, which are nevertheless well known and have efficient numerical implementations in common software) which is clearly just as useful for analytical or numerical computations of expectations for quite large classes of functions of paths. This contrasts with some results in the literature where expectations of specific functions are expressed in terms of convergent function series, which generally converge quite slowly (see [7]).

Comparing to the results in [8] where expectations in QNV models are expressed in terms of martingale changes of measure that depend on the terminal value of stopped Brownian motion, a theoretical novelty here lies in the observation that expectations can also be expressed in terms of martingale measure changes of reflected Brownian motion. These expressions for expectations of the form $\mathbb{E}\left(f\left(X_t\right)\right)$ lead to similarly general expressions for path dependent functionals to those in [8] (see Section 6) and also efficient Monte-Carlo numerical methods for evaluating functionals where an analytical expression may not be available. Perhaps not too surprisingly, under close inspection it turns out that the specific measure change processes obtained in each case happen to be basically identical to those in [8] with the main difference being, as mentioned, that here the Brownian motion is reflected on the boundary as opposed to stopped. A surprising additional feature of several of these measure changes is that the martingale that performs the measure change is not strictly positive in contrast to the usual Girsanov-type measure changes.

Thus, we believe these results complement those in [8] nicely and contribute to better understanding of QNV models, especially since they are obtained from a largely independent perspective.

Interestingly, these same measure changes were also identified, in a somewhat heuristic fashion, in [9] from yet another point of view, namely Doobs h-transform (although certain assertions made in [9] are incorrect).

A minor novel result is also the observation that in the imaginary root configuration the solution to (1) has simple limit distribution, which seems to have gone unnoticed in the literature. Although this appears at first as counterintuitive (indeed, there are claims in the literature that this diffusion even explodes in finite time), and it is not at all obvious even from the results in [8], this fact becomes immediately evident here.

3. Preliminaries

We are seeking closed form expressions for weak solutions to Equation (1) with

$$\sigma \left(x\right)={(ax+b)}^2+c$$

The nature of the solution $X_t$ clearly depends on the sign of c and in the case of two real roots also the position of the initial value $X_0=x$.

By linearly rescaling of $X_t$ and t, it is straightforward to see that any Equation (1) with quadratic volatility term $\sigma$ can be reduced to one of the following four cases.

• Two real roots

  $$\begin{array}{ccc}\hfill d{X}_t& =& (1-X_t^2)d{W}_t\hfill \end{array}$$

  (3)

  with two choices for the initial condition $X_0=x$

  $$\begin{array}{ccc}\hfill X_0& =& x\in (-1,1)\hfill \\ & \mathrm{o}\mathrm{r}& \hfill \end{array}$$

  (4)

  $$\begin{array}{ccc}\hfill X_0& =& x\in (1,\infty )\hfill \end{array}$$

  (5)
  The case $x<-1$ is essentially the same as $x>1$ by symmetry $X\mapsto -X$ and $W\mapsto -W$.
• Imaginary roots

  $$\begin{array}{ccc}\hfill d{X}_t& =& (1+X_t^2)d{W}_t\hfill \\ \hfill X_0& =& x\in \mathbb{R}\hfill \end{array}$$

  (6)
  In many ways, this is the example with the most interesting analysis.
• Double roots

  $$\begin{array}{ccc}\hfill d{X}_t& =& X_t^{2}d{W}_t\hfill \\ \hfill X_0& =& x\in (0,\infty )\hfill \end{array}$$

  (7)
  Again, the case $x<0$ is the same by symmetry.

Our reference for properties of the general Equation (1) is [10], specifically Section 5.5. These references are of fundamental importance in the following since they not only guarantee the unique existence of the solutions to our equations, but also determine their precise domains and asymptotic properties.

Theorem 5.4. in [10] shows that a driftless SDE has a non-exploding weak solution for any initial distribution at least for continuous functions $\sigma \left(x\right)$. Using Theorem 5.7. in [10], we can also directly establish uniqueness in law for the solutions to all our Equations (3), (6) and (7). Furthermore, Theorem 5.9 in [10] gives sufficient conditions for pathwise uniqueness and strong uniqueness. These conditions can also be verified to hold for all our equations.

In the following, we will therefore rely on the fact that continuous global non-exploding unique weak and strong solutions exist for all the equations we will be dealing with.

Instead of dealing with Equations (3), (6) and (7) directly, we will actually prefer to work with transformed versions of these equations.

By the Itô formula, we have for $\phi \in C^2$

$$d\phi \left(X_t\right)=\sigma \left(X_t\right){\phi }^{\prime }\left(X_t\right)d{W}_t+\frac{1}{2}{\sigma }^2\left(X_t\right){\phi }^{″}\left(X_t\right)dt$$

(8)

The transformations we have in mind are those that reduce the local martingale term above to Brownian motion (at the expense of introducing drift). To be explicit, we will be working with equations for $Y_t=g\left(X_t\right)$ where

$$g\left(x\right)=\int \frac{1}{\sigma \left(x\right)}dx$$

(9)

We will refer to such a function as a natural rescaling of $X_t$, although it should be noted that in standard terminology $X_t$ would be regarded as the natural scale of the process $Y_t$ since $X_t$ is driftless.

Of course the fact that the processes (3), (6) and (7) do not explode as processes on $\mathbb{R}$ does not mean they could not explode when viewed as processes on subintervals $I\subset \mathbb{R}$. For instance, it is not immediately obvious that the solution to (7) does not reach the value 0 in finite time.

In order to determine the asymptotic behaviour of our processes $X_t$ and $Y_t$, we apply Feller’s test for explosions as described in [10], Section 5.5.C.

By carrying out the calculations, case by case, required for the application of Proposition 5.22 and Theorem 5.29 in [10], we can verify that

• The solution to (3) with initial value (4) remains on the interval $(-1,1)$ for all $t<\infty$ a.s., and has an a.s. limit ${lim}_{t\to \infty }{X}_t=Z$, where Z is a random variable with distribution $\mathbb{P}(Z=1)=1-\mathbb{P}(Z=-1)=(x+1)/2$.
• The solution to (3) with initial value (5) remains on the interval $(1,\infty )$ for all $t<\infty$ a.s. and has an a.s. limit ${lim}_{t\to \infty }{X}_t=1$.
• The solution to (6) remains on $(-\infty ,\infty )$, i.e., $X_t$ has no explosions (which we already know).
• The solution to (7) remains on the interval $(0,\infty )$ for all $t<\infty$ a.s. and has an a.s. limit ${lim}_{t\to \infty }{X}_t=0$.

Another thing worth mentioning is that we will be referring to the differential operator implicit in the bounded variation term of (8)

$${\mathcal{L}}^X=\frac{1}{2}{\sigma }^2\left(x\right)\frac{{\partial }^2}{\partial x^2}$$

(10)

as the ‘generator’ of $X_t$ and likewise for the operators found in the equations for $Y_t$. However, we will not be interpreting ${\mathcal{L}}^X$ as an infinitesimal generator of a diffusion defined by

$${\mathcal{L}}^X\phi =\underset{t\to 0}{lim}{\dot{P}}_t\phi$$

for a Markov semigroup $P_t$ unless or until we have a specific semigroup $P_t$ in mind. We will be viewing such ‘generators’ primarily as linear operators acting on $C^2$ functions.

4. Eigenfunction Analysis and Rescaling

A known property of our equations is that their solutions are local martingales, but not martingales (except for the bounded interval case (3) with initial value (4)), as can be verified using the results from [11]. This clearly complicates taking expectations of the solutions or functions of the solutions.

As mentioned, here, the general idea is to use eigenfunctions of the generators associated to our SDEs to construct sufficiently large families of bounded local martingales (and hence true martingales), which then almost directly produce Fourier transforms, or Fourier series, of the transition densities of the natural rescalings of our processes simply by the defining property of martingales.

The results of this eigenfunction analysis are listed bellow. The eigenfunctions can be found using standard techniques for solving second order ordinary differential equations (our reference is [12]), but the fact that the listed functions are indeed eigenfunctions can of course be verified directly.

For each of our equations for $X_t$, we list the generators ${\mathcal{L}}^X$ from (10) together with their eigenfunctions and the natural rescalings of the equations for the convenient choice of primitive function in (9).

Since we are imposing no boundary conditions we obtain a two-dimensional eigenspace of ${\mathcal{L}}^X$ with eigenvalue $\lambda \left(w\right)=\lambda (-w)$ for each $w\in \mathbb{C}$ (see (11), (20) and (25)). The forms of these eigenfunctions in each case are chosen so that we obtain a real-valued eigenfunction for $w\in \mathbb{R}$.

For the rescaled processes $Y_t$ we also list their generators and eigenfunctions for the corresponding eigenvalues, which are the same as for the generators of $X_t$.

An important thing to note are the supports of the random variables $X_t$ and $Y_t$ (for $t<\infty$). For $X_t$, the fact that our equations have no explosions in finite time with regards to the supports listed below was verified using the Fellers test in the previous section. For $Y_t$ the support for each case can be obtained by considering the image of the specific scale mapping (9), or independently using Fellers test.

• For Equation (3) with initial value (4), we have

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}X& :& {\mathcal{L}}^X\phi =\frac{1}{2}{(1-x^2)}^2{\phi }^{\prime \prime }\hfill \\ \hfill \mathrm{support}\mathrm{of}{X}_t& :& (-1,1)\hfill \\ \hfill \mathrm{eigenvalues}& :& \lambda \left(w\right)=\frac{w^2-1}{2}\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^X& :& {\psi }_X(x,w)=\sqrt{1-x^2}{e}^{\pm watanh\left(x\right)}\hfill \\ \hfill \mathrm{natural}\mathrm{rescaling}& :& g\left(x\right)=\frac{1}{2}log\frac{1+x}{1-x}=atanh\left(x\right)\hfill \\ & & Y_t=atanh\left(X_t\right)\hfill \\ \hfill & & X_t=tanh\left(Y_t\right)\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \mathrm{equation}\mathrm{for}Y& :& d{Y}_t=d{W}_t+tanh\left(Y_t\right)dt\hfill \\ \hfill \mathrm{support}\mathrm{of}{Y}_t& :& (-\infty ,\infty )\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}Y& :& {\mathcal{L}}^Y\phi =\frac{1}{2}{\phi }^{\prime \prime }+tanh\left(y\right){\phi }^{\prime }\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^Y& :& {\psi }_Y(y,w)=\frac{e^{\pm wy}}{cosh\left(y\right)}\hfill \end{array}$$

(15)

• For Equation (3) and initial value (5), we have

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}X& :& {\mathcal{L}}^X\phi =\frac{1}{2}{(1-x^2)}^2{\phi }^{\prime \prime }\hfill \\ \hfill \mathrm{support}\mathrm{of}{X}_t& :& (1,\infty )\hfill \\ \hfill \mathrm{eigenvalues}& :& \lambda \left(w\right)=\frac{w^2-1}{2}\hfill \\ \hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^X& :& {\psi }_X(x,w)=\sqrt{x^2-1}{e}^{\pm wacotanh\left(x\right)}\hfill \\ \hfill \mathrm{natural}\mathrm{rescaling}& :& g\left(x\right)=\frac{1}{2}log\frac{x+1}{x-1}=acotanh\left(x\right)\hfill \\ & & Y_t=acotanh\left(X_t\right)\hfill \\ & & X_t=cotanh\left(Y_t\right)\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \mathrm{equation}\mathrm{for}Y& :& d{Y}_t=d{W}_t+cotanh\left(Y_t\right)dt\hfill \\ \hfill \mathrm{support}\mathrm{of}{Y}_t& :& (0,\infty )\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}Y& :& {\mathcal{L}}^Y\phi =\frac{1}{2}{\phi }^{\prime \prime }+cotanh\left(y\right){\phi }^{\prime }\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^Y& :& {\psi }_Y(y,w)=\frac{e^{\pm wy}}{sinh\left(y\right)}\hfill \end{array}$$

(19)

• For Equation (6), we have

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}X& :& {\mathcal{L}}^X\phi =\frac{1}{2}{(1+x^2)}^2{\phi }^{\prime \prime }\hfill \\ \hfill \mathrm{support}\mathrm{of}{X}_t& :& (-\infty ,\infty )\hfill \\ \hfill \mathrm{eigenvalues}& :& \lambda \left(w\right)=\frac{1+w^2}{2}\hfill \\ \hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^X& :& {\psi }_X(x,w)=\sqrt{1+x^2}{e}^{\pm warctan\left(x\right)}\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill \mathrm{natural}\mathrm{rescaling}& :& g\left(x\right)=arctan\left(x\right)\hfill \\ & & Y_t=arctan\left(X_t\right)\hfill \\ & & X_t=tan\left(Y_t\right)\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill \mathrm{equation}\mathrm{for}Y& :& d{Y}_t=d{W}_t-tan\left(Y_t\right)dt\hfill \\ \hfill \mathrm{support}\mathrm{of}{Y}_t& :& (-\pi /2,\pi /2)\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}Y& :& {\mathcal{L}}^Y\phi =\frac{1}{2}{\phi }^{\prime \prime }-tan\left(y\right){\phi }^{\prime }\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^Y& :& {\psi }_Y(y,w)=\frac{e^{\pm wy}}{cos\left(y\right)}\hfill \end{array}$$

(24)

• For Equation (7), we have

$$\begin{array}{ccc}\hfill \mathrm{generator}\mathrm{of}X& :& {\mathcal{L}}^X\phi =\frac{1}{2}{x}^4{\phi }^{\prime \prime }\hfill \\ \hfill \mathrm{support}\mathrm{of}{X}_t& :& (0,\infty )\hfill \\ \hfill \mathrm{eigenvalues}& :& \lambda \left(w\right)=\frac{w^2}{2}\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^X& :& {\psi }_X(x,w)=xexp(\pm w/x)\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill \mathrm{natural}\mathrm{rescaling}& :& g\left(x\right)=\frac{1}{x}\hfill \\ & & Y_t=\frac{1}{X_t}\hfill \\ & & X_t=\frac{1}{Y_t}\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill \mathrm{equation}\mathrm{for}Y& :& d{Y}_t=-d{W}_t+\frac{1}{Y_t}dt\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill \mathrm{support}\mathrm{of}{Y}_t& :& (0,\infty )\hfill \\ \hfill \mathrm{generator}\mathrm{of}Y& :& {\mathcal{L}}^Y\phi =\frac{1}{2}{\phi }^{\prime \prime }+\frac{1}{y}{\phi }^{\prime }\hfill \\ \hfill \mathrm{eigenspace}\mathrm{of}{\mathcal{L}}^Y& :& {\psi }_Y(y,w)=exp(\pm wy)/y\hfill \end{array}$$

(29)

For each of our Equations (3) with (4) and (5), (6) and (7), corresponding eigenfunctions ${\psi }_X$ and ${\psi }_Y$ listed above, and a fixed value of $w\in \mathbb{C}$, applying the Itô formula shows that

$$L_t^w:=e^{-\lambda \left(w\right)t}{\psi }_X(X_t,w)=e^{-\lambda \left(w\right)t}{\psi }_Y(Y_t,w)$$

(30)

is (in general) a complex-valued local martingale. In general, the eigenfunctions ${\psi }_X$ and ${\psi }_Y$ listed above, defined on their respective supports for $X_t$ and $Y_t$, are unbounded. But it is possible to find linear subspaces of functions spanned by these eigenfunctions that are bounded. The corresponding subspaces spanned by the local martingales $L_t^w$ in (30) will then contain only bounded martingales. As we show, taking expectations of these martingales will yield enough information about the processes $X_t$ and $Y_t$ to be able to explicitly determine their probability measures.

At this point, our analysis of the equations diverges, as we need to consider each case separately.

5. Weak Solutions

5.1. Two Real Roots, Bounded Interval $(-1,1)$

The first case we examine is the simplest, Equation (3) with initial value (4). This is also the only case where the solution is a true martingale. For illustration, a few sample path are simulated in Figure 1.

For the functions ${\psi }_Y$ from (15), we notice that for $Re\left(w\right)=0$ and $y\in \mathbb{R}$ we have

$$|{\psi }_Y(y,w)|=\left|\frac{e^{wy}}{cosh\left(y\right)}\right|\le 1$$

and the same goes for $|{\psi }_X(x,w)|$ for $x\in [-1,1]$. Hence, the local martingales from (30)

$$L_t^{ui}=e^{-\lambda \left(ui\right)t}{\psi }_Y(Y_t,ui)=e^{\frac{1+u^2}{2}t}\frac{e^{ui{Y}_t}}{cosh\left(Y_t\right)}$$

are bounded martingales for all $u\in \mathbb{R}$. By taking expectations, we then establish that

$$\mathbb{E}\left(\frac{e^{ui{Y}_t}}{cosh\left(Y_t\right)}|Y_0=x\right)=e^{-\frac{1+u^2}{2}t}\frac{e^{uix}}{cosh\left(x\right)}$$

or

$${\int }_{-\infty }^{\infty }\frac{e^{uiy}}{cosh\left(y\right)}\mathbb{P}(Y_t\in dy|Y_0=x)=e^{-\frac{1+u^2}{2}t}\frac{e^{uix}}{cosh\left(x\right)}$$

(31)

which means we have found the Fourier transform of the expression under the integral sign. Without too much work we can invert the right-hand side of (31). To arrive at the exact form stated below requires an amount of algebraic manipulation, but the details are omitted since they are all elementary.

Proposition 1.

Let $Y_t$ be the solution to Equation (13) with initial condition $Y_0=x\in \mathbb{R}$. Then,

$$\mathbb{P}(Y_t\in dy|Y_0=x)=\frac{1}{\sqrt{2\pi t}}\frac{cosh\left(y\right)e^{-\frac{t}{2}}}{cosh\left(x\right)}{e}^{-\frac{{(y-x)}^2}{2t}}dy$$

(32)

An interesting way of rewriting (32) is

$$\mathbb{P}(Y_t\in dy|Y_0=x)=\frac{1}{\sqrt{2\pi t}}\left(\frac{1}{1+e^{-2x}}{e}^{-\frac{{(y-x-t)}^2}{2t}}+\frac{1}{1+e^{2x}}{e}^{-\frac{{(y-x+t)}^2}{2t}}\right)dy$$

This shows the probability measure of $Y_t$ is actually a weighted average of diffusing Gaussian measures moving in opposite direction with speed 1 and also produces a neat expression for the weak solution to (13).

Corollary 1.

Let $Y_t$ be the solution to (13) and let ${\tilde{W}}_t$ be a Brownian motion with ${\tilde{W}}_0=x$. For Wiener integrable functions, f we have

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\frac{1}{1+e^{-2x}}\mathbb{E}(f({\tilde{W}}_t+t))+\frac{1}{1+e^{2x}}\mathbb{E}(f({\tilde{W}}_t-t))$$

Another way of interpreting (32) is as a measure change of Brownian motion.

Corollary 2.

Let $Y_t$ and ${\tilde{W}}_t$ be as before. Define the martingale

$$M_t=\frac{cosh(\widetilde{W_t})e^{-\frac{t}{2}}}{cosh(x)}$$

with $M_0=1$. For Wiener integrable functions f, we have

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\mathbb{E}(M_{t}f({\tilde{W}}_t))$$

(33)

5.2. Two Real Roots, Half-Line $(1,\infty )$

Here, we turn to the two roots case (3) with initial value to the right of the roots. It is well-known the solution is a supermartingale. Some simulations are shown in Figure 2.

An important difference compared to the finite interval case is that the functions (19) have singularities at 0 for any $w\in \mathbb{C}$. But this is only a minor difficulty because there are plenty of linear combinations of functions from (19) that are bounded and continuous on $[0,\infty )$. For $u\in \mathbb{R}$, we define

$${\varphi }_Y(y,u)=\frac{1}{2i}\left({\psi }_Y(y,ui)-{\psi }_Y(y,-ui)\right)=\frac{sin\left(uy\right)}{sinh\left(y\right)}$$

and confirm that ${\varphi }_Y(y,u)$ has a limit when $y\to 0$ for any $u\in \mathbb{R}$ and is still an eigenfunction of (18) with $\lambda \left(ui\right)=\lambda (-ui)=-\frac{1+u^2}{2}$. The Itô formula then shows that

$$M_t^u=e^{-\lambda \left(ui\right)t}{\varphi }_Y(Y_t,u)=e^{\frac{1+u^2}{2}t}\frac{sin\left(u{Y}_t\right)}{sinh\left(Y_t\right)}$$

(34)

are bounded martingales for all $u\in \mathbb{R}$, which allows us to take expectations at any time $t<\infty$. Since we know the distribution of $Y_t$ is supported on $(0,\infty )$, this yields the transition densities via the sine Fourier transform.

Proposition 2.

Let $Y_t$ be the solution to (17) with $Y_0=x\in (0,\infty )$. Then,

$$\mathbb{P}(Y_t\in dy|Y_0=x)=\frac{1}{\sqrt{2\pi t}}\frac{sinh\left(y\right)e^{-t/2}}{sinh\left(x\right)}\left(e^{-\frac{{(y-x)}^2}{2t}}-e^{-\frac{{(y+x)}^2}{2t}}\right)dy$$

(35)

for $y\ge 0$ (and $=0$ for $y<0$).

Proof. 

Taking expectation conditional on $Y_0=x$ of the bounded martingale (34) we obtain

$$\mathbb{E}\left(\frac{sin\left(u{Y}_t\right)}{sinh\left(Y_t\right)}|Y_0=x\right)=e^{-\frac{1+u^2}{2}t}\frac{sin\left(ux\right)}{sinh\left(x\right)}$$

or

$${\int }_0^{\infty }\frac{sin\left(uy\right)}{sinh\left(y\right)}\mathbb{P}(Y_t\in dy|Y_0=x)=e^{-\frac{1+u^2}{2}t}\frac{sin\left(ux\right)}{sinh\left(x\right)}$$

Computing the inverse sine Fourier transform of the right-hand side (as a function of u), we obtain the expression (35) for the density of $Y_t$ conditional on $Y_0=x$. □

By elementary algebraic manipulation of (35), it is also possible to express the weak solution to (17) in terms of Brownian motion with drift reflected at the origin.

Corollary 3.

Let $Y_t$ be the solution to (17) and ${\tilde{W}}_t$ be a Brownian motion with ${\tilde{W}}_0=x>0$. For Wiener integrable functions f on $(0,\infty )$, we have

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\frac{1}{1-e^{-2x}}\mathbb{E}\left(f\right(|{\tilde{W}}_t+t\left|\right))-\frac{1}{e^{2x}-1}\mathbb{E}\left(f\right(|{\tilde{W}}_t-t\left|\right))$$

(36)

Proof. 

We denote the weights above by $u_1\left(x\right)=1/(1-e^{-2x})$ and $u_2\left(x\right)=1/(e^{2x}-1)$ (by the way, we see that $u_1\left(x\right)-u_2\left(x\right)=1$). With a little reworking of (35), we can write

$$\begin{array}{ccc}\hfill \mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)& =& \frac{1}{\sqrt{2\pi t}}\left(u_1\left(x\right){\int }_0^{\infty }f\left(y\right)e^{-\frac{{(y-x-t)}^2}{2t}}dy-u_2\left(x\right){\int }_0^{\infty }f\left(y\right)e^{-\frac{{(y-x+t)}^2}{2t}}dy\right.\hfill \\ & & \left.+u_1\left(x\right){\int }_0^{\infty }f\left(y\right)e^{-\frac{{(y+x+t)}^2}{2t}}dy-u_2\left(x\right){\int }_0^{\infty }f\left(y\right)e^{-\frac{{(y+x-t)}^2}{2t}}dy\right)\hfill \end{array}$$

Clearly, we need f to be defined only on $(0,\infty )$. But we can extend the domain of f to $(-\infty ,\infty )$ by defining $f(-y)=f\left(y\right)$ and substitute $y\mapsto -y$ in the third and fourth integrals above to rewrite the expression as

$$=\frac{1}{\sqrt{2\pi t}}\left(u_1\left(x\right){\int }_{-\infty }^{\infty }f\left(y\right)e^{-\frac{{(y-x-t)}^2}{2t}}dy-u_2\left(x\right){\int }_{-\infty }^{\infty }f\left(y\right)e^{-\frac{{(y-x+t)}^2}{2t}}dy\right)$$

which is exactly what one obtains while computing the right-hand side of (36). □

Yet another way of interpreting the weak solution to (17) given by (35) is as a martingale measure change of Brownian motion reflected at the origin, albeit a rather odd kind of measure change.

Corollary 4.

Let $Y_t$ and ${\tilde{W}}_t$ be as above. Define the martingale

$$M_t=\frac{sinh({\tilde{W}}_t)e^{-t/2}}{sinh\left(x\right)}$$

with $M_0=1$. Then, for Weiner integrable functions f on $(0,\infty )$, we have

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\mathbb{E}(M_{t}f\left(\right|{\tilde{W}}_t\left|\right))$$

(37)

Proof. 

Similarly to the proof of the previous Corollary, we write out the expression for $\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)$ according to (35), define $f(-y)=f\left(y\right)$, substitute $y\mapsto -y$ in the second integral and remember that $sinh(-y)=-sinh\left(y\right)$ to find that the result is the same as computing the right-hand side of (37). □

The ‘measure change’ expressed in (37) is of course strange considering the martingale $M_t$ attains negative values. Other examples of curios ‘negative measure changes’ akin to (37) will be found in the next subsections. We attempt a heuristic explanation of such phenomena in Section 6.

5.3. Imaginary Roots

The most mathematically interesting case is (6). For illustration, some sample path simulations are shown in Figure 3.

Here, we find only countably many linear combinations of the eigenfunctions in (24), which are bounded on $(-\frac{\pi }{2},\frac{\pi }{2})$, but this will suffice.

For each $n\in \mathbb{N}$, we can define

$$\begin{array}{ccc}\hfill {\varphi }_{2n-1}\left(x\right)& =& \frac{1}{2}\left({\psi }_Y(x,(2n-1)i)+{\psi }_Y(x,-(2n-1)i)\right)=\frac{cos\left(\right(2n-1\left)x\right)}{cos\left(x\right)}\hfill \\ \hfill {\varphi }_{2n}\left(x\right)& =& \frac{1}{2i}\left({\psi }_Y(x,ni)-{\psi }_Y(x,-2ni)\right)=\frac{sin\left(2nx\right)}{cos\left(x\right)}\hfill \end{array}$$

(38)

and confirm these functions have well defined limits when $x\to \pm \frac{\pi }{2}$ and can thus be considered $C([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ functions. Note also that since $\lambda (-w)=\lambda \left(w\right)$ the functions (38) remain eigenfunctions of ${\mathcal{L}}^Y$.

Clearly, these functions are related to the usual Fourier basis. In particular, the odd functions ${\varphi }_{2n-1}\left(x\right)$ are a rescaled variant of the Dirichlet kernel. They are not quite the same as the usual basis though, and it seems necessary to prove that (38) actually provide a worthwhile basis of functions. This is probably most easily accomplished by appealing to the Stone–Weierstrass theorem.

Proposition 3.

The linear hull of ${\left\{{\varphi }_k\left(x\right)\right\}}_{k=1,\dots }$ from (38) is dense in $C([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ with regard to the topology of uniform convergence.

Proof. 

Using basic trigonometric identities one finds that

$$\begin{array}{ccc}\hfill cos\left(2nx\right)& =& \frac{1}{2}\left({\varphi }_{2n-1}\left(x\right)+{\varphi }_{2n+1}\left(x\right)\right)\hfill \\ \hfill sin\left(\right(2n+1\left)x\right)& =& \frac{1}{2}\left({\varphi }_{2n}\left(x\right)+{\varphi }_{2n+2}\left(x\right)\right)\hfill \end{array}$$

(39)

The usual Fourier basis on $[-\frac{\pi }{2},\frac{\pi }{2}]$ consists of $cos\left(2nx\right)$ and $sin\left(2nx\right)$ and the difference here is that the sine functions are shifted. Clearly then, (38) and (39) span the same space (in the sense that any finite linear combination of the functions in (38) can be expressed as a finite linear combination of (39) and vice-versa), so it suffices to prove the claim for (39).

According to the Stone–Weierstrass theorem we need to see that the functions (39) form a unital algebra (for point-wise products) that separates points.

The former can be verified using basic trigonometric identities. Indeed, one can see this algebra is generated by $\{1,sin(x\left)\right\}$.

The function $sin\left(x\right)$ is a one-to-one map on $[-\frac{\pi }{2},\frac{\pi }{2}]$ so it separates any distinct pair $x,y\in [-\frac{\pi }{2},\frac{\pi }{2}]$, including the edge points. So, in contrast to the usual Fourier basis, there is even no need to identify the endpoints $-\frac{\pi }{2}$ and $\frac{\pi }{2}$. □

We now view (23) as an operator acting upon the $C^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ subspace of an $L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ space for an appropriate weighted inner product

$${\langle f,g\rangle }_{\rho }={\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(x\right)g\left(x\right)\rho \left(x\right)dx$$

The weight $\rho$ by which (23) becomes a symmetric operator (and hence with orthogonal eigenspaces) on the $L_{\rho }^2$ space is determined by

$$\rho \left(x\right)=A{e}^{-2\int tan\left(x\right)dx}=A{cos}^2\left(x\right)$$

for any non-zero constant A (see [12]). We choose the constant A to make $\rho$ a probability density which means

$$\rho \left(x\right)=\frac{2}{\pi }{cos}^2\left(x\right)$$

(40)

With the choice (40) for the weight the eigenfunctions (38) define an orthonormal basis for a Hilbert space ${\mathcal{H}}_{\rho }\subseteq L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$, meaning no additional normalizing constants are required.

Note that the Hilbert space ${\mathcal{H}}_{\rho }$ spanned by (38) contains at least $C([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ by Proposition 3 since uniform convergence implies convergence in the $L_{\rho }^2$ norm.

Since indicator functions of Borel subsets of $[-\frac{\pi }{2},\frac{\pi }{2}]$ can be approximated by $C([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ functions and $L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ functions can approximated by linear combinations of indicators we can assert that we actually have equality of Hilbert spaces

$${\mathcal{H}}_{\rho }=L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$$

The point of using the eigenfunctions (38) is of course that the local martingales

$$L_t^k=e^{\frac{k^2-1}{2}t}{\varphi }_k\left(Y_t\right)$$

(41)

are bounded martingales and we can thus take expectations to obtain

$$\mathbb{E}\left({\varphi }_k\left(Y_t\right)\right|Y_0=x)=e^{\frac{1-k^2}{2}t}{\varphi }_k\left(x\right)$$

This means the operator $T_t$ defined by

$$T_{t}f\left(x\right)=\mathbb{E}\left(f\left(Y_t\right)\right|X_0=x)$$

(42)

acts on the ${\mathcal{H}}_{\rho }$ space spanned by (38) as a contraction along the diagonal:

$$T_t{\varphi }_k=e^{\frac{1-k^2}{2}t}{\varphi }_k$$

(43)

This enables us to compute the expectation $\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)$ for any ${\mathcal{H}}_{\rho }$ function and also to express the conditional probability density of $Y_t$.

In the following, we use auxiliary theta functions as defined in [13].

$$\begin{array}{ccc}\hfill {\vartheta }_3(u,q)& =& \sum _{n=-\infty }^{\infty }{q}^{n^2}exp\left(2niu\right)=1+2\sum _{n=1}^{\infty }{q}^{n^2}cos\left(2niu\right)\hfill \\ \hfill {\vartheta }_4(u,q)& =& \sum _{n=-\infty }^{\infty }{q}^{n^2}{(-1)}^{n}exp\left(2niu\right)=1+2\sum _{n=1}^{\infty }{q}^{n^2}{(-1)}^{n}cos\left(2niu\right)\hfill \end{array}$$

Proposition 4.

Let

$$f\in {\mathcal{H}}_{\rho }=L_{\rho }^2\left([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R}\right)$$

(44)

where the weight ρ is defined by (40) and let $Y_t$ be the solution to (22) with $Y_0=x\in (-\frac{\pi }{2},\frac{\pi }{2})$. Then, we have

$$\begin{array}{ccc}\hfill \mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)& =& {\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)\rho \left(y\right)\sum _{k=1}^{\infty }{e}^{\frac{1-k^2}{2}t}{\varphi }_k\left(x\right){\varphi }_k\left(y\right)dy\hfill \\ & =& {\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)\rho \left(y\right)\left(1+\sum _{k=2}^{\infty }{e}^{\frac{1-k^2}{2}t}{\varphi }_k\left(x\right){\varphi }_k\left(y\right)\right)dy\hfill \end{array}$$

(45)

where the functions ${\varphi }_k$ are defined by (38).

The Markov kernel above can be expressed as

$$\mathbb{P}(Y_t\in dy|Y_0=x)=\frac{1}{2\pi }\frac{cos\left(y\right)e^{\frac{t}{2}}}{cos\left(x\right)}\left({\vartheta }_3\left(\frac{y-x}{2},e^{-\frac{t}{2}}\right)-{\vartheta }_4\left(\frac{y+x}{2},e^{-\frac{t}{2}}\right)\right)dy$$

(46)

for $y\in [-\frac{\pi }{2},\frac{\pi }{2}]$.

Proof. 

By assumption we have

$$f=\sum _{k=1}^{\infty }f,{\varphi }_k{\rangle }_{\rho }{\varphi }_k$$

Applying the operator (42) and noting (43), which works because of the fact that (41) are bounded martingales, we have

$$\begin{array}{ccc}\hfill \mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)& =& T_t\sum _{k=1}^{\infty }\langle f,{\varphi }_k{\rangle }_{\rho }{\varphi }_k\left(x\right)\hfill \\ & =& \sum _{k=1}^{\infty }{e}^{\frac{1-k^2}{2}t}\langle f,{\varphi }_k{\rangle }_{\rho }{\varphi }_k\left(x\right)\hfill \\ & =& \sum _{k=1}^{\infty }{e}^{\frac{1-k^2}{2}t}\left({\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right){\varphi }_k\left(y\right)\rho \left(y\right)dy\right){\varphi }_k\left(x\right)\hfill \end{array}$$

We can interchange the sum and integral due to $L^2$ convergence which gives (45). Obtaining (46) is just a matter of manipulating the expression for the integration kernel. Here are the details:

$$\begin{array}{ccc}\hfill \rho \left(y\right)\sum _{k=1}^{\infty }{e}^{\frac{1-k^2}{2}t}{\varphi }_k\left(x\right){\varphi }_k\left(y\right)& =& \frac{2}{\pi }{cos}^2\left(y\right)\left(\sum _{n=1}^{\infty }{e}^{\frac{1-{(2n-1)}^2}{2}t}\frac{cos\left(\right(2n-1\left)x\right)}{cos\left(x\right)}\frac{cos\left(\right(2n-1\left)y\right)}{cos\left(y\right)}\right.\hfill \\ & & +\left.\sum _{n=1}^{\infty }{e}^{\frac{1-{\left(2n\right)}^2}{2}t}\frac{sin\left(2nx\right)}{cos\left(x\right)}\frac{sin\left(2ny\right)}{cos\left(y\right)}\right)\hfill \\ & =& \frac{1}{\pi }\frac{cos\left(y\right)e^{\frac{t}{2}}}{cos\left(x\right)}\left(\sum _{n=1}^{\infty }{e}^{\frac{-{(2n-1)}^2}{2}t}\left(cos\left(\right(2n-1\left)\right(y-x\left)\right)+cos\left(\right(2n-1\left)\right(y+x\left)\right)\right)\right.\hfill \\ & & +\left.\sum _{n=1}^{\infty }{e}^{\frac{-{\left(2n\right)}^2}{2}t}\left(cos\left(2n\right(y-x\left)\right)-cos\left(2n\right(y+x\left)\right)\right)\right)\hfill \\ & =& \frac{1}{\pi }\frac{cos\left(y\right)e^{\frac{t}{2}}}{cos\left(x\right)}\sum _{k=1}^{\infty }{e}^{\frac{-k^2}{2}t}\left(cos\left(k(y-x)\right)-{(-1)}^{k}cos\left(k(y+x)\right)\right)\hfill \\ & =& \frac{1}{2\pi }\frac{cos\left(y\right)e^{\frac{t}{2}}}{cos\left(x\right)}\left({\vartheta }_3\left(\frac{y-x}{2},e^{-\frac{t}{2}}\right)-{\vartheta }_4\left(\frac{y+x}{2},e^{-\frac{t}{2}}\right)\right)\hfill \end{array}$$

(47)

□

Of course, wherever there is one expression involving theta functions, there are bound to be many more. We use the Poisson summation formula to derive only one such representation, which can then be related directly to reflected Brownian motion.

Proposition 5.

Let ${\tilde{W}}_t$ be a Brownian motion with ${\tilde{W}}_0=x$ and $Y_t$ the solution to (22) with $Y_0=x$ where $x\in (-\frac{\pi }{2},\frac{\pi }{2})$.

Denote the functions

$$m(y,x,t)=\frac{cos\left(y\right)e^{\frac{t}{2}}}{cos\left(x\right)}$$

and

$$r\left(y\right)=\sum _{k=-\infty }^{\infty }(y-2k\pi )·1_{\{[-\frac{\pi }{2},\frac{\pi }{2}]+2k\pi \}}\left(y\right)+((2k+1)\pi -y)·1_{\{[\frac{\pi }{2},\frac{3\pi }{2}]+2k\pi \}}\left(y\right)=arcsin(sin\left(y\right))$$

This implies

$$M_t=m({\tilde{W}}_t,x,t)$$

is a martingale with $M_0=1$ and that

$$R_t=r({\tilde{W}}_t)$$

represents Brownian motion reflected on the edges of $[-\frac{\pi }{2},\frac{\pi }{2}]$.

For $f\in L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ we have

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\mathbb{E}\left(M_{t}f\left(R_t\right)\right)$$

(48)

Proof. 

We proceed from the expression for the integration kernel (47) to write $\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)$ as

$$\begin{array}{ccc}& =& \frac{1}{\pi }{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)\sum _{k=1}^{\infty }{e}^{\frac{-k^2}{2}t}\left(cos\left(k(y-x)\right)-{(-1)}^{k}cos\left(k(x+y)\right)\right)dy\hfill \\ & =& \frac{1}{2\pi }{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)\sum _{k=-\infty }^{\infty }{e}^{\frac{-k^2}{2}t}\left(exp\left(k(y-x)i\right)-{(-1)}^{k}exp\left(k(x+y)i\right)\right)dy\hfill \end{array}$$

(49)

since the term for $k=0$ in the sum above equals zero. We apply Poisson summation to each of the two sums:

$$\begin{array}{ccc}\hfill \sum _{k=-\infty }^{\infty }g\left(k\right)& =& \sum _{k=-\infty }^{\infty }{e}^{-\frac{k^{2}t}{2}}exp\left(k(y-x)i\right)\hfill \\ \hfill \sum _{k=-\infty }^{\infty }h\left(k\right)& =& \sum _{k=-\infty }^{\infty }{e}^{-\frac{k^{2}t}{2}}{(-1)}^{k}exp\left(k(y+x)i\right)=\sum _{k=-\infty }^{\infty }{e}^{-\frac{k^{2}t}{2}}exp\left(k(y+x+\pi )i\right)\hfill \end{array}$$

Computing the Fourier transforms, we obtain

$$\widehat{g}\left(w\right)={\int }_{-\infty }^{\infty }g\left(u\right)e^{-2\pi iwu}du=\sqrt{\frac{2\pi }{t}}exp\left(-\frac{{(y-x-2\pi w)}^2}{2t}\right)$$

and similarly

$$\widehat{h}\left(w\right)=\sqrt{\frac{2\pi }{t}}exp\left(-\frac{{(y+x+\pi -2\pi w)}^2}{2t}\right)$$

Then (49) can be expressed as

$$\begin{array}{ccc}& =& \frac{1}{\sqrt{2\pi t}}{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)\sum _{k=-\infty }^{\infty }exp\left(-\frac{{(x-y+2\pi k)}^2}{2t}\right)-exp\left(-\frac{{(x+y+\pi +2k\pi )}^2}{2t}\right)dy\hfill \end{array}$$

Noting that $m(y+\pi ,x,t)=-m(y,x,t)$, we can write

$$\begin{array}{ccc}& =& \frac{1}{\sqrt{2\pi t}}\sum _{k=-\infty }^{\infty }{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)e^{-\frac{{(y-x+2\pi k)}^2}{2t}}dy-{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)e^{-\frac{{(-y-x-\pi +2k\pi )}^2}{2t}}dy\hfill \\ & =& \frac{1}{\sqrt{2\pi t}}\sum _{k=-\infty }^{\infty }{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)e^{-\frac{{(y-x+2\pi k)}^2}{2t}}dy+{\int }_{\frac{\pi }{2}}^{\frac{3\pi }{2}}f(y-\pi )m(y,x,t)e^{-\frac{{(-y-x+2k\pi )}^2}{2t}}dy\hfill \end{array}$$

by substituting $y+\pi \mapsto y$ in the second integral. We expand the domain of f to a periodic function with period $\pi$, note that $m(y,x,t)$ is even with period $2\pi$ and substitute $y\mapsto -y$ in the second integral above to obtain

$$\begin{array}{ccc}& =& \frac{1}{\sqrt{2\pi t}}\sum _{k=-\infty }^{\infty }{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)m(y,x,t)e^{-\frac{{(y-x+2\pi k)}^2}{2t}}dy+{\int }_{\frac{\pi }{2}}^{\frac{3\pi }{2}}f(-y)m(y,x,t)e^{-\frac{{(y-x+2k\pi )}^2}{2t}}dy\hfill \\ & =& \frac{1}{\sqrt{2\pi t}}{\int }_{{\bigcup }_k\{[-\frac{\pi }{2},\frac{\pi }{2}]+2k\pi \}}f\left(y\right)m(y,x,t)e^{-\frac{{(y-x)}^2}{2t}}dy+{\int }_{{\bigcup }_k\{[\frac{\pi }{2},\frac{3\pi }{2}]+2k\pi \}}f(-y)m(y,x,t)e^{-\frac{{(y-x)}^2}{2t}}dy\hfill \end{array}$$

Specializing to $f\left(y\right)=1_A\left(y\right)$ for Borel subset $A\subset [-\frac{\pi }{2},\frac{\pi }{2}]$ (so that we have $1_A(-y)=1_{-A}\left(y\right)$) we see that the last expression is what one could write in order to compute

$$\mathbb{E}\left(1_A\left(R_t\right)m(\widetilde{W_t},x,t)\right|\widetilde{W_0}=x)$$

with $R_t$ being the Brownian motion ${\tilde{W}}_t$ reflected at $\pm \pi /2$. □

From (46) (or more obviously from (45)) we can see that our process $Y_t$ even has an elementary limit distribution, at least in a weak sense.

Corollary 5.

For f and $Y_t$, as before we have

$$\mathbb{E}\left(f\left(Y_{\infty }\right)\right|Y_0=x)={\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right)\rho \left(y\right)dy=\frac{2}{\pi }{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}f\left(y\right){cos}^2\left(y\right)dy$$

This implies the limit distribution is independent from the initial value $Y_0=x$.

The results for $Y_t$ from (22) of course carry over to the solutions of our original Equation (6) via the map $X_t=tan\left(Y_t\right)$ (recalling (21)). In particular, it follows that $X_t$ also has a limit distribution.

Corollary 6.

The solution $X_t$ to Equation (6) has a limit distribution in the sense that for f such that $f\circ tan\in L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$

$$\mathbb{E}\left(f\left(X_{\infty }\right)\right|X_0=x)=\frac{2}{\pi }{\int }_{-\infty }^{\infty }\frac{f\left(x\right)}{{(1+x^2)}^2}dx$$

5.4. Double Roots

The remaining case is the double roots Equation (7). Like the case (3) with initial value to the right of the roots, it is well-known that the solution is a supermartingale. Some sample path simulations are shown in Figure 4.

The functions (29) are not bounded when $y\to 0$ (which corresponds to $x\to \infty$ in (26)), but we can form linear combinations that are bounded:

$${\varphi }_Y(y,u)=\frac{1}{2i}\left({\psi }_Y(y,ui)-{\psi }_Y(y,-ui)\right)=\frac{sin\left(uy\right)}{y}$$

We then have that

$$M_t^u=e^{-\lambda \left(ui\right)t}{\varphi }_Y(Y_t,u)=e^{\frac{u^2}{2}t}\frac{sin\left(u{Y}_t\right)}{Y_t}$$

(50)

is a bounded martingale. Taking expectation produces the sine Fourier transform of the transition density.

Proposition 6.

Let $Y_t$ be the solution to (28) with initial condition $Y_0=x>0$. Then,

$$\mathbb{P}\left(Y_t\in dy|Y_0=x\right)=\frac{1}{\sqrt{2\pi t}}\frac{y}{x}\left(e^{-\frac{{(y-x)}^2}{2t}}-e^{-\frac{{(y+x)}^2}{2t}}\right)$$

(51)

for $y\ge 0$.

Proof. 

Since (50) is a martingale, we can take expectations at any fixed time $t<\infty$ conditional on $Y_0=x>0$ to obtain

$$\mathbb{E}\left(\frac{sin\left(u{Y}_t\right)}{Y_t}|Y_0=x\right)=e^{-\frac{u^2}{2}t}\frac{sin\left(ux\right)}{x}$$

Considering the distribution of $Y_t$ is supported on $(0,\infty )$, this can be rewritten as

$${\int }_0^{\infty }\frac{sin\left(uy\right)}{y}\mathbb{P}\left(Y_t\in dy|Y_0=x\right)=e^{-\frac{u^2}{2}t}\frac{sin\left(ux\right)}{x}$$

By computing the inverse sine Fourier transform of the right-hand side and after some elementary manipulations, we can express the transition density in the form (51). □

Again, a little reworking of the expression for $\mathbb{E}\left(f\left(Y_t\right)\right)$ using (51) yields another instance of a ‘negative’ measure change.

Corollary 7.

Let $Y_t$ and $Y_0=x>0$ be as above. Let ${\tilde{W}}_t$ be a Brownian motion with ${\tilde{W}}_0=x$ and define

$$M_t=\frac{{\tilde{W}}_t}{x}$$

which is clearly a martingale with $M_0=1$. Then, for Wiener integrable functions f on $(0,\infty )$, we have

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\mathbb{E}(M_{t}f\left(\right|{\tilde{W}}_t\left|\right))$$

(52)

Proof. 

Using (51), we can write

$$\begin{array}{ccc}\hfill \mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)& =& \frac{1}{\sqrt{2\pi t}}\left({\int }_0^{\infty }f\left(y\right)\frac{y}{x}{e}^{-\frac{{(y-x)}^2}{2t}}dy\right.\hfill \\ & & \left.-{\int }_0^{\infty }f\left(y\right)\frac{y}{x}{e}^{-\frac{{(y+x)}^2}{2t}}dy\right)\hfill \end{array}$$

By extending the definition of f to negative values by $f(-y)=f\left(y\right)$ and substituting $y\mapsto -y$ in the second integral we can rewrite this as

$$\mathbb{E}\left(f\left(Y_t\right)\right|Y_0=x)=\frac{1}{\sqrt{2\pi t}}{\int }_{-\infty }^{\infty }f\left(y\right)\frac{y}{x}{e}^{-\frac{{(y-x)}^2}{2t}}dy$$

which is the same as computing the right-hand side of (52). □

6. Summary of Results and Some Remarks

All the expressions for $Y_t$ obviously carry over to $X_t$. We summarize only the results obtained by Formulas (33), (37), (48) and (52) by substitutions (12), (16), (21) and (27). Explicit expressions for transition kernels can be found by substitution in the integrals.

• For the solution to (3) with initial value (4) and function f on $(-1,1)$ such that $f\circ tanh$ is Wiener integrable, we have

  $$\mathbb{E}\left(f\left(X_t\right)\right|X_0=x)=\mathbb{E}(M_{t}f(tanh({\tilde{W}}_t))|{\tilde{W}}_0=atanh(x))$$

  (53)

  where

  $$M_t=\frac{cosh({\tilde{W}}_t)e^{-\frac{t}{2}}}{cosh({\tilde{W}}_0)}$$
• For the solution to (3) with initial value (5) and function f on $(1,\infty )$ such that $f\circ cotanh$ is Wiener integrable, we have

  $$\mathbb{E}\left(f\left(X_t\right)\right|X_0=x)=\mathbb{E}(M_{t}f(cotanh(|{\tilde{W}}_t\left|\right)\left)\right|{\tilde{W}}_0=acotanh\left(x\right))$$

  (54)

  where

  $$M_t=\frac{sinh({\tilde{W}}_t)e^{-\frac{t}{2}}}{sinh({\tilde{W}}_0)}$$
• For the solution to (6) and f on $\mathbb{R}$ such that $f\circ tan\in L_{\rho }^2([-\frac{\pi }{2},\frac{\pi }{2}],\mathbb{R})$ (with $\rho$ defined by (40)), we have

  $$\mathbb{E}\left(f\left(X_t\right)\right|X_0=x)=\mathbb{E}\left(M_{t}f(tan\left(R_t\right))|{\tilde{W}}_0=arctan\left(x\right)\right)$$

  (55)

  where

  $$M_t=\frac{cos({\tilde{W}}_t)e^{\frac{t}{2}}}{cos({\tilde{W}}_0)}$$

  and $R_t=arcsin(sin({\tilde{W}}_t))$ is Brownian motion reflected on the edges of $[-\pi /2,\pi /2]$.
• For the solution to (7) and function f on $(0,\infty )$ such that $xf(1/x)$ is Wiener integrable, we have

  $$\mathbb{E}\left(f\left(X_t\right)\right|X_0=x)=\mathbb{E}\left(M_{t}f(1/|{\tilde{W}}_t\left|\right)|{\tilde{W}}_0=1/x\right)$$

  (56)

  where

  $$M_t=\frac{{\tilde{W}}_t}{{\tilde{W}}_0}$$

Perhaps the most interesting aspect of the computed weak solutions are their interpretation of martingale measure changes expressed in (37), (48) and (52).

The martingales $M_t$ from each of these formulas all have the the property $M_0=1$ and the results certainly look like usual martingale measure changes of a reflected Brownian motion $R_t$ (reflected at the origin in the cases (37 and (52) or at $\pm \frac{\pi }{2}$ in the case of (48)). But from a probabilistic viewpoint, this should be nonsense since the martingales $M_t$ are all negative much of the time.

However, while the reflected Brownian motion $R_t$ and the associated martingale $M_t$ are both $\sigma \left({\tilde{W}}_t\right)$-measurable, $M_t$ is not $\sigma \left(R_t\right)$-measurable: for fixed t and $\omega \in \Omega$, the weight $M_t\left(\omega \right)$ assigned to the value $R_t\left(\omega \right)$ can be either positive or negative depending the trajectory of ${\tilde{W}}_t\left(\omega \right)$. It all still works out because the negative weights have a smaller probability than positive weights, even for a fixed value of $R_t$.

So, although (37), (52) and (48) cannot be regarded as a changes of measure in the usual sense, it shows it might be possible to broaden the concept of a change of measure to include the possibility of weights attaining negative values (or even complex values) without compromising fundamental probabilistic properties of the resulting semigroup, as long as the weights are not simply functions of the process undergoing the measure-change.

In our specific cases, the expressions for these measure changes were obtained by elementary means via simple algebraic manipulations of the transition kernels, as can be carried out to obtain special cases of the Girsanov transform for example, but it might be interesting to see if one could develop a more general framework of possibly negative martingale measure changes.

The fact that functionals of the type $\mathbb{E}\left(f\left(X_t\right)\right|X_0=x)$ can be expressed using these measure changes also implies that expressions for more general functionals of paths can also be obtained in the same way and for the same reasons as for familiar Girsanov-type measure changes. For instance, for a functional of the path ${\left\{X_s\right\}}_{t\in [0,t]}$ of the type

$${\int }_0^{t}g\left(X_s\right)ds$$

where g is a Riemann integrable function, one can approximate this integral with a finite sum and express the expectation of this sum as

$$\begin{array}{ccc}\hfill \mathbb{E}(\sum _{i}g(X_{t_i})\Delta t_i)& =& \sum _i\mathbb{E}(g(X_{t_i}))\Delta t_i\hfill \\ & =& \sum _i\mathbb{E}(M_{t_i}g(\widetilde{W_{t_i}}))\Delta t_i\hfill \\ & =& \sum _i\mathbb{E}(\mathbb{E}(M_t|M_{t_i})g(\widetilde{W_{t_i}}))\Delta t_i\hfill \\ & =& \sum _i\mathbb{E}(M_{t}g(\widetilde{W_{t_i}}))\Delta t_i\hfill \\ & =& \mathbb{E}(M_t\sum _{i}g(\widetilde{W_{t_i}})\Delta t_i)\hfill \end{array}$$

which implies $\mathbb{E}({\int }_0^{t}g(X_s)ds)=\mathbb{E}(M_t{\int }_0^{t}g({\tilde{W}}_s)ds)$ when these expectations exist. The point is that in order to write expectations for functionals of the paths ${\left\{X_s\right\}}_{t\in [0,t]}$ using a measure change, the assumption that $M_t$ is non-negative is not needed, only that $M_t$ is a martingale is required.

It seems worth mentioning that a heuristic derivation of the ‘measure changes’ (37), (52) and (48), or at least something consistent with these formulas, is also described in [9] and uses Doobs h-transform. We give only a brief description.

With the semigroup of Brownian motion $P_t$ we associate the generator

$$\mathcal{L}=\underset{t\to 0}{lim}\frac{d}{dt}{P}_t=\frac{1}{2}\frac{{\partial }^2}{\partial x^2}$$

as usual. Take h to be an eigenfunction of $\mathcal{L}$

$$\mathcal{L}h=\lambda h$$

and define a transform of $P_t$ by

$$P_t^{h}f\left(x\right)=\frac{e^{-\lambda t}}{h\left(x\right)}\mathbb{E}\left(h\left(W_t\right)f\left(W_t\right)\right|W_0=x)$$

We can then explicitly compute the expression for the transformed generator for a given eigenfunction h

$${\mathcal{L}}^{h}f=\underset{t\to 0}{lim}\frac{d}{dt}{P}_t^{h}f=\underset{t\to 0}{lim}\frac{d}{dt}\left(\frac{e^{-\lambda t}}{h}{P}_t\left(hf\right)\right)$$

For instance, the choice $h\left(x\right)=cosh\left(x\right)$ yields

$${\mathcal{L}}^h=\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+tanh\left(x\right)\frac{\partial }{\partial x}$$

which is the generator for $Y_t$ from (14). This gives a link between the generator of Brownian motion $\mathcal{L}$ and (14), from which we can deduce the expression for the weak solution to (13) as the measure change (33).

On the other hand, the choice $h\left(x\right)=cos\left(x\right)$ gives

$${\mathcal{L}}^h=\frac{1}{2}\frac{{\partial }^2}{\partial x^2}-tan\left(x\right)\frac{\partial }{\partial x}$$

which is the generator for $Y_t$ from (23). Again, this gives us a direct link between the generator of Brownian motion and the generator (23) via the ‘change of measure’

$$\frac{d\mathbb{Q}}{d\mathbb{P}}=M_t=\frac{cos\left(W_t\right)e^{t/2}}{cos\left(x\right)}$$

One objection to this can be raised in the name of principle since $M_t$ takes negative values, which is strange for a probability measure change (although as mentioned this can actually be meaningful). A more practical problem in this case is that although we obtain the correct generator $\mathcal{L}$ associated with (23), without knowing its domain this does not uniquely determine the process generated by $\mathcal{L}$. Namely, the generator of Brownian motion $\mathcal{L}$ also happens to be the generator of reflected Brownian motion (and also of $W_t mod 1$ for instance). Equation (48) reveals that (23) is actually related to reflected Brownian motion. But when working with transformations only on the level of generators, it is not necessarily immediately clear what the underlying Markov semigroups are.

Another point that becomes apparent due to Proposition 4, is that the SDE (6) actually has a limit distribution (Corollary 6). This can even be explained intuitively if we recall that any local martingale $X_t$ can be expressed as a time-change of some Brownian motion ${\tilde{W}}_t$, i.e., $X_t={\tilde{W}}_{\tau }$, where the time change process $\tau$ is related to the quadratic variation $\langle X_t\rangle$. In our case, the time increment can be expressed as $d\tau =d\langle X_t\rangle =(1+X_t^2)dt$. This means $X_t={\tilde{W}}_{\tau }$ moves at an accelerated rate the further it is from 0. But since Brownian motion is recurrent, this means it will also return to 0 more quickly and will thus spend more time in its vicinity where it moves more slowly. The fact that it has a limit distribution (unlike standard Brownian motion) in simply due to the fact it spends much of the time near 0.

7. Results of Numerical Simulations

In Figure 5, Figure 6, Figure 7 and Figure 8, we present some graphs, which are the results of numerical simulations for each of our four SDEs (two graphs for each of the four equations).

The blue lines in the graphs show the theoretical probability densities of the variable $X_t$ at the specified times t and initial values $X_0$, which are computed according to the formulas in Section 5.

The histograms show the results of numerical simulations of sample paths of $X_t$ taken at the same specified times and initial value (examples of such simulated sample paths are as shown in the illustrations in Section 5). Each histogram is the result for $n={10}^4$ simulations. These histograms match the theoretical values to around the expected $\approx 1\%$ error for a sample size of $n={10}^4$. This clearly demonstrates the accuracy of our results.

The red dots represent the Monte-Carlo evaluation of $\mathbb{E}\left(f\left(X_t\right)\right)$ for a range of indicator functions $f\left(x\right)=1_{[x_i,x_{i+1}]}\left(x\right)$ of given small intervals $[x_i,x_{i+1}]$ using the Formulas (53), (54), (55) and (56) from Section 6. For each red dot in each of the graphs, we were able to take a sample size of $n={10}^8$, since instead of simulating the entire sample path $X_t$ for each simulation (as is the case with the histograms), it is enough to choose just one normal random number that simulates the value of ${\tilde{W}}_t$, compute $M_{t}f({\tilde{W}}_t)$ according to the Formulas (53)–(56), and take the average of these values to approximate the expected value $\mathbb{E}(M_{t}f({\tilde{W}}_t))$. All the dots agree with the theoretical values up to an error, which is expected to be $\approx 0.01\%$ for a sample size of $n={10}^8$. This additionally confirms the correctness of the results and also demonstrates the efficiency of the proposed Monte Carlo schemes for evaluation of $\mathbb{E}\left(f\left(X_t\right)\right)$ without the need to simulate the trajectories of $X_t$.

8. Conclusions

We have shown that the described method provides a successful and rigorous way of treating QNV models that apparently also clarifies some confusing aspects and results in the existing literature.

From a theoretical viewpoint, this ’spectral method’, i.e., the blending of Itô calculus, classic ODE eigensystem analysis and Fourier analysis to solve an SDE, is of course far from completely original. Indeed, Fourier and related transforms are of fundamental importance throughout probability in general, including Itô calculus. These connections, however, do not appear to have been used much in this specific way for SDEs, most likely due to the fact that there are other ways of treating SDEs, which are often based on more general and sophisticated theory. And although a more sophisticated theory clearly may offer advantages, it may also be an advantage if it is possible to keep things simple. The more traditional martingale theory treatment avoids technical and practical issues of convergence and integrability sometimes present in the treatment of SDEs by PDE or Markov process methods, by using localisation techniques. Here, we avoid the pitfalls at a basic level by exploiting the direct link between SDEs and their ’generators’ (without referring to Markov semigroups) provided by the Itô formula, and then by choosing an expedient complete eigensystem of bounded functions.

This works well for QNV models due to their special analytic properties, which are discussed in the literature (for instance in [9]), and is not as widely applicable as the mentioned traditional and general martingale theory methods. But this method should be applicable to at least Sturm–Liouville-type SDEs, which include a wide range of classic diffusions (like Bessel, hypergeometric, Jacobi diffusions and so on). Such classic SDEs have of course been widely studied by various means and a ’spectral method’ may not necessarily produce significant new results, but we believe it can be valuable, particularly in instances where the strict local martingale or semimartingale nature of the SDEs must be accounted for and can produce exact quantitative results where applicable, as in the case of QNV models.

Another theoretical aspect that may be worth exploring is the question of whether a general framework for non-strictly positive martingale measure changes, which spontaneously appeared in our treatment, can be developed, as touched upon in Section 6.

From a more application-oriented viewpoint, we have shown that the derived expectation formulas provide for efficient numerical and also Monte Carlo methods for evaluating expectations that can be applied to for instance numerically efficient pricing of arbitrary derivatives based on QNV modelled assets. These results can be useful in other fields where a stochastic process of interest can be modelled by, or transformed to, a QNV type SDE. We believe that QNV models and related SDEs may find more applications, in particular if the fact that these SDEs are essentially solved was more widely known outside mathematical finance.

We also think that similarly efficient exact and numerical ways of computing quantities associated to other SDEs could be provided by such a method where it is applicable.