%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% %% This is the mc_template.tex file for the mc document class. %% %% It is used to prepare s manuscript for Mathematical %% %% Communications journal. %% %% %% %% The mc.cls class works only with a pdflatex engine. %% %% The file newmc.cls should be placed where LaTeX %% %% can find it, e.g. in the current working directory. %% %% %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \documentclass{mc} %%===============================================================%% %% Please add here your own packages, macros and enviroments. %% %% It is not necessary to include ams* and graphicx packages %% %% since they are automatically included by the mc class. %% %% Avoid defining your own environments and use the already %% %% defined ones (e.g.~theorem, lemma etc.) %% %%===============================================================%% \newcommand{\BS}{\operatorname{BS}} \newcommand{\Cov}{\operatorname{Cov}} \renewcommand{\d}{\,\mathrm{d}} \newcommand{\dx}{\,\mathrm{d}x} \newcommand{\dv}{\,\mathrm{d}v} \newcommand{\e}{\mathrm{e}} \newcommand{\D}{\mathbb{D}} \newcommand{\E}{\mathbb{E}} \newcommand{\F}{\mathcal{F}} \renewcommand{\H}{\mathcal{H}} \newcommand{\K}{\mathfrak{K}} \newcommand{\N}{\mathbb{N}} \renewcommand{\P}{\mathbb{P}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\R}{\mathbb{R}} \newcommand{\Var}{\operatorname{Var}} \newcommand{\1}{\mathbf{1}} \newcommand{\phy}{\varphi} \newcommand{\eps}{\varepsilon} % \usepackage{refcheck} %%===============================================================%% %%===============================================================%% %% Journal info will be edited by the typesetter %% %% DO NOT CHANGE THIS PART %% %%===============================================================%% \setcounter{page}{1} \renewcommand\thisnumber{x} \renewcommand\thisyear {2026} \renewcommand\thismonth{xxx} \renewcommand\thisvolume{31} \renewcommand\datereceived{December 24, 2024} \renewcommand\dateaccepted{October 6, 2025} \renewcommand\doinum{10.1000/100} %%===============================================================%% % 24.12.2024. % 6.10.2025. \begin{document} % \begin{linenumbers} %%===============================================================%% %% TITLE %% %% Please add the title with \title[Short title]{Title} %% %% Short title is a running head apearing in the header. %% %%===============================================================%% \title[Lower bound of density of SDDE driven by fBm] % at most 50 characters including spaces {Gaussian lower bound and positivity of the density of stochastic delay differential equations driven by a fractional Brownian motion} % at most 150 characters including spaces () %%===============================================================%% %%===============================================================%% %% AUTHOR(S) %% %% %% %% Add author's details in the following format. For each author %% %% provide the affiliation, address and Orcid identifier. %% %% Mark the corresponding author with \comma\corrauth. %% %%===============================================================%% \author[Ò.~Burés and C.~Rovira] % put here short author names for header {Òscar Burés\affil1\comma\corrauth\orcidnumber{0009-0003-1741-1308} and Carles Rovira\affil2\orcidnumber{0000-0001-9021-9804} } \address{\affilnum1 Departament de Matemàtica Econòmica, Financera i Actuarial, Facultat d'Economia i Empresa, Universitat de Barcelona, 08\,034 Barcelona, Spain \\ \affilnum2 Departament de Matemàtiques i Informàtica, Facultat de Matemàtiques i Informàtica, Universitat de Barcelona, 08\,007 Barcelona, Spain } \emails{% \email{oscar.bures@ub.edu} (Ò.~Burés), \email{carles.rovira@ub.edu} (C.~Rovira) } %% For single author use the following format. %% %\author[F.~Author]% % {First Author\orcidnumber{0000-0000-0000-0000} % } % %\address{Affiliation 1} % %\emailsingle{% % \email{fauthor@mathos.hr} % } %%===============================================================%% %%===============================================================%% %% ABSTRACT %% %%===============================================================%% \begin{abstract} In this paper, we prove that the density of a stochastic delay differential equation driven by a fBm with Hurst parameter $H > 1/2$ is strictly positive combining Nourdin-Viens' and Kohatsu-Higa's method. \end{abstract} %%===============================================================%% %%===============================================================%% %% KEYWORDS %% %%===============================================================%% \keywords{Stochastic delay equation; Malliavin calculus; fractional Brownian motion; estimates of the density} %%===============================================================%% %%===============================================================%% %% AMS subject classification %% %%===============================================================%% \ams{60H10; 60H07; 60H05} %%===============================================================%% %%===============================================================%% \maketitle %%===============================================================%% %%===============================================================%% %% MAIN BODY %% %%===============================================================%% \section{Introduction} Let $\Omega = \mathcal{C}_0([0,T])$ represent the space of continuous functions on the interval $[0,T]$ that vanish at zero, equipped with the supremum norm, thus making it a Banach space. We consider the Borel $\sigma$-field $\mathcal{F} := \mathcal{B}(\Omega)$ and a probability measure $P$ under which the canonical process defined by $B^H_t(\omega) = \omega(t)$ is a fractional Brownian motion (fBm) with Hurst parameter $H > 1/2$. Throughout this paper, we will work in this underlying probability space $(\Omega, \mathcal{F}, P)$. A fractional Brownian motion \(B^H = \{B^H_t; t \in [0,T]\}\) is a centered Gaussian process with covariance function given by \[ R_H(t,s) := \E[B_t^H B_s^H] = \frac{|t|^{2H} + |s|^{2H} - |t - s|^{2H}}{2}. \] Although this process does not exhibit independent increments, its increments are stationary. This feature, together with the Gaussianity of $B^H$, allows us to establish the inequality \[ \E(|B^H_t - B^H_s|^p) \leq C |t - s|^{pH}, \] which, by Kolmogorov's continuity criterion, indicates that $B^H$ has $\gamma$-Hölder continuous paths for any $\gamma \in (0,H)$. This paper focuses on stochastic delay differential equations (SDDEs) driven by a fractional Brownian motion with $H > 1/2$, specifically following the form introduced by Ferrante and Rovira in \cite{ferrante2006stochastic}: \begin{equation} \label{eq: main equation} \begin{cases} X_t = \eta_0 + \int_0^t \sigma(X_{s - r}) \, dB^H_s + \int_0^t b(X_s) \, ds, \quad t \in (0,T], \\ X_t = \eta_t, \quad t \in [-r,0], \end{cases} \end{equation} where \(\eta\) is a smooth function, \(r > 0\) is a constant delay parameter and $B^H$ is a fBm of Hurst parameter $H > 1/2$. Given the regularity properties of $B^H$ for $H > 1/2$, the stochastic integral in this equation can be interpreted as a Riemann-Stieltjes integral, leveraging Young's results in \cite{Young}. Furthermore, Zähle's framework for fractional calculus \cite{Zhle1998IntegrationWR} allows us to represent this integral as a Lebesgue integral via an integration by parts formula. In 2002, Nualart and Rascanu established in their work \cite{rascanu2002differential} the existence and uniqueness of solutions for general stochastic differential equations (SDEs) driven by a fBm with Hurst parameter $H > 1/2$. Their proof can be easily adapted to prove the existence and uniqueness of solutions for stochastic delay differential equations of the form \eqref{eq: main equation}. In 2006, Ferrante and Rovira presented in \cite{ferrante2006stochastic} an alternative approach to proving the existence and uniqueness of solutions to equation \eqref{eq: main equation}. They utilized an inductive strategy that specifically leveraged the delay property, thereby providing a different perspective from that of Nualart and Rascanu. Additionally, Ferrante and Rovira established the existence and smoothness of the solution's density using methods tailored for the delayed framework. A significant advancement in understanding the law of general SDEs was made by Nualart and Saussereau in 2009. In \cite{nualart2009malliavin}, they proved the existence of the density by connecting Malliavin differentiability with Fréchet differentiability. Again, the proof presented for the general framework in \cite{nualart2009malliavin} can also be straightforwardly adapted to SDDEs. In 2012, León and Tindel proved in \cite{leontindel} the smoothness of the density of solutions to SDDEs using rough path techniques, arriving at the same conclusion as Rovira and Ferrante in \cite{ferrante2006stochastic} for $H > 1/2$ using a different method. In this paper, as a starting point we take the results of \cite{ferrante2006stochastic} and \cite{leontindel}, which ensure that the density of the solution to \eqref{eq: main equation} is smooth, and we aim to prove that under smoothness conditions on the coefficients this density function is strictly positive. In order to prove this, we derive a Gaussian-type lower bound for the density of the solution to a stochastic delay differential equation combining two methods widely used in the literature: the Nourdin and Viens method (see \cite{NourdinViens}), which is the most frequently used method for bounding the density of solutions to stochastic differential equations with additive noise and Kohatsu-Higa's method (see \cite{kohatsu2003lower}), which is widely used for bounding density functions of solutions to stochastic differential equations where the noise is not additive but the conditional law of the solution with respect to the $\sigma$-algebras generated by the noise driving the equation can be compared to a Gaussian distribution. The advantage of delay equations lies in their structure, which provides greater flexibility when applying different methods to bound the density of their solutions. In a general setting of stochastic differential equations (SDEs), such flexibility is typically not available, and in some cases, it may not even be possible to find a method that yields satisfactory results. SDDEs of this nature are significant in various fields, including mathematical finance and biological modeling, where delays in data and external noise—modeled as fractional noise—play a crucial role. For instance, Arriojas et al. \cite{arriojas2006delayed} examined financial models driven by SDDEs, while similar frameworks are utilized in biology to account for delayed feedback influenced by noisy environments. Our main result, stated formally in Theorem \ref{th: main result}, provides a lower bound for the density function of the solution to the SDDE. Specifically, we prove that under certain regularity and ellipticity conditions on \(\sigma\) and \(b\), the density \(p_t(x)\) of the solution to equation \eqref{eq: main equation} satisfies: \begin{equation} \label{eq: bound introduction 1} p_t(x) \geq \frac{\E(|X_t-m_t|)}{c_1t^{2H}}\exp\left( -\frac{c_2(x-m_t)^2}{t^{2H}}\right), \end{equation} for $t \in (0,r]$ and \begin{equation} \label{eq: bound introduction 2} p_t(x) \geq \frac{c_3}{t^H} \exp \left( -\frac{c_4(x - \eta_0)^2}{t^{2H}} \right), \end{equation} for $t \in (r,T]$, where \(m_t = \E[X_t]\). This result implies the strict positivity of the density as an immediate corollary. The paper is organized as follows: Section \ref{sec: preliminaries} introduces the Malliavin calculus tools and techniques used to get to the conclusion of this work. In section \ref{sec: lower bound} we prove that the density function of the solution to \eqref{eq: main equation} satisfies bounds \eqref{eq: bound introduction 1} and \eqref{eq: bound introduction 2}. In the same section, we also state the corollary that concludes the work presented in this paper. Finally, there is an appendix containing two auxiliary results of real analysis that are key to derive the Gaussian-type lower bound of the density $p_t(x)$ of the solution to \eqref{eq: main equation}. \section{Preliminaries} \label{sec: preliminaries} \subsection{Wiener space associated to the fractional Brownian motion} We define $\mathcal{E}$ as the space of step functions of the form \[ s(t) = \sum_{j=1}^n a_j \1_{[0,t_j]}(t), \quad a_j \in \R, t_j \in [0,T] \text{ for all } j = 1, \dots, n \] with respect to the inner product defined by \begin{equation} \label{eq: inner product} \langle \1_{[0,t]}, \1_{[0,s]} \rangle_{\mathcal{H}} = R_H(t,s), \end{equation} where $R_H$ is the covariance function of the fBm. This inner product extends to $\mathcal{E}$ by linearity. We define $\mathcal{H}$ as the closure of $\mathcal{E}$ with respect to the inner product \eqref{eq: inner product}. Moreover, for $\phy, \psi \in \mathcal{H}$, the inner product \eqref{eq: inner product} can be written as \[ \langle \phy, \psi\rangle_{\mathcal{H}} = \alpha_H \int_0^T \int_0^T \phy_u \psi_v|u-v|^{2H-2} du dv, \] where \[ \alpha_H = H(2H-1). \] Consider the kernel $K_H : [0,T]^2 \to \R$ defined by \[ \begin{split} K_H(t,s) &= c_H s^{1/2 - H}\int_s^t (u-s)^{H-3/2}u^{H-1/2} du, \\ c_H &= \sqrt{\frac{H(2H-1)}{\beta(2-2H, H-1/2)}} \end{split} \] for $t > s$ and $K_H(t,s) = 0$ in the case $t \leq s$. This kernel induces an isometry between $\mathcal{H}$ and $L^2([0,T])$. Indeed, the operator $K^*_{H}$ defined by \[ (K_H^* \phy)_s = \int_s^T \phy_t\partial_t K_H(s,t) dt \] satisfies that, for every $\phy, \psi \in \mathcal{H}$, \[ \langle \phy, \psi \rangle_{\mathcal{H}} = \langle K^*_H \phy, K_H^* \psi \rangle_{L^2}. \] Moreover, if $W = \{W_t; t \in [0,T]\}$ denotes a standard Brownian motion and we fix $H > 1/2$, then the Volterra process \[ B^H_t := \int_0^t K_H(t,s) dW_s \] is a fractional Brownian motion of Hurst index $H$. Abusing the notation, we also define $B^H := \{ B^H(h); h \in \mathcal{H}\}$ the isonormal process associated to $\mathcal{H}$, that is, a centered Gaussian process indexed by $\H$ with covariance function $\E[B^H(\phy) B^H(\psi)] = \langle \phy, \psi \rangle_{\H}$. Notice that the process $B^H_t := B^H(\1_{[0,t]})$ is, again, a fractional Brownian motion of Hurst index $H$. We can represent $B^H(h)$ as \[ B^H(h) = \int_0^T h_t dB^H_t, \] where this last integral can be understood in the Young sense (see \cite{Young}). \subsection{Malliavin calculus tools} \begin{definition} We define $\mathcal{S}$ as the space of random variables $F$ of the form \[ F = f(B^H(h_1), \dots, B^H(h_n)), \] where $n \geq 1$ and $f \in \mathcal{C}^{\infty}_p(\R^n)$ (i.e. $f$ is a real-valued smooth function and $f$, together with all its partial derivatives, has at most polynomial growth). We say that $\mathcal{S}$ is the space of cylindrical random variables. \end{definition} We can define the Malliavin derivative in this class of random variables $\mathcal{S}$. \begin{definition} For $F \in \mathcal{S}$, we define $DF$ as the $\mathcal{H}$-valued random variable of the form \[ DF = \sum_{j=1}^n \partial_j f (B^H(h_1), \dots, B^H(h_n))h_j. \] Given $h \in \mathcal{H}$, we can also define the derivative of $F$ in the direction of $h$ as \[ D_h F = \langle DF, h\rangle_{\mathcal{H}}. \] We also define the process $\{D_tF; t \in [0,T]\}$ as \[ D_tF = \sum_{j=1}^n \partial_j f (B^H(h_1), \dots, B^H(h_n))h_j(t). \] \end{definition} In the same way in which we have introduced first order derivatives, we can define derivatives of order $k \geq 2$ as follows: for $t_1, \dots t_k \in [0,T]$, we define $D^{(k)}_{t_1, \dots, t_k}$ as the the random variable \[ D^{(k)}_{t_1, \dots, t_k} F = D_{t_1} \cdots D_{t_k}F. \] Notice that the random variable $D^{(k)}F$ is therefore an $\H ^{\otimes k}$-valued random variable. For every $k \geq 1$ and $p \geq 1$, the operators $D^{(k)}$ are closable from $\mathcal{S}$ to $L^p(\Omega; \mathcal{H}^{\otimes k}$). We denote by $\D^{k,p}$ the closure of $\mathcal{S}$ with respect to the norm \[ ||F||_{k,p} = \left( \E[|F|^p] + \sum_{j=1}^k \E[||D^{(j)}F||^p_{\mathcal{H}^{\otimes j}}]\right)^{1/p}. \] Notice that $\D^{k, p} \subset \D^{l, p}$ if $l\leq k$ and $\D^{k,p} \subset \D^{k,q}$ if $q \leq p$. We also define the space $\D^{\infty}:= \cap_{k \geq 1} \cap_{p \geq 1} \D^{k,p}$. The same Malliavin-Sobolev spaces can be defined for an isonormal process based on a standard Brownian motion $W$. In such case, we will write $\D^{k,p}_{W}$ to stress the underlying isonormal process. For a random variable $F \in \D^{1,2}$, we define $\Gamma_F$ and $\Gamma_F^{-1}$ (provided $||DF||^2_{\mathcal{H}} > 0$ a.s.) as \[ \Gamma_F = ||DF||^2_{\mathcal{H}}, \quad \Gamma_F^{-1} = \left( ||DF||^2_{\mathcal{H}}\right)^{-1}. \] In the same way in which we have introduced the Malliavin derivative in the classical Malliavin calculus setting, we need to introduce the conditional Malliavin calculus. To this end, let $\mathbb{F} =\{ \mathcal{F}_t; t \in [0,T]\}$ be the natural filtration associated to $B^H$. For a given $t \in [0,T]$ and a random variable $F \in L^1(\Omega)$, we set \[ E_t[F] = E[F |\mathcal{F}_t]. \] As in the classical Malliavin calculus, we denote by $||F||_{k,p,t}$ and $\Gamma_{F,t}$ the following objects \begin{equation} \label{eq: def conditional quantities BH} ||F||_{k,p,t} = \left( E_t[|F|^p] + \sum_{j=1}^k E_t[||D^{(j)}F||^p_{\mathcal{H}[t,T]^{\otimes j}}]\right)^{1/p}_{,} \quad \Gamma_{F,t} = ||DF||^2_{\mathcal{H}[t,T]}, \end{equation} for $t \in [0,T]$, where for $s, t \in [0,T]$ with $s\leq t$ and $\phy \in \mathcal{H}$, \[ ||\phy||^2_{\mathcal{H}[s, t]} = \int_s^t \int_s^t \phy_u \phy_v |u-v|^{2H-2} dudv. \] We are now interested in an integration by parts formula similar to the one found in \cite[Proposition 2.1.4]{nualart2006malliavin}. For our purposes, we will use a conditional version of this formula. The result is stated in the Wiener process $W$ in the $[0,T]$ framework. Since our aim is a conditional integration by parts formula for the fBm case, the first thing we shall show in order to make sense of this result is that the Malliavin derivative with respect to $W$ and the Malliavin derivative with respect to $B^H$ are related in some way. The relation between these two operators can be stated as follows: \begin{proposition} \label{relation derivatives} Let $\D^{1,2}_W$ be the Malliavin-Sobolev space associated to $W$. Then $\D^{1,2} = (K_H^*)^{-1} \D^{1.2}_W$ and for every $F \in \D_W^{1,2}$, we have \[ D^WF = K_H^* DF, \] whenever both sides of the equality are well defined. \end{proposition} \begin{proof} See \cite[Proposition 5.2.1]{nualart2006malliavin}. \end{proof} In the same way in which we have defined the objects \eqref{eq: def conditional quantities BH} in the fractional Brownian motion case, we can define them in the standard Brownian motion case. When working in this framework, we define: \begin{equation} ||F||_{k,p,t, W} = \left( E_t[|F|^p] + \sum_{j=1}^k E_t[||(D^W)^jF||^p_{L^2[t,T]^{\otimes j}}]\right)^{1/p}, \quad \Gamma_{F,t}^{W} = ||D^W F||^2_{L^2[t,T]}. \end{equation} Now that we know that there exists a relation between the derivative operators with respect to $W$ and $B^H$, we can state the conditional integration by parts formula for $W$ and the result for $B^H$ will be a direct consequence of applying Proposition \ref{relation derivatives}. The result for the Wiener process $W$ is stated in \cite[Proposition 2.1.4]{nualart2006malliavin}. Its conditional version (see, for instance, \cite{kohatsu2003lower} or \cite{besalu2016gaussian}) can be formulated as follows: \begin{proposition} \label{p: IBP} Fix $n \geq 1$. Suppose $F$ is a non-degenerate random variable; let $Z \in \D_{W}^{\infty}$ be an $\mathcal{F}_s$-measurable random variable such that $F +Z$ is a non-degenerate random variable and let $G \in \D_{W}^{\infty}$. We denote by $\1^n:=(1,\dots, 1) \in \R^n$. For any function $g \in \mathcal{C}_p^{\infty}(\R)$ there exists a random variable $H^s_{(\1^n)}(F,G)$ such that \begin{equation} \label{eq: IBP} E_{s}[g^{(n)}(F + Z) G] = E_s[g(F+Z) H^s_{(\1^n)}(F,G)], \end{equation} where $H^s_{(\1^n)}(F,G)$ is defined recursively by \[ H^s_{(1)}(F,G) = \delta^W_s(G (\Gamma^{-1}_{F,s})DF) \] and \[ H^s_{(\1^n)}(F,G) = H^s_{(1)}(F, H^s_{(\1^{n-1})}(F,G)). \] Here, $\delta^W_s$ denotes the Skorohod integral in the interval $[s,T]$. Moreover, for $q_1, q_2, q_3$ such that $\frac{1}{p} = \frac{1}{q_1} + \frac{1}{q_2} + \frac{1}{q_3}$, the following estimate holds: \begin{equation} \label{eq: IBP estimate} ||H^s_{(\1^n)}(F,G)||_{p,s, W} \leq c ||(\Gamma^W_{F, s})^{-1}||_{2^{n-1}q_1, s, W}^n||F||^{2(n+1)}_{n+2, 2^n q_2, s, W} ||G||_{n, q_3, s, W}. \end{equation} \end{proposition} Some remarks can be deduced from this proposition. \begin{remark} \label{remark 1} Notice that using an approximation argument, the conclusion of Proposition \ref{p: IBP} holds for $g(F+Z) = \1_{\{F+Z > x\}}$. Indeed, consider a sequence $\{g_k; k \geq 1\}$ of smooth compactly supported functions that converge to de Dirac's delta centered at a point $x \in \R$. Then, using the dominated convergence theorem, it is clear that the right-hand side of \eqref{eq: IBP} converges to \[ E_s[\1_{\{F+Z > x\}}H_{(\1^n)}^s(F,G)], \] while, concerning the left-hand side, the limit is \[ E_s[g^{(n)}(F+Z)G], \] where $g^{(n)}$ is the $n$-th order distributional derivative of $g$. \end{remark} \begin{remark} \label{remark 2} It is interesting to study (and will be useful in the future) how the integration by parts formula changes with rescalings of $F$. More precisely, given $\mu \in \R$, we want to see how estimate \eqref{eq: IBP estimate} changes when we consider $\mu F$ instead of $F$. We will prove by using induction on $n$ that \begin{equation*}% \label{eq: IBP reescaling} H^s_{(\1^{n})}(\mu F, G) = \frac{1}{\mu^n}H^s_{(\1^n)}(F,G). \end{equation*} First, we cover the case $n = 1$. Observe that \[ (\Gamma^{W}_{\mu F,s})^{-1} = (||\mu D^W F||^2_{L^2[s,T]})^{-1} = \mu^{-2} (\Gamma^W_{F,s})^{-1} \] and obviously $D^W\mu F = \mu D^W F$. Then using the linearity of $\delta^W_s$ we have \[ H^s_{(1)}(\mu F, G) = \frac{1}{\mu} H^s_{(1)}(F,G). \] Observe also that $H^s_{(1)}(F, \mu G) = \mu H^s_{(1)}(F, G)$. Assume that it holds for $n-1$. For $n$ we have \[ \begin{split} H^s_{(\1^{n})}(\mu F, G) = &H_{(1)}(\mu F, H^s_{(\1^{n-1})}(\mu F, G)) \\ = & \frac{1}{\mu}H_{(1)}\left( F, \frac{1}{\mu^{n-1}} H^s_{(\1^{n-1})}(F, G)\right) \\ = & \frac{1}{\mu^n}H_{(1)}\left( F, H^s_{(\1^{n-1})}(F, G)\right) \\ = & \frac{1}{\mu^n}H^s_{(\1^n)}(F,G). \end{split} \] This proves the desired scaling property of $H^s_{(\1^n)}$. \end{remark} Notice that from the relation given by Proposition \ref{relation derivatives}, the formula still holds for the Malliavin spaces associated to the fractional Brownian motion with the change of the underlying Hilbert space and shifting the random variables via $K_H^*$. Hence, Proposition \ref{p: IBP} also holds in the fractional Brownian motion framework. \section{Lower bound for the density} \label{sec: lower bound} The objective of this section is to give a proof of the lower bound of the density of $X_t$. The main result concerning this bound is encapsulated in the following theorem: \begin{theorem} \label{th: main result} Let $X_t$ be the solution to equation \eqref{eq: main equation} with $\sigma, b \in \mathcal{C}^{\infty}_b(\R)$, $\eta \in \mathcal{C}^{\infty}_b((-r,0))$ (i.e. $\sigma$, $b$ and $\eta$ are real-valued smooth functions and these functions are, together with all their partial derivatives, bounded), and moreover, there exist two constants $0 < \lambda < \Lambda$ such that \[ \lambda \leq \inf_{x \in \R} \sigma(x) \leq \sup_{x \in \R} \sigma(x) \leq \Lambda.\] Then, for every $t \in (0,T]$, the density function $p_t(x)$ of $X_t$ satisfies \begin{equation} \label{eq: bound main result 0r} p_t(x) \geq \frac{\E(|X_t-m_t|)}{c_1 t^{2H}}\exp \left( - \frac{c_2(x-m_t)^2}{t^{2H}}\right) \end{equation} if $t \in (0,r]$, and \begin{equation} \label{eq: bound main result rT} p_t(x) \geq \frac{c_3}{t^H}\exp \left( - \frac{c_4(x-\eta_0)^2}{t^{2H}}\right) \end{equation} if $t \in (r,T]$ for every $x \in \R$, where $c_1, c_2, c_3, c_4 > 0$ are real constants and $m_t = \E[X_t]$. \end{theorem} \begin{remark} Notice that since $\sigma$ and $b$ are assumed to be bounded, then one can show that $m_t = \E[X_t] < \infty$ using direct estimates. \end{remark} \begin{remark} The result is stated for one-dimensional stochastic delay differential equations. The generalization to the multidimensional case is straightforward (see, for instance, \cite{besalu2016gaussian} for an example of a multidimensional case). \end{remark} \begin{remark} Notice that this result is highly dependent on the fact that the delay $r > 0$ is constant. For another types of delay, another strategy has to be considered in order to prove the lower bound for the density. \end{remark} \begin{remark} The derivation of bound \eqref{eq: bound main result 0r} relies heavily on the restriction $t \in (0,r]$ and the proof of \eqref{eq: bound main result rT} is simplified by assuming $t > r$. Nevertheless, the method used for \eqref{eq: bound main result rT} can in fact be extended to all $t > 0$ (see \cite{besalu2016gaussian} and \cite{kohatsu2003lower} for some applications of the method in this scenario), thereby providing results valid on the whole interval $(0,T]$. Consequently, for the region $t \in (0,r]$, one can obtain two different bounds: the one specific to this interval, as in \eqref{eq: bound main result 0r}, and the one arising from the extension of the second method. \end{remark} Throughout the section, we will provide the auxiliary results needed to prove Theorem \ref{th: main result}. As mentioned, in order to prove this result, we will make a distinction in the cases $(0,r]$ and $(r, T]$. For the case $(0,r]$, we will make use of the techniques developed in \cite{NourdinViens}, as it is the most comfortable method to study the density of the solution $X_t$ when $t \in (0,r]$ because of the structure of the equation in this interval of time. \subsection{The case $t \in (0,r]$} \subsubsection{A general bounding technique} In order to illustrate how to obtain a bound for the density in this case, we will briefly recall the method developed by Nourdin and Viens in \cite{NourdinViens}. The bound relies on the following results: \begin{theorem} \label{th: g_F} Let $F \in \D^{1,2}$ with zero mean, let $g_F(x)$ be the function defined as \[ g_F(x) = E\left[\langle DF, -DL^{-1}F\rangle_{\mathcal{H}} | F = x\right], \] where $L$ denotes the Ornstein-Uhlenbeck operator associated to $B^H$ (see \cite[Section 1.4]{nualart2006malliavin} for a full discussion about the operator $L$). The law of $F$ has a density $p_F(x)$ if and only if the random variable $g_F(F)$ is strictly positive almost surely. In this case, the support of $p_F$, $supp(p_F$), is a closed interval of $\R$ of the form $[A,B]$ with $-\infty \leq A < B \leq \infty$ containing zero and, for almost $x \in supp(p_F)$, \[ p_F(x) = \frac{\E[|F|]}{2g_F(x)} \exp \left(-\int_0^x \frac{z}{g_F(z)}dz\right). \] \end{theorem} This result has the following consequence. \begin{corollary} \label{c: corol g_F} If there exist $\sigma_{min}, \sigma_{max} > 0$ such that \[ \sigma_{min}^2 \leq g_F(F) \leq \sigma^2_{max}, \] then $F$ has density function satisfying \[ \frac{\E[|F|]}{2\sigma^2_{max}}\exp\left( -\frac{x^2}{2\sigma^2_{min}} \right) \leq p_F(x) \leq \frac{\E[|F|]}{2\sigma^2_{min}}\exp\left( -\frac{x^2}{2\sigma^2_{max}} \right). \] \end{corollary} Even though this result gives us the lower and the upper bound, it is not clear from the definition of $g_F$ how to analyze this function. Indeed, we have to keep in mind that in our case, $F$ will be related to the solution of an SDDE, so computing $DL^{-1}F$ seems, a priori, a difficult task. Concerning our case, let $B^H$ be a fractional Brownian motion with Hurst parameter $H> 1/2$. As an abuse of notation, we will denote by $B^H$ the isonormal process asociated to $\mathcal{H}$. In \cite{NourdinViens}, the authors obtain a formula for $g_F(x)$ which avoids the computation of $-DL^{-1}F$. The result from \cite{NourdinViens} concerning the computation of $g_F$ is the following. \begin{proposition} \label{p: reformulation g_F} Assume $DF = \Phi_F(B^H)$ for a measurable function $\Phi_{F} : \R^{\mathcal{H}} \to \mathcal{H}$. Then, we have \[ \langle DF, -DL^{-1}F \rangle_{\mathcal{H}} = \int_{0}^{\infty} e^{-\theta} \langle \Phi_F(B^H), \E' \left(\Phi_F(e^{-\theta}B^H + \sqrt{1-e^{-2\theta}}B^{H'}) \right)\rangle_{\mathcal{H}} d\theta \] and, therefore, \[ g_F(F) = \int_{0}^{\infty} e^{-\theta} \mathbf{E}\left(\langle \Phi_F(B^H), \Phi_F(e^{-\theta}B^H + \sqrt{1-e^{-2\theta}}B^{H'}) \rangle_{\mathcal{H}}|F\right) d\theta, \] where $B^{H'}$ stands for an independent copy of $B^H$ such that $B^H$ and $B^{H'}$ are defined in the probability space $(\Omega \times \Omega', \mathcal{F} \otimes \mathcal{F}', P \times P')$, $\E'$ denotes the expectation with respect to $P'$ and $\mathbf{E}$ denotes the expectation with respect to $P \times P'$. \end{proposition} We will rely heavily on this way of computing $g_F$ in order to prove the following bounds for $t\in (0,r]$: \begin{proposition} \label{p: density in (0,r]} Let $X$ be the solution \eqref{eq: main equation} with $\sigma, b$ satisfying the hypothesis of Theorem \ref{th: main result}. If $\sigma$ is $(\lambda, \Lambda)$-elliptic, then for every $t \in (0,r]$, there exist constants $c_1 < c_2$ such that the density $p_t(x)$ of $X_t$ satisfies: \[ \begin{split} \frac{c_1\E[|X_t-m_t|]}{2\Lambda^2 t^{2H}}\exp\left( -\frac{(x-m_t)^2}{2\lambda^2t^{2H}}\right) \leq &p_t(x) \\ \leq &\frac{c_2\E[|X_t - m_t|]}{2\lambda^2 t^{2H}}\exp\left( -\frac{(x-m_t)^2}{2\Lambda^2 t^{2H}}\right), \end{split} \] where $m_t = \E[X_t]$. \end{proposition} \begin{proof} We will apply Corollary \ref{c: corol g_F}. Notice that for $t \in (0,r]$, $X_t$ solves the equation \[ X_t = \eta_0 + \int_0^t \sigma(\eta(s-r))dB^H_s + \int_0^t b(X_s) ds. \] Hence, the Malliavin derivative of $X_t$ in the direction $s < t$ satisfies the following equation: \[ D_sX_t = \sigma(\eta(s-r)) + \int_0^t \sigma'(\eta(u-r)) D_s\eta(u-r) dB^H_u + \int_0^t b'(X_u) D_sX_u du. \] Since $\eta$ is deterministic, $D_s\eta(u-r) = 0$, so the equation satisfied by $D_sX_t$ is reduced to \[ D_sX_t = \sigma(\eta(s-r)) + \int_0^t b'(X_u)D_sX_u du. \] This is an ODE with initial condition $\sigma(\eta(s-r))$, so the explicit solution to this equation is \[ D_sX_t = \sigma(\eta(s-r)) \exp\left( \int_s^t b'(X_u) du \right). \] Moreover, from the fact that $|b'(X_u)| \leq ||b'||_{\infty}$ and the ellipticity condition on $\sigma$ we deduce that there exists $M > 0$ such that \begin{equation} \label{eq: der0r} \lambda e^{-Mr} \leq D_sX_t \leq \Lambda e^{Mr}. \end{equation} The important conclusion of this bound is that \eqref{eq: der0r} holds uniformly with respect to $\omega \in \Omega$. Notice that one has that $X_t \in \D^{1,2}$, but Theorem \ref{th: g_F} (and its reformulation in Proposition \ref{p: reformulation g_F}) works for centered random variables. Hence, in order to apply the bounding technique, we will prove the result for the centered random variable $F = X_t - \E(X_t)$ and we will deduce from there the result for $X_t$. Since $\E(X_t)$ is a real number, $DF = DX_t$. Moreover, recall that if $\phy, \psi \in \mathcal{H}$, then \[ \langle \phy, \psi\rangle_{\mathcal{H}} = \int_0^T \int_0^T \phy_u\psi_v|u-v|^{2H-2} du dv. \] Hence, if we define $F^{\theta} := F(e^{-\theta}\omega + \sqrt{1-e^{-2\theta}}\omega')$, we can write $g_F$ as \[ g_F(F) = \int_{0}^{\infty} e^{-\theta} \E \left[\E' \left(\int_0^t \int_0^t D_uX_t (D_vX_t)^{\theta} |u-v|^{2H-2} du dv | F \right)\right] d\theta. \] Indeed, since we are working under the canonical space of the fractional Brownian motion, one has that \[ \Phi_F(e^{-\theta} B^H(\omega) + \sqrt{1-e^{-2\theta}}B^{H'}(\omega')) = \Phi_F(e^{-\theta} \omega + \sqrt{1-e^{-2\theta}}\omega'). \] Hence, if $DF(\omega) = \Phi_F(\omega)$, then $$\Phi_F(e^{-\theta} \omega + \sqrt{1-e^{-2\theta}}\omega') = DF(e^{-\theta} \omega + \sqrt{1-e^{-2\theta}}\omega') =: DF^{\theta}(\omega, \omega').$$ Using the bounds \eqref{eq: der0r} we can easily find that there exist $c_1, c_2 > 0$ such that \[ c_1 \lambda^2 t^{2H} \leq g_F(F) \leq c_2 \Lambda^2 t^{2H}, \] which finishes the proof. \end{proof} Notice that renaming the constants, we derive bound (\ref{eq: bound main result 0r}). Some comments about this method are that, as the reader can observe, it is extremely comfortable to study stochastic differential equations driven by an additive noise and, under a suitable hypothesis on $\sigma$, the same arguments work for equations of type \eqref{eq: main equation} for $t \in (0,r]$. The fact that for $t > r$ the diffusion coefficient is random makes it impossible to keep applying this same method. \subsection{The case $t > r$} As mentioned, a natural approach to address the case $t > r$ is continuing with the same approach as in the case $t \in (0,r]$. However, one can see that this method turns out to be difficult to apply. Indeed, consider now $t > r$. Then, the Malliavin derivative of $X_t$ in the direction $s < r 0$. The strategy found in \cite{kohatsu2003lower} has been utilized in several contexts for proving lower bounds for the density of solution to stochastic differential equations. For instance, an application of the method to general equations driven by a fractional Brownian motion with $H > 1/2$ can be found in \cite{besalu2016gaussian}. One big difference between our adaptation to the delay case of the method presented in \cite{kohatsu2003lower} and the one used for the case $t \in (0,r]$ based on the results in \cite{NourdinViens} is that the latter method also produces l an upper bound for $p_t(x)$. The method that we will use will only give us a lower bound, but it can be complemented by an upper bound using similar arguments as in \cite{baudoin}, adapting them properly to the delay scenario. In general, proving Gaussian type lower bounds for density functions of the solutions to stochastic differential equations is generally a long and technically demanding task. One of the main objectives of this work is to illustrate how the delays are actually helpful in order to get simpler proofs. The first step to work with this method is to represent the density function $p_t(x)$ as $\E[\delta_x(X_t)]$. To do so, we rely on the following result. \begin{theorem} \label{th: representation} Under the hypothesis of Theorem \ref{th: main result}, the unique solution to \eqref{eq: main equation} is a non-degenerate random variable in the sense of Malliavin for all $t \in (0,T]$, that is, \begin{itemize} \item [(i)] $X_t \in \D^{\infty}$. \item [(ii)] $\Gamma^{-1}_{X_t} > 0$ almost surely, and $\Gamma^{-1}_{X_t}\in \bigcap_{p\geq 1}L^p(\Omega)$. \end{itemize} Moreover, the density $p_t(x)$ of $X_t$ admits the representation $p_t(x) = \E[\delta_x(X_t)]$, where $\delta_x$ denotes Dirac's delta measure at $x$. \end{theorem} \begin{proof} Items $(i)$ and $(ii)$ are proved in \cite{ferrante2006stochastic}. Hence, by the same argument as in \cite{kohatsu2003lower}, we obtain that the density $p_t(x)$ of $X_t$ can be expressed as \[ p_t(x) = \E[\delta_{x}(X_t)], \] as desired. \end{proof} In order to analyze $p_t(x) = \E[\delta_x(X_t)]$, we will construct an approximating sequence $\{F_n; 0 \leq n \leq N\}$ such that $F_N = X_t$, and we will evaluate $p_t(x)$ via evaluating the conditional densities of every $F_n$. Notice that this method considers $N$ as large as needed but fixed. The value of $N$ will be chosen appropriately later on, since several constraints on $N$ will need to be simultaneously satisfied. By construction, analyzing $p_t(x)$ is equivalent to analyzing $\E[\delta_x(F_N)]$ thanks to the representation given by Theorem \ref{th: representation}. In order to construct such an approximating sequence, we will construct a partition $\pi = \{0 = t_0 < t_1 < \cdots < t_N = t\}$ such that $|t_n - t_{n-1}| < \eps$ for a small enough $\eps > 0$ that will be chosen conveniently in the future. Then, $F_n$ will be an $\mathcal{F}_{t_n}$-measurable random variable. In order to exploit the property of the delay, we will consider $F_{n}$ of the form \[ F_n = F_{n-1} + I_n + R_n, \] where $I_n$ is a stochastic integral of a predictable process and $R_n$ is the remainder term. In the case of a general stochastic differential equation, the choice of $I_n$ is not direct. In fact, choosing $I_n$ as the integral of a predictable process makes the term $R_n$ difficult to study. Thanks to the delays, our choice of $I_n$ is straightforward and natural. We therefore define \begin{align*} F_{n-1} =& \ \eta_0 + \int_0^{t_{n-1}} \sigma(X_{s-r})dB^H_{s} + \int_0^{t_{n-1}}b(X_s)ds, \\ I_n =& \int_{t_{n-1}}^{t_n} \sigma(X_{s-r})dB^H_s, \end{align*} and \[ R_n = \int_{t_{n-1}}^{t_n} b(X_s) ds. \] This allows us to write \[ p_t(x) = \E[\delta_x(F_N)] = \E[\delta_x(F_{N-1} + I_{N} + R_N)], \] which, by the properties of the conditional expectation, can also be written as \[ p_t(x) = \E[ E_{t_{N-1}}(\delta_x(F_{N-1} + I_{N} + R_N))]. \] Therefore, in order to get information about $p_t(x)$ we first need to get as much information as we can about $E_{t_{N-1}}[\delta_x(F_{N-1} + I_{N} + R_N)]$. Using a Taylor expansion, this term can be written as \[ \begin{split} E_{t_{N-1}}[\delta_x(F_{N-1} + I_{N} + R_N)] = & E_{t_{N-1}}[\delta_x(F_{N-1} + I_N)] \\ &+E_{t_{N-1}}\left[\int_0^1 \delta_x'(F_{N-1} + I_N + \rho R_N) R_N d\rho \right], \end{split} \] where the derivative in the third term has to be understood as the second order derivative in the distributional sense of $\1_{\{F_{N-1} + I_N + \rho R_N > x\}}$. Inspired by this decomposition, we write for all $1 \leq n \leq N$, \begin{equation} \label{eq: Js} J_{1,n} = E_{t_{n-1}}[\delta_x(F_{n-1} + I_n)], \quad J_{2,n} = E_{t_{n-1}}\left[\int_0^1 \delta_x'(F_{n-1} + I_n + \rho R_n) R_n d\rho \right]. \end{equation} So, a first (and natural) lower bound is \[ E_{t_{n-1}}[\delta_x(F_{n-1}+I_n + R_n)] = J_{1,n} + J_{2,n} \geq J_{1.n} - |J_{2,n}|. \] The following sections are devoted to the derivation of a lower bound for $J_{1,n}$ and an upper bound for $|J_{2,n}|$. From now on, we will use $C, c$ as constants that may depend on $\sigma$,$b$, $T$, $H$, $r$ or the exponents appearing in the integration by parts estimate \eqref{eq: IBP estimate}. However, these constants are independent of $N$ and $t$, so we will not track their precise dependencies as it not relevant to the proof of the Gaussian lower bound. It is important to note, nonetheless, that the dependency on $r$ prevents us from obtaining a similar lower bound for the density of general stochastic differential equations driven by a fBm with $H > 1/2$. \subsubsection{Lower bound for $J_{1,n}$} The first constraint we will impose on the constant $N$ to make the lower bound possible is that $N$ is large enough so that \[ \alpha_H\frac{t^{2H}}{N} < r^{2H}. \] The key to give a lower bound for $J_{1,n}$ is the following result: \begin{proposition} \label{p: lower bound J_1} Assume there exist two real constants $0 < \lambda < \Lambda$ such that $\lambda \leq \inf_{x \in \R} \sigma(x) \leq \sup_{x \in \R} \sigma(x) \leq \Lambda$. Let $J_{1,n}$ be defined as in \eqref{eq: Js}. Then, if the partition $\pi$ is such that $|t_n - t_{n-1}|^{2H} = \sigma^2_N := \alpha_H\frac{t^{2H}}{N} $, then \[ J_{1,n}\geq \frac{1}{\sqrt{2\pi \Lambda^2 \sigma^2_N}} \exp \left( - \frac{(x-F_{n-1})^2}{2 \lambda^2 \sigma^2_N} \right). \] \end{proposition} \begin{proof} Notice that the term $J_{1,n}$ is the conditional density of $F_{n-1} + I_n$ with respect to $\mathcal{F}_{t_n}$. Hence, in order to get a lower bound for $J_{1,n}$, we need to find and bound this conditional density. On the one hand, $F_{n-1}$ is $\mathcal{F}_{t_{n-1}}$-measurable, so it behaves like a constant when we condition it with respect to $\mathcal{F}_{t_{n-1}}$. From the fact that $t_{n} - t_{n-1} < r$, we have that $\sigma(X_{s-r})$ is $\mathcal{F}_{t_{n-1}}$ measurable for all $s \in [t_{n-1}, t_n]$ and therefore, the law of $F_{n-1} + I_n$, when conditioned to $\mathcal{F}_{t_{n-1}}$, follows a Gaussian distribution with mean $F_{n-1}$ and variance $||\sigma(X_{\cdot - r})||^2_{\mathcal{H}[t_{n-1}, t_n]}$. Moreover, from the ellipticity condition on $\sigma$ and the definition of $\sigma^2_N$ we readily find that \[ \lambda^2 \sigma^2_N \leq ||\sigma(X_{\cdot - r})||^2_{\mathcal{H}[t_{n-1}, t_n]} \leq \Lambda^2 \sigma^2_N. \] This allows us to conclude that \[ J_{1,n}\geq \frac{1}{\sqrt{2\pi \Lambda^2 \sigma^2_N}} \exp \left( - \frac{(x-F_{n-1})^2}{2 \lambda^2 \sigma^2_N} \right), \] as claimed. \end{proof} \begin{remark} It remains to check the fact that such a partition defined as in Proposition \ref{p: lower bound J_1} exists. However, proving that there exists a unique partition $\pi = \{0 = t_0 < t_1 < \cdots < t_N = t\}$ of the interval $[0,t]$ with the property that $|t_n - t_{n-1}|^{2H} = \sigma^2_N := \alpha_H \frac{t^{2H}}{N}$ is direct. \end{remark} Observe that the proof of Proposition \ref{p: lower bound J_1} holds true under the assumption of the existence of such a partition $\pi$, which is equivalent to assuming that $N$ can be chosen arbitrarily large. In the upcoming sections we will see that the constraints on $N$ allow us to choose this parameter as large as desired. \subsubsection{Upper bound for $|J_{2,n}|$} The result concerning the bound for $|J_{2,n}|$ as follows. \begin{proposition} \label{p: bound J_2n a priori} Let $J_{2,n}$ be defined as in \eqref{eq: Js}. If the partition $\pi$ is the same as in Proposition \ref{p: lower bound J_1}, there exist two constants $C, \gamma > 0$ such that \begin{equation} \label{eq: bound J_2n a priori} |J_{2,n}| \leq C\frac{N^{-\gamma}}{\sigma_N}. \end{equation} \end{proposition} The proof of this proposition is not immediate. In order to conclude \eqref{eq: bound J_2n a priori}, we need to derive a first estimate using the integration by parts formula and then we will bound each of the terms involved using some technical lemmas. First, recall that \[ J_{2,n} = E_{t_{n-1}}\left[\int_0^1 \delta_x'(F_{n-1} + I_n + \rho R_n) R_n d\rho \right]. \] We introduce a new random variable $U^{\rho}_n$ defined by \begin{equation*}% \label{eq: definition Un} \sigma_N U^{\rho}_n = I_n + \rho R_n. \end{equation*} With this random variable, $J_{2,n}$ can now be written as \[ J_{2,n} = E_{t_{n-1}}\left[ \int_0^1\delta_x' (F_{n-1} + \sigma_N U^{\rho}_n)R_n d\rho \right]. \] Using Fubini's theorem, we can rewrite $J_{2,n}$ as \[ J_{2,n} = \int_0^1 E_{t_{n-1}} \left[\delta_x'(F_{n-1} + \sigma_N U^{\rho}_n) R_n \right]d\rho. \] Integration by parts formula \eqref{eq: IBP} for second order derivatives and Remark \ref{remark 1} now yields \[ J_{2,n} = \int_0^1 E_{t_{n-1}} \left[ \1_{\{I_n + \rho R_n > x- F_{n-1}}\} H^{t_{n-1}}_{(\1^2)}( \sigma_N U^{\rho}_n, R_n)\right] d\rho. \] Then, using Remark \ref{remark 2}, this last expression can be expressed as \[ J_{2,n} = \sigma_N^{-2} \int_0^1 E_{t_{n-1}} \left[ \1_{\{I_n + \rho R_n > x- F_{n-1}}\} H^{t_{n-1}}_{(\1^2)}( U^{\rho}_n, R_n)\right] d\rho. \] Now, from the fact that $\1_{\{I_n + \rho R_n > x-F_{n-1}\}} \leq 1$, Hölder's inequality and the estimate \eqref{eq: IBP estimate}, $J_{2,n}$ is bounded in the following way: \[ |J_{2,n}| \leq c \sigma_N^{-2} A_1 \int_0^1 A_2(\rho)A_3(\rho) d\rho, \] where \[ A_1 = ||R_n||_{k_1, p_1, t_{n-1}}, \] \[ A_2(\rho) = || \Gamma^{-1}_{U^{\rho}_n, t_{n-1}}||_{p_3, t_{n-1}}^{k_3}, \] and \[ A_3(\rho) = ||U^{\rho}_n ||^{k_4}_{k_2, p_2, t_{n-1}} \] for some constants $k_1, k_2, k_3, k_4, p_1, p_2, p_3 > 0$. The exact values of the constants are not extremely important, but they can be determined according to the integration by parts formula estimate \eqref{eq: IBP estimate}. Indeed, for the second order derivative case, we know that $k_1 = 2$, $k_2 = 4$, $k_3 = 2$, $k_4 = 6$ and $p_1 = q_3$, $p_2 = 4q_2$, $p_3 = 2q_1$, where $q_1, q_2, q_3$ satisfy $\frac{1}{q_1} + \frac{1}{q_2} + \frac{1}{q_3} = 1$. The bound on the quantities $A_j$ relies on the following lemma: \begin{lemma} \label{l: X uniformly bounded} Let $\pi$ be the partition defined as in Proposition \ref{p: lower bound J_1}. If $s_1, \dots, s_j, \tau \in [t_{n-1},t_n]$ with $s_1 \vee \dots \vee s_j <\tau$, then $D^{(j)}_{s_1, \dots, s_j} X_{\tau}$ is uniformly bounded in $\omega \in \Omega$. \end{lemma} \begin{proof} The result for the first order derivative is straightforward. Indeed, differentiating equation \eqref{eq: main equation} in the direction $s \in [t_{n-1}, t_n]$ we obtain \[ D_sX_{\tau} = \sigma(X_{s-r}) + \int_{0}^{\tau} D_sX_{u-r}\sigma'(X_{u-r})dB^H_u + \int_0^{\tau} b'(X_u)D_sX_u du. \] Now, since $D_sX_{u-r} = 0$ for all $u \in [0,\tau]$ due to the choice of the partition, the last expression can be written as \[ D_sX_{\tau} = \sigma(X_{s-r}) + \int_0^{\tau} b'(X_u) D_sX_u du. \] This equation has an explicit solution \[ D_sX_{\tau} = \sigma(X_{s-r}) \exp\left(\int_s^{\tau} b'(X_u) du \right). \] Finally, from the hypothesis on $\sigma$ and $b'$, we have that there exists a constant $C > 0$ such that \[ |D_sX_{\tau}| \leq C. \] In order to prove the result for higher order derivatives, we will use an induction argument. Assume that all derivatives up to order $j-1$ are bounded. The derivative of order $j$ satisfies the following equation: \begin{equation} \label{eq: Higher Order Derivatives} \begin{split} D^{(j)}_{s_1, \dots , s_j} X_{\tau} = &\sum_{q = 1}^j D_{s_1} \cdots \check{D_{s_q}} \cdots D_{s_j}\sigma(X_{s-r}) \\ & + \int_{0}^{\tau} D_{s_1} \dots D_{s_j} \sigma(X_{u-r}) dB^H_u \\ & + \int_{0}^{\tau} D_{s_1} \dots D_{s_j} b(X_{u}) du. \end{split} \end{equation} Now, since $s_1\vee \cdots \vee s_j - r < t_{n-1}$ and $\tau \in [t_{n-1}, t_n]$, then equation \eqref{eq: Higher Order Derivatives} can be written as \begin{equation} \label{eq: induction} D^{(j)}_{s_1,\dots, s_j} X_{\tau} = \sum_{q = 1}^j D_{s_1} \cdots \check{D_{s_q}} \cdots D_{s_j}\sigma(X_{s_{q}-r}) + \int_{s_1 \vee \cdots \vee s_j}^{\tau} D_{s_1} \dots D_{s_j} b(X_{u}) du, \end{equation} where $\check{D_{s_q}}$ means that the term $D_{s_q}$ is omitted (for instance, $D_{s_1}\check{D_{s_2}}D_{s_3} = D_{s_1}D_{s_3}$). On the one hand, each term of the form \[ D_{s_1} \cdots \check{D_{s_q}} \cdots D_{s_j}\sigma(X_{s_{q}-r}) \] involves derivatives of orders $k = 1, \dots, j-1$ of $X_t$ and derivatives of orders $l = 1, \dots j$ of $\sigma$. Since $\sigma$ has bounded derivatives of all orders and, by induction hypothesis, all derivatives up to order $j-1$ of $X_t$ are bounded, we get that there exists $C_1 > 0$ such that \[ \sum_{q = 1}^{j} |D_{s_1} \cdots \check{D_{s_q}} \cdots D_{s_j}\sigma(X_{s_{q}-r})| \leq C_1. \] Concerning the term $D_{s_1} \dots D_{s_j} b(X_{u})$, we use the product rule and the chain rule for Malliavin derivatives $j$ times and we find that \[ D_{s_1} \dots D_{s_j} b(X_{u}) = \sum_{k=1}^j b^{(k)}(X_u) \sum_{\Pi \in \mathcal{P}(j,k)} \prod_{l=1}^kD^{(|\Pi_l|)}_{s_{\Pi_l}}X_u, \] where $\mathcal{P}(j,k)$ denotes the set of all partitions of $\{s_1, \dots, s_j\}$ into $k$ subsets and $\Pi_l$ denotes the $l$-th subset of this partition. $|\Pi_l|$ denotes its length and $s_{\Pi_l} = (s_{i_1}, \dots, s_{i_{|\Pi_l}})$ if $\Pi_l = \{s_{i_1}, \dots, s_{i_{|\Pi_l|}}\}$. Even though it is a sum with a lot of terms, the only one that involves a $j$-th order derivative of $X_u$ is \[ b'(X_u) D^{(j)}_{s_1, \dots, s_j}X_u. \] Taking that all derivatives of $b$ are bounded and all derivatives of $X_u$ up to order $j-1$ are bounded, we find that there exist $C_2, C_3 > 0$ such that \begin{equation} \label{eq: higher order derivatives of b} D_{s_1} \dots D_{s_j} b(X_{u}) \leq C_2 + C_3 D^{(j)}_{s_1, \dots, s_j}X_u. \end{equation} Inserting this into \eqref{eq: induction} leads us to \[ D^{(j)}_{s_1, \dots, s_j}X_{\tau} \leq C_1 + \int_0^{\tau} \left( C_2 + C_3 D^{(j)}_{s_1, \dots, s_j}X_u \right) du. \] Using Gronwall's lemma, we conclude that $|D^{(j)}_{s_1, \dots, s_j}X_{\tau}|$ is uniformly bounded. \end{proof} We can now proceed to derive the estimates. Concerning the following 3 lemmas, in the proofs we will use $C$ as a universal constant that may switch from line to line. If several different constants are involved, we will name them $C_1, C_2, \dots$. To streamline the reading of the paper, we omit the dependency of the constant on the parameters. \begin{lemma}[Bound for $A_1$] \label{l: bound for A_1} Recall that $A_1 = ||R_n||_{k_1, p_1, t_{n-1}}$. Under the same hypothesis as in Lemma \ref{l: X uniformly bounded}, there exists a constant $C > 0$ depending on the maximum of the derivatives of the function $b$ (which is finite because $b \in \mathcal{C}^{\infty}_b(\R)$), the Hurst parameter $H$ and $T$ such that \[ A_1 \leq C N^{-1/(2H)}. \] In particular, there exists $\gamma > 0$ depending on the Hurst parameter $H$ and $C > 0$ depending on $b$ (and its derivatives), $T$, and $H$ such that \[ \sigma_N^{-2}A_1 \leq C \frac{N^{-\gamma}}{\sigma_N}. \] \end{lemma} \begin{proof} By definition, \[ A_1^{p_1} = ||R_n||^{p_1}_{k_1, p_1, t_{n-1}} = E_{t_{n-1}}(|R_n|^{p_1}) + \sum_{j=1}^{k_1} E_{t_{n-1}}\left(||D^{(j)}R_n||^{p_1}_{\mathcal{H}^{\otimes j}[t_{n-1},t_n]}\right). \] The strategy of the proof is to verify that $|R_n|$ is of the order of $N^{-1/(2H)}$, and the Malliavin derivatives are $o(N^{-1})$, so that the asymptotic behaviour of $A_1$ is dominated by the first term of $A_1$, and the terms involving the derivatives of $R_n$ are negligible. For the first term, notice that $|R_n| \leq \int_{t_{n-1}}^{t_n}| b(X_s)| ds \leq ||b||_{L^{\infty}}|t_n-t_{n-1}|$. Since the partition is chosen so that \[|t_n - t_{n-1}|^{2H} =\alpha_H \frac{t^{2H}}{N} \leq \alpha_H\frac{T^{2H}}{N},\] we can conclude that there exists a constant $C > 0$ depending on $b$, $T$ and $H$ such that $|R_n| \leq C N^{-1/(2H)}$, so we can find $C > 0$ such that $E_ {t_{n-1}}(|R_n|^{p_1}) \leq C N^{-p_1/(2H)}$, which is a consistent estimate with what we want to prove. The rest of the proof is done by bounding all the derivatives of order $j$ from $1$ up to $k_1$. We will show the bound for the first derivative since its computation will be useful in future lemmas, and the conclusion for the higher order derivatives can be easily derived by using estimate \eqref{eq: higher order derivatives of b} and the fact that the Malliavin derivatives of $X_t$ are bounded thanks to Lemma \ref{l: X uniformly bounded}. For the first derivative, \[ ||DR_n||^{p_1}_{\mathcal{H}[t_{n-1}, t_n]} = \left(\int_{t_{n-1}}^{t_n} \int_{t_{n-1}}^{t_n} D_uR_nD_vR_n|u-v|^{2H-2}du dv\right )^{p_1/2}. \] Now, using the fact that \[ |D_sR_n| \leq \int_{t_{n-1}}^{t_n} |b'(X_u) D_sX_u| du \leq C|t_n - t_{n-1}| \leq CN^{-1/(2H)}, \] we find that \begin{equation} \label{eq: bound norm of DR_n} ||DR_n||^{p_1}_{\mathcal{H}[t_{n-1}, t_n]} \leq C \left( N^{-1/H} |t_n - t_{n-1}|^{2H}\right)^{p_1/2} \leq CN^{-\frac{p_1(1+\frac{1}{H})}{2}}. \end{equation} Hence, \[ E_{t_{n-1}}\left( ||DR_n||^{p_1}_{\mathcal{H}[t_{n-1}, t_n]}\right) \leq CN^{-\frac{p_1(1+\frac{1}{H})}{2}}, \] as desired. \end{proof} \begin{lemma}[Bound for $A_2(\rho)$] \label{l: bound for A_2} Recall that $A_2(\rho) = ||\Gamma^{-1}_{U^{\rho}_n, t_{n-1}}||^{k_3}_{p_3, t_{n-1}}$. Under the same hypothesis as in Lemma \ref{l: X uniformly bounded}, the quantity $A_2(\rho)$ is uniformly bounded. \end{lemma} \begin{proof} Recall that \[ \sigma_N U^{\rho}_n = I_n + \rho R_n, \] so \[ \sigma^2_N ||DU^{\rho}_n||^2_{\mathcal{H}[t_{n-1}, t_n]} = ||DI_n + \rho DR_n||^2_{\mathcal{H}[t_{n-1}, t_n]}. \] Now, using Lemma \ref{l: Ap 1} in the appendix and the fact that $-\rho \geq -1$, we have \[ \sigma^2_N ||DU^{\rho}_n||^2_{\mathcal{H}[t_{n-1}, t_n]} \geq \frac{||DI_n||^2_{\mathcal{H}[t_{n-1}, t_n]}}{2} - ||DR_n||^2_{\mathcal{H}[t_{n-1}, t_n]}. \] Since \[ I_n = \int_{t_{n-1}}^{t_n} \sigma(X_{s-r})dB^H_s, \] it is clear that for $s \in [t_{n-1}, t_n]$ we have \[ D_sI_n = \sigma(X_{s-r}). \] From the fact that $\sigma(X_{s-r}) > \lambda$ and $\int_{t_{n-1}}^{t_n}\int_{t_{n-1}}^{t_n} |u-v|^{2H-2} dudv = C \sigma^2_N$ (see Lemma \ref{l: Ap 2} in the appendix) we derive the following estimate: \[ ||DI_n||^2_{\mathcal{H}[t_{n-1},t_n]} = \int_{t_{n-1}}^{t_n} \int_{t_{n-1}}^{t_n} \sigma(X_{u-r})\sigma(X_{v-r}) \cdot |u-v|^{2H-2} du dv \geq C \sigma^2_N. \] Moreover, by inequality \eqref{eq: bound norm of DR_n} in the proof of the previous lemma for $p_1 = 2$, we have that $$||DR_n||^2_{\mathcal{H}[t_{n-1}, t_n]} \leq CN^{-1-\frac{1}{H}}.$$ Since $\sigma^2_N = O(N^{-1})$, for $N$ big enough we deduce that there exist $C_1, C_2 > 0$ not depending on $t$ such that \[ \sigma^2_N ||DU^{\rho}_n||^2_{\mathcal{H}[t_{n-1},t_n]} \geq C_1\sigma^2_N - C_2 N^{-1-\frac{1}{H}}. \] Now, since \begin{align*} C_1\sigma^2_N - C_2 N^{-1-\frac{1}{H}} = &\sigma^2_N\left(C_1 - \frac{C_2}{\alpha_H t^{2H}N^{\frac{1}{H}}} \right) \\ \geq &\sigma^2_N\left(C_1 - \frac{C_2}{\alpha_H r^{2H}N^{\frac{1}{H}}} \right), \end{align*} because $t > r$, we can choose $N$ large enough so that there exists $C_3 > 0$ with \[ C_1\sigma^2_N - C_2 N^{-1-\frac{1}{H}} \geq \sigma^2_N\left(C_1 - \frac{C_2}{\alpha_H r^{2H}N^{\frac{1}{H}}} \right) \geq C_3 \sigma_N^2, \] so \[ ||DU^{\rho}_n||_{\mathcal{H}[t_{n-1}, t_n]}^2 \geq C_3 \frac{\sigma^2_N}{\sigma^2_N} = C_3 > 0. \] This implies that \[ \Gamma^{-1}_{U^{\rho}_n, t_{n-1}} \leq \frac{1}{C_3}. \] Taking conditional norms, we deduce that there exists $C > 0$ satisfying $||\Gamma^{-1}_{U^{\rho}_n, t_{n-1}}||^{k_3}_{p_3, t_{n-1}} \leq C$, which proves the uniform bound of $A_2(\rho)$. \end{proof} \begin{lemma}[Bound for $A_3(\rho)$] \label{l: bound for A_3} Recall that $A_3(\rho) = ||U^{\rho}_n||^{k_4}_{k_2, p_2, t_{n-1}}$. Under the same hypothesis as in Lemma \ref{l: X uniformly bounded}, the quantity $A_3(\rho)$ is uniformly bounded in $\rho$ and in $\omega \in \Omega$. \end{lemma} \begin{proof} Applying the relation $\sigma_N U^{\rho}_n = I_n + \rho R_n$, we readily see by applying norms that \[ \sigma^{k_4}_N||U^{\rho}_n||^{k_4}_{k_2,p_2, t_{n-1}} = ||I_n + \rho R_n||^{k_4}_{k_2,p_2, t_{n-1}}. \] By virtue of Minkowski's inequality, the fact that $(|a|+|b|)^{k_4} \leq C(|a|^{k_4} + |b|^{k_4})$ and $\rho \leq 1$, we have \begin{equation} \label{eq: bound Un} ||U^{\rho}_n||^{k_4}_{k_2,p_2, t_{n-1}} \leq C \sigma_N^{-k_4} \left( ||I_n||^{k_4}_{k_2,p_2, t_{n-1}} + ||R_n||^{k_4}_{k_2,p_2, t_{n-1}} \right), \end{equation} so it is enough to bound the quantities $||I_n||^{k_4}_{k_2,p_2, t_{n-1}}$ and $||R_n||^{k_4}_{k_2,p_2, t_{n-1}}$ separately. Lemma \ref{l: bound for A_1} gives us the bound \begin{equation*}% \label{eq: estimate 1} ||R_n||^{k_4}_{k_2, p_2, t_{n-1}} \leq C N^{-k_4/2H}. \end{equation*} In particular, this bound implies that \begin{equation} \label{eq: Bound Rn} \sigma^{-k_4}_N||R_n||^{k_4}_{k_2,p_2,t_{n-1}} \leq C N^{\frac{-k_4}{2}\left( \frac{1}{H}-1\right)}. \end{equation} For the $I_n$ term, we know that $D_sI_n = \sigma(X_{s-r})$, and therefore, the $j$-th order derivatives of $I_n$ in directions belonging to $[t_n, t_{n-1}]^j$ vanish if $j \geq 2$ due to the fact that the partition is chosen so that $s-r < t_{n-1}$ for all $s \in [t_{n-1}, t_n]$. This implies that \[ ||I_n||^{k_4}_{k_2, p_2, t_{n-1}} = \left(E_{t_{n-1}}(|I_n|^{p_2}) + E_{t_{n-1}}(||DI_n||_{\mathcal{H}[t_{n-1},t_n]}^{p_2}) \right)^{k_4/p_2}. \] On the one hand, since $I_n = \int_{t_{n-1}}^{t_n} \sigma(X_{s-r})dB^H_s$ and $s-r < t_{n-1}$, the law of $I_n$ conditioned to $\mathcal{F}_{t_{n-1}}$ is conditionally Gaussian with zero mean and variance $||\sigma(X_{\cdot -r})||_{\mathcal{H}[t_{n-1},t_n]}^2$. Hence, it is then well-known that by the properties of the moments of Gaussian random variables we have \[ E_{t_{n-1}}(|I_n|^{p_2}) = C||\sigma(X_{\cdot -r})||^{p_2}_{\mathcal{H}[t_{n-1},t_n]}. \] Finally, taking that $\sigma(X_{s-r}) \leq \Lambda$ it is easy to find that \begin{equation} \label{eq: bound variance} ||\sigma(X_{\cdot-r})||^{p_2}_{\mathcal{H}[t_{n-1},t_n]} \leq \Lambda^{p_2} \sigma_N^{p_2}, \end{equation} from which we conclude that \[ E_{t_{n-1}}(|I_n|^{p_2}) \leq C\sigma^{p_2}_N \] for a certain constant $C > 0$. Now, since $||DI_n||^{p_2}_{\mathcal{H}[t_{n-1},t_n]} = ||\sigma(X_{\cdot-r})||^{p_2}_{\mathcal{H}[t_{n-1},t_n]}$, we resort to inequality \eqref{eq: bound variance} to conclude that there exists $C > 0$ such that \[ E_{t_{n-1}}(||DI_n||^{p_2}_{\mathcal{H}[t_{n-1},t_n]}) \leq C\sigma^{p_2}_N. \] All in all, the conclusion of the previous estimates is that \begin{equation*} ||I_n||^{k_4}_{k_2, p_2, t_{n-1}} \leq C \sigma_N^{k_4}, \end{equation*} and in particular, \begin{equation}\label{eq: Bound In} \sigma_N^{-k_4}||I_n||^{k_4}_{k_2, p_2, t_{n-1}} \leq C. \end{equation} To end the proof of this lemma, we insert in estimates \eqref{eq: Bound Rn} and \eqref{eq: Bound In} into \eqref{eq: bound Un} and get \[ ||U^{\rho}_n||^{k_4}_{k_2, p_2, t_{n-1}} \leq C \] for a constant $C > 0$ depending on $p_2, \Lambda, k_4$ and $H$. \end{proof} Let us now illustrate how these 3 lemmas prove Proposition \ref{p: bound J_2n a priori}. Recall that, using the integration by parts formula, we derived the estimate \[ |J_{2,n}| \leq \sigma^{-2}_N A_1 \int_0^1 A_2(\rho) A_3(\rho) d \rho. \] First, the conclusion of Lemma \ref{l: bound for A_1} is that \[ \sigma_N^{-2}A_1 \leq \frac{C N^{-\gamma}}{\sigma_N} \] for real constants $C, \gamma > 0$. Moreover, the conclusion of Lemmas \ref{l: bound for A_2} and \ref{l: bound for A_3} is that there exists a constant $C > 0$ such that $A_2(\rho) \leq C$ and $A_3(\rho) \leq C$ uniformly in $\rho \in [0,1]$ and in $\omega \in \Omega$, respectively. This information leads us to the following estimate: \begin{equation} \label{eq: Final bound J2n} |J_{2,n}| \leq \sigma^{-2}_N A_1 \int_0^1 A_2(\rho) A_3(\rho) d \rho \leq \frac{C N^{-\gamma}}{\sigma_N}, \end{equation} which proves Proposition \ref{p: bound J_2n a priori}. With these bounds on $J_{1,n}$ and $J_{2,n}$, we have all the tools needed in order to prove the lower bound for $t > r$. \subsection{Proof of the lower bound} \label{ss: Final proof} In Section 3.1, we have already proved the result when $t \in (0,r]$, and in Subsections 3.2.1 and 3.2.2 we have proved all the auxiliary results that allow us to prove the lower bound. This section is devoted to putting all the information together and concluding that the density of the solution $X_t$ to \eqref{eq: main equation} is strictly positive. \begin{proof}[Proof of Theorem \ref{th: main result}] Recall that, thanks to Proposition \ref{p: lower bound J_1} and the bound we have concluded in equation \eqref{eq: Final bound J2n}, we can bound $E_{t_{n-1}}[\delta_x(F_n)]$ as \begin{equation} \label{eq: lower bound conditioned} E_{t_{n-1}}[\delta_x(F_n)] \geq \frac{1}{\sqrt{2\pi \Lambda^2 \sigma^2_N}} \exp \left( - \frac{(x-F_{n-1})^2}{2 \lambda^2 \sigma^2_N} \right) - \frac{C N^{-\gamma}}{\sigma_N}. \end{equation} We define the intervals $I_n = I(y_n, c_1 \sigma_N) := \{z \in \R; |z-y_n|< c_1\sigma_N\}$ for some constant $c_1 > 0$ to be determined and $y_n = \eta_0 + \frac{n}{N}(x-\eta_0)$, where $y_0 = \eta_0 = F_0$ is the initial condition of the SDDE. We also define $\{x_n; n = 0, \dots, N\}$, where $x_0 = \eta_0 = F_0$, $x_j \in I_j$ for $j = 1, \dots, N-1$ and $x_N = x$. Moreover, we will assume that $|y_n - y_{n-1}| \leq c_1 \sigma_N$ (later in the proof we will see that the constants $c_1$ and $N$ can be chosen so that it is satisfied). We first rewrite the density of $F_N$ as \[ \E[\delta_x(F_N)]= \int_{\R} \E[\delta_x(F_N) \delta_{x_{N-1}}(F_{N-1})] dx_{N-1}. \] Then, due to the positivity of $\E[\delta_x(F_N) \delta_{x_{N-1}}(F_{N-1})]$ (see \cite[Theorem 1]{kohatsu2003lower}) we have the bound \[ \int_{\R} \E[\delta_x(F_N) \delta_{x_{N-1}}(F_{N-1})] dx_{N-1} \geq \int_{I_{N-1}}\E[\delta_x(F_N) \delta_{x_{N-1}}(F_{N-1})] dx_{N-1}. \] Using the properties of the conditional expectation and estimate \eqref{eq: lower bound conditioned} we find that \[ \begin{split} &\E[\delta_x(F_N)] \\ \geq & \frac{c}{\sigma_N} \int_{I_{N-1}}\E\left[\left(\exp \left( - \frac{(x-F_{N-1})^2}{2 \lambda^2 \sigma^2_N} \right) - C N^{-\gamma}\right)\delta_{x_{N-1}}(F_{N-1})\right]dx_{N-1}. \\ \end{split} \] We want this integral to give a strictly positive contribution (otherwise, the resulting bound is $p_t(x) \geq 0$, which gives us no information). In order to achieve this objective, since $x_{N-1} \in I_{N-1}$, we can observe by Fubini's theorem that the expectation is taken in the set $\{\omega \in \Omega; F_{N-1}(\omega) \in I_{N-1}\}$ since the term $\delta_{x_{N-1}}(F_{N-1})$ vanishes if $F_{N-1} \notin I_{N-1}$. Moreover, if $F_{N-1} \in I_{N-1}$ we have \begin{equation} \label{eq: x-F(N-1)} |x-F_{N-1}|\leq |x - y_N| + |y_N - y_{N-1}| + |y_{N-1} + F_{N-1}| \leq 3c_1 \sigma_N \leq 4c_1 \sigma_N, \end{equation} where we choose bound $4c_1 \sigma_N$ instead of $3c_1\sigma_N$ to make the computations nicer. In order to obtain a Gaussian type bound, we will choose the constants $c_1$ and $N$ appropriately. First, we choose $c_1$ small enough so that \[ \exp\left( - \frac{8c_1^2}{\lambda^2}\right) \geq \frac{1}{2}. \] Regarding the choice of $N$, we consider a constant $c_2 \geq \frac{1}{4c_1^2}$ and choose $N$ of the form \[ N = \frac{c_2|x-\eta_0|^2}{t^{2H}}. \] Notice that under this definition of $c_2$ and $N$, condition \eqref{eq: x-F(N-1)} is satisfied. The idea is that $c_2$ is a constant that controls how big we want $N$ to be, so the larger $c_2$, the larger $N$. Since $c_2 \geq \frac{1}{4c_1^2}$ can be chosen arbitrarily large, we will demand $c_2$ to be large enough so that $N = \frac{c_2 |x-\eta_0|^2}{t^{2H}}$ satisfies $C N^{-\gamma} \leq 1/4$. Under this first restriction we have that \begin{equation*}% \label{eq: lower bound exponential} \exp \left(- \frac{(x-F_{N-1})^2}{2\lambda^2 \sigma^2_N} \right) \geq \exp \left(- \frac{16c_1^2 \sigma_N^2}{2\lambda^2 \sigma^2_N} \right) =\exp \left(- \frac{8c_1^2}{\lambda^2} \right) \geq \frac{1}{2}. \end{equation*} On the other hand, taking that $CN^{-\gamma} \leq 1/4$, we find that the integrand is stricly positive. We can therefore proceed to bound the density as follows: \begin{align*} \E[\delta_x(F_N)] \geq &\frac{c}{\sigma_N} \int_{I_{N-1}}\E\left[\left(\exp \left( - \frac{(x-F_{N-1})^2}{2 \lambda^2 \sigma^2_N} \right) - C N^{-\gamma}\right)\delta_{x_{N-1}}(F_{N-1})\right]dx_{N-1}\\ \geq &\frac{c}{\sigma_N} \int_{I_{N-1}} \E \left[ \left( \frac{1}{2} - \frac{1}{4}\right)\delta_{x_{N-1}}(F_{N-1}) \right] dx_{N-1}, \end{align*} so we are led to \[ \E[\delta_x(F_N)] \geq \frac{c}{4\sigma_N} \int_{I_{N-1}} \E[\delta_{x_{N-1}}(F_{N-1})] dx_{N-1}. \] In particular, these constraints imply that $E_{t_{n-1}}[\delta_{x_{n}}(F_n)] \geq \frac{c}{4\sigma_N}$ for all $n = 1, \dots, N$. Indeed, by equation \eqref{eq: lower bound conditioned} we have \[ \begin{split} E_{t_{n-1}}[\delta_{x_{n}}(F_n)] \geq &\frac{c}{\sigma_N} \left[\exp \left( - \frac{(x_{n}-F_{n-1})^2}{2 \lambda^2 \sigma^2_N} \right) - C N^{-\gamma}\right]\\ \geq &\frac{c}{\sigma_N} \left[ \exp \left( - \frac{(x_{n}-F_{n-1})^2}{2 \lambda^2 \sigma^2_N} \right) - \frac{1}{4} \right]. \end{split} \] Now, using the fact that \[ |x_n - F_{n-1}| \leq |x_n - y_n| + |y_n - y_{n-1}| + |y_{n-1} + F_{n-1}| \leq 3c_1 \sigma_N \leq 4c_1 \sigma_N, \] if $|y_n - y_{n-1}| \leq c_1 \sigma_N$, we have that \[ \exp \left( - \frac{(x_{n}-F_{n-1})^2}{2 \lambda^2 \sigma^2_N} \right) \geq \exp \left( - \frac{16c_1^2 \sigma_N^2}{2 \lambda^2 \sigma^2_N} \right) \geq \frac{1}{2}. \] We can proceed iterating this same procedure one step from $N-1$ to $N-2$ and taking that $\delta_{x_{N-2}}(F_{N-2})$ is $\mathcal{F}_{t_{N-2}}$-measurable to get \[ \begin{split} \E[\delta_x(F_N)] \geq &\frac{c}{4\sigma_N} \int_{I_{N-1}}\E[\delta_{x_{N-1}}(F_{N-1})] dx_{N-1}\\ \geq &\frac{c}{4\sigma_N}\int_{I_{N-1}} \int_{I_{N-2}} \E[\delta_{x_{N-1}}(F_{N-1}) \delta_{x_{N-2}}(F_{N-2})] dx_{N-2} dx_{N-1} \\ = & \frac{c}{4\sigma_N}\int_{I_{N-1}} \int_{I_{N-2}} \E\left[ E_{t_{N-2}}[\delta_{x_{N-1}}(F_{N-1})] \delta_{x_{N-2}}(F_{N-2})\right] dx_{N-2} dx_{N-1} \\ \geq & \left(\frac{c}{4\sigma_N}\right)^2 |I_{N-1}|\int_{I_{N-2}} \E[\delta_{x_{N-2}}(F_{N-2})]dx_{N-2}, \end{split} \] as long as $y_{N-1}, y_{N-2} \in I_{N-1}\cap I_{N-2} \neq \emptyset$, which holds because we are working under the assumption that $|y_{N-1} - y_{N-2}| \leq c_1 \sigma_N$. One can check by iterating this procedure from $N-1$ all the way to $1$ and taking that all intervals $I_j$ have the same length as the interval centered at the origin with radius $c_1 \sigma_N$ that \[ \begin{split} \E[\delta_x(F_N)] \geq &\left(\frac{c}{4\sigma_N}\right)^N |I(0,c_1\sigma_N)|^{N-1} \\ = &\left( \frac{c}{4}\right)^N \left(\frac{N^{1/2}}{t^{H}} \right)^N \left(\frac{t^{H}} {N^{1/2}} \right)^{N-1} (2c_1)^{N-1} \\ = & \frac{N^{1/2}}{t^{H}} \left( \frac{c}{4}\right)^N (2c_1)^{N-1} \\ = &\frac{1}{2c_1 t^H}\exp\left( N\log\left( \frac{c \cdot c_1}{2}\right)+\frac{1}{2}\log(N)\right). \end{split} \] To the previous restriction on $c_1$, i.e. $\exp\left(\frac{-8c_1^2}{\lambda^2} \right) \geq \frac{1}{2}$, we also impose that \[ L := -\log\left( \frac{c \cdot c_1}{4}\right) > 0, \] and \[ N L^N \geq 1. \] Now, taking into account that $N = \frac{c_2 |x-\eta_0|^2}{t^{2H}}$, and assuming that $N \geq 1$ (because we want $N$ to be a natural number) we get \[ \E\left[ \delta_x(F_N)\right] \geq \frac{1}{c_1 t^H}\exp \left( -\frac{L c_2 |x-\eta_0|^2}{t^{2H}}\right), \] which is a bound consistent with the one in \eqref{eq: bound main result rT}. The only task left to do is to adjust the constant $c_2$ so that it is compatible with all of the previous constraints, i.e.: \begin{itemize} \item $|y_n - y_{n-1}| \leq c_1 \sigma_N$. \item $c_2 \geq \frac{1}{4c_1^2}$. \item $N = \frac{c_2 |x-\eta_0|^2}{t^{2H}}$. \end{itemize}On the one hand, we have to take $c_2$ so that \[ |y_{n+1}-y_n| \leq c_1 \sigma_N. \] Now, from the definition of $y_n$'s, we have \[ |y_{n+1}-y_n| = \frac{|x-\eta_0|}{N} = \frac{|x-\eta_0|}{\sqrt{N}}\frac{1}{\sqrt{N}} \] and taking that $N = \frac{c_2|x-\eta_0|^2}{t^{2H}}$ we get \[ |y_{n+1}-y_n| = \frac{|x-\eta_0|}{\sqrt{N}}\frac{1}{\sqrt{c_2}} \frac{t^H}{|x-\eta_0|} = \frac{\sigma_N}{\sqrt{c_2}}. \] Hence, it is enough to consider $c_2$ in such a way that \[ \frac{1}{\sqrt{c_2}} \leq c_1, \quad c_2 \geq \frac{1}{4c_1^2}, \] where the first constraint follows from the fact that we need $|y_{n}-y_{n-1}| \leq c_1 \sigma_N$ and the second constraint follows from relation \eqref{eq: x-F(N-1)}. Choosing therefore an arbitrarily large $c_2$ satisfying \[ \frac{1}{\sqrt{c_2}}\leq c_1 \] is enough to conclude the bound. Finally, renaming the constants, we obtain the desired bound \eqref{eq: bound main result rT}. \end{proof} As a corollary of this result, we can state the conclusion of this work. \begin{corollary} Let $X$ be the solution to \eqref{eq: main equation}. Under the same hypothesis as in Theorem \ref{th: main result}, the density function $p_t(x)$ of $X_t$ is strictly positive. \end{corollary} \begin{proof} A direct consequence of Theorem \ref{th: main result}. \end{proof} %%===============================================================%% %% ACKNOWLEDGEMENTS %% %% Acknowledgements can be added here. %% %%===============================================================%% \section*{Acknowledgements} %%===============================================================%% Author Òscar Burés is supported by the predoctoral program AGAUR-FI ajuts (2025 FI-1 00580). Author Carles Rovira is supported by the grant PID2021-123733NB-I00 from SEIDI, Ministerio de Economia y Competividad. %%===============================================================%% %% REFERENCES %% %% References should be provided in bibtex file. %% %% We suggest using MR Lookup for finding bibtex entries. %% %%===============================================================%% \bibliography{references} %%===============================================================%% %% APPENDICES %% %% Appendices can be added here. %% %%===============================================================%% \normalsize \begin{appendices} \section{Real analysis tools} \begin{lemma} \label{l: Ap 1} Let $\mathcal{H}$ be a real Hilbert space. Denote by $\langle \cdot, \cdot \rangle$ the scalar product in $\mathcal{H}$ and $||\cdot||$ its norm. Then, for any $f,g \in \mathcal{H}$, we have \[ ||f+g||^2 \geq \frac{||f||^2}{2}-||g||^2. \] \end{lemma} \begin{proof} Notice that \begin{equation} \label{a: 1} \langle \frac{1}{\sqrt{2}}f + \sqrt{2}g, \frac{1}{\sqrt{2}}f + \sqrt{2}g \rangle = \frac{1}{2}||f||^2 + 2||g||^2 + 2\langle f, g \rangle \geq 0. \end{equation} Then, taking that \[ ||f + g||^2 = ||f||^2 + ||g||^2 + 2\langle f, g\rangle, \] equation \eqref{a: 1} turns into \[ ||f+g||^2 - \frac{||f||^2}{2} + ||g||^2 \geq 0, \] which proves the inequality. \end{proof} \begin{lemma} \label{l: Ap 2} Let $0