# Finite Time Explosion of Stochastic Differential Equations: A survey into
Khasminskii’s Lyapunov Method and its Consistency with the Osgood Criterion

Seungsoo Lee

###### Abstract

Solutions of Stochastic Differential Equations can have three types of
explosive behaviors: almost-sure non-explosive, explosion with positive
probability, and almost sure explosion. In this paper, we will provide a
survey of Khasminskii’s Lyapunov method for classifying explosive behaviors of
solutions of stochastic differential equations. We will embark our expedition
by examining the renowned Feller’s test for explosion and observing its
shortfalls. Afterwards, we will present Khasminskii’s Lyapunov method for
almost-sure non-explosion, explosion with positive probability, and almost-
sure explosion. Ample examples will be provided to illuminate the power of
Khasminskii’s Lyapunov methods. Furthermore, quick layovers will be made to
extend Khasminskii’s Lyapunov method for almost-sure non-explosion and
explosion with positive probability for jump processes with constant Poisson
intensities.

## 1 Background

In this paper, we will go over the established methods to analyze explosive
behaviors for stochastic differential equations and their advantages and
disadvantages. Namely, we will look at the renowned Feller’s test for
explosion, and the lesser known Lyapunov function method. First show the
explosive behavior for the deterministic ordinary differential equations (ODE)
to provide a motivation for studying the explosive behavior of stochastic
differential equations (SDE).

It has already been shown in [8] that finite time explosion in ODE has already
been fully characterized by the Osgood condition which states that for an ODE
of the form $\frac{dx}{dt}=b(x(t))>0$ with $x(0)=\xi>0$ and $t>0$ explodes if
and only if

$\int_{\xi}^{\infty}\frac{1}{b(s)}ds<\infty$

However, when our differential equation is perturbed with noise, our
characterization becomes more complicated as the solution path can explode
with positive probability, probability 1, and not explode with probability 1.
It has already been shown in [5] that when we have noise with constant
intensity, that the Osgood condition of the drift parameter dictates the
almost sure explosion and almost sure non-explosion of solution paths. Yet, we
can see that the Osgood condition is no longer applicable when we have noise
with high powers as shown in Example 1.2, where the drift indicates explosion,
but the addition of large noise prevents explosion. Thus, further analysis is
required to characterize the explosive behavior of SDEs.

Throughout the paper, we shall define a solution of our SDE to be explosive
when the hitting time for our solution, $X_{t}$, is finite,
$\tau_{\infty}<\infty$, and non-explosive if $\tau_{\infty}=\infty$. Thus,
there are three possible cases involving explosions: positive probability of
explosion $\mathbb{P}(\tau_{\infty}<\infty)>0$, almost sure non-explosion
$\mathbb{P}(\tau_{\infty}=\infty)=1$, and almost sure explosion
$\mathbb{P}(\tau_{\infty}<\infty)=1$. Of the three methods for testing
explosion, only the Lyapunov function method will fully explore all three
types of explosions.

One prominent test for finite time explosion is Feller’s test for explosion as
shown by Karatzas and Shreve [2]

######  Theorem 1.1.

Let our SDE be of the form

$dX_{t}=b(X_{t})dt+\sigma(X_{t})dW_{t}\>(pg.329).$

Assume we have $\sigma^{2}(x)>0\>\forall x\in\mathbb{R}$ (non-degeneracy), and
$\forall x\in\mathbb{R}\exists\epsilon>0$ st.
$\int_{x-\epsilon}^{x+\epsilon}\frac{\lvert b(y)\rvert
dy}{\sigma^{2}(y)}<\infty$ (local integrability). Now, let the scale function
$p(x)$ be defined as:

$p(x):=\int_{c}^{x}\exp{\left(-2\int_{c}^{\xi}\frac{b(\zeta)d\zeta}{\sigma^{2}(\zeta)}\right)}d\xi\>\>\>(pg.339).$

Now, we define
$v(x):=\int^{x}_{c}p^{\prime}(y)\int_{c}^{y}\frac{2dz}{p^{\prime}(z)\sigma^{2}(z)}dy\>\>(pg.347).$
Now let $I:=(l,r)$ denote the interval of interest and $\tau$ be the hitting
time for our end points $l$ or $r$. Then, Feller’s Test states that when
$v(l^{+})=v(r^{-})=\infty$, $\mathbb{P}(\tau=\infty)=1$ and when
$v(l^{+})=v(r^{-})=\infty$ is not true, $\mathbb{P}(\tau=\infty)<1$.

Now, let us show two examples of applying Feller’s test to show explosion of a
process and non-explosion of a process. We will start with an explosive
example with constant intensity noise.

######  Example 1.2.

Explosiveness of Constant Intensity Noise. Consider the following SDE

$dX_{t}=X_{t}^{2}dt+\sigma dW_{t}$

where $\sigma\in\mathbb{R}$ is a constant. It has been shown in Leon-Villa’s
Theorem 3.1, that such process explodes in finite time if and only if the
Osgood criterion ($\int_{a}^{\infty}\frac{1}{b(t)}dt<\infty$) is satisfied
[5]. Since $\int_{a}^{\infty}\frac{1}{t^{2}}dt<\infty$, we expect Feller’s
test to state $\mathbb{P}(\tau=\infty)<1$ or at least one of
$v(\infty),\>v(-\infty)$ to be bounded. Let $\l^{+}=-\infty$ and
$r^{-}=\infty$, so we have $I=(-\infty,\infty)=\mathbb{R}$. Note our
$p(x)=\int_{c}^{x}\exp{(-2\int_{c}^{\xi}\frac{\zeta^{2}d\zeta}{\sigma^{2}})}d\xi$
and so we know our $v(x)$ is of the form

$v(x)=\int^{x}_{c}e^{-\frac{2}{3\sigma^{2}}(y^{3}-c^{3})}\left(\frac{2}{\sigma^{2}}\int_{c}^{y}e^{\frac{2}{3\sigma^{2}}(\zeta^{3}-c^{3})}d\zeta\right)dy$
$=\int^{x}_{c}e^{-\frac{2}{3\sigma^{2}}(y^{3})}\left(\frac{2}{\sigma^{2}}\int_{c}^{y}e^{\frac{2}{3\sigma^{2}}(\zeta^{3})}d\zeta\right)dy$

Now, when we have $v(\infty)$, we have

$\displaystyle v(\infty)$
$\displaystyle=\int^{\infty}_{c}e^{-\frac{2}{3\sigma^{2}}(y^{3})}\left(\frac{2}{\sigma^{2}}\int_{c}^{y}e^{\frac{2}{3\sigma^{2}}(\zeta^{3})}d\zeta\right)dy$

Let us first address the integral inside the integral. Note that we can
establish the following inequality for sufficiently large $y$

$\int_{c}^{y}e^{\zeta^{3}}d\zeta=\int_{c}^{y-1}e^{\zeta^{3}}d\zeta+\int_{y-1}^{y}\frac{3\zeta^{2}}{3\zeta^{2}}e^{\zeta^{3}}d\zeta$

Now, we apply integration by parts by letting $u=\frac{1}{\zeta^{2}}$ and
$dv=3\zeta^{2}e^{\zeta^{3}}d\zeta$. Thus, we have

$\displaystyle\int_{c}^{y}\frac{3\zeta^{2}e^{\zeta^{3}}}{3\zeta^{2}}d\zeta$
$\displaystyle=\frac{e^{y^{3}}}{y^{2}}-\frac{e^{c^{3}}}{c^{2}}+\int_{c}^{y}\frac{2e^{\zeta^{3}}}{\zeta^{3}}d\zeta$
$\displaystyle\leq k\frac{e^{y^{3}}}{y^{2}}$

Thus we can state by reparameterizing
$\tilde{y}^{3}=\frac{2}{3\sigma^{2}}y^{3}$

$\displaystyle v(\infty)$
$\displaystyle=\int^{\infty}_{c}e^{-\frac{2}{3\sigma^{2}}(y^{3})}\left(\frac{2}{\sigma^{2}}\int_{c}^{y}e^{\frac{2}{3\sigma^{2}}(\zeta^{3})}d\zeta\right)dy$
$\displaystyle\leq\int^{\infty}_{c}e^{-\tilde{y}^{3}}k^{\prime}\frac{e^{\tilde{y}^{3}}}{\tilde{y}^{2}}d\tilde{y}$
$\displaystyle=\int_{c}^{\infty}\frac{k}{\tilde{y}^{2}}d\tilde{y}\>\>\>\><\infty$
(1)

Thus, since $v(\infty)<\infty$ we know from Feller’s test that
$\mathbb{P}(\tau_{\infty}<\infty)>0$ as desired.

######  Example 1.3.

Non-explosiveness of High Power Noise. Consider the following SDE

$dX_{t}=X_{t}^{2}dt+X_{t}^{5}dW_{t}.$

Then we know
$p(x)=\int_{c}^{x}\exp{\left(-2\int_{c}^{\xi}\frac{\zeta^{2}}{\zeta^{10}}\right)d\zeta}d\xi$.
Then we know

$v(x)=k\int_{c}^{x}e^{k_{2}y^{-7}}\int_{c}^{y}\zeta^{-10}e^{-k_{2}\zeta^{-7}}d\zeta
dy$

We want to establish a lower estimate for $v(x)$ and show that they go to
infinity as $x\rightarrow\pm\infty$. Let $c>1$. Then we know

$\int_{c}^{y}\zeta^{-10}e^{-k_{2}\zeta^{-7}}d\zeta\geq\int_{c}^{y}\zeta^{-15}e^{-k_{2}\zeta^{-7}}d\zeta$

Now we apply u-substitution for the integral on the right with $u=x^{-7}$ to
attain

$\int_{c}^{y}\zeta^{-15}e^{-k_{2}\zeta^{-7}}d\zeta=-\frac{1}{7k_{2}^{2}}\int
e^{u}udu$

Now we apply integration by parts and substitute to find

$\int_{c}^{y}\zeta^{-15}e^{-k_{2}\zeta^{-7}}d\zeta=-\frac{-k_{2}c^{7}e^{-\frac{k_{2}}{y^{7}}}-c^{7}y^{7}e^{-\frac{k_{2}}{y^{7}}}+k_{2}y^{7}e^{-\frac{k_{2}}{c^{7}}}+c^{7}y^{7}e^{-\frac{k_{2}}{c^{7}}}}{7k_{2}^{2}c^{7}y^{7}}$
$=\frac{(c^{7}y^{7}+k_{2}c^{7})e^{-k_{2}y^{-7}}}{7k_{2}^{2}c^{7}y^{7}}+k_{3}$

Thus, we know

$\displaystyle v(x)$ $\displaystyle\geq
k\int_{c}^{x}e^{k_{2}y^{-7}}\int_{c}^{y}\zeta^{-15}e^{-k_{2}\zeta^{-7}}d\zeta
dy$
$\displaystyle=k\int_{c}^{x}e^{k_{2}y^{-7}}\left(\frac{(c^{7}y^{7}+k_{2}c^{7})e^{-k_{2}y^{-7}}}{7k_{2}^{2}c^{7}y^{7}}+k_{3}\right)dy$

Thus as $x\rightarrow\pm\infty$ we can see that we are integrating an
approximate constant over infinite length and thus
$v(\infty)=v(-\infty)=\infty$ as desired.

######  Remark 1.4.

One can see that by the Osgood criterion, the ODE $dX=X^{2}dt$ explodes as
$\int\frac{1}{X^{2}}dX<\infty$. When a small amount of noise is added through
$dX_{t}=X_{t}^{2}+\sigma dW_{t}$, we observe explosion. However, when the
noise is sufficiently large enough, like $dX_{t}=X_{t}^{2}+X_{t}^{5}dW_{t}$,
we see that our SDE does not explode. In fact, as $X_{t}>>0$, we see that the
diffusion component becomes more significant than the force component.
However, for dimension less than 2, the noise component cannot cause explosion
as a drunk man will always find his way home while a drunk bird may not.

######  Remark 1.5.

One benefit of Feller’s test is that one can apply it in a systematic manner.
However, one downside of Feller’s test is that it is difficult to extend to
higher dimensions. Note that local integrability, which is used to construct
$p(x)$, provides the following constraint $\frac{|b(x)|}{\sigma^{2}(x)}$.
Thus, in higher dimensions, where $b(x)$ is a vector and $\sigma(x)$ is a
matrix, it is not obvious how one can extend this test. Furthermore, Feller’s
test does not explore all three cases of explosion as it provides no
information about almost sure explosion.

Another method to test for positive probability of explosion is through the
use of Lyapunov function as described by Chow and Khasminskii [1]. Unlike
Feller’s test for explosion, which was methodical but involved constructing
multiple integral functions that could be difficult to compute and was limited
to 1-Dimensional processes, the Lyapunov method involves finding a function
$V$ that satisfies certain criterion involving its end behavior and generator,
which is computed through its partial derivatives, which are often easier to
evaluate than integral functions.

## 2 Almost Sure Non-explosion of SDE through Lyapunov Functions

In the following proofs for Khasminskii’s theorems using Lyapunov functions,
we use the concept of generators to establish a connection between the
expectation and the initial condition as used in Dynkin’s formula [6]. Thus,
in this section we establish how to find the generator of an SDE with and
without jumps for the sake of completeness. For a stochastic process
$X_{t}\in\mathbb{R}^{d}$, we define its generator $L$ acting on $f$ as

$Lf=\lim_{h\downarrow
0}\frac{1}{h}\mathbb{E}(f(X_{t+h})-f(X_{t})|X_{t})\cite[cite]{[\@@bibref{}{oksendal2019stochastic}{}{}]}$

Now, for the first case without jumps, let $X_{t}$ be a continuous process
that is a solution of $dX_{t}=b(t,X_{t})dt+\sigma(t,X_{t})dW_{t}$ and
$f\in\mathcal{C}^{1,2}$. From Ito’s lemma, we know

$df(X_{t})=\left(b(t,X_{t})\frac{\partial f}{\partial
X_{t}}+\frac{1}{2}\sigma(t,X_{t})^{2}\frac{\partial^{2}f}{\partial
X_{t}^{2}}\right)dt+\sigma(t,X_{t})\frac{\partial f}{\partial X_{t}}dW_{t}$

Thus, it follows that

$f(X_{t})-f(X_{0})=\int_{0}^{t}b(u,X_{u})\frac{\partial f}{\partial
X_{u}}+\frac{1}{2}\sigma(u,X_{u})^{2}\frac{\partial^{2}f}{\partial
X_{u}^{2}}du+\int_{0}^{t}\sigma(u,X_{u})\frac{\partial f}{\partial
X_{u}}dW_{u}$

Now, by taking the expectation and recalling that Ito integrals are
martingales, we observe

$\displaystyle\lim_{t\downarrow 0}\frac{1}{t}\mathbb{E}(f(X_{t})-f(X_{0}))$
$\displaystyle=\lim_{t\downarrow
0}\frac{1}{t}\mathbb{E}\left(\int_{0}^{t}b(u,X_{u})\frac{\partial f}{\partial
X_{u}}+\frac{1}{2}\sigma(u,X_{u})^{2}\frac{\partial^{2}f}{\partial
X_{u}^{2}}du\right)$ $\displaystyle=\lim_{t\downarrow
0}\frac{1}{t}\int_{0}^{t}\mathbb{E}\left(b(u,X_{u})\frac{\partial f}{\partial
X_{u}}+\frac{1}{2}\sigma(u,X_{u})^{2}\frac{\partial^{2}f}{\partial
X_{u}^{2}}\right)du$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})$

Now, for the case with jumps, let $X_{t}$ be a solution of the SDE of the form
$dX_{t}=b(t,X_{t})dt+\sigma(t,X_{t})dW_{t}+J(t,X_{t})dN_{t}$. Note from Ito’s
lemma, we know

$df(X_{t})=\left(b(t,X_{t})\frac{\partial f}{\partial
X_{t}}+\frac{1}{2}\sigma(t,X_{t})^{2}\frac{\partial^{2}f}{\partial
X_{t}^{2}}\right)dt+\sigma(t,X_{t})\frac{\partial f}{\partial
X_{t}}dW_{t}+(f(X_{t^{-}}+Y_{t})-f(X_{t^{-}}))dN_{t}$

Thus, it follows from the integral form that

$f(X_{t})-f(X_{0})=\int_{0}^{t}\left(b(t,X_{t})\frac{\partial f}{\partial
X_{u}}+\frac{1}{2}\sigma(u,X_{u})^{2}\frac{\partial^{2}f}{\partial
X_{u}^{2}}\right)du+\int_{0}^{t}\sigma(u,X_{u})\frac{\partial f}{\partial
X_{u}}dW_{u}$ $+\int_{0}^{t}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}}))dN_{s}$

Now, note that by taking the expectation, we have

$\displaystyle\lim_{t\downarrow 0}\frac{1}{t}\mathbb{E}(f(X_{t})-f(X_{0}))$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})+\mathbb{E}(\int_{0}^{t}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}}))dN_{s})$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})+\mathbb{E}(\sum^{N_{t}}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}})))$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})+\mathbb{E}^{N_{t}}(\mathbb{E}(\sum^{N_{t}}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}}))|N_{t}))$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})+\mathbb{E}^{N_{t}}(N_{t}\mathbb{E}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}})|N_{t}))$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})+\mathbb{E}^{N_{t}}(N_{t}\mathbb{E}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}})))$
$\displaystyle=b(t,X_{t})f^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}f^{\prime\prime}(X_{t})+\lambda\mathbb{E}^{Y_{t}}(f(X_{s^{-}}+Y_{s})-f(X_{s^{-}})))$

Thus, when $f$ is a polynomial, the jump component of the generator becomes
the moments for $Y$, with respect to the jump size distribution.

######  Remark 2.1.

Note that when we have an $f$ depending explicitly on time, $f(t,X_{t})$, we
have a $\frac{\partial f}{\partial t}$ term in our force from our Ito formula.
In the proofs for the stability theorems, we will choose our Lyapunov
functions to be of the form $f(X_{t})=V(X_{t})$, but we will multiply it by
$e^{-C(t-t_{0})}$ so we will end up with a
$f(t,X_{t})=V(X_{t})e^{-C(t-t_{0})}$ and utilize the $\frac{\partial}{\partial
t}$ present in the Ito formula.

Now, we will show how Khasminskii’s Lyapunov function can be used to show
almost sure non-explosion $\mathbb{P}(\tau_{\infty}=\infty)=1$

######  Theorem 2.2.

Khasminskii’s theorem 3.5 for almost sure non-explosion in [3]. Let
$\mathcal{D}\subset\mathbb{R}^{d}$ be a bounded open set with regular boundary
$\partial\mathcal{D}$ and let $\mathcal{D}^{c}$ denote the complement of
$\mathcal{D}$. Suppose our SDE parameters satisfy the Lipschitz condition on
any compact subset of $\mathbb{R}^{d}$ for any $t\geq t_{0}$. Now, assume
there exist a positive function $V(x)$ belonging to the class
$\mathcal{C}^{2}$ and positive constants $C$ such that the following
conditions hold

$LV\leq CV$ $V=\inf_{|x|>R}V(x)\rightarrow\infty\>\text{as
}R\rightarrow\infty$

Then it follows that $\mathbb{P}(\tau_{\infty}<\infty)=0$.

###### Proof.

Let $\tau_{R}:=\tau_{R}\wedge t$. Note, with the chain rule argument, we see

$\displaystyle L\\{V(x)e^{c(t-t_{0})}\\}$
$\displaystyle=e^{-c(t-t_{0})}LV(x)-CV(x)e^{-c(t-t_{0})}$
$\displaystyle=e^{-c(t-t_{0})}(LV(x)-CV(x))$ $\displaystyle\leq
0\>\>\>\text{since $e^{-c(t-t_{0})}>0$ and $LV(x)\leq CV(x)$ for
$x\in\mathcal{D}^{c}$}$

Now note from Dynkin’s formula that we have

$\mathbb{E}(V(X_{\tau_{R})}e^{-C(\tau_{R}-t_{0})})=V(X_{t_{0}})+\mathbb{E}\left(\int_{0}^{\tau_{R}}LV(x)e^{-C(s-t_{0})}ds\right)$
$\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\leq\mathbb{E}(V(X_{\tau_{R}}))e^{-C(t-t_{0})}\leq
V(X_{t_{0}})$ $\mathbb{E}(V(X_{\tau_{R}}))e^{-C(t-t_{0})}\leq V(X_{t_{0}})$
$\mathbb{E}(V(X_{\tau_{R}}))\leq e^{C(t-t_{0})}V(X_{t_{0}})$

Now, let us split our expectation into two cases where $t\leq\tau_{R}$ and
$t>\tau_{R}$. When $t>\tau_{R}$, we know from definition of hitting time that
$X_{t}$ has reached $D^{c}$. Thus, our expectation on the left hand side
becomes

$\mathbb{E}(V(X_{\tau_{R}}))=\mathbb{E}(V(X_{\tau_{R}}))\mathbbm{1}_{\\{\tau_{R}<t\\}}+\mathbb{E}(V(X_{\tau_{R}}))\mathbbm{1}_{\\{t<\tau_{R}\\}}$

Now note that $\tau_{R}<t$, we have $V(X_{\tau_{R}})\geq\inf_{|x|>R}V(x)$ by
definition of hitting time. For $t<\tau_{R}$, we know $V\geq 0$ from our
positivity assumption. Thus, we have the following string of inequality

$V(X_{0})e^{C(t-t_{0})}\geq\mathbb{E}(V({X_{\tau_{R}}}))\geq\inf_{|x|>R}V(x)\mathbb{P}(\tau_{R}\leq
t)+0\mathbb{P}(\tau_{R}>t)$

Thus, it follows that

$\frac{V(X_{0})e^{C(t-t_{0})}}{\inf_{|x|>R}V(x)}\geq\mathbb{P}(\tau_{R}\leq
t)$

Thus, as $R\rightarrow\infty$, we know from our condition $V\rightarrow\infty$
and thus $\mathbb{P}(\tau_{R}\leq t)=0$ as desired. ∎

######  Remark 2.3.

In [1], Khasminskii and Chow states how the following condition serves as a
sufficient condition for non-explosion for a $d-$dimensional SDE with
$x\in\mathbb{R}^{d}$ and
$d\mathbb{X}_{t}=b(\mathbb{X}_{t},t)dt+\sigma(\mathbb{X}_{t},t)d\mathbb{W}_{t}$
with $A(t,x):=\sigma(\mathbb{X}_{t},t)\sigma^{T}(\mathbb{X}_{t},t)$:

$\langle b(t,x),x\rangle+\frac{1}{2}tr(A(t,x))-\lambda_{\max}(A(t,x))\leq
C\lVert x\rVert^{2}\ln(\lVert x\rVert^{2})$

Although Khasminskii indicates that the condition was derived based on
$V(x)=\ln(\lVert x\rVert^{2})$, the actual gritty calculation is omitted. We
will show the calculation in this section. From [6], it is well known that the
multi variable generator is of the form

$LV=\sum_{i=1}^{d}b_{i}(t,x)\frac{\partial V}{\partial
x_{i}}+\frac{1}{2}\sum_{i,j=1}^{d}a_{i,j}(t,x)\frac{\partial^{2}V}{\partial
x_{i}\partial x_{j}}$

Let $V(x)=\ln(\lVert x\rVert^{2})=\ln(x_{1}^{2}+\dots+x_{d}^{2})$. Thus, it
follows

$\frac{\partial V}{\partial x_{i}}=\frac{2x_{i}}{\lVert x\rVert^{2}}$
$\frac{\partial^{2}V}{\partial x_{i}\partial x_{j}}=\frac{-2x_{i}x_{j}}{\lVert
x\rVert^{4}}\>\>\>\>\text{for $i\neq j$}$ $\frac{\partial^{2}V}{\partial
x_{i}^{2}}=\frac{-4x_{i}^{2}}{\lVert x\rVert^{4}}+\frac{2}{\lVert
x\rVert^{2}}$

Thus, we have

$\displaystyle LV$ $\displaystyle=\frac{2\langle b(t,x),x\rangle}{\lVert
x\rVert^{2}}+\frac{1}{2}\sum_{i,j=1}^{d}a_{i,j}(t,x)\frac{-4x_{i}x_{j}}{\lVert
x\rVert^{4}}+\frac{1}{2}\sum_{i=1}^{d}a_{ii}(t,x)\frac{2}{\lVert x\rVert^{2}}$
$\displaystyle=\frac{2\langle b(t,x),x\rangle}{\lVert
x\rVert^{2}}-2\frac{\langle Ax,x\rangle}{\lVert
x\rVert^{4}}+\frac{\text{Tr}(A)}{\lVert x\rVert^{2}}$

Now, note that $A=\sigma^{T}(x,t)\sigma(x,t)$ is a positive symmetric definite
matrix. Thus, we know

$0\leq\lambda_{\text{min}}\leq\dots\leq\lambda_{\text{max}}$

It then follows that

$\lambda_{\text{min}}\lVert x\rVert^{2}\leq\langle
Ax,x\rangle\leq\lambda_{\text{max}}\lVert x\rVert^{2}$

Thus, since $\langle Ax,x\rangle\leq\lambda_{\text{max}}\lVert x\rVert^{2}$,
we know

$\displaystyle LV$ $\displaystyle=\frac{2\langle b(t,x),x\rangle}{\lVert
x\rVert^{2}}-2\frac{\langle Ax,x\rangle}{\lVert
x\rVert^{4}}+\frac{\text{Tr}(A)}{\lVert x\rVert^{2}}$
$\displaystyle\leq\frac{2\langle b(t,x),x\rangle}{\lVert
x\rVert^{2}}-2\frac{\lambda_{\text{max}}}{\lVert
x\rVert^{2}}+\frac{\text{Tr}(A)}{\lVert x\rVert^{2}}$

Thus, when $\langle
b(t,x),x\rangle+\frac{1}{2}tr(A(t,x))-\lambda_{\max}(A(t,x))\leq C\lVert
x\rVert^{2}\ln(\lVert x\rVert^{2})$ holds, it follows that $LV\leq CV$ as
desired.

######  Example 2.4.

Consider the following one dimensional SDE $dX_{t}=X_{t}^{p}dt+\sigma dW_{t}$
where $\sigma\in\mathbb{R}$ is a constant. When $0<p<1$ is true,
$\int_{1}^{\infty}\frac{1}{x^{p}}dx=\infty$ and our SDE remains non-explosive
almost surely according to [5]. We can also see from Chow and Khasminskii’s
condition that we have

$x^{p}x+\frac{1}{2}\sigma^{2}-\sigma^{2}<x^{p+1}\leq x^{2}\ln(x^{2})$

which is true since $\frac{x^{p+1}}{x^{2}\ln(x^{2})}\rightarrow 0$ as
$x\rightarrow\infty$ by applying L’Hospital and noting $p<1$.

######  Example 2.5.

Consider the following 2 dimensional SDE with $p>2$
$d\mathbb{X}_{t}=\begin{pmatrix}\lVert X_{t}\rVert^{p}&0\\\ 0&\lVert
X_{t}\rVert^{p}\end{pmatrix}d\mathbb{W}_{t}$. Thus, $A=\begin{pmatrix}\lVert
X_{t}\rVert^{2p}&0\\\ 0&\lVert X_{t}\rVert^{2p}\end{pmatrix}$. Then, since we
have a diagonal matrix, we know $\mathbb{X}_{t}$ is non-explosive from
Khasminskii’s condition

$\frac{2\lVert x\rVert^{2p}}{2}-\lVert x\rVert^{2p}=0\leq C\lVert
x\rVert^{2}\ln(\lVert x\rVert^{2})$

However, we can see that from dimension 3, our condition no longer holds
despite being a martingale.

Now, we will show the extension of Khasminskii’s Theorem for jump processes.
The extension for jump processes follows in an analogous manner to
Khasminskii’s theorem for non-explosion. The most notable difference is that
our condition for Lyapunov function must also hold for the interior of $D$.

######  Theorem 2.6.

Let our SDE with jumps be of the form
$dX_{t}=b(X_{t})dt+\sigma(X_{t})dW_{t}+Y_{t}dN_{t}$ where $dN_{t}$ has
constant poisson arrival. Let $\mathcal{D}\subset\mathbb{R}^{d}$ be a bounded
open set with regular boundary $\partial\mathcal{D}$ and let $\mathcal{D}^{c}$
denote the complement of $\mathcal{D}$. Suppose our SDE parameters satisfy the
Lipschitz condition on any compact subset of $\mathbb{R}^{d}$ for any $t\geq
t_{0}$. Now, assume there exist a positive function $V(x)$ belonging to the
class $\mathcal{C}^{2}$ and positive constants $C$ such that the following
conditions hold

$LV\leq CV$ $V=\inf_{|x|>R}V(x)\rightarrow\infty\>\text{as
}R\rightarrow\infty$

Then it follows that $\mathbb{P}(\tau_{\infty}<\infty)=0$.

###### Proof.

Let $\tau_{R}:=\tau_{R}\wedge t$. We can use the same chain rule argument to
note

$\displaystyle L\\{V(x)e^{c(t-t_{0})}\\}$
$\displaystyle=e^{-c(t-t_{0})}LV(x)-CV(x)e^{-c(t-t_{0})}$
$\displaystyle=e^{-c(t-t_{0})}(LV(x)-CV(x))$ $\displaystyle\leq
0\>\>\>\text{since $e^{-c(t-t_{0})}<0$ and $LV(x)\leq CV(x)$ for
$x\in\mathcal{D}^{c}$}$

Then, we apply Dynkin’s formula again to find

$\mathbb{E}(V(X_{\tau_{R})}e^{-C(\tau_{R}-t_{0})})=V(X_{t_{0}})+\mathbb{E}\left(\int_{0}^{\tau_{R}}LV(x)e^{-C(s-t_{0})}ds\right)$
$\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\leq\mathbb{E}(V(X_{\tau_{R}}))e^{-C(t-t_{0})}\leq
V(X_{t_{0}})$ $\mathbb{E}(V(X_{\tau_{R}}))e^{-C(t-t_{0})}\leq V(X_{t_{0}})$
$\mathbb{E}(V(X_{\tau_{R}}))\leq e^{C(t-t_{0})}V(X_{t_{0}})$

Now, let us split our expectation into two cases where $t\leq\tau_{R}$ and
$t>\tau_{R}$. When $t>\tau_{R}$, we know from definition of hitting time that
$X_{t}$ has reached $D^{c}$. Thus, our expectation on the left hand side
becomes

$\mathbb{E}(V(X_{\tau_{R}}))=\mathbb{E}(V(X_{\tau_{R}}))\mathbbm{1}_{\\{\tau_{R}<t\\}}+\mathbb{E}(V(X_{\tau_{R}}))\mathbbm{1}_{\\{t<\tau_{R}\\}}$

Now note that $\tau_{R}<t$, we have $V(X_{\tau_{R}})\geq\inf_{|x|>R}V(x)$ by
definition of hitting time. For $t<\tau_{R}$, we know $V\geq 0$ from our
positivity assumption. Thus, we have the following string of inequality

$V(X_{0})e^{C(t-t_{0})}\geq\mathbb{E}(V({X_{\tau_{R}}}))\geq\inf_{|x|>R}V(x)\mathbb{P}(\tau_{R}\leq
t)+0\mathbb{P}(\tau_{R}>t)$

Thus, it follows that

$\frac{V(X_{0})e^{C(t-t_{0})}}{\inf_{|x|>R}V(x)}\geq\mathbb{P}(\tau_{R}\leq
t)$

Thus, as $R\rightarrow\infty$, we know from our condition $V\rightarrow\infty$
and thus $\mathbb{P}(\tau_{R}\leq t)=0$ as desired. ∎

######  Remark 2.7.

Note that our Lyapunov function cannot have points of singularities anymore.
Thus, a potential candidate for our Lyapunov function is $V(x)=x^{2}$. Then,
when we consider a Merton Jump Diffusion model with lognormal jump size as
described by the following SDE

$dX_{t}=\mu X_{t}dt+\sigma X_{t}dW_{t}+X_{t}(Y_{t}-1)dN_{t}$

We can show stability with $V(x)=x^{2}$. Note $V(x)\rightarrow\infty$ as
$|x|\rightarrow\infty$. Furthermore

$\displaystyle LV$
$\displaystyle=b(t,X_{t})V^{\prime}(X_{t})+\frac{1}{2}\sigma(t,X_{t})^{2}V^{\prime\prime}(X_{t})+\lambda(X_{t}^{2}\mathbb{E}(Y_{t}^{2})-X_{t}^{2})$
$\displaystyle=2\mu
X_{t}^{2}+\sigma^{2}X_{t}^{2}+\lambda(e^{2\mu_{J}+2\sigma^{2}_{J}}-1)X_{t}^{2}$
$\displaystyle=(2\mu+\sigma^{2}+\lambda(e^{2\mu_{J}+2\sigma^{2}_{J}}-1))X_{t}^{2}$
$\displaystyle\leq(2\mu+\sigma^{2}+\lambda(e^{2\mu_{J}+2\sigma^{2}_{J}}-1))V(X_{t})$

Thus, we have found our $C$ as desired.

######  Theorem 2.8.

Corollary 3.1 from Khasminskii[3]. Assume the conditions from theorem 3.1
hold. Now, let $\mathcal{D}$ be a set with regular boundaries
$\partial\mathcal{D}$. Let $D_{n}$ be an increasing sequence such that
$\cup^{\infty}\mathcal{D}_{n}=\mathcal{D}$. Assume $V(t,x)$ be twice
differentiable with respect to $x$ and once with respect to $t$. Let
$X_{0}\in\mathcal{D}$.

###### Proof.

The proof follows analogously as before but with slight modifications to our
hitting time conditions. We will provide an example to illuminate the details,
but for the sake of brevity the proof will be omitted. ∎

######  Example 2.9.

Consider an SDE of the form: $dX_{t}=X_{t}^{-\alpha}dt+dW_{t}$ where
$\alpha>1$. Note that when $x=0$, $b(x)=x^{-\alpha}=\infty$. Now, let
$D=(0,\infty)$, and $D_{n}=(\frac{1}{n},\infty)$. Let $V(x)=\frac{1}{x}$. Then
$V^{\prime}(x)=-x^{-2}$ and $V^{\prime\prime}(x)=2x^{-3}$. Note:

$\displaystyle LV$ $\displaystyle=-x^{-2-\alpha}+x^{-3}$
$\displaystyle=x^{-1}(x^{-2}-x^{1-\alpha})$

Note since $\alpha>1$, we know $1+\alpha>2$. Thus:

$\lim_{x\rightarrow 0}(x^{-2}-x^{-1-\alpha})=-\infty$
$\lim_{x\rightarrow\infty}(x^{-2}-x^{-1-\alpha})=0$
$\sup_{x>0}(x^{-2}-x^{-1-\alpha})=C<\infty$

Note, since we have $LV\leq CV$, we can use the chain rule argument to note

$\mathbb{E}(X_{t})^{-1}\leq X_{0}^{-1}e^{Ct}$

Now let $\tau_{n}:=\inf\\{t>0:X_{t}\leq\frac{1}{n}\\}$. Then it follows:

$\mathbb{E}(X_{t\wedge\tau_{n}})^{-1}\leq X_{0}^{-1}e^{Ct\wedge\tau_{n}}$

Now when $\tau_{n}\leq t$, we know $X_{t\wedge\tau_{n}}=\frac{1}{n}$. Thus, we
have

$n\mathbb{P}(\tau_{n}\leq t)\leq X_{0}^{-1}e^{Ct}$ $\mathbb{P}(\tau_{n}\leq
t)\leq\frac{e^{Ct}}{nX_{0}}$

Now, as $n\rightarrow\infty$, we know $\mathbb{P}(\tau_{n}\leq t)=0$ as
desired. Thus, we have shown that $X_{t}$ cannot become negative or explode to
infinity and will always reside within $D=(0,\infty)$.

## 3 Positive Probability of Explosion Through Lyapunov Functions

######  Theorem 3.1.

Chow and Khasminskii’s Theorem 1. Let $\mathcal{D}\subset\mathbb{R}^{d}$ be a
bounded open set with regular boundary $\partial\mathcal{D}$ and let
$\mathcal{D}^{c}$ denote the complement of $\mathcal{D}$. Suppose our SDE
parameters satisfy the Lipschitz condition on any compact subset of
$\mathbb{R}^{d}$ for any $t\geq t_{0}$. Now, assume there exist a positive
function $V(t,x)$ belonging to the class
$\mathcal{C}^{1,2}([t_{0},\infty]\times\mathcal{D}^{c})$ and positive
constants $K_{1},\>K_{2},\>K_{3}$ and $C$ such that the following conditions
hold

$\sup_{t\geq t_{0},x\in\mathcal{D}^{c}}V(x)=K_{1}<\infty$ $\sup_{t\geq
t_{0},x\in\mathcal{D}}V(x)=K_{2}<\inf_{t\geq t_{0}x\in\Gamma}V(x)=K_{3}$

for some set $\Gamma\subset\mathcal{D}^{c}$, and

$LV(x)\geq CV(x)$

for $t\geq t_{0},\>x\in\mathcal{D}^{c}$. Then the solution $X_{t}$ with the
initial condition $X_{t_{1}}=x\in\Gamma$ for some $t_{1}\geq t_{0}$ has an
explosion with positive probability.

###### Proof.

Let $\tau_{R}:=\min\\{\tau_{R(t)},t,\tau_{\partial\mathcal{D}}\\}$ (the
process stops once it hits the boundary of $\mathcal{D}$ or leaves the disk
$R$). The following is intended to serve as a more intuitive explanation of
the proof provided by Khasminskii and Chow. Consider a disk $R$ of radius $r$,
a disk inside $D$ in $R$ and a disk $\Gamma$ encircling $D$ (for our
construction, $\Gamma\cap D=\emptyset$) as shown below.

$R$$\Gamma$$D$

Let our process start from $\Gamma$ and our goal is to show there is a
positive probability of reaching the boundary of $R$. Note, for
$x\in\mathcal{D}^{c}$, we have by applying the product rule

$\displaystyle L\\{V(x)e^{c(t-t_{0})}\\}$
$\displaystyle=e^{-c(t-t_{0})}LV(x)-CV(x)e^{-c(t-t_{0})}$
$\displaystyle=e^{-c(t-t_{0})}(LV(x)-CV(x))$ $\displaystyle\geq
0\>\>\>\text{since $e^{-c(t-t_{0})}>0$ and $LV(x)\geq CV(x)$ for
$x\in\mathcal{D}^{c}$}$

Now, from Dynkin’s formula, we find

$\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})=V(X_{t_{0}})+\mathbb{E}\left(\int_{0}^{\tau_{R}}LV(x)e^{-C(t-t_{0})}ds\right)$

Thus

$\mathbb{E}(V(\tau_{R})e^{-C(\tau_{R}-t_{0})})\geq V(X_{t_{0}}).$

Now, we can divide our expectation into three cases: the case where we have
hit $R$, the case where we have hit $\partial D$, and the case where we have
hit neither $R$ or $\partial D$. This splitting of expectation can be denoted
in the following manner

$\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\mathbbm{1}_{\\{\tau_{R}<t\wedge\tau_{\partial
D}\\}}+\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\mathbbm{1}_{\\{\tau_{\partial
D}<t\wedge\tau_{R}\\}}$
$+\mathbb{E}(V(X_{t})e^{-C(t-t_{0})})\mathbbm{1}_{\\{t<\tau_{R}\wedge\tau_{\partial
D}\\}}\geq V(X_{t_{0}})$

Now from the conditions on our Lyapunov function, we note

$\displaystyle\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\mathbbm{1}_{\\{\tau_{\partial
R}<t\wedge\tau_{\partial D}\\}}$ $\displaystyle\geq
V(X_{t_{0}})-K_{2}-K_{1}e^{-C(t-t_{0})}$ $\displaystyle\geq
K_{3}-K_{2}-K_{1}e^{-C(t-t_{0})}$

Now, as $R\rightarrow\infty,\>t\rightarrow\infty$, we note

$\mathbb{E}(K_{1}e^{-C(\tau_{\infty}-t_{0})})\geq K_{3}-K_{2}>0$

Thus, it follows that $\mathbb{P}(\tau_{\infty}<\infty)>0$ since
$\mathbb{P}(\tau_{\infty}=\infty)=1$ would imply that our expectation on the
left is zero. ∎

######  Remark 3.2.

Analogous to Remark 2.3, Khasminskii and Chow provides the condition for
positive probability for explosion in [1]

$\langle
b(t,x),x\rangle+\frac{1}{2}\text{Tr}(A(t,x))-\lambda_{\text{max}}\left(1+\frac{1+\epsilon}{\ln(\lVert
x\rVert^{2})}\right)\geq C\lVert x\rVert^{2}(\ln(\lVert
x\rVert^{2}))^{1+\epsilon}$

The derivation follows the exact same steps as shown in Remark 2.3 but with
$V(x)=K-\frac{1}{\ln(\lVert x\rVert)^{\epsilon}}$ and will thus be skipped.

######  Example 3.3.

Consider our 3 dimensional extension of example 2.4. Let our SDE be of the
form: $d\mathbb{X}_{t}=\begin{pmatrix}\lVert X_{t}\rVert^{p}&0&0\\\ 0&\lVert
X_{t}\rVert^{p}&0\\\ 0&0&\lVert X_{t}\rVert^{p}\end{pmatrix}d\mathbb{W}_{t}$,
where $p>1$. Note $A=\begin{pmatrix}\lVert X_{t}\rVert^{2p}&0&0\\\ 0&\lVert
X_{t}\rVert^{2p}&0\\\ 0&0&\lVert X_{t}\rVert^{2p}\end{pmatrix}$. Then, from
remark 3.2, we need to show there exists $C$ and $\epsilon$ such that

$\frac{3}{2}\lVert x\rVert^{2p}-\lVert
x\rVert^{2p}\left(1+\frac{1+\epsilon}{\ln(\lVert x\rVert^{2})}\right)\geq
C\lVert x\rVert^{2}(\ln(\lVert x\rVert^{2})^{1+\epsilon}$

holds for all $\lVert x\rVert>R$. Note that it suffices to check

$\frac{\ln(\lVert x\rVert^{2})^{1+\epsilon}}{\lVert x\rVert^{p-1}}\rightarrow
0\>\text{as $\lVert x\rVert\rightarrow\infty$}$

Note by properties of exponents and L’Hospital, we have

$\frac{2^{1+\epsilon}\frac{(1+\epsilon)\ln(\lVert x\rVert)}{\lVert
x\rVert}}{\frac{p-1}{\lVert x\rVert^{2-p}}}\equiv\frac{\ln(\lVert
x\rVert)}{\lVert x\rVert^{p-1}}\rightarrow 0$

Thus, despite being a martingale, our process without drift can explode with
positive probability from the noise alone in dimension 3.

Now, let us extend theorem 3.1 to jump processes

######  Theorem 3.4.

Extending Theorem 3.1 to Jump Processes. Let our SDE be of the form
$dX_{t}=b(X_{t},t)dt+\sigma(X_{t},t)dW_{t}+Y_{t}dN_{t}$ where $Y_{t}$ is the
jump size and $N_{t}$ is poisson with constant intensity. Now, let
$\mathcal{D}\subset\mathbb{R}^{d}$ be a bounded open set with regular boundary
$\partial\mathcal{D}$ and let $\mathcal{D}^{c}$ denote the complement of
$\mathcal{D}$. Suppose our SDE parameters satisfy the Lipschitz condition on
any compact subset of $\mathbb{R}^{d}$ for any $t\geq t_{0}$. Now, assume
there exist a positive function $V(t,x)$ belonging to the class
$\mathcal{C}^{1,2}([t_{0},\infty]\times\mathcal{D}^{c})$ and positive
constants $K_{1},\>K_{2},\>K_{3}$ and $C$ such that the following conditions
hold

$\sup_{t\geq t_{0},x\in\mathcal{D}^{c}}V(x)=K_{1}<\infty$ $\sup_{t\geq
t_{0},x\in\mathcal{D}}V(x)=K_{2}<\inf_{t\geq t_{0}x\in\Gamma}V(x)=K_{3}$

for some set $\Gamma\subset\mathcal{D}^{c}$, and

$LV(x)\geq CV(x)$

for $t\geq t_{0},\>x\in\mathcal{D}^{c}$.

###### Proof.

Note that even with the presence of jumps, since our choice of $V$ depends
solely on $x$, we can still apply the product rule argument used in theorem
1.1

$\displaystyle L\\{V(x)e^{c(t-t_{0})}\\}$
$\displaystyle=e^{-c(t-t_{0})}LV(x)-CV(x)e^{-c(t-t_{0})}$
$\displaystyle=e^{-c(t-t_{0})}(LV(x)-CV(x))$ $\displaystyle\geq
0\>\>\>\text{since $e^{-c(t-t_{0})}>0$ and $LV(x)\geq CV(x)$ for
$x\in\mathcal{D}^{c}$}$

However, since our process can jump over the boundary of $D$, we simply change
our hitting time condition to be the whole set $D$ and the proof follows
analogously. From Dynkin’s formula, we find

$\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})=V(X_{t_{0}})+\mathbb{E}\left(\int_{0}^{\tau_{R}}LV(x)e^{-C(t-t_{0})}ds\right)$

Thus

$\mathbb{E}(V(\tau_{R})e^{-C(\tau_{R}-t_{0})})\geq V(X_{t_{0}}).$

Now, we can divide our expectation into three cases: the case where we have
hit $R$, the case where we have hit $\partial D$, and the case where we have
entered neither $R$ or $D$. This splitting of expectation can be denoted in
the following manner

$\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\mathbbm{1}_{\\{\tau_{R}<t\wedge\tau_{D}\\}}+\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\mathbbm{1}_{\\{\tau_{D}<t\wedge\tau_{R}\\}}$
$+\mathbb{E}(V(X_{t})e^{-C(t-t_{0})})\mathbbm{1}_{\\{t<\tau_{R}\wedge\tau_{D}\\}}\geq
V(X_{t_{0}})$

Now from the conditions on our Lyapunov function, we note

$\displaystyle\mathbb{E}(V(X_{\tau_{R}})e^{-C(\tau_{R}-t_{0})})\mathbbm{1}_{\\{\tau_{R}<t\wedge\tau_{D}\\}}$
$\displaystyle\geq V(X_{t_{0}})-K_{2}-K_{1}e^{-C(t-t_{0})}$ $\displaystyle\geq
K_{3}-K_{2}-K_{1}e^{-C(t-t_{0})}$

Now, as $R\rightarrow\infty,\>t\rightarrow\infty$, we note

$\mathbb{E}(K_{1}e^{-C(\tau_{\infty}-t_{0})})\geq K_{3}-K_{2}>0$

Thus, it follows that $\mathbb{P}(\tau_{\infty}<\infty)>0$ since
$\mathbb{P}(\tau_{\infty}=\infty)=1$ would imply that our expectation on the
left is zero. ∎

## 4 Almost Sure Explosion

In Khasminskii and Chow, a theorem containing the sufficient condition for
almost sure explosion using Lyapunov functions is provided. In this section,
we will provide an intuitive explanation of their proof and provide a
sufficient condition for almost sure explosion in 1 dimension, namely
ellipticity.

######  Theorem 4.1.

Khasminskii and Chow’s Almost Sure Explosion [1]. Suppose the conditions for
theorem 3.1 hold. Assume

$\inf_{t\geq t_{0},x\in\mathcal{D}^{c}}V(t,x)=K_{0}>0$

and for $x\in\partial D$

$\mathbb{P}(\tau_{\Gamma}<\infty)=1$

Then for any $t\geq t_{0},x\in\mathbb{R}^{n}$,
$\mathbb{P}(\tau_{\infty}<\infty)=1$.

###### Proof.

The proof can be found in [1]. ∎

######  Remark 4.2.

The main additional constraint is having $\mathbb{P}(\tau_{\Gamma}<\infty)=1$.
Note, because we are assuming the conditions for theorem 2.2, we know that
when our process reaches $\Gamma$, there is a positive probability of reaching
the outer disk $R$ in finite time, or $\mathbb{P}(\tau_{\infty}<\infty)>0$.
Thus, our process starts from $\Gamma$ and either returns to $\Gamma$ in a
finite amount of time or hits $\partial\mathcal{D}$ and returns to $\Gamma$ in
a finite amount of time with probability 1. Thus, every time our process
returns to $\Gamma$ we have a positive probability of explosion, and so
eventually our process must explode as you cannot have an infinite coin toss
without observing a heads.

A sufficient criterion for $\mathbb{P}(\tau_{\Gamma}<\infty)=1$ is
ellipticity[4], which for one dimensional SDE of the form
$dX_{t}=b(x,t)dt+\sigma(x,t)dW_{t}$, has $\sigma(x,t)>c_{0}>0$. Or in short,
the noise cannot be killed. It is worth mentioning that ellipticity in
multiple dimension is defined by the following condition

$v^{T}\sigma(x)\sigma(x)^{T}v\geq c_{0}\lVert v\rVert^{2}$

where $\sigma$ is a matrix and $v$ is a vector. In short, as long as the noise
stays bounded above zero, our process cannot get stuck inside $D$ and must
reach $\Gamma$ in a finite amount of time as desired.

######  Example 4.3.

Leon and Villa has already shown that for $b(X_{t})>0$ and increasing with
$\sigma(t,X_{t})=\sigma$, where $\sigma$ is a constant, we have almost sure
explosion when $\int_{1}^{\infty}\frac{1}{b(y)}dy<\infty$ and almost sure no-
explosion when $\int_{1}^{\infty}\frac{1}{b(y)}dy=\infty$[5]. Let us further
assume that $\frac{b^{\prime}(x)}{b(x)^{2}}\rightarrow 0$ as
$x\rightarrow\infty$. Note that for the almost sure explosion case, we can
choose $V(x)=\int_{1}^{x}\frac{1}{b(y)}dy$. Then
$\lim_{x\rightarrow\infty}V(x)=K_{1}<\infty$ and thus our supremum and infimum
conditions are satisfied. Furthermore, since $b(x)$ is positive monotonic
increasing, we know $\sup_{x\in\partial D}V(x)<\inf_{x\in\Gamma}V(x)$. Now, we
need to check that there exists a positive $C$ such that $LV\geq CV$. Note

$LV=\frac{b(x)}{b(x)}-\frac{\sigma}{2}\frac{b^{\prime}(x)}{b(x)^{2}}$

Thus, from our assumption, we know for a sufficiently large $x$, we know our
$LV>0$. Furthermore, since $V(x)=\int_{1}^{x}\frac{1}{b(y)}dy<K_{1}$, we can
choose always choose $C$ to be sufficiently smaller than our $LV$, so our
condition $LV\geq CK_{1}\geq CV$ also holds.

######  Remark 4.4.

For the cases when $b(x)$ is a simple polynomial like $b(x)=x^{p}$, we know
$\int^{\infty}_{1}\frac{1}{x^{p}}dx<\infty$ when $p>1$. Then, we know
$\frac{b^{\prime}(x)}{b(x)^{2}}=\frac{(p-1)x^{p-1}}{x^{2p}}\rightarrow 0$ as
$x\rightarrow\infty$. Thus, it is worth noticing that the added assumption
$\frac{b^{\prime}(x)}{b(x)^{2}}$ does not provide much hindrance for many of
the cases one may find in nature.

## 5 Conclusion

Thus, in this paper, we started our expedition of finite time explosion for
SDEs through ODE comparisons as presented by Leon and Villa, which stated that
explosion for SDEs with constant intensity noise follows analogously with the
explosion of ODEs through the Osgood criterion. However, from the prominent
Feller’s test for explosion, we were able to observe that for different types
of noise functions, we could have almost sure non-explosion and positive
probability of explosion regardless of the Osgood criterion. Then, we explored
Khasminskii’s Lyapunov function methods to classify explosion for SDE in
multiple dimensions and provide an extension to jump processes. Nevertheless,
because of the one-dimensional aspect of Feller’s test for explosion and the
lengthy integration tests, we adopted Khasminskii’s method of Lyapunov
functions to characterize explosions for SDEs.

We first establish Khasminskii’s conditions for almost sure non-explosion
using Khasminskii’s Lyapunov method. We then show the almost sure non-
explosion for constant intensity is consistent with Lyon Villa’s finding with
the Osgood criterion as shown in Example 2.4. We show further examples on a
2-dimensional Ito integral that will be compared with the 3-dimensional Ito
integral to illuminate the dependence on dimension with explosions.

We then use Khasminskii’s theorem on almost sure explosion and show that a
sufficient condition for explosion is ellipticity, which in one dimension is
consistent with the findings of Leon and Villa with the Osgood criterion.
Afterwards, we show through Khasminskii’s criterion for positive probability
of explosion which is derived through Khasminskii’s Lyapunov method for
positive probability of explosion. We then show that the 3-dimensional Ito
integral, which has zero mean, can explode with positive probability.

## Acknowledgement

I am grateful to Professor M. Salins for introducing me to the problem of
finite time explosion in SDEs and the Osgood criterion. His super power of
glancing at a paper once and explaining it in simple terms helped me choose
the most illuminating examples to include.

## 6 Appendix

We show that an Ito integral is a martingale with zero mean. Note from
definition of Ito integral $I_{n}$

$\displaystyle\mathbb{E}(I_{n+1}|\mathcal{F}_{n})$
$\displaystyle=\mathbb{E}\left(\sum_{k=0}^{n}\Delta_{k}(W_{n+1}-W_{n})\bigg{|}\mathcal{F}_{n}\right)$
$\displaystyle=\mathbb{E}\left(\sum_{k=0}^{n-1}\Delta_{k}(W_{n}-W_{n-1})+\Delta_{n}(W_{n+1}-W_{n})\bigg{|}\mathcal{F}_{n}\right)$
$\displaystyle=I_{n}+\mathbb{E}(\Delta_{n}(W_{n+1}-W_{n})|\mathcal{F}_{n})$
$\displaystyle=I_{n}\>\>\>\>\text{as
$\mathbb{E}(W_{n+1}-W_{n}|\mathcal{F}_{n})=0$}$

Thus, we find the Ito integral is a martingale as desired.

## References

  * [1] Pao-Liu Chow and Rafail Khasminskii.  Almost sure explosion of solutions to stochastic differential equations.  Stochastic Processes and their Applications, 124(1):639–645, 2014.
  * [2] Ioannis Karatzas and Steven Shreve.  Brownian motion and stochastic calculus, volume 113.  springer, 2014.
  * [3] Rafail Khasminskii.  Stochastic stability of differential equations.  Springer, 2012.
  * [4] Leonid Koralov and Yakov G Sinai.  Theory of probability and random processes.  Springer Science & Business Media, 2007.
  * [5] Jorge A León and José Villa.  An osgood criterion for integral equations with applications to stochastic differential equations with an additive noise.  Statistics & probability letters, 81(4):470–477, 2011.
  * [6] Bernt Oksendal.  Stochastic differential equations: an introduction with applications.  Springer Science & Business Media, 2013.
  * [7] Bernt Øksendal and Agnes Sulem.  Stochastic control of jump diffusions.  In Applied Stochastic Control of Jump Diffusions, pages 93–155. Springer, 2019.
  * [8] William Fogg Osgood.  Beweis der existenz einer lösung der differentialgleichung $\frac{dy}{dx}=f\left(x,y\right)$ ohne hinzunahme der cauchy-lipschitz’schen bedingung.  Monatshefte für Mathematik und Physik, 9:331–345, 1898.


# Comparison of the LASSO and Integrative LASSO with Penalty Factors (IPF-
LASSO) methods for multi-omics data: Variable selection with Type I error
control

Charlotte Castel, Zhi Zhao, Magne Thoresen C. Castel, Postdoctoral Researcher,
Z. Zhao, Researcher, and M. Thoresen, Professor, are with the Centre for
Biostatistics and Epidemiology (OCBE), Institute of Basic Medical Sciences,
Faculty of Medicine, University of Oslo, Oslo, Norway (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>magne.thoresen@medisin.uio.no).

###### Abstract

Variable selection in relation to regression modeling has constituted a
methodological problem for more than 60 years. Especially in the context of
high-dimensional regression, developing stable and reliable methods,
algorithms, and computational tools for variable selection has become an
important research topic. Omics data is one source of such high-dimensional
data, characterized by diverse genomic layers, and an additional analytical
challenge is how to integrate these layers into various types of analyses.
While the IPF-LASSO model has previously explored the integration of multiple
omics modalities for feature selection and prediction by introducing distinct
penalty parameters for each modality, the challenge of incorporating
heterogeneous data layers into variable selection with Type I error control
remains an open problem. To address this problem, we applied stability
selection as a method for variable selection with false positives control in
both IPF-LASSO and regular LASSO. The objective of this study was to compare
the LASSO algorithm with IPF-LASSO, investigating whether introducing
different penalty parameters per omics modality could improve statistical
power while controlling false positives. Two high-dimensional data structures
were investigated, one with independent data and the other with correlated
data. The different models were also illustrated using data from a study on
breast cancer treatment, where the IPF-LASSO model was able to select some
highly relevant clinical variables.

###### Index Terms:

High-dimensional data, penalized regression, multi-omics integration, variable
selection, error control

##  I Introduction

Variable selection in regression models has been an important methodological
challenge for more than 60 years [1]. Developing stable and dependable
methods, algorithms, and computational tools for variable selection stands as
a crucial task, especially in the challenging context of high-dimensional
regression, where the number of potential explanatory variables exceeds the
number of observations (p $\gg$ n). The conventional approach assumes
sparsity, implying that only a small subset of the explanatory variables
genuinely influences the outcome. Methods like the LASSO[2] have become very
popular for integrating variable selection into the estimation processes.

A hot topic in omics research during recent years is integration of different
layers of genomic information. Due to technological developments, we often
have available diverse high-dimensional measurements such as DNA methylation,
copy number variation, and mRNA expression for the same individual, in
addition to more traditional clinical data. A structured integration of such
different layers gives an opportunity to gain a holistic understanding of the
biological system in question. For prediction purposes, the IPF-LASSO [3], the
Priority-Lasso [4], the IPFStructPenalty [5], the mix-lasso [6] and the
cooperative learning[7] methods have been suggested. They all perform variable
selection, but none of them incorporate any type of error control. Thus, the
challenge of integrating heterogeneous data layers into methods for variable
selection with Type I error control (false positive) remains an open problem.

The IPF-LASSO model was initially developed in 2017 to address the integration
of multiple omics modalities for variable selection in predictive modeling.
Boulesteix et al. [3] proposed a penalized regression method to address this
challenge, applying distinct penalty factors for different data modalities.
The results in this article indicate that the IPF-LASSO method and regular
LASSO are comparable in terms of prediction accuracy, but with fewer variables
selected for the IPF-LASSO. These results would indicate that the IPF-LASSO
could select variables with fewer false positives compared to the regular
LASSO, which again begs the question if the IPF-LASSO would achieve higher
power under proper error control, compared to the regular LASSO.

Managing error control in statistical analyses can be done in various ways; by
controlling the false discovery rate (FDR), controlling the family-wise error
rate (FWR), or controlling the number of false positives. To our knowledge,
there are four main approaches for applying error control in variable
selection in high-dimensional regression models.(i) Stability selection [8, 9]
aggregates results from applying a selection procedure to subsamples,
controlling the expected number of false positives. (ii) Data splitting
methods [10, 11, 12, 13, 14], controlling FDR or FWR, dividing the data
randomly into two parts for variable selection and inference, respectively.
(iii) Selective inference [15, 16, 17] also controls FDR or FWR, by formally
conditioning on the selection event. (iv) Model-X knockoff and its like [18]
generates knockoff predictors to mimic the covariance pattern of the original
predictors, controlling FDR in variable selection. The work presented here
focuses on the control of the number of false positives.
To our knowledge, there is no existing work on variable selection with error
control in high-dimensional models integrating multiple omics modalities.
Stability selection is a widely used method for variable selection with
control over the expected number of false positives. We applied stability
selection to the IPF-LASSO. The objective was to compare the popular LASSO
algorithm with the IPF-LASSO to evaluate if introducing different penalty
factors per data modality improves power while controlling the number of false
positives, compared to the standard LASSO. Two high-dimensional data settings
were investigated: one setting with independent variables and one setting with
correlated variables. The simulation design is the same as in the IPF-LASSO
article [3], determining the best penalty parameters through cross-validation
(CV). The focus in this article was on binary outcomes and logistic models,
and we will keep the same focus, although the methods are of course equally
relevant also for other models.
In our simulation studies, we observe that both the IPF-LASSO and the LASSO
control the number of false positives well; however, the IPF-LASSO increases
the power compared to standard LASSO, particularly in scenarios with a small
true model size, a high absolute difference between proportions of relevant
variables in the two modalities, and a small ratio between the size of the
smaller modality and the larger modality. For a real data application, we
chose to test our models for variable selection with error control on the same
breast cancer dataset from the study by Hatzis et al. [19], as utilized in the
article by Boulesteix et al. [3] This dataset includes a relatively small set
of clinical variables in addition to a high-dimensional set of gene
expressions, which constitutes our two data modalities.
This article is divided as follows. After this introduction, Section 2
describes the different methods, Section 3 presents the results from the
simulation studies, and Section 4 presents the results from the analysis of
the breast cancer data. In Section 5 we discuss our findings and conclude.

##  II Methods

###  II-A Integrative LASSO with Penalty Factors (IPF-LASSO)

Let $x_{ij}$ represent the predictor variable $j$ measured on subject $i$, and
$y_{i}$ represent the (continuous) response values associated, where
$i=1,...,n$ and $j=1,...,p$. The idea of the standard LASSO method for the
linear model involves solving the $\ell_{1}$-penalized regression problem by
determining $\boldsymbol{\beta}$ = $\left\\{\beta_{j}\right\\}$ that minimizes
the loss function:

$\sum_{i=1}^{n}\left(y_{i}-\sum_{j=1}^{p}x_{ij}\beta_{j}\right)^{2}+\lambda\lVert\boldsymbol{\beta}\rVert_{1},$

where the $\ell_{1}$-norm is denoted by

$\lVert.\rVert_{1}=\sum_{j=1}^{p}|\beta_{j}|$

and $\lambda$ is the penalty parameter.

Readers can explore Tibshirani’s work for more information about LASSO
regression [2]. The LASSO framework has been extended to other types of
regression models, e.g. logistic regression [20] which we will take advantage
of in this work

Boulesteix et al. [3] proposed to deal with the case of multiple data
modalities by introducing a small modification to the classical LASSO. Let us
denote by $X_{1}^{(m)}$, … , $X_{p_{m}}^{(m)}$ the variables belonging to
modality $m$ (for $m=1$,…,$M$), with their respective observed values
$x_{i1}^{(m)}$, …, $x_{ip_{m}}^{(m)}$ for subject $i$ (for $i=1,...,n$). The
value $p_{m}$ corresponds to the number of variables within modality $m$. In
the same way, $\beta_{j}^{(m)}$ $(j=1,...,p_{m})$ represents the regression
parameters associated with variable $X_{j}^{(m)}$ from modality $m$. The idea
of the IPF-LASSO is to use a weighted sum of the $l_{1}$-norm of the
coefficient vector of each data modality as the penalty term to take into
account the differences between data modalities. The new $l_{1}$-penalized
regression problem is solved by minimizing

$\sum_{i=1}^{n}\left(y_{i}-\sum_{m=1}^{M}\sum_{j=1}^{p_{m}}x_{ij}^{(m)}\beta_{j}^{(m)}\right)^{2}+\sum_{m=1}^{M}\lambda_{m}\lVert\boldsymbol{\beta}^{(m)}\rVert_{1}$

with $\lambda_{m}>0$ being the penalty parameter applied to the variables from
modality $m$. How to choose $\lambda_{m}$ for $m=1,...,M$ becomes a
challenging task. Boulesteix et al. introduce the penalty factors to solve
this problem. As explained in Boulesteix et al., the term “penalty factors”
refers to the multiplicative factors applied to the penalty term. For
simplicity, the reference modality is the first modality with penalty
$\lambda_{1}$ and associated penalty factor 1. The penalty factor for any
given modality $\lambda_{m}$ is then determined by the ratio
$\lambda_{m}/\lambda_{1}$. The estimation, the connection with the other LASSO
variations and the choice of the penalty parameters are well explained in the
work of Boulesteix et al.[3].

###  II-B Stability selection

Meinshausen and Bühlmann introduced stability selection in 2010 [8, 12],
establishing it as a widely used method for stable variable selection and for
controlling false positives. Stability selection improves the performance of a
variable selection algorithm by applying this algorithm multiple times to
random subsamples of the data, each of size $\lfloor n/2\rfloor$. The
variables selected are those chosen most frequently across subsamples.
The main attractive aspect of stability selection comes from its control of
false positives through an upper bound on the expected number of falsely
selected variables, termed V for the rest of the article. However, this
control holds under specific assumptions about the selection procedure,
requiring it to surpass random guessing and a robust exchangeability condition
concerning the selection of noise variables.
To address the issues of these assumptions, a variation called Complementary
Pairs Stability Selection (CPSS) was proposed by Rajen et al. [9] in 2013.
CPSS divides variables into two groups based on their probability of
selection: one with a low probability and the other with a high probability.
Subsamples as complementary pairs are drawn, and the selection procedure is
applied to each. The selection probability for each variable is estimated,
representing the proportion of fitted models where the associated coefficient
estimate differs from zero. The final selected variables are those with a high
selection probability, exceeding a user-defined cut-off. Stability selection
can be used for Type I error control through the upper bound V. However, in
practice, stability selection is often used with the aim of stabilizing the
obtained results in terms of selected variables, simply by choosing some (non-
optimal) selection threshold. We will introduce stability selection to the
IPF-lasso for variable selection with error control.

##  III Simulation

###  III-A Simulation design

In this simulation study, we choose to explore three stability selection
thresholds: manually chosen at 70%, manually chosen at 80%, and finally the
formally optimal threshold, determined by the r-concave bound proposed by
Rajen et al. [9]. This bound is a function of the stability selection
threshold. Thus, by selecting an appropriate upper bound for the expected
number of false positive selections, one can calculate the optimal threshold.
All parameters that go into this calculation are selected as described in
Rajen et al. [9].

The value of B, corresponding to the number of complementary pairs of
subsamples is set to 50. Stability selection with these three different
thresholds will be applied both to the regular LASSO and to the IPF-LASSO.
This will create six different procedures to compare. The error control, by
the number of false positives, and the power of these six models will be
evaluated on simulated data. Given the multitude of potential simulation
configurations, a lot of settings could be tested out. To make the comparison
more comprehensible, we chose to design our simulation study very close to the
one in the original IPF-LASSO article[3].
We simulate situations with two data modalities. Each modality is
characterized by two parameters: the number of variables in each modality
($p_{1}$ and $p_{2}$) and the number of non-zero coefficients in each modality
($b_{1}$ and $b_{2}$). We chose a fixed signal strength of $\beta=1$ for both
modalities across all simulations. Each design had a sample size of $n=100$,
and we generated 100 simulated datasets for each configuration. For all
designs, we used a binary outcome. All considered designs in this simulation
study are detailed in Table I.

| $p_{1}$ | $p_{2}$ | $b_{1}$ | $b_{2}$
---|---|---|---|---
Design A | 1000 | 1000 | 10 | 10
Design B | 100 | 1000 | 3 | 30
Design C | 100 | 1000 | 10 | 10
Design D | 100 | 1000 | 20 | 0
Design E | 20 | 1000 | 3 | 10
Design F | 20 | 1000 | 15 | 3
Design G | 20 | 1000 | 10 | 10
Design H | 20 | 1000 | 3 | 3
Design I | 20 | 2000 | 20 | 0
TABLE I:  Configurations of the design used in the simulations

These particular designs were considered to be more or less realistic, as in
the situation of “useless omics data” where one of the modalities does not
contain non-zero variables.
We considered two data settings: the first, called “independent data”, with no
correlation within or between the two modalities, and the second, called
“correlated data”, with a correlation both within and across modalities in the
same dataset. This second setting with correlated data is obviously considered
to be closer to real data.

The “independent data” are simulated as follows:

  1. 1.

The binary outcome is sampled from the Bernoulli distribution
${\mathcal{B}}(0.5)$.

  2. 2.

The variables of the datasets are sampled from the multivariate normal
distribution:
$X_{1},X_{2},\ldots,X_{p_{1}+p_{2}}|Y=0\sim\mathcal{M}\mathcal{N}(0_{p_{1}+p_{2}},\Sigma)$
$X_{1},X_{2},\ldots,X_{p_{1}+p_{2}}|Y=1\sim\mathcal{M}\mathcal{N}(\mu,\Sigma)$
with $\mu^{T}=\left(\overbrace{\underbrace{\beta,\ldots,\beta}_{\text{length
}b_{1}},0,\ldots,0,}^{\text{length
}p_{1}}\overbrace{\underbrace{\beta,\ldots,\beta}_{\text{length
}b_{2}},0,\ldots,0}^{\text{length }p_{2}}\right)$
and $\Sigma=\mathrm{I}_{p_{1}+p_{2}}={\begin{pmatrix}1&0&\cdots&0\\\
0&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&0\\\
0&\cdots&0&1\end{pmatrix}}$

The “correlated data” are simulated in the same way as the “independent data”
but the $\Sigma$ matrix is modified. To take into account both the correlation
between variables of the same data modality and between variables of different
modalities, the Identity matrix is replaced by a covariance matrix, assuming
that there are b = 10 sets of variables that are mutually correlated within
each modality. The correlation between variables belonging to the same block
group is fixed at $\rho=0.4$, and defined by the matrix $\Sigma$ as in Figure
1.

Figure 1: Covariance matrix $\Sigma$ simulated for the “correlated data”
$\Sigma=\left(\begin{array}[]{@{}c|c@{}}\begin{matrix}A_{p_{1}/b}(\rho)&&&{\text{0}}\\\
&A_{p_{1}/b}(\rho)&&\\\ &&\ddots&\\\ {\text{0}}&&&A_{p_{1}/b}(\rho)\\\
\end{matrix}&\begin{matrix}B_{p_{1}/b,p_{2}/b}(\rho)&&{\text{0}}\\\
&\ddots&\\\ {\text{0}}&&B_{p_{1}/b,p_{2}/b}(\rho)\\\ \end{matrix}\\\
\cline{1-2}\cr\begin{matrix}B_{p_{2}/b,p_{1}/b}(\rho)&&{\text{0}}\\\
&\ddots&\\\ {\text{0}}&&B_{p_{2}/b,p_{1}/b}(\rho)\\\
\end{matrix}&\begin{matrix}A_{p_{2}/b}(\rho)&&&{\text{0}}\\\ &\ddots&\\\
{\text{0}}&&A_{p_{2}/b}(\rho)\\\ \end{matrix}\\\ \end{array}\right)$

with $A_{q}(\rho)$ and $B_{q,q}(\rho)$ are $(q\times q)$ matrix defined as
follow: $A_{q}(\rho)=\begin{pmatrix}1&&\rho\\\ &\ddots&\\\ \rho&&1\\\
\end{pmatrix}$ and $B_{q,q}(\rho)=\begin{pmatrix}\rho&&\rho\\\ &\ddots&\\\
\rho&&\rho\\\ \end{pmatrix}$.

In order to randomly distribute the non-zero variables over the correlated
blocks, the columns are randomly permuted after generating the data. These two
simulation settings are the ones proposed in the original article about the
IPF-LASSO.

We applied the six procedures presented previously to these simulation
designs. For the IPF-LASSO method, the 11 following combinations of penalty
factors were tested:

$(1,1),(1,2),(1,4),(1,8),(1,16),(1,32),(1,64),$
$(1,\frac{1}{2}),(1,\frac{1}{4}),(1,\frac{1}{8}),(1,\frac{1}{16})$.

The selection of candidate values for the tuning parameters is often
arbitrary; however, the chosen combinations are based on the experience of
Boulesteix et al. With the availability of bigger computational resources, the
grid of combinations could be significantly refined. To choose the reference
parameter $\lambda_{1}$ and the best combination of penalty factors, 5-fold
cross-validation repeated 10 times is performed [21], and because of the
binary outcome, the error measure computed in the cross-validation is the
classification error. No constraint is set on the maximal number of predictors
allowed in the final model, i.e., the best model is selected based on CV
without restriction. For the standard LASSO, the parameter $\lambda$ is
determined using 5 folds cross validation repeated 10 times.
In stability selection, the two first thresholds are chosen at 0.7 and 0.8
(70% and 80% selection probability). In the third stability selection
configuration, we determine the threshold using the r-concave bound with V = 2
for the number of false positives, denoted as the “optimal threshold” for the
rest of the article. The value of this “optimal threshold” can vary across
simulation configurations but remains consistently between 0.5 and 0.55 in our
case. This also means that the two first thresholds lead to very strict error
control.
To compare the performance of the six procedures, we compared the number of
false positives in relation to the upper bound V = 2 . To evaluate the
effectiveness of the selection procedures, we calculated the true positive
proportion (TPP), also known as power.

###  III-B Simulation results for “independent data”

For the “independent data”, the results of the nine simulated designs applied
to the six procedures with the two criteria described in the previous section
are shown respectively in Figure 2 and in Figure 3. Figure 2 and Figure 3
represent the boxplot of all models.

Figure 2:  Boxplot of the true positive proportion for all 9 designs and all 6
models for the independent data. The models are presented from left to right
in the same order as in the legend (LASSO70: Lasso with threshold at 70%,
IPF_LASSO80: IPF-Lasso with threshold at 70%, LASSO80: Lasso with threshold at
80%, IPF_LASSO80: IPF-Lasso with threshold at 80%, LASSO_OPTI: Lasso with
optimal threshold, IPF_LASSO_OPTI: IPF-Lasso with optimal threshold) Figure 3:
Boxplot of the number of false positives for all 9 designs and all 6 models
for the independent data. The models are presented from left to right in the
same order as in the legend (LASSO70: Lasso with threshold at 70%,
IPF_LASSO80: IPF-Lasso with threshold at 70%, LASSO80: Lasso with threshold at
80%, IPF_LASSO80: IPF-Lasso with threshold at 80%, LASSO_OPTI: Lasso with
optimal threshold, IPF_LASSO_OPTI: IPF-Lasso with optimal threshold)

In all nine designs, all six procedures performed well in terms of controlling
the number of false positives that we set to V = 2. Overall, very few false
positives were selected at all in any of the experiments. This is no surprise,
as the bound is known to be conservative. Since there is no difference in
terms of false positives between the design and the procedures, it is
interesting to look at the power (TPP). Here, the results are varying. In
Design D, F and I, where the two datasets are very different in terms of the
proportions of non-zero variables, the IPF-LASSO demonstrates superior power
compared to regular LASSO. It is reasonable that when one of the modalities
contains more informative variables, the different penalties introduced by the
IPF-LASSO lead to higher power while maintaining error control. On the other
hand, regular LASSO seems to outperform the IPF-LASSO in Designs A, B and C,
where respectively, the two modalities, the proportion of truly relevant
variables, and the number of truly relevant variables are the same. These
results correspond to the results regarding prediction performance in the
original IPF-LASSO paper.

###  III-C Simulation results for “correlated data”

For the “correlated data”, the results of the nine simulated designs applied
to the six procedures with the two criteria described in Section 3.1 are shown
in Figure 4 and Figure 5. Figure 4 and Figure 5 represent the boxplot of all
models.

Figure 4:  Boxplot of the true positive proportion for all 9 designs and all 6
models for the correlated data. The models are presented from left to right in
the same order as in the legend (LASSO70: Lasso with threshold at 70%,
IPF_LASSO80: IPF-Lasso with threshold at 70%, LASSO80: Lasso with threshold at
80%, IPF_LASSO80: IPF-Lasso with threshold at 80%, LASSO_OPTI: Lasso with
optimal threshold, IPF_LASSO_OPTI: IPF-Lasso with optimal threshold) Figure 5:
Boxplot of the number of false positives for all 9 designs and all 6 models
for the correlated data. The models are presented from left to right in the
same order as in the legend (LASSO70: Lasso with threshold at 70%,
IPF_LASSO80: IPF-Lasso with threshold at 70%, LASSO80: Lasso with threshold at
80%, IPF_LASSO80: IPF-Lasso with threshold at 80%, LASSO_OPTI: Lasso with
optimal threshold, IPF_LASSO_OPTI: IPF-Lasso with optimal threshold)

The results for correlated data are very similar to the results obtained for
independent data in terms of variable selection. These six procedures do not
seem to be disrupted by the correlated data and they still perform well in
terms of error control. The power is generally lower compared to the situation
with independent data, as expected. Similar to the independent data scenario,
IPF-LASSO outperforms regular LASSO in Designs D, F, and I, exhibiting higher
power. In the correlated data context, IPF-LASSO also demonstrates superior
power compared to regular LASSO in two additional designs: G and H. In Design
H, there is little no variation in the results for the IPF-Lasso model with
the optimal threshold. The boxplot is concentrated on 0.5, surrounded by some
extremes values. In these two settings, the particular characteristics include
a low number of true active variables (less than 0.5%) in one of the two data
modalities (in Design H) and a large size difference between the two
modalities. It seems logical that when the size of the modalities differ
significantly, the different penalties introduced by IPF-LASSO will balance
the information given by both modalities. Giving them the same penalty could
cause the larger modality to be favored at the expense of the smaller one,
which could still contain relevant information.

In summary, in this simulation study, we discovered three patterns in which
the IPF-LASSO algorithm performs better than the regular LASSO. The first
pattern is when the difference in proportions of non-zero variables between
the two datasets is large. The bigger the difference, the better is the IPF-
LASSO as compared to the regular LASSO. The second pattern which benefits from
the IPF-LASSO is when the number of non-zero variables is small in all of the
modalities. The last pattern in favour of the IPF-LASSO is the one with two
modalities of very different sizes. The larger the difference between the the
sizes, the better is the performance of the IPF-LASSO as compared to the
regular LASSO. These three patterns are the same as those highlighted in the
original article by Boulesteix et al. for prediction on multi-omics data.

##  IV Real Data Application

###  IV-A Breast Cancer Data

Aligning with our idea to stay close to the work of Boulesteix et al. [3], we
use the same breast cancer data as in their article to illustrate our
suggested variable selection procedure. The data has two modalities: one low-
dimensional modality containing clinical data and one high-dimensional
modality with gene expression data. We are interested in classification of
patients with a minimal residual disease (RCB-I) vs. patients with a moderate
or extensive residual disease (RCB-II and RCB-III) after chemotherapy. The
data were initially analyzed by Hatzis et al. [19], and are publicly available
through the Gene Expression Omnibus repository (GSE25066). We decided not to
apply any preselection step, for two reasons: firstly, to make use all the
information provided by all variables, and secondly, to avoid introducing bias
by choosing among various preselection methods. The clinical set contains the
same variables as being used in the work of De Bin et al. [22]. The clinical
variables are:

  * •

Age (continuous)

  * •

Nodal status before treatment (4 categories): N0 (no lymph node involved), N1
(1-3), N2 (4–9), and N3 (10 and more).

  * •

Tumor size (4 categories): T1 (less than 2cm), T2 (between 2-5cm), T3 (more
than 5cm), and T4 (any size in the chest)

  * •

Grade (3 categories): Grade 1 (score 3,4,5), Grade 2 (score 6,7), and Grade 3
(score 8,9) [23]

  * •

Estrogen receptor (positive / negative)

  * •

Progesterone receptor (positive / negative)

For the analysis, the variables nodal status, tumor size and grade are
represented by dummy variables in the model. The categories N0 (nodal status
0), T1 (tumor size 1) and Grade 1 are used as reference categories.

The high-dimensional modality contains 22,283 microarray gene expression
measurements (Affymetrix-U133A GeneChip). More details about the data are
available in the original paper [19]. After removing patients with missing
data, the final dataset contains 385 patients, of which 277 are classified as
RCB-I.

###  IV-B Analysis results

We apply the six procedures presented previously with logistic regression on
the breast cancer data. As for the simulated data, the following combinations
of penalty factors were used for the IPF-LASSO:

$(1,1),(1,2),(1,4),(1,8),(1,16),(1,32),(1,64),$
$(1,\frac{1}{2}),(1,\frac{1}{4}),(1,\frac{1}{8}),(1,\frac{1}{16})$

and the best combination was chosen by cross-validation.
Among the six procedures, only the IPF model with the optimal stability
selection threshold (c=0.50) was able to select any variables. The upper bound
was in this case set to V = 5. The only variable selected was age, which is
not particularly interesting. However, it is well known that with
$\ell_{1}$-penalty, the presence of highly correlated variables can cause
problems for stability selection [24]. Adding a small $\ell_{2}$-penalty,
leading to the Elastic Net penalty[25], can alleviate this issue. Thus, our
six procedures were applied again to the real data with the Elastic Net
penalty (the mixing parameter alpha set to 0.7). Again, only the IPF model
with the optimal threshold for stability selection (c=0.5) was able to select
any variables. Only clinical variables were selected. This is of course not
very surprising, as the clinical variables are much “closer” to the outcome
(treatment response) compared to the gene expressions[19]. The variables
selected were Age, Estrogen Receptor, Nodal status, and Grade.

##  V Conclusion

This article focuses on an investigation of whether the IPF-LASSO method can
achieve increased statistical power while controlling false positives compared
to the standard LASSO approach in a multi-omics data settings. To our
knowledge, this work is the first to treat variable selection with error
control while integrating several different correlated data sources. Stability
selection was implemented in both IPF-LASSO and regular LASSO on a binary
outcome with the aim of controlling the number of false positives.

Simulation studies highlight that, in scenarios where there are differences
between the two data sources in terms of size of the data, quantity, or ratio
of true variables, the IPF-LASSO method shows increased statistical power
while controlling false positives compared to the conventional LASSO approach.
In these scenarios, the IPF-LASSO should be preferred over the regular LASSO.

On the cancer data, one clinical data modality and one omics data modality
were used to predict the binary outcome and only the IPF-LASSO model was able
to select any variables. All the selected variables were clinical variables.
In this particular situation, this shows the superiority of the IPF-LASSO,
which was capable of selecting variables even among the smallest dataset. The
fact that no gene expressions were selected could be explained by our choice
of response variable (treatment response), which is obviously closely related
to the clinical variables. However, we should be aware that our upper bound
for false positives selection was higher than the number of variables that was
actually selected.

While the IPF-LASSO appears advantageous for variable selection with false
positive control in multi-omics data, like all LASSO-type methods, this method
may not perform well when there are highly correlated predictors. An upgrade
of the IPF-LASSO could be to substitute the LASSO penalty with an Elastic Net
penalty, like we did in our applied example, see also Djordjilovic et al.
[26]. Comparing this modified “IPF-elastic-net” model to the standard Elastic-
net model for variable selection with false positive control in multi-omics
data settings holds potential interest and merit for further exploration. This
work is focused on controlling the number of false positives, but other types
of error control exist, like the False Discovery Rate (FDR). Looking into ways
of controlling FDR in the IPF-LASSO would be a natural extension of this work.
Finally, continued progress in computational and statistical methods is
crucial to achieve biologically meaningful integration of high-dimensional
multi-omics data [27].

## Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the
research, authorship, and/or publication of this article.

## Funding

The author(s) disclosed receipt of the following financial support for the
research, authorship, and/or publication of this article: This research is
funded by Norwegian Cancer Society, Grant no. 216137.

## References

  * [1] E. George, “The variable selection problem,” _Journal of the American Statistical Association_ , vol. 95, no. 452, pp. 1304–1308, 2000.
  * [2] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [3] A.-L. Boulesteix, R. De Bin, X. Jiang, M. Fuchs _et al._ , “Ipf-lasso: integrative-penalized regression with penalty factors for prediction based on multi-omics data,” _Computational and mathematical methods in medicine_ , vol. 2017, 2017.
  * [4] S. Klau, V. Jurinovic, R. Hornung, T. Herold, and A.-L. Boulesteix, “Priority-lasso: a simple hierarchical approach to the prediction of clinical outcome using multi-omics data,” _BMC bioinformatics_ , vol. 19, pp. 1–14, 2018.
  * [5] Z. Zhao and M. Zucknick, “Structured penalized regression for drug sensitivity prediction,” _Journal of the Royal Statistical Society Series C: Applied Statistics_ , vol. 69, no. 3, pp. 525–545, 2020.
  * [6] Z. Zhao, S. Wang, M. Zucknick, and T. Aittokallio, “Tissue-specific identification of multi-omics features for pan-cancer drug response prediction,” _Iscience_ , vol. 25, no. 8, 2022.
  * [7] D. Y. Ding, S. Li, B. Narasimhan, and R. Tibshirani, “Cooperative learning for multiview analysis,” _Proceedings of the National Academy of Sciences_ , vol. 119, no. 38, p. e2202113119, 2022.
  * [8] N. Meinshausen and P. Buhlmann, “Stability selection,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , vol. 72, no. 4, pp. 417–473, 2010.
  * [9] R. D. Shah and R. J. Samworth, “Variable selection with error control: another look at stability selection,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , vol. 75, no. 1, pp. 55–80, 2013.
  * [10] D. R. Cox, “A note on data-splitting for the evaluation of significance levels,” _Biometrika_ , vol. 62, no. 2, pp. 441–444, 1975.
  * [11] L. Wasserman and K. Roeder, “High dimensional variable selection,” _The Annals of Statistics_ , vol. 37, no. 5A, p. 2178, 2009.
  * [12] N. Meinshausen, L. Meier, and P. Buhlmann, “P-values for high-dimensional regression,” _Journal of the American Statistical Association_ , vol. 104, no. 488, pp. 1671–1681, 2009.
  * [13] C. Dai, B. Lin, X. Xing, and J. S. Liu, “False discovery rate control via data splitting,” _Journal of the American Statistical Association_ , pp. 1–18, 2022.
  * [14] D. G. Rasines and G. A. Young, “Splitting strategies for post-selection inference,” _Biometrika_ , vol. 110, no. 3, pp. 597–614, 2023.
  * [15] W. Fithian, D. Sun, and J. Taylor, “Optimal inference after model selection,” _arXiv preprint arXiv:1410.2597_ , 2014.
  * [16] R. Lockhart, J. Taylor, R. J. Tibshirani, and R. Tibshirani, “A significance test for the lasso,” _The Annals of Statistics_ , vol. 42, no. 2, p. 413, 2014.
  * [17] J. D. Lee, D. L. Sun, Y. Sun, and J. E. Taylor, “Exact post-selection inference, with application to the lasso,” _The Annals of Statistics_ , vol. 44, no. 3, pp. 907 – 927, 2016.
  * [18] E. Candes, Y. Fan, L. Janson, and J. Lv, “Panning for gold:‘model-x’knockoffs for high dimensional controlled variable selection,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , vol. 80, no. 3, pp. 551–577, 2018.
  * [19] C. Hatzis, L. Pusztai, V. Valero, D. J. Booser, L. Esserman, A. Lluch, T. Vidaurre, F. Holmes, E. Souchon, H. Wang _et al._ , “A genomic predictor of response and survival following taxane-anthracycline chemotherapy for invasive breast cancer,” _Journal of the American Medical Association_ , vol. 305, no. 18, pp. 1873–1881, 2011.
  * [20] J. Friedman, T. Hastie, and R. Tibshirani, “Regularization paths for generalized linear models via coordinate descent,” _Journal of statistical software_ , vol. 33, no. 1, p. 1, 2010.
  * [21] A.-L. Boulesteix, A. Richter, and C. Bernau, “Complexity selection with cross-validation for lasso and sparse partial least squares using high-dimensional data,” in _Algorithms from and for Nature and Life: Classification and Data Analysis_. Springer, 2013, pp. 261–268.
  * [22] R. De Bin, W. Sauerbrei, and A.-L. Boulesteix, “Investigating the prediction ability of survival models based on both clinical and omics data: two case studies,” _Statistics in medicine_ , vol. 33, no. 30, pp. 5310–5329, 2014\.
  * [23] A. C. Society, “About breast cancer,” https://www.cancer.org/content/dam/CRC/PDF/Public/8577.00.pdf, 2021, accessed: 2024.
  * [24] P. Kirk, A. Witkover, C. R. Bangham, S. Richardson, A. M. Lewin, and M. P. Stumpf, “Balancing the robustness and predictive performance of biomarkers,” _Journal of computational biology_ , vol. 20, no. 12, pp. 979–989, 2013.
  * [25] H. Zou and T. Hastie, “Regularization and variable selection via the elastic net,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , vol. 67, no. 2, pp. 301–320, 2005.
  * [26] V. Djordjilovic, E. Ponzi, T. H. Nost, and M. Thoresen, “penalizedclr: an r package for penalized conditional logistic regression for integration of multiple omics layers,” _arXiv preprint: 240201398_ , 2024.
  * [27] D. Wissel, D. Rowson, and V. Boeva, “Systematic comparison of multi-omics survival models reveals a widespread lack of noise resistance,” _Cell Reports Methods_ , vol. 3, no. 4, 2023.