risk parameters, which can create severe losses and subsequent liquidity and
solvency issues, both for institutions and their customers. As risk
classification is widely considered a Markov process both in theory and by
practitioners, considering transiton probabilities for loans can allow us to
forececast PD values and estimate worst-case scenario provisions for loan
exposures. We note that, in practice, estimating the parameters of the asset
prices under each regime may be difficult. However, many financial
institutions consider such models and its mathematical framework is well
established, see e.g., [20] and [17]. Another approach is to use historical
parameters from Stage 1 and Stage 2 loans to estimate the changes that occur
when a loan transitions between Stages. For this example, we will be
estimating the provisions under the regime switching model developed above. To
this end, we consider an IFRS 9 compliant transition matrix:

$P=\left(\begin{array}[]{c|ccc}\text{IFRS 9 Rating}&\text{Stage 1}&\text{Stage
2}&\text{Stage 3}\\\ \hline\cr\text{Stage 1}&p_{11}&p_{12}&p_{12}\\\
\text{Stage 2}&p_{12}&p_{22}&p_{13}\\\ \text{Stage
3}&0&0&1\end{array}\right).$

For a loan originating in Stage 1, we can now forecast credit losses by taking
into account the probability of a SICR (significant increase in credit risk)
event. Under the regime switching model, in the case of a transition to
another Stage, we will need to estimate the PD values for the asset process
governed by the new parameters. Using the straightforward notation $\Phi^{i}$
or $PD^{i}$ to emphasize the Stage (regime) under which the specific PD value
is estimated, we can then define the ” _Stage-weighted provisions_ ”, given
by:

$\displaystyle
WP_{t}:=p_{11}EAD_{t}LGD_{t}PD^{1}+p_{12}\sum_{i=t}^{T}\big{(}\Phi^{2}(x,i-1)-\Phi^{2}(x,i)\big{)}\mathbb{E}\Big{[}\frac{LGD_{i}EAD_{i}}{(1+r_{i})^{i-t}}|\mathcal{F}_{t}\Big{]}+p_{13}EAD_{t}LGD_{t},$
(144)

where the third term occurs in the case of default (i.e. transition to Stege
3), we have that $PD=100\%$. This calculation holds for the case where we
consider that the transition to Stage 2 occurs one period (e.g., year) after.
However, we can also consider the cases where the deterioration occurs at any
point $k>t$. For this calculation we require the $k$-step transition matrix of
the underlying rating process, which is known to be $P^{k}$, whose elements
will be symbolized as below:

$P^{k}=\left(\begin{array}[]{c|ccc}\text{IFRS 9 Rating}&\text{Stage
1}&\text{Stage 2}&\text{Stage 3}\\\ \hline\cr\text{Stage
1}&p^{k}_{11}&p^{k}_{12}&p^{k}_{12}\\\ \text{Stage
2}&p^{k}_{12}&p^{k}_{22}&p^{k}_{13}\\\ \text{Stage
3}&0&0&1\end{array}\right).$

with the understanding that $p^{k}_{ij}$ represents the $k$-th step transition
probability. We then have:

$\displaystyle\mathbb{E}[WP_{k}|\mathcal{F}_{t}]=p^{k}_{11}EAD_{k}LGD_{k}PD^{1}_{k}+p^{k}_{12}\sum_{i=k}^{T}\big{(}\Phi^{2}(x,i-1)$
$\displaystyle-\Phi^{2}(x,i)\big{)}\mathbb{E}\Big{[}\frac{LGD_{i}EAD_{i}}{(1+r_{i})^{i-t}}|\mathcal{F}_{t}\Big{]}$
$\displaystyle+p^{k}_{13}EAD_{t}LGD_{t}.$ (145)

At any point, with the dynamics of the underlying Markov process, we can
obtain the corresponding $WP_{t}$ values and calculate the above _Expected
Stage-weighted provisions_. This calculation takes the future evolution of the
loan, as well as the regime into consideration to provide an estimation that
incorporates all scenarios. For illustrative purposes, we consider the example
below.

######  Example 7.3.

Consider an asset process governed by the regime-switching model below:

$\displaystyle
dG_{t}=\begin{cases}k_{1}(\theta_{1}-G_{t})dt+\sigma_{1}dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\
G_{0}=x,\text{ if $R_{t}=$Stage 1, }\\\
k_{2}(\theta_{2}-G_{t})dt+\sigma_{2}dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\
G_{0}=x,\text{ if $R_{t}=$Stage 2, }\\\
k_{3}(\theta_{3}-G_{t})dt+\sigma_{3}dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\
G_{0}=x,\text{ if $R_{t}=$Stage 3, }\\\ \end{cases}$ (146)

with regime-specific parameters given in Fig. 2 and $\mathcal{D}=[-6.0,6.0]$,
with $S=4.0$ (we can also consider a different limit value $S$ for each
regime. However, in this example the Monte Carlo estimates indicate that the
same value suffices). We consider normally distributed jumps with $Z\sim
N(0.0,0.5)$, rate $\lambda=1.0$. We have set $N=1001,T=1001$ for the space and
time grids, respectively. Furthermore, the generator matrix of the underlying
Markov process given by:

$Q=\left(\begin{array}[]{ccc}-0.5&0.3&0.2\\\ 0.3&-0.6&0.3\\\ 0.0&0.0&0.0\\\
\end{array}\right).$

The graphs in Fig. 2 display the estimated survival probability in each regime
(Stage), resulting from solving scheme $(\ref{matrix-form})$.

Figure 2: Survival probability for (146) with
$(k_{1},k_{2},k_{3})=(0.3,0.2,0.0),(\theta_{1},\theta_{2},\theta_{3})=(0.8,0.5,0.0),(\sigma_{1},\sigma_{2},\sigma_{3})=(0.3,0.5,0.0)$,
with $R_{0}=$ Stage 1 (top left), $R_{0}=$ Stage 2 (top right). The average
survival probability across all regimes is also shown (bottom), which also
account for the survival probability when $R_{0}=$ Stage 3, for which we have
$\Phi(x,u)\equiv 0$.

We consider an initial position of $x=0.3$ and maturity $T=10$. The transition
matrix of the underlying Markov process is obtained by calculating
$P=\exp(Q)$:

$\displaystyle P=\left(\begin{array}[]{c|ccc}\text{IFRS 9 Rating}&\text{Stage
1}&\text{Stage 2}&\text{Stage 3}\\\ \hline\cr\text{Stage 1}&0.63&0.18&0.19\\\
\text{Stage 2}&0.18&0.57&0.25\\\ \text{Stage 3}&0.00&0.00&1.00\\\
\end{array}\right).$ (151)

We will perform the provisioning scenario analysis by forecasting the Stage-
weighted provisions, given by (7.1.3) for the next four years. We first
calculate the $k-$step transition matrices:

$\displaystyle P^{2}=\left(\begin{array}[]{ccc}0.43&0.21&0.36\\\
0.21&0.36&0.43\\\ 0.00&0.00&1.00\\\
\end{array}\right),\,\,\,\,P^{3}=\left(\begin{array}[]{ccc}0.31&0.20&0.50\\\
0.20&0.24&0.56\\\ 0.00&0.00&1.00\\\
\end{array}\right),\,\,\,\,P^{4}=\left(\begin{array}[]{ccc}0.23&0.17&0.60\\\
0.17&0.18&0.66\\\ 0.00&0.00&1.00\\\ \end{array}\right).$ (161)

At time $t=0$ we consider the forward looking scenarios and can calculate the
Stage 1 and Stage 2 provisions. Recall that Stage 1 provisions are given by
(138). Stage 2 (expected lifetime) provisions are given in column $ECL_{t}$ in
Table 3 below, which also contains the Point-in-Time Stage 1 and Stage 2 PD
required to calculate the provisions.

|  Time until maturity ($u$)
---
$EAD_{u}$ | $\Psi^{1}(x,u)(\%)$ | $\Psi^{2}(x,u)(\%)$ |  Stage 1 $PD_{u}^{PiT}(\%)$ |  Stage 2 $PD_{u}^{PiT}(\%)$ | $LGD(\%)$ | $EL_{u}$
$10.0$ | 100 | 4.76 | 10.55 | 0.52 | 1.53 | 75% | 1.147
$9.0$ | 90 | 4.24 | 9.02 | 0.50 | 1.50 | 75% | 1.012
$8.0$ | 80 | 3.74 | 7.52 | 0.48 | 1.44 | 75% | 0.864
$7.0$ | 70 | 3.26 | 6.08 | 0.48 | 1.33 | 75% | 0.698
$6.0$ | 60 | 2.78 | 4.75 | 0.47 | 1.17 | 75% | 0.527
$5.0$ | 50 | 2.31 | 3.58 | 0.47 | 0.99 | 75% | 0.371
$4.0$ | 40 | 1.84 | 2.59 | 0.46 | 0.80 | 75% | 0.240
$3.0$ | 30 | 1.38 | 1.79 | 0.46 | 0.65 | 75% | 0.146
$2.0$ | 20 | 0.92 | 1.14 | 0.47 | 0.58 | 75% | 0.087
$1.0$ | 10 | 0.45 | 0.56 | 0.45 | 0.56 | 75% | 0.042

Table 3: Stage 1 and 2 PDs and expected losses.

As shown in Example 7.2, the Lifetime Provisions can be calculated by the sum
of the expected losses column, $EL_{u}$,. We can now calculate the Stage 1 and
Stage 2 provisions at each subsequent time period, which we will use to
calculate the Stage-weighted provisions for the next four years, calculated by
(7.1.3). The results are shown in Table 4, where the final column below
contains the Stage-weighted provisions.

|  Time until maturity ($u$)
---
$EAD_{u}$ | Stage 1 Provisions | $ECL$ | Stage 3 Provisions | $WP_{u}$
$10.0$ | 100 | 0.390 | 5.13 | 75.00 | 15.42
$9.0$ | 90 | 0.338 | 3.99 | 67.50 | 25.28
$8.0$ | 80 | 0.288 | 2.98 | 60.00 | 30.69
$7.0$ | 70 | 0.247 | 2.11 | 52.50 | 31.92

Table 4: Stage-weighted provision calculations for the next 4 years.

The large difference observed even between the Lifetime and Stage-weighted
provisions is evidence of the importance of such scenario analysis in
provision calculations. Particularly in cases similar to that examined in this
example, where the probability of transitioning to a default state is quite
high the results can have an extremely large effect, which risk managers must
account for in risk and provisioning policies. $\triangleleft$

###  7.2 Further Applications in credit risk modelling

####  7.2.1 Pricing of Credit Default Swaps

Another financial field in which the PD function plays a paramount role is
credit derivatives pricing. In particular, we consider the fair price of
Credit Default Swap (CDS). A default swap is a contract that protects the
holder of an underlying swap from the losses caused by the default to the
obligation’s issuer. Therefore, the evolution of the PD values can be used for
the pricing, hedging and managing of such options. Extensive work has been
done on modeling and pricing CDSs, such as in [18] and [35]. Specifically, it
can be shown that the price of the CDS is given by:

$\displaystyle
CDS=(1-R)\left(-\int_{0}^{T}e^{-rs}d\Phi(x,s)\right)-c\int_{0}^{T}e^{-rs}\Phi(x,s)ds,$

and the corresponding par spread:

$\displaystyle
c^{*}=\frac{(1-R)\Big{(}-\int_{0}^{T}e^{-rs}d\Phi(x,s)\Big{)}}{\int_{0}^{T}e^{-rs}\Phi(x,s)ds},$

where $R$ is the specific recovery rate and $r$ is the risk-free rate. The
above expression can be discetized as follows:

$\displaystyle
c^{*}=\frac{(1-R)\sum_{i=1}^{n}e^{-rt_{i}}(\Phi(x,t_{i-1})-\Phi(x,t_{i}))}{\frac{1}{2}\sum_{i=1}^{n}e^{-rt_{i}}(\Phi(x,t_{i-1})+\Phi(x,t_{i}))\Delta
t_{i}},$ (162)

where the Trapezoidal rule has been used for the discretization of the
denominator. Estimating the price and par rate of CDS therefore requires the
term structure of the underlying risk-free and survival probability processes.
We present a simplified example, whereby the interest rate is again considered
to be zero.

######  Example 7.4.

Consider a CDS with maturity $T=10$ years and recovery rate $R=0.5$, where the
asset process evolves according to the following stochastic volatility model:

$\displaystyle\begin{cases}dG_{t}=k(\theta-
G_{t})dt+\sigma(Y_{t})dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\,\,\ G_{0}=x\\\
dY_{t}=\kappa(\mu-Y_{t})dt+\xi\sqrt{Y_{t}}dW_{t},\,\,\,\ Y_{0}=y.\\\
\end{cases}$

The parameters of the stochastic model are given in Fig. 3, we let
$\mathcal{D}=[-5.0,75.0]$ and consider a spatial and temporal discretization
with 200 and 1000 steps respectively. Jumps are again normally distributed,
with size $Z\sim N(0.3,0.5)$, rate $\lambda=1.0$ and $J=90$. Furthermore, we
set a grid with 200 steps for the volatility $\mathcal{V}=[0.0,200.0]$. The
parameters and resulting graph of the PD function can be seen in Figure 3.

Figure 3: (Left) Evolution of the PD under the stochastic volatility model
with $(k,\theta,\kappa,\mu,\xi)=(2.0,2.0,0.05,0.1,0.07)$, for various values
of the starting volatility $Y_{0}$. (Right) The average survival probability
across all volatility values.

We assume an initial position of $x=3.0$ and plot the evolution of the average
survival probability in Fig. 4.

Figure 4: Evolution of the average survival probability under the stochastic
volatility model described in Example 7.4, with starting point $x=0.3$.

The resulting par spread is calculated using (162) to obtain $c^{*}=0.33$.

$\triangleleft$

####  7.2.2 Credit Portfolio Optimization

For many financial institutions, one of the most important tasks is
securitization of credit exposures. Ultimately, this can be formulated as an
optimization problem. The PD function affects the risk of each exposure and,
by extension, the corresponding return as well. To this end, we present a
simple example to show how such an optimization problem can be solved under
the stochastic volatility PD model.

######  Example 7.5.

Suppose a securitization agency creates a portfolio consisting of loans (or
credit derivatives), each with different underlying asset process. The
resulting PD functions will differ depending on the loan’s (or derivative’s)
characteristics and asset value processes. Slightly abusing notation, we
suppose that, for a portfolio of three loans, the corresponding PD functions
are given by $PD_{i},i=1,2,3$, estimated using the methodology developed
above.

The agency aims to select the investment allocated to each of the credit
exposures. Specifically, it poses the following portfolio optimization
problem: suppose $w_{i},i=1,2,3$ and $r_{i},i=1,2,3$ represent the weight of
total investment allocated to each institution’s set of loans and their
average return, respectively. Consider, furthermore, that the required
portfolio rate of return is set to be $R^{*}$. For the credit exposure $i$, at
time $t$, the expected loss is given by
$EL^{i}_{t}=EAD^{i}_{t}LGD^{i}_{t}PD^{i}_{t}$, and we can then define the
total loss function for the agency by $L(t)=\sum_{1}^{3}w_{i}EL^{i}_{t}$. In
order to rebalance the portfolio at each period, the securitization agency is
then interested in solving a portfolio optimization problem (we present a very
simple such problem, which can be solved analytically to illustrate the use of
the method). The optimization we consider is the following:

$\displaystyle\min_{\begin{subarray}{c}\bf
w\end{subarray}}\mathbb{E}[U(L(t))],\text{ subject to }$ $\displaystyle
w_{1}r_{1}+w_{2}r_{2}+w_{3}r_{3}=R,$ $\displaystyle w_{1}+w_{2}+w_{3}=1,$

for an appropriate loss function $U$. While the following analysis can be
extended to any convex loss function, for the sake of simplicity, we
illustrate the calculation selecting the quadratic loss function
$U(L)=bL^{2}-L$ (in the sense of a negative utility function). To standardize
the optimization problem, we consider that $EAD_{i}$ is given as a percentage
of the original loan value and for simplicity we consider a constant
$LGD^{i}_{t}=1$ for $i=1,2,3$ and all $t$. At any point in time $t$, the
expected loss utility is then:

$\sum_{i=1}^{3}U(w_{i}EAD_{i})PD^{PiT}_{i},$

where $PD^{PiT}_{i}$ is the current Point-in-Time default probability
corresponding to exposure $i$. The agency must optimize the portfolio by
solving:

$\displaystyle\text{minimize
}f(w_{1},w_{2},w_{3}):=\sum_{i=1}^{3}b(w_{i}EAD_{i})^{2}-w_{i}EAD_{i},\text{
subject to}$ $\displaystyle w_{1}r_{1}+w_{2}r_{2}+w_{3}r_{3}=R,$
$\displaystyle w_{1}+w_{2}+w_{3}=1,$

This simple quadratic optimization problem can now be solved either
analytically or numerically. It is straightforward to calculate:

$w_{3}^{*}=\frac{PD^{PiT}_{1}(2bEAD_{1}^{2}\delta\epsilon-\epsilon)+PD^{PiT}_{2}(\gamma-2baEAD_{2}^{2}\beta\gamma)-PD^{PiT}_{3}}{2(bPD^{PiT}_{1}EAD_{1}^{2}\epsilon^{2}+bPD^{PiT}_{2}EAD_{2}^{2}\gamma^{2}+bPD^{PiT}_{3}EAD_{3}^{2})},$

where
$\beta=\frac{R-r_{1}}{r_{2}-r_{1}},\gamma=\frac{r_{3}-r_{1}}{r_{2}-r_{1}},\delta=\frac{r_{2}-R}{r_{2}-r_{1}}$
and $\epsilon=\frac{r_{3}-r_{2}}{r_{2}-r_{1}}$. A straightforward substitution
using the two conditions will result in the corresponding values $w_{1}^{*}$
and $w_{2}^{*}$. We consider the above setting with average returns from each
institution’s instruments $r=(r_{1},r_{2},r_{3})^{T}$ and current exposures
$EAD=(EAD_{1},EAD_{2},EAD_{3})^{T}$ given by:

$\displaystyle r=(0.1\,\ 0.3\,\ 0.1)^{T},\,\,\,\,EAD=(0.9\,\ 0.8\,\ 0.7)^{T}.$

In order to obtain the vector containing the PD values, we consider the three
asset classes described by the processes below:

$\displaystyle
dG^{1}_{t}=k_{1}(\theta_{1}-G^{1}_{t})dt+\sigma_{1}dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\,\
x_{0}=1.00$ $\displaystyle
dG^{2}_{t}=k_{2}(\theta_{2}-G^{2}_{t})dt+\sigma_{2}dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\,\
x_{0}=0.20$ $\displaystyle
dG^{3}_{t}=k_{3}(\theta_{3}-G^{3}_{t})dt+\sigma_{3}dB_{t}+\int_{\mathbb{R}}zN(dt,dz),\,\,\
x_{0}=0.50,$

with
$(k_{1},k_{2},k_{3})=(0.5,0.8,0.5),(\theta_{1},\theta_{2},\theta_{3})=(3.5,3.0,2.5),(\sigma_{1},\sigma_{2},\sigma_{3})=(2.0,1.5,2.5)$,
$J=150,\lambda=1.0$ and jump distributions $Z\sim N(0.0,0.2)$ for all three.
Solving PIDE $(\ref{pide2})$, we obtain the $PD^{PiT}$ values:

$PD^{PiT}=(0.0783\,\ 0.1447\,\ 0.0447)^{T},$

and fixing the expected total return to be $R=25.00\%$, the resulting optimal
weights $w^{*}=(w_{1}^{*}\,\ w_{2}^{*}\,\ w_{3}^{*})^{T}$ are:

$w^{*}=(0.163\,\ 0.750\,\ 0.087)^{T}.$

$\triangleleft$

For extensive work on portfolio optimization problems with defaultable assets,
we refer the interested reader to e.g., [4]. Furthermore, empirical studies of
the applicability of standard ruin probabilities in practice can be found in
[15]. In the example above, we focus on a case where an agency must assess and
optimize a portfolio of loan exposures with varying characteristics. Such
cases could be loans originating in different sectors; e.g., in [48], the
authors consider a portfolio of risky bonds originating from the Industry and
Service sector.

##  8 Conclusion

In this paper we have focused on a generalized approach of estimating PD
values, accounting for both the cases of variable starting times and
maturities. We show that under certain conditions imposed on the models
representing the asset processes, these two cases can be dealt with
equivalently and lead to important novel representations of the PD function.
Specifically, with the integral equation approach, we can construct a robust
mathematical framework that allows us to develop both theoretical and
numerical tools for the calculation of the PD values. This methodology has
important advantages over standard Monte Carlo methods, as well as over
existing approaches using PIDEs, as we are able to consider sophisticated
models that incorporate multiple latent variables, without sacrificing
mathematical rigor for required regularity assumptions.

In terms of practicality and applications, the proposed framework covers many
of the difficulties financial institutions face due to the new regulatory
requirements for provision calculations, as well as continuous credit risk
monitoring for SICR events. We hypothesize that this approach could be useful
for practitioners, given that it constitutes a complete and efficient
modelling framework, with which one can calculate Point-in-Time and Lifetime
PD values, each of which are used extensively in credit risk management.
Specifically, this framework is motivated by the needs created by the IFRS 9
regulations, under which forecasting credit losses accurately and efficiently
is of paramount importance. We show how the PD estimations can be used to
calculate Stage 2 provisions, as well as more advanced, scenario-based
provisions, and of extensive further applications in credit risk modelling.

Finally, we note that this approach most likely is best fit for corporate and
small business loans, where the estimation of asset processes has been
documented in well-established work. Of course, it is possible that with new
developments in payment services and Open-Banking solutions (in accordance to
the Payment Services Directive 2), such methods could be applied to individual
consumers, given sufficient historical data. An example of recent work done in
this direction is [51]. To conclude, it is important to mention that the LGD
parameter is also of great importance for provision calculations; in this work
we considered a constant LGD, however, in practice, LGD values require
separate model development, often related to current macroeconomic variables,
as shown in [13]. Future research could focus on considering appropriate
models for the evolution of the LGD, in combination with the PD function.

##  Appendix A Stochastic processes

###  A.1 The continuous Ornstein Uhlenbeck process

In its simplest form, the OU process $X_{t}$ is defined as the stochastic
process satisfying the SDE:

$\displaystyle dX_{t}=k(\theta-X_{t})dt+\sigma dB_{t},\,\,X_{s}=x,$ (163)

for some known $x$, where, as above, $B_{t}$ represents the standard Brownian
motion and $k,\theta$ and $\sigma$ are positive real constants. The OU process
is a mean-reverting, Gaussian and Markov process, which is also temporally
homogeneous. We can therefore equivalently write (163) as:

$\displaystyle dX_{u}=k(\theta-X_{u})du+\sigma dB_{u},\,\,X_{0}=x,$ (164)

where $u=t-s$. For simplicity, we write $X_{t}^{x}$ to indicate the OU process
with $X_{0}=x$. We adopt this convention for all stochastic processes in the
remainder of this work. Employing Itô’s formula we can obtain the solution to
the above SDE:

$X_{t}=xe^{-kt}+\theta(1-e^{-kt})+\sigma\int_{0}^{t}e^{-k(t-u)}dB_{u},$

from which is it easy to see that $X_{t}\sim
N(\theta+(x-\theta)e^{-kt},\sigma^{2}(1-e^{-2kt})/2k).$ These properties are
what make this particular family of processes widely used in many
applications. We will also need the following regarding the transition density
and hitting time for the OU process.

######  Theorem A.1.

The transition density of the OU process, with initial condition $X_{0}=x$ is
given by:

$\displaystyle
p(y,x,t)\equiv\mathbb{P}\big{(}X_{t}=y|X_{0}=x\big{)}=\sqrt{\frac{k}{\pi\sigma^{2}(1-e^{-2kt})}}\exp\Big{(}-\frac{k(y-\theta-(x-\theta)e^{-kt})^{2}}{\sigma^{2}(1-e^{-2kt})}\Big{)}.$
(165)

Furthermore, for the OU process as defined in (163), we define the
corresponding survival probability distribution, given by
$Q(x,t):=\mathbb{P}\Big{(}\inf_{\begin{subarray}{c}r\leq
t\end{subarray}}G_{r}^{x}>0\Big{)}$. The distribution can be obtained via
appropriate Volterra equations, for which we refer the reader to [40]. We will
not employ this representation explicitly, however we will use the fact that
both $p(y,x,t)$ and $Q(x,t)$ are continuous in $(x,t)$. More details can be
found in Appendix B.

###  A.2 Continuous Time Markov Chain

In this section we outline, for the reader’s convenience, the background and
important results pertaining to Continuous Time Markov Chains (CTMC), which
are used for the regime-switching models.

######  Definition A.2.

A continuous time Markov chain is a continuous stochastic process $X_{t},t\geq
0$, with a discrete state space $\mathcal{R}$, of cardinality
$|\mathcal{R}|<\infty$, satisfying the Markov property, and such that:

$\displaystyle\mathbb{P}(X_{t+\delta}=j|X_{t}=i)=\begin{cases}q_{ij}\delta+o(\delta),\,\
i\neq j\\\ 1+q_{ii}\delta+o(\delta),\end{cases}$ (166)

as $\delta\downarrow 0$. In the above $q_{ij}$ are known as the transition
rates, for which we have $\sum_{i\in S}q_{ij}=0$, $q_{ij}\geq 0$ for $i\neq
j$.

The matrix $Q$ with entries $(Q)_{ij}=q_{ij}$, for $i,j\in\mathcal{R}$ is
known as the generator matrix of the Markov process (also referred to as the
transition rate matrix). Similar to the discrete time Markov chains, we can
define the transition matrix for a CTMS, $P(t),t\geq 0$, with entries:

$\displaystyle p_{ij}=\mathbb{P}(X_{t}=j|X_{0}=i),$ (167)

for $i,j\in\mathcal{R}$. The following result holds for the transition matrix,
from which we are also able to obtain a connection between the transition and
generator matrices.

######  Theorem A.3.

The transition matrix $P(t)$ satisfies the Kolmogorov forward equation:

$\displaystyle P^{\prime}(t)=P(t)Q,$

and hence:

$\displaystyle P(t)=e^{tQ}.$ (168)

Finally, we remind the reader that we say that a state $i$ is transient if,
given that the chain starts at $i$, it is possible, but not certain, that the
chain will return to $i$. Equivalently, there exists a non-zero probability
that the chain will never return to $i$. On the other hand, a state $i$ is
defined as absorbing if $\mathbb{P}(X_{t}=i|X_{0}=i)=1$, for all $t\geq 0$,
i.e. the probability of transitioning from $i$ to any other state is zero.

###  A.3 Lévy processes

Throughout this paper, we have abopted the notation used in [46]. We start by
defining a Lévy process:

######  Definition A.4 (Lévy process).

A Lévy process $\\{L_{t}\\}_{t\geq 0}$ is a stochastic process for which the
following conditions hold:

  * •

$L_{0}=0$.

  * •

$L$ has independent and stationary increments, i.e., if $t>s$ then
$L_{t}-L_{s}$ is independent from $L_{s}$ and
$L_{t}-L_{s}\stackrel{{\scriptstyle D}}{{=}}L_{t-s}$.

  * •

$L$ is stochastically continuous, i.e for all $\epsilon>0$ and all $s>0$ we
have

$\lim_{t\rightarrow s}\mathbb{P}(|X(t)-X(s)|>\epsilon)=0.$

A consequence of the above definition is the celebrated Itô - Lévy
decomposition. First, we define the following required quantities:

######  Definition A.5.

Let $L_{t}$ be a Lévy process, whose jump is defined as $\Delta
L_{t}=L_{t}-L_{t_{-}}$. Furthermore, let $\textbf{B}_{0}$ be the family of
Borel sets $U\subset\mathbb{R}$, whose closure does not contain 0. Then, for
$U\in\textbf{B}_{0}$, define the Poisson random measure of the Lévy process
$L_{t}$ by:

$N(t,U)=\sum_{0<s\leq t}\mathds{1}(\Delta L_{s}).$

The Poisson random measure represents the number of jumps of size $\Delta
L_{s}\in U$, which occur up to time $t$. We can therefore define the intensity
of the jumps as follows:

######  Definition A.6.

The intensity of a Lévy jump process $L_{t}$, known as the Lévy measure of
$L_{t}$ is defined as:

$\nu(U)=\mathbb{E}[N(1,U)],$

where, as above, $U\in\textbf{B}_{0}$.

A useful consequence of the above definitions is that if $\nu$ is the Lévy
measure of a simple Compound Poisson Process with rate $\lambda$ and jump size
density $f(z)$, then we have that

$\nu(U)=\lambda f(U).$

To this end, we will employ the following result in the subsequent sections,
due to [39]:

######  Theorem A.7.

Consider the Poisson random measure $N(t,U)$, with $U\in\textbf{B}_{0}$, and
corresponding Lévy measure $\nu(U)$. Then the process:

$X_{t}=\int_{0}^{t}\int_{B}zN(ds,dz),$

where $B\in\mathcal{B}(\mathbb{R})$, is a Compound Poisson Process with rate
$\nu(B)$ and jump distribution $\frac{\nu(dx)|_{B}}{\nu(B)}$.

In this paper, we will focus on jump terms in the form above. We can now
present the following theorem:

######  Theorem A.8 (Itô - Lévy decomposition).

Let $\\{L_{t}\\}_{t\geq 0}$ be a Lévy process. Then, we have

$\displaystyle L_{t}=bt+\sigma
B_{t}+\int_{|z|<1}z\tilde{N}(t,dz)+\int_{|z|\geq 1}zN(t,dz),$ (169)

for $t\geq 0$, where $b,\sigma\in\mathbb{R}$, $B_{t}$ is a Brownian motion and
$\tilde{N}(t,dz):=N(t,dz)-\nu(dz)t$ is the compensated Poisson measure.

More generally, we can define the stochastic process $X(t)$, as:

$\displaystyle
dX_{t}=a(t)dt+\sigma(t)dB(t)+\int_{|z|<1}H(t,z)\tilde{N}(dt,dz)+\int_{|z|\geq
1}H(t,z)N(dt,dz),$ (170)

known as Lévy - Itô processes. Moreover, by combining the compensator with the
drift term the above can be written as:

$\displaystyle
dX_{t}=a(t)dt+\sigma(t)dB(t)+\int_{z\in\mathbb{R}}H(t,z)N(dt,dz).$ (171)

We will adopt this formulation throughout the remainder of this work. For such
processes, we have the following results, which are extension of the standard
Itô and Generator formulas.

######  Theorem A.9 (Itô formula).

Let $X_{t}\in\mathbb{R}$ be an Itô- Lévy process and consider a function
$f(x,t)$, with $f\in C^{2}(\mathbb{R}\times[0,T])$. Then, the dynamics of the
process $f(X_{t},t)$ are given by the following version of the Itô formula:

$\displaystyle df(X_{t},t)=\frac{\partial f}{\partial
t}(X_{t},t)dt+\frac{\partial f}{\partial
x}(X_{t},t)\big{(}a(t)dt+\sigma(t)dB_{t}\big{)}+\frac{1}{2}\frac{\partial^{2}f}{\partial
x^{2}}(X_{t},t)\sigma^{2}(t)dt$
$\displaystyle+\int_{\mathbb{R}}\big{(}f\big{(}X_{t-}+H(t,z),t\big{)}-f(X_{t-},t)\big{)}N(dt,dz)$

######  Definition A.10 (Generator).

For a Lévy-Itô process, given by (171), and function
$f:\mathbb{R}\times[0,T]\rightarrow\mathbb{R}$ we define the generator
$\mathcal{L}$ by:

$\displaystyle\mathcal{L}f(x,t)=\lim_{t\downarrow
0}\frac{\mathbb{E}^{x}[f(X_{t},t)]-f(x,t)}{t},$ (172)

where $\mathbb{E}^{x}[f(X_{t},t)]=\mathbb{E}[f(X_{t},t)|X_{0}=x].$

It particular, it can be shown that the generator admits the following form:

$\displaystyle\mathcal{L}f(x,t)=\frac{\partial f}{\partial
t}+a(t)\frac{\partial f}{\partial
x}+\frac{1}{2}\sigma(t)^{2}\frac{\partial^{2}f}{\partial
x^{2}}+\int_{\mathbb{R}}\big{(}f(x+z,t)-f(x,t)\big{)}\nu(dz),$ (173)

##  Appendix B Infinitesimal generators and PDEs for the continuous Ornstein -
Uhlenbeck models

In this Appendix we detail the generators for the continuous counterparts of
the processes considered in this paper, built upon the continuous OU process
given by (164) in Appendix A.1. Firstly, we recall that the survival
probability $Q(x,t):=\mathbb{P}\Big{(}\inf_{\begin{subarray}{c}r\leq
t\end{subarray}}X_{r}^{x}>0\Big{)}$, as well as the corresponding transition
density $p(\cdot,x,t)$ satisfy the equation:

$\displaystyle\frac{\partial f}{\partial t}(x,t)=\mathcal{L}f(x,t),$ (174)

where $\mathcal{L}$ represents the operator:

$\displaystyle\mathcal{L}f(x,t)=k(\theta-x)\frac{\partial f}{\partial
t}(x,t)+\frac{1}{2}\sigma^{2}\frac{\partial f^{2}}{\partial x^{2}}(x,t).$
(175)

This is known as the Kolmogorov backward equation. In general, when
considering the survival probabilities and transition densities under the
continuous versions of the regime switching, stochastic volatility and
generalized models, analogous equations to (174) are produced. These depend on
the generators of the processes, which now include terms to capture the
evolution of the regime and/or volatility processes. For more details on the
generators of regime switching and stochastic volatility models see e.g., [31,
58]. Below we display the Kolmogorov backward equations under the continuous
regime switching, stochastic volatility and generalized models, respectively:

$\displaystyle\frac{\partial f}{\partial t}(x,\rho,t)$
$\displaystyle=\mathcal{L}_{1}Q(x,\rho,t)$
$\displaystyle:=k_{\rho}(\theta_{\rho}-x)\frac{\partial f}{\partial
x}(x,\rho,t)+\frac{1}{2}\sigma_{\rho}^{2}\frac{\partial^{2}f}{\partial
x^{2}}(x,\rho,t)+\sum_{j\neq\rho}q_{\rho j}\Big{(}Q(x,j,t)-f(x,\rho,t)\Big{)}$
(176) $\displaystyle\frac{\partial f}{\partial t}(x,y,t)$
$\displaystyle=\mathcal{L}_{2}f(x,y,t)$
$\displaystyle:=k(\theta-x)\frac{\partial f}{\partial
x}(x,y,t)+\kappa(\mu-y)\frac{\partial f}{\partial
y}(x,y,t)+\frac{1}{2}y\frac{\partial^{2}f}{\partial
x^{2}}(x,y,t)+\frac{1}{2}\xi^{2}y\frac{\partial^{2}Q}{\partial y^{2}}(x,y,t)$
(177) $\displaystyle\frac{\partial f}{\partial t}(x,\rho,y,t)$
$\displaystyle=\mathcal{L}_{3}f(x,\rho,y,u):=k_{\rho}(\theta_{\rho}-x)\frac{\partial
f}{\partial x}(x,\rho,y,t)+\kappa(\mu-y)\frac{\partial f}{\partial
y}(x,\rho,y,t)$
$\displaystyle+\frac{1}{2}\sigma_{\rho}^{2}y\frac{\partial^{2}f}{\partial
x^{2}}(x,\rho,y,t)+\frac{1}{2}\xi^{2}y\frac{\partial^{2}f}{\partial
y^{2}}(x,\rho,y,t)+\sum_{j\neq\rho}q_{\rho
j}\Big{(}f(x,j,y,t)-f(x,\rho,y,t)\Big{)},$ (178)

where we define separately the operators $\mathcal{L}_{1},\mathcal{L}_{2}$ and
$\mathcal{L}_{3}$ for notational convenience. Note that we write the
dependency on the regime as a subscript on the right hand side of the
equations above.

###  B.1 Regularity of solutions to parabolic PDEs

We will require the following results pertaining to the regularity of
solutions of the second order parabolic PDE (174). The following are due to
[27]. We first define some relevant function spaces that will be required for
the subsequent regularity results.

######  Definition B.1.

Consider $\Omega\subset\mathbb{R}^{n}$ an open set, with closure
$\bar{\Omega}$. Furthermore, consider a fixed time horizon $T>0$ and define
$Q_{T}=\Omega\times[0,T]$, with closure $\bar{Q}_{T}$. We then define the
following spaces, for $0<\alpha<1$:

  * •

$C^{0}(\bar{\Omega})$ is the Banach space of bounded continuous functions in
$\bar{\Omega}$, with the natural supremum norm:

$\|\cdot\|_{C^{0}(\bar{\Omega})}\equiv\|\cdot\|_{0,\bar{\Omega}}=\sup_{\Omega}|\cdot|$

  * •

$C^{2,1}\left(\bar{Q}_{T}\right)$ is the Banach space of functions
$\varphi(x,t)$ belonging to $C^{0}\left(\bar{Q}_{T}\right)$ together their
derivatives $\frac{\partial f}{\partial x},\frac{\partial^{2}f}{\partial
x^{2}},\frac{\partial f}{\partial t}$ in $\bar{Q}_{T}$ with natural norm.

  * •

$C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)$ is the Banach space of
function $\varphi$ in $C^{0}\left(\bar{Q}_{T}\right)$ which are Hölder
continuous in $\bar{Q}_{T}$ with exponent $\alpha$ in $x$ and
$\frac{\alpha}{2}$ in $t$ i.e. having a finite value for the seminorm

$\langle f\rangle_{\bar{Q}_{T}}^{(\alpha)}\equiv\langle
f\rangle_{x,\bar{Q}_{T}}^{(\alpha)}+\langle
f\rangle_{t,\bar{Q}_{T}}^{\left(\frac{\alpha}{2}\right)}$

where

$\displaystyle\langle f\rangle_{x,\bar{Q}_{T}}^{(\alpha)}=\inf\left\\{C\geq
0:\left|f(x,t)-f\left(x^{\prime},t\right)\right|\leq
C\left|x-x^{\prime}\right|^{\alpha},\forall x,x^{\prime},t\right\\}$
$\displaystyle\langle
f\rangle_{t,\bar{Q}_{T}}^{\left(\frac{\alpha}{2}\right)}=\inf\left\\{C\geq
0:\left|f(x,t)-f\left(x,t^{\prime}\right)\right|\leq
C\left|t-t^{\prime}\right|^{\frac{\alpha}{2}},\forall x,t,t^{\prime}\right\\}$

The quantity

$\|f\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}\equiv\|f\|_{\alpha,\bar{Q}_{T}}=\|f\|_{0,\bar{Q}_{T}}+\langle
f\rangle_{\bar{Q}_{T}}^{(\alpha)}$

defines a norm.

  * •

$C^{2+\alpha,\frac{2+\alpha}{2}}\left(\bar{Q}_{T}\right)$ is the Banach space
of functions $f(x,t)$ in $C^{2,1}\left(\bar{Q}_{T}\right)$ having a finite
value for the seminorm:

$\langle
f\rangle_{\bar{Q}_{T}}^{(2+\alpha)}=\left\langle\partial_{t}f\right\rangle_{\bar{Q}_{T}}^{(\alpha)}+\sum_{i,j=1}^{d}\left\langle\partial_{ij}f\right\rangle_{\bar{Q}_{T}}^{(\alpha)}+\sum_{i=1}^{d}\left\langle\partial_{i}f\right\rangle_{t,\bar{Q}_{T}}^{\frac{1+\alpha}{2}}.$

Then, the quantity

$\|f\|_{C^{2+\alpha,\frac{2+\alpha}{2}}\left(\bar{Q}_{T}\right)}\equiv\|f\|_{2+\alpha,\bar{Q}_{T}}=\sum_{2r+s\leq
2}\left\|\partial_{t}^{r}\partial_{x}^{s}f\right\|_{0,\bar{Q}_{T}}+\langle
f\rangle_{\bar{Q}_{T}}^{(2+\alpha)}$

defines a norm.

######  Theorem B.2.

Consider a bounded domain $\Omega$, the operator $L:=\frac{\partial
f}{\partial t}(x,t)-\mathcal{L}f(x,t)$ and the PDE:

$\displaystyle\begin{cases}Lf=g(x,t)&\text{ for }(x,t)\in Q_{T}\\\
f(x,0)=\varphi(x)&\text{ for }x\in\Omega\\\ f(x,t)=\psi(x,t)&\text{ for
}x\in\Sigma_{T}:=\partial\Omega\times[0,T].\end{cases}$ (179)

Then, for any $g\in
C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right),\varphi\in
C^{2+\alpha}(\bar{\Omega}),\psi\in
C^{2+\alpha,\frac{2+\alpha}{2}}\left(\Sigma_{T}\right)$, with $0<a<1$, (188)
has a unique solution from the class
$C^{2+\alpha,\frac{2+\alpha}{2}}\left(\bar{Q}_{T}\right)$ and satisfies the
inequality:

$\|f\|_{2+\alpha,\bar{Q}_{T}}\leq
C\left(\|g\|_{\alpha,\bar{Q}_{T}}+\|\varphi\|_{2+\alpha,\bar{\Omega}}+\|\psi\|_{2+\alpha,\Sigma_{T}}\right),$

with the constant $C$ independent of $f,\varphi$ and $\psi$.

When extending to Lévy models and the corresponding integro-differential
equations, we will need the following result.

######  Lemma B.3.

Consider $f\in C^{\alpha+2,\frac{2+\alpha}{2}}\left(\bar{Q}_{T}\right)$ and
the differential operator:

$If(x,t)=\int_{\Omega}[f(x+z,t)-f(x,t)]\nu(dz).$

Then, for $0<a<1$, we have that:

$\|If\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}\leq\varepsilon\|\nabla
f\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}+C(\varepsilon)\|f\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}.$

Note that Lemma F.2 is a simplified version of the corresponding results in
[27], where the authors consider additional integral operators of higher
orders.

##  Appendix C Kolmogorov equations for regime switching and stochastic
volatility models

Below we recall the Kolmogorov backward equations under the continuous regime
switching, stochastic volatility and generalized models, respectively. These
depend on the generators of the processes, which now include terms to capture
the evolution of the regime and/or volatility processes.

$\displaystyle\frac{\partial f}{\partial t}(x,\rho,t)$
$\displaystyle=\mathcal{L}_{1}Q(x,\rho,t):=$ $\displaystyle
k_{\rho}(\theta_{\rho}-x)\frac{\partial f}{\partial
x}(x,\rho,t)+\frac{1}{2}\sigma_{\rho}^{2}\frac{\partial^{2}f}{\partial
x^{2}}(x,\rho,t)+\sum_{j\neq\rho}q_{\rho j}\Big{(}Q(x,j,t)-f(x,\rho,t)\Big{)}$
(180) $\displaystyle\frac{\partial f}{\partial t}(x,y,t)$
$\displaystyle=\mathcal{L}_{2}f(x,y,t):=$ $\displaystyle k(\theta-$
$\displaystyle x)\frac{\partial f}{\partial
x}(x,y,t)+\kappa(\mu-y)\frac{\partial f}{\partial
y}(x,y,t)+\frac{1}{2}y\frac{\partial^{2}f}{\partial
x^{2}}(x,y,t)+\frac{1}{2}\xi^{2}y\frac{\partial^{2}Q}{\partial y^{2}}(x,y,t)$
(181) $\displaystyle\frac{\partial f}{\partial t}(x,\rho,y,t)$
$\displaystyle=\mathcal{L}_{3}f(x,\rho,y,u):=k_{\rho}(\theta_{\rho}-x)\frac{\partial
f}{\partial x}(x,\rho,y,t)+\kappa(\mu-y)\frac{\partial f}{\partial
y}(x,\rho,y,t)$ $\displaystyle+\frac{1}{2}\sigma_{\rho}^{2}$ $\displaystyle
y\frac{\partial^{2}f}{\partial
x^{2}}(x,\rho,y,t)+\frac{1}{2}\xi^{2}y\frac{\partial^{2}f}{\partial
y^{2}}(x,\rho,y,t)+\sum_{j\neq\rho}q_{\rho
j}\Big{(}f(x,j,y,t)-f(x,\rho,y,t)\Big{)},$ (182)

where we define separately the generator operators
$\mathcal{L}_{1},\mathcal{L}_{2}$ and $\mathcal{L}_{3}$, for notational
convenience. For more details on the generators of regime switching and
stochastic volatility models see e.g., [31, 58].

##  Appendix D PIDEs for the PD functions in Sobolev spaces

For completeness, we first recall some basic definitions pertaining to the
theory of weak derivatives in Sobolev spaces.

######  Definition D.1.

(Weak derivative) Consider an open subset $\Omega\subset\mathbb{R}^{n}$ and
the space of continuous functions which are $k$ times continuously
differentiable, for $k=1,2,\dots,$ denoted by $C^{k}(\Omega)$ and
$L^{1}_{loc}(\Omega)$ the space of locally integrable functions. Furthermore,
let $\alpha=(\alpha_{1},\dots,\alpha_{n})$ be a multi-index, with order
$|\alpha|:=\sum_{i}\alpha_{i}$, and denote $D^{a}u$ by:

$D^{\alpha}u=\frac{\partial^{|\alpha|}u}{\partial
x_{1}^{\alpha_{1}}\ldots\partial
x_{n}^{\alpha_{n}}}=\frac{\partial^{\alpha_{1}}}{\partial
x_{1}^{\alpha_{1}}}\ldots\frac{\partial^{\alpha_{n}}}{\partial
x_{n}^{\alpha_{n}}}u.$

Then, for $f\in L^{1}_{\text{loc}}(\Omega)$ we define $u\in
L^{1}_{\text{loc}}(\Omega)$ to be the $\alpha$th weak derivative of $f$,
$D^{\alpha}f=u$, if:

$\int_{\Omega}fD^{\alpha}\varphi dx=(-1)^{|\alpha|}\int_{\Omega}u\varphi dx,$

for every smooth test functon with compact support, $\varphi$.

With $\mathcal{Q}:=\mathcal{D}\times\mathcal{V}\times[0,T]$, we are interested
in functions which are twice weakly differentiable with respect to the initial
condition and once with respect to the time until maturity. Hence, we can
therefore work in the Sobolev space containing all such functions
$W^{2,1}(\mathcal{Q})=\left\\{f\in L^{p}(\mathcal{Q}):D^{\alpha}f\in
L^{1}(\mathcal{Q}),|\alpha|\leqslant 2\right\\}$.

To work in this space, we also need an appropriate weak version of the Itô
formula, pertaining to cases where the underlying function may not enjoy the
required regularity properties; these results are given by Theorems D.2 and
D.3, due to [38] and [45], respectively. We include the results below, for
completeness:

######  Theorem D.2.

Consider the stochastic process

$dX_{t}=a(x,t)dt+\sigma(x,t)dB_{t}$

and a region $\mathcal{Q}$, where $B_{t}$ is a standard Brownian motion, with
function $f$ such that function $f\in W^{2,1}(\mathcal{Q})$. Moreover, let
$\tau$ be some Markov time such that $\tau<\tau_{\mathcal{Q}}$, where $t_{Q}$
is the exit time of the process from the region $\mathcal{Q}$. Then, if there
exists some constant $K$ such that $|\sigma(x,t)|+|a(x,t)|\leq K$, for some
fixed time $s$ we have that:

$\displaystyle f(X_{\tau},s+\tau)$
$\displaystyle-f(X_{t},s+t)=\int_{t}^{\tau}\frac{\partial f}{\partial
u}(X_{u},s+u)du$ $\displaystyle+\int_{t}^{\tau}\frac{\partial f}{\partial
x}(X_{u},s+u)\sigma(X_{u},u)dB(t)+\frac{1}{2}\int_{t}^{\tau}\frac{\partial^{2}f(X_{u},s+u)}{\partial
x^{2}}\sigma^{2}(X_{u},u)du,$ (183)

almost surely.

######  Theorem D.3.

Consider the stochastic process with representation

$X_{t}=\gamma t+\int_{0}^{t}\int_{\mathbb{R}}zN(dz,du),$

where $\gamma\in\mathbb{R}$. Assume $f:\mathcal{Q}\rightarrow\mathbb{R}$ is a
continuous function on $U\mathcal{Q}$ such that $f\in
L_{loc}^{1}(\mathcal{Q})$, i.e., $f$ is locally integrable. Furthermore,
assume the existence of locally bounded weak first order derivatives, as
defined in D.1. Then:

$\displaystyle f(X_{t},t)=f(X_{0},0)+\int_{0}^{t}\frac{\partial f}{\partial
s}(X_{u},u)du+\gamma\int_{0}^{t}\frac{\partial f}{\partial x}(X_{u},u)du+$
$\displaystyle\quad+\int_{0}^{t}\int_{\mathbb{R}}f(X_{u-}+z,u)-f(X_{u-},u)N(du,dz),$
(184)

where all derivatives are understood in the weak sense.

In the Sobolev setting, we can now derive a PIDE using the common approach of
appropriate martingale arguments (similar analysis has been given in e.g.,
[44]).

######  Lemma D.4.

The survival probability with a variable starting time and fixed maturity $T$,
$\Phi(x,\rho,y,s;T)$ satisfies the partial integro-differential equation,
almost surely:

$\displaystyle\frac{\partial\Phi}{\partial
s}(x,\rho,y,s;T)+\mathcal{L}_{3}\Phi(x,\rho,y,s;T)+\int_{\mathbb{R}}\Big{(}\Phi(x+z,\rho,y,s;T)-\Phi(x,\rho,y,s;T)\Big{)}\nu(dz)=0,$
(185)

for $(x,\rho,y,s)\in\mathcal{D}\times\mathcal{R}\times\mathcal{V}\times[0,T]$,
with initial and boundary conditions:

$\displaystyle\Phi(x,\rho,y,T;T)=\mathds{1}_{\\{x>0\\}},\,\,\
(x,\rho,y)\in\mathcal{D}\times\mathcal{R}\times\mathcal{V},$
$\displaystyle\Phi(0,\rho,y,s;T)=0,\,\,\
(\rho,y,s)\in\mathcal{R}\times\mathcal{V}\times[0,T],$
$\displaystyle\Phi(x,\rho,y,s;T)\rightarrow 1,\text{ as
}x\rightarrow\infty,\,\,\
(\rho,y,s)\in\mathcal{R}\times\mathcal{V}\times[0,T],$
$\displaystyle\frac{\partial\Phi}{\partial y}(x,\rho,y,s;T)=0,\text{ as
}y\rightarrow\infty\,\,\
(x,\rho,s)\in\mathcal{D}\times\mathcal{R}\times[0,T],$ (186)

with the generator operator $\mathcal{L}_{3}$ as given in (182).

###### Proof.

We begin by considering the dynamics of the survival probability. As $\Phi$ is
differentiable in $W^{2,1}(Q)$, we will employ Theorems D.2 and D.3 above. We
then obtain:

$\displaystyle\Phi(G_{w},R_{w},$ $\displaystyle
Y_{w},w)-\Phi(G_{s},\rho,y,s)=\int_{s}^{w}\Big{(}\frac{\partial\Phi}{\partial
r}(G_{r},R_{r},Y_{r},r)+\mathcal{L}_{3}\Phi(x,\rho,y,s)\Big{)}dr$
$\displaystyle+\int_{s}^{w}\sigma(R_{r})\sqrt{Y_{r}}\frac{\partial\Phi}{\partial
x}(G_{r},R_{r},Y_{r},r)dB_{r}+\int_{s}^{w}\xi\sqrt{Y_{r}}\frac{\partial\Phi}{\partial
x}(G_{r},R_{r},Y_{r},r)dW_{r}$
$\displaystyle+\int_{s}^{w}\int_{\mathbb{R}}\Big{(}\Phi(G_{r}+z,R_{r},Y_{r},r)-\Phi(G_{r},R_{r},Y_{r},r)\Big{)}N(dr,dz),\quad\quad$
(187)

where the derivatives are understood in the weak sense in accordance to
definition D.1. Also note that we omit the dependence on the $t$ parameter,
for brevity. We are now able to formulate the following result regarding the
survival probability. We write the dynamics above in terms of the compensated
Poisson measure $\tilde{N}(dt,dz)=N(dt,dz)-\nu(dz)dt$. The last term then
becomes:

$\int_{s}^{w}\int_{\mathbb{R}}\Big{(}\Phi(G_{r}+z,R_{r},Y_{r},r)-\Phi(G_{r},R_{r},Y_{r},r)\Big{)}(\tilde{N}(dr,dz)+\nu(dz)dr).$

Combining with the dynamics of $\Phi$ above and using the fact that the sum of
the non-martingale quantities must be identically zero, we obtain PIDE (185),
as required. The boundary conditions follow by definition of the survival
probability. ∎

##  Appendix E Existence and continuity of the PD function

######  Theorem E.1.

Arzelà-Ascoli. Let $(X,d)$ be a compact metric space and $C(X)$ the space of
continuous functions on $X$. Then, if a sequence of continuous functions
$\\{f\\}_{n=1}^{\infty}$ in $C(X)$ is bounded and equicontinuous it has a
uniformly convergent subsequence.

######  Theorem E.2.

Schauder Fixed Point Let $(X,\|\cdot\|)$ be a Banach space and $S\subset X$ is
compact, convex, and nonempty. Any continuous operator $A:S\rightarrow S$ has
at least one fixed point.

##  Appendix F Regularity of solutions to parabolic PDEs

We will require the following results pertaining to the regularity of
solutions of the second order parabolic PDE (174). The following are due to
[27].

######  Theorem F.1.

Consider a bounded domain $\Omega$, the operator $L:=\frac{\partial
f}{\partial t}(x,t)-\mathcal{L}f(x,t)$, where $\mathcal{L}$ is the generator
operator, and the PDE:

$\displaystyle\begin{cases}Lf=g(x,t)&\text{ for }(x,t)\in Q_{T}\\\
f(x,0)=\varphi(x)&\text{ for }x\in\Omega\\\ f(x,t)=\psi(x,t)&\text{ for
}x\in\Sigma_{T}:=\partial\Omega\times[0,T].\end{cases}$ (188)

Then, for any $g\in
C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right),\varphi\in
C^{2+\alpha}(\bar{\Omega}),\psi\in
C^{2+\alpha,\frac{2+\alpha}{2}}\left(\Sigma_{T}\right)$, with $0<a<1$, (188)
has a unique solution from the class
$C^{2+\alpha,\frac{2+\alpha}{2}}\left(\bar{Q}_{T}\right)$ and satisfies the
inequality:

$\|f\|_{2+\alpha,\bar{Q}_{T}}\leq
C\left(\|g\|_{\alpha,\bar{Q}_{T}}+\|\varphi\|_{2+\alpha,\bar{\Omega}}+\|\psi\|_{2+\alpha,\Sigma_{T}}\right),$

with the constant $C$ independent of $f,\varphi$ and $\psi$.

When extending to Lévy models and the corresponding integro-differential
equations, we will need the following result.

######  Lemma F.2.

Consider $f\in C^{\alpha+2,\frac{2+\alpha}{2}}\left(\bar{Q}_{T}\right)$ and
the differential operator:

$If(x,t)=\int_{\Omega}[f(x+z,t)-f(x,t)]\nu(dz).$

Then, for $0<a<1$, we have that:

$\|If\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}\leq\varepsilon\|\nabla
f\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}+C(\varepsilon)\|f\|_{C^{\alpha,\frac{\alpha}{2}}\left(\bar{Q}_{T}\right)}.$

Note that Lemma F.2 is a simplified version of the corresponding results in
[27], where the authors consider additional integral operators of higher
orders.

## References

  * [1] Damien Ackerer, Damir Filipović, and Sergio Pulido.  The jacobi stochastic volatility model.  Finance and Stochastics, 22:667–700, 2018.
  * [2] David Applebaum.  Lévy processes and stochastic calculus.  Cambridge university press, 2009.
  * [3] Chrystalleni Aristidou.  The meta-phillips curve: Modelling us inflation in the presence of regime change.  Journal of Macroeconomics, 57:367–379, 2018.
  * [4] Sujith Asanga, Alexandru Asimit, Alexandru Badescu, and Steven Haberman.  Portfolio optimization under solvency constraints: a dynamical approach.  North American Actuarial Journal, 18(3):394–416, 2014.
  * [5] Laura Ballotta and Efrem Bonfiglioli.  Multivariate asset models using lévy processes and applications.  The European Journal of Finance, 22(13):1320–1350, 2016.
  * [6] Laura Ballotta and Gianluca Fusai.  Counterparty credit risk in a multivariate structural model with jumps.  Finance, 36(1):39–74, 2015.
  * [7] Guy Barles and Cyril Imbert.  Second-order elliptic integro-differential equations: viscosity solutions’ theory revisited.  In Annales de l’IHP Analyse non linéaire, volume 25, pages 567–585, 2008.
  * [8] Ole E Barndorff-Nielsen and Neil Shephard.  Modelling by lévy processess for financial econometrics.  Lévy processes, pages 283–318, 2001.
  * [9] Ole E Barndorff-Nielsen and Neil Shephard.  Non-gaussian ornstein–uhlenbeck-based models and some of their uses in financial economics.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(2):167–241, 2001.
  * [10] David S Bates.  Jumps and stochastic volatility: Exchange rate processes implicit in deutsche mark options.  The Review of Financial Studies, 9(1):69–107, 1996.
  * [11] Dirk Beerbaum.  Significant increase in credit risk according to ifrs 9: Implications for financial institutions.  International Journal of Economics & Management Sciences, 4(9):1–3, 2015.
  * [12] T Belkina and Yu Kabanov.  Viscosity solutions of integro-differential equations for nonruin probabilities.  Theory of Probability & Its Applications, 60(4):671–679, 2016\.
  * [13] Tony Bellotti and Jonathan Crook.  Loss given default models incorporating macroeconomic variables for credit cards.  International Journal of Forecasting, 28(1):171–182, 2012.
  * [14] Eric Benhamou, Emmanuel Gobet, and Mohammed Miri.  Time dependent heston model.  SIAM Journal on Financial Mathematics, 1(1):289–325, 2010.
  * [15] Alexander Braun, Hato Schmeiser, and Florian Schreiber.  Solvency ii’s market risk standard formula: How credible is the proclaimed ruin probability?  Journal of insurance issues, pages 1–30, 2015.
  * [16] Alexandre Brouste, Masaaki Fukasawa, Hideitsu Hino, Stefano M Iacus, Kengo Kamatani, Yuta Koike, Ryosuke Nomura, Y Shimizu, M Uchida, and N Yoshida.  The yuima project: A computational framework for simulation and inference of stochastic differential equations.  American statistical association, 2014.
  * [17] Max Bruche.  Estimating structural bond pricing models via simulated maximum likelihood.  2005\.
  * [18] Jessica Cariboni and Wim Schoutens.  Pricing credit default swaps under lévy models.  Journal of Computational Finance, 10(4):71, 2007.
  * [19] Peter Carr and Anita Mayo.  On the numerical evaluation of option prices in jump diffusion processes.  European Journal of Finance, 13(4):353–372, 2007.
  * [20] Somnath Chatterjee.  Centre for central banking studies.  Modelling credit risk. Bank of England, 2015.
  * [21] Rama Cont and Ekaterina Voltchkova.  A finite difference scheme for option pricing in jump diffusion and exponential lévy models.  SIAM Journal on Numerical Analysis, 43(4):1596–1626, 2005.
  * [22] Rama Cont and Ekaterina Voltchkova.  Integro-differential equations for option prices in exponential lévy models.  Finance and Stochastics, 9(3):299–325, 2005.
  * [23] John C Cox, Jonathan E Ingersoll Jr, and Stephen A Ross.  A theory of the term structure of interest rates.  In Theory of valuation, pages 129–164. World Scientific, 2005.
  * [24] Michael G Crandall and Pierre-Louis Lions.  Viscosity solutions of hamilton-jacobi equations.  Transactions of the American mathematical society, 277(1):1–42, 1983.
  * [25] Yann d’Halluin, Peter A Forsyth, and Kenneth R Vetzal.  Robust numerical methods for contingent claims under jump diffusion processes.  IMA Journal of Numerical Analysis, 25(1):87–112, 2005.
  * [26] Jin-Chuan Duan, Ivilina Popova, and Peter Ritchken.  Option pricing under regime switching.  Quantitative Finance, 2(2):116, 2002.
  * [27] Maria Giovanna Garroni, José Luis Menaldi, et al.  Green functions for second order parabolic integro-differential problems, volume 275.  Chapman and Hall/CRC, 1992.
  * [28] K Georgiou, GN Domazakis, D Pappas, and AN Yannacopoulos.  Markov chain lumpability and applications to credit risk modelling in compliance with the international financial reporting standard 9 framework.  European Journal of Operational Research, 292(3):1146–1164, 2021\.
  * [29] Stéphane Goutte.  Pricing and hedging in stochastic volatility regime switching models.  2013\.
  * [30] Stéphane Goutte and Armand Ngoupeyou.  Conditional laplace formula in regime switching model: Application to defaultable bond.  Preprint, 2011.
  * [31] Donatien Hainaut.  Financial modeling with switching lévy processes.  ESC, 2011.
  * [32] Donatien Hainaut and Olivier Le Courtois.  An intensity model for credit risk with switching lévy processes.  Quantitative Finance, 14(8):1453–1465, 2014.
  * [33] Said Hamadène and Marie-Amelie Morlais.  Viscosity solutions for second order integro-differential equations without monotonicity condition: the probabilistic approach.  Stochastics, 88(4):632–649, 2016.
  * [34] James D Hamilton.  Regime switching models.  In Macroeconometrics and time series analysis, pages 202–209. Springer, 2010.
  * [35] Patrick Houweling and Ton Vorst.  Pricing default swaps: Empirical evidence.  Journal of international Money and Finance, 24(8):1200–1225, 2005\.
  * [36] John Hull, Izzy Nelken, and Alan White.  Merton’s model, credit risk, and volatility skews.  Journal of Credit Risk Volume, 1(1):05, 2004.
  * [37] Dah-Jing Jwo.  Investigation of phase-locked loop statistics via numerical implementation of the fokker–planck equation.  Applied Sciences, 10(7):2625, 2020.
  * [38] Nikolaj Vladimirovič Krylov.  Controlled diffusion processes, volume 14.  Springer Science & Business Media, 2008.
  * [39] Andreas E Kyprianou.  Introductory lectures on fluctuations of Lévy processes with applications.  Springer Science & Business Media, 2006.
  * [40] Alexander Lipton and Vadim Kaushansky.  On the first hitting time density of an ornstein-uhlenbeck process.  arXiv preprint arXiv:1810.02390, 2018.
  * [41] Elisa Luciano and Wim Schoutens.  A multivariate jump-driven financial asset model.  Quantitative finance, 6(5):385–402, 2006.
  * [42] Yuliya Mishura and Olena Ragulina.  Ruin Probabilities: Smoothness, Bounds, Supermartingale Approach.  Elsevier, 2016.
  * [43] Peter Miu and Bogie Ozdemir.  Basel requirement of downturn lgd: Modeling and estimating pd & lgd correlations.  Journal of Credit Risk, 2(2):43–68, 2006.
  * [44] Christian Max Møller.  Stochastic differential equations for ruin probabilities.  Journal of applied probability, pages 74–89, 1995.
  * [45] Ramin Okhrati and Uwe Schmock.  Ito’s formula for finite variation levy processes: The case of non-smooth functions.  Journal of Mathematical Analysis and Applications, 430(2):1163–1174, 2015.
  * [46] Bernt Karsten Øksendal and Agnes Sulem.  Applied stocstic control of jump diffusions, volume 498.  Springer, 2007.
  * [47] O Onalan.  Financial modelling with ornstein-uhlenbeck processes driven by lévy process.  In Proceedings of the World Congress on Engineering, volume 2, pages 1–3, 2009.
  * [48] Puneet Pasricha, Dharmaraja Selvamuthu, Guglielmo D’Amico, and Raimondo Manca.  Portfolio optimization of credit risky bonds: a semi-markov process approach.  Financial Innovation, 6:1–14, 2020.
  * [49] Halis Sak and Wolfgang Hörmann.  Fast simulations in credit risk.  Quantitative Finance, 12(10):1557–1569, 2012.
  * [50] Wim Schoutens and Jessica Cariboni.  Lévy processes in credit risk, volume 519.  John Wiley & Sons, 2010.
  * [51] Ellen Tobback and David Martens.  Retail credit scoring using fine-grained payment data.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 182(4):1227–1246, 2019.
  * [52] Jari Toivanen.  A componentwise splitting method for pricing american options under the bates model.  In Applied and numerical partial differential equations, pages 213–227. Springer, 2010.
  * [53] Oldrich Vasicek.  An equilibrium characterization of the term structure.  Journal of financial economics, 5(2):177–188, 1977.
  * [54] Kimmo Virolainen.  Macro stress testing with a macroeconomic credit risk model for finland.  Bank of Finland Research Discussion Paper, (18), 2004.
  * [55] Jiri Witzany.  A two-factor model for pd and lgd correlation.  Available at SSRN 1476305, 2011.
  * [56] Xin Xu.  Estimating lifetime expected credit losses under ifrs 9.  Available at SSRN 2758513, 2016.
  * [57] Chunsheng Zhou.  A jump-diffusion approach to modeling credit risk and valuing defaultable securities.  Available at SSRN 39800, 1997.
  * [58] Chao Zhu, George Yin, and Nicholas A Baran.  Feynman–kac formulas for regime-switching jump diffusions and their applications.  Stochastics An International Journal of Probability and Stochastic Processes, 87(6):1000–1032, 2015.


$Z_{\mathrm{eff}}$ factor. The re-scaling ensures that the overlap integrals
of the wavefunctions are consistent with the TB parameters that reproduce the
band structure of graphene listed in Tab. 2. For completeness, we list these
wavefunctions in App. D.1.

In contrast, we model the carbon wavefunctions using Roothaan-Hartree-Fock
(RHF) wavefunctions Bunge  _et al._ (1993), which we summarize in App. D.2. By
comparison, we find that the re-scaled hydrogenic wavefunctions are a poor
approximation for the required ground state wavefunctions of atomic electrons
in carbon, as can be seen in Fig. 11.

While the RHF wavefunctions are a better description of carbon electrons, it
should be noted that their overlap integrals are not consistent with the
overlap parameters for graphene in Tab. 2. A re-scaling of the RHF
wavefunctions similar to the approach by Hochberg et al. is possible, but
spoils the accurate description of electrons in atomic carbon. However, this
seems to be a general feature of evaluating electron wavefunctions in the TB
approximation.

###  E.3 Electron ejection rate

Figure 12: Comparison of the energy spectrum between our TB results and the
results by Hochberg et al. (spectrum digitized from Fig. 1 of Hochberg  _et
al._ (2017)). The solid lines show the total spectrum, whereas the dashed
(dotted) lines show the contributions of the $\pi$\- ($\sigma$-) electrons.
For this figure, we used the SHM parameters from Hochberg  _et al._ (2017) in
order to facilitate the comparison (Hochberg et al. use $v_{0}=220\text{
km/s}$ as opposed to our choice of $v_{0}=238\text{ km/s}$.)

In Hochberg  _et al._ (2017), the total rate of DM-induced electron ejections
in graphene is given as

$\displaystyle R$
$\displaystyle=2\frac{\rho_{\chi}}{m_{\chi}}N_{C}A_{\mathrm{uc}}\sum_{i}\int\frac{\mathrm{d}^{2}\boldsymbol{\ell}}{(2\pi)^{2}}\int\mathrm{d}^{3}\mathbf{v}g(\mathbf{v})v\sigma_{i}(\boldsymbol{\ell})\,,$
(120) where $\displaystyle v\sigma_{i}(\boldsymbol{\ell})$
$\displaystyle=\frac{\overline{\sigma}_{e}}{\mu^{2}}\int\frac{\mathrm{d}^{3}\mathbf{k}_{f}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}\mathbf{q}}{4\pi}\left|F_{\mathrm{DM}}(q)\right|^{2}\left|\widetilde{\Psi}_{i}\left(\boldsymbol{\ell},\mathbf{q}-\mathbf{k}_{f}\right)\right|^{2}$
$\displaystyle\times\delta\left(\frac{k_{f}^{2}}{2m_{e}}-E_{i}(\boldsymbol{\ell})+\Phi+\frac{q^{2}}{2m_{\chi}}-\mathbf{q}\cdot\mathbf{v}\right)\,.$
(121)

Here, we use the notation of Hochberg et al. This needs to be compared to our
Eqs. (33) and (34). Here, we use the replacement Eq. (46) for two-dimensional
targets, and further identify

$\displaystyle R_{\mathrm{free}}\rightarrow\frac{16\pi
m_{e}^{2}m_{\chi}^{2}}{\mu^{2}_{e\chi}}\overline{\sigma}_{e}|F_{\mathrm{DM}}(q)|^{2}$
(122)

to facilitate the comparison. We find agreement between our expressions with
one exception. Instead of the number of carbon atoms $N_{C}$, we find that the
electron ejection rate is proportional to the number of unit cells
$N_{\mathrm{cell}}$. Hence our expressions for the electron ejection rates
differ by a factor of 2.

In Fig. 12, we compare the energy spectrum for a DM particle of 100 MeV mass
using the TB approach presented by Hochberg et al. and compare it to the
improved version presented in this work. Note that the interaction model used
by Hochberg et al. corresponds to $\mathcal{O}_{1}$ interactions in our
general framework. We find that our predicted spectrum is about one order of
magnitude lower than the one presented in Hochberg  _et al._ (2017).

## References

  * Bertone and Hooper (2018) Gianfranco Bertone and Dan Hooper, “History of dark matter,” Rev. Mod. Phys. 90, 045002 (2018), arXiv:1605.04909 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) Nabila Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Clowe  _et al._ (2006) Douglas Clowe, Marusa Bradac, Anthony H. Gonzalez, Maxim Markevitch, Scott W. Randall, Christine Jones, and Dennis Zaritsky, “A direct empirical proof of the existence of dark matter,” Astrophys. J. 648, L109–L113 (2006), arXiv:astro-ph/0608407 [astro-ph] .
  * Persic  _et al._ (1996) Massimo Persic, Paolo Salucci, and Fulvio Stel, “The Universal rotation curve of spiral galaxies: 1. The Dark matter connection,” Mon. Not. Roy. Astron. Soc. 281, 27 (1996), arXiv:astro-ph/9506004 .
  * Drukier and Stodolsky (1984) Andrzej Drukier and Leo Stodolsky, “Principles and Applications of a Neutral Current Detector for Neutrino Physics and Astronomy,” Phys. Rev. D30, 2295 (1984).
  * Goodman and Witten (1985) Mark W. Goodman and Edward Witten, “Detectability of Certain Dark Matter Candidates,” Phys. Rev. D31, 3059 (1985).
  * Schumann (2019) Marc Schumann, “Direct Detection of WIMP Dark Matter: Concepts and Status,” J. Phys. G 46, 103003 (2019), arXiv:1903.03026 [astro-ph.CO] .
  * Battaglieri  _et al._ (2017) Marco Battaglieri _et al._ , “US Cosmic Visions: New Ideas in Dark Matter 2017; Community Report,” FERMILAB-CONF-17-282-AE-PPD-T (2017), arXiv:1707.04591 [hep-ph] .
  * Agnes  _et al._ (2018) P. Agnes _et al._ (DarkSide), “Constraints on Sub-GeV Dark-Matter–Electron Scattering from the DarkSide-50 Experiment,” Phys. Rev. Lett. 121, 111303 (2018), arXiv:1802.06998 [astro-ph.CO] .
  * Essig  _et al._ (2012a) Rouven Essig, Aaron Manalaysay, Jeremy Mardon, Peter Sorensen, and Tomer Volansky, “First Direct Detection Limits on sub-GeV Dark Matter from XENON10,” Phys. Rev. Lett. 109, 021301 (2012a), arXiv:1206.2644 [astro-ph.CO] .
  * Essig  _et al._ (2017) Rouven Essig, Tomer Volansky, and Tien-Tien Yu, “New Constraints and Prospects for sub-GeV Dark Matter Scattering off Electrons in Xenon,” Phys. Rev. D96, 043017 (2017), arXiv:1703.00910 [hep-ph] .
  * Aprile  _et al._ (2019) E. Aprile _et al._ (XENON), “Light Dark Matter Search with Ionization Signals in XENON1T,” Phys. Rev. Lett. 123, 251801 (2019), arXiv:1907.11485 [hep-ex] .
  * Essig  _et al._ (2012b) Rouven Essig, Jeremy Mardon, and Tomer Volansky, “Direct Detection of Sub-GeV Dark Matter,” Phys. Rev. D85, 076007 (2012b), arXiv:1108.5383 [hep-ph] .
  * Graham  _et al._ (2012) Peter W. Graham, David E. Kaplan, Surjeet Rajendran, and Matthew T. Walters, “Semiconductor Probes of Light Dark Matter,” Phys. Dark Univ. 1, 32–49 (2012), arXiv:1203.2531 [hep-ph] .
  * Lee  _et al._ (2015) Samuel K. Lee, Mariangela Lisanti, Siddharth Mishra-Sharma, and Benjamin R. Safdi, “Modulation Effects in Dark Matter-Electron Scattering Experiments,” Phys. Rev. D92, 083517 (2015), arXiv:1508.07361 [hep-ph] .
  * Essig  _et al._ (2016) Rouven Essig, Marivi Fernandez-Serra, Jeremy Mardon, Adrian Soto, Tomer Volansky, and Tien-Tien Yu, “Direct Detection of sub-GeV Dark Matter with Semiconductor Targets,” JHEP 05, 046 (2016), arXiv:1509.01598 [hep-ph] .
  * Crisler  _et al._ (2018) Michael Crisler, Rouven Essig, Juan Estrada, Guillermo Fernandez, Javier Tiffenberg, Miguel Sofo haro, Tomer Volansky, and Tien-Tien Yu (SENSEI), “SENSEI: First Direct-Detection Constraints on sub-GeV Dark Matter from a Surface Run,” Phys. Rev. Lett. 121, 061803 (2018), arXiv:1804.00088 [hep-ex] .
  * Agnese  _et al._ (2018) R. Agnese _et al._ (SuperCDMS), “First Dark Matter Constraints from a SuperCDMS Single-Charge Sensitive Detector,” Phys. Rev. Lett. 121, 051301 (2018), [erratum: Phys. Rev. Lett.122,no.6,069901(2019)], arXiv:1804.10697 [hep-ex] .
  * Abramoff  _et al._ (2019) Orr Abramoff _et al._ (SENSEI), “SENSEI: Direct-Detection Constraints on Sub-GeV Dark Matter from a Shallow Underground Run Using a Prototype Skipper-CCD,” Phys. Rev. Lett. 122, 161801 (2019), arXiv:1901.10478 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2019) A. Aguilar-Arevalo _et al._ (DAMIC), “Constraints on Light Dark Matter Particles Interacting with Electrons from DAMIC at SNOLAB,” Phys. Rev. Lett. 123, 181802 (2019), arXiv:1907.12628 [astro-ph.CO] .
  * Amaral  _et al._ (2020) D. W. Amaral _et al._ (SuperCDMS), “Constraints on low-mass, relic dark matter candidates from a surface-operated SuperCDMS single-charge sensitive detector,” Phys. Rev. D 102, 091101 (2020), arXiv:2005.14067 [hep-ex] .
  * Andersson  _et al._ (2020) Erik Andersson, Alex Bökmark, Riccardo Catena, Timon Emken, Henrik Klein Moberg, and Emil Åstrand, “Projected sensitivity to sub-GeV dark matter of next-generation semiconductor detectors,” JCAP 05, 036 (2020), arXiv:2001.08910 [hep-ph] .
  * Hochberg  _et al._ (2018) Yonit Hochberg, Yonatan Kahn, Mariangela Lisanti, Kathryn M. Zurek, Adolfo G. Grushin, Roni Ilan, Sinéad M. Griffin, Zhen-Fei Liu, Sophie F. Weber, and Jeffrey B. Neaton, “Detection of sub-MeV Dark Matter with Three-Dimensional Dirac Materials,” Phys. Rev. D97, 015004 (2018), arXiv:1708.08929 [hep-ph] .
  * Geilhufe  _et al._ (2020) R. Matthias Geilhufe, Felix Kahlhoefer, and Martin Wolfgang Winkler, “Dirac Materials for Sub-MeV Dark Matter Detection: New Targets and Improved Formalism,” Phys. Rev. D 101, 055005 (2020), arXiv:1910.02091 [hep-ph] .
  * Hochberg  _et al._ (2017) Yonit Hochberg, Yonatan Kahn, Mariangela Lisanti, Christopher G. Tully, and Kathryn M. Zurek, “Directional detection of dark matter with two-dimensional targets,” Phys. Lett. B772, 239–246 (2017), arXiv:1606.08849 [hep-ph] .
  * Geilhufe  _et al._ (2018) R. Matthias Geilhufe, Bart Olsthoorn, Alfredo Ferella, Timo Koski, Felix Kahlhoefer, Jan Conrad, and Alexander V. Balatsky, “Materials Informatics for Dark Matter Detection,” Phys. Status Solidi RRL12, 0293 (2018), arXiv:1806.06040 [cond-mat.mtrl-sci] .
  * Capparelli  _et al._ (2015) L. M. Capparelli, G. Cavoto, D. Mazzilli, and A. D. Polosa, “Directional Dark Matter Searches with Carbon Nanotubes,” Phys. Dark Univ. 9-10, 24–30 (2015), [Erratum: Phys.Dark Univ. 11, 79–80 (2016)], arXiv:1412.8213 [physics.ins-det] .
  * Cavoto  _et al._ (2016) G. Cavoto, E. N. M. Cirillo, F. Cocina, J. Ferretti, and A. D. Polosa, “WIMP detection and slow ion dynamics in carbon nanotube arrays,” Eur. Phys. J. C 76, 349 (2016), arXiv:1602.03216 [physics.ins-det] .
  * Cavoto  _et al._ (2018) G. Cavoto, F. Luchetta, and A. D. Polosa, “Sub-GeV Dark Matter Detection with Electron Recoils in Carbon Nanotubes,” Phys. Lett. B 776, 338–344 (2018), arXiv:1706.02487 [hep-ph] .
  * Cavoto  _et al._ (2020) G. Cavoto, M. G. Betti, C. Mariani, F. Pandolfi, A. D. Polosa, I. Rago, and A. Ruocco, “Carbon nanotubes as anisotropic target for dark matter,” J. Phys. Conf. Ser. 1468, 012232 (2020), arXiv:1911.01122 [physics.ins-det] .
  * Betts  _et al._ (2013) S. Betts _et al._ , “Development of a Relic Neutrino Detection Experiment at PTOLEMY: Princeton Tritium Observatory for Light, Early-Universe, Massive-Neutrino Yield,” in  _Community Summer Study 2013: Snowmass on the Mississippi_ (2013) arXiv:1307.4738 [astro-ph.IM] .
  * Baracchini  _et al._ (2018) E. Baracchini _et al._ (PTOLEMY), “PTOLEMY: A Proposal for Thermal Relic Detection of Massive Neutrinos and Directional Detection of MeV Dark Matter,” (2018), arXiv:1808.01892 [physics.ins-det] .
  * Apponi  _et al._ (2022) A. Apponi _et al._ (PTOLEMY), “Heisenberg’s uncertainty principle in the PTOLEMY project: A theory update,” Phys. Rev. D 106, 053002 (2022), arXiv:2203.11228 [hep-ph] .
  * Apponi  _et al._ (2021) Alice Apponi, Gianluca Cavoto, Carlo Mariani, Francesco Pandolfi, Ilaria Rago, Alessandro Ruocco, Christopher Tully, and Fang Zhao, “Carbon Nanostructures for Directional Light Dark Matter Detection,” PoS ICHEP2020, 648 (2021).
  * Catena  _et al._ (2020) Riccardo Catena, Timon Emken, Nicola A. Spaldin, and Walter Tarantino, “Atomic responses to general dark matter-electron interactions,” Phys. Rev. Res. 2, 033195 (2020), arXiv:1912.08204 [hep-ph] .
  * Catena  _et al._ (2022) Riccardo Catena, Daniel Cole, Timon Emken, Marek Matas, Nicola Spaldin, Walter Tarantino, and Einar Urdshals, “Dark matter - electron interactions in materials beyond the dark photon model,” (2022), arXiv:2210.07305 [hep-ph] .
  * Aprile  _et al._ (2022) E. Aprile _et al._ (XENON), “Emission of single and few electrons in XENON1T and limits on light dark matter,” Phys. Rev. D 106, 022001 (2022), arXiv:2112.12116 [hep-ex] .
  * Catena  _et al._ (2023) Riccardo Catena, Timon Emken, Marek Matas, Nicola Spaldin, and Einar Urdshals, “Direct searches for general dark matter-electron interactions with graphene detectors: Part II. Sensitivity studies,” (2023), arXiv:2303.xxxxx [hep-ph] .
  * Emken (2023) Timon Emken, “Darphene [Code, v0.1.0],” The code can be found under https://github.com/temken/Darphene. Version 0.1.0 is archived as DOI:10.5281/zenodo.7774375 (2023).
  * Urdshals and Matas (2021) Einar Urdshals and Marek Matas, “Qedark-eft,”  (2021).
  * Catena  _et al._ (2021) Riccardo Catena, Timon Emken, Marek Matas, Nicola A. Spaldin, and Einar Urdshals, “Crystal responses to general dark matter-electron interactions,” Phys. Rev. Res. 3, 033149 (2021), arXiv:2105.02233 [hep-ph] .
  * Baxter  _et al._ (2021) D. Baxter _et al._ , “Recommended conventions for reporting results from direct dark matter searches,” Eur. Phys. J. C 81, 907 (2021), arXiv:2105.00599 [hep-ex] .
  * Puschnig and Lüftner (2015) Peter Puschnig and Daniel Lüftner, “Simulation of angle-resolved photoemission spectra by approximating the final state by a plane wave: From graphene to polycyclic aromatic hydrocarbon molecules,” Journal of Electron Spectroscopy and Related Phenomena 200, 193–208 (2015), special Anniversary Issue: Volume 200.
  * Fayet (1980) Pierre Fayet, “Effects of the Spin 1 Partner of the Goldstino (Gravitino) on Neutral Current Phenomenology,” Phys. Lett. B 95, 285–289 (1980).
  * Fayet (1981) Pierre Fayet, “On the Search for a New Spin 1 Boson,” Nucl. Phys. B 187, 184–204 (1981).
  * Holdom (1986) Bob Holdom, “Two U(1)’s and Epsilon Charge Shifts,” Phys. Lett. 166B, 196–198 (1986).
  * Boehm and Fayet (2004) C. Boehm and Pierre Fayet, “Scalar dark matter candidates,” Nucl. Phys. B 683, 219–263 (2004), arXiv:hep-ph/0305261 .
  * Fayet (2004) Pierre Fayet, “Light spin 1/2 or spin 0 dark matter particles,” Phys. Rev. D 70, 023514 (2004), arXiv:hep-ph/0403226 .
  * Nijboer and Wette (1957) B. R. A. Nijboer and F. W. de Wette, “On the calculation of lattice sums,” Physica 23, 309–321 (1957).
  * Guennebaud  _et al._ (2010) Gaël Guennebaud, Benoît Jacob, _et al._ , “Eigen v3,” http://eigen.tuxfamily.org (2010).
  * Bunge  _et al._ (1993) C. F. Bunge, J. A. Barrientos, and A. V. Bunge, “Roothaan-Hartree-Fock Ground-State Atomic Wave Functions: Slater-Type Orbital Expansions and Expectation Values for Z = 2-54,” Atom. Data Nucl. Data Tabl. 53, 113–162 (1993).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864–B871 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * Martin (2004) Richard M. Martin,  _Electronic Structure: Basic Theory and Practical Methods_ (Cambridge University Press, 2004).
  * Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, “Ground state of the electron gas by a stochastic method,” Phys. Rev. Lett. 45, 566 (1980).
  * Trevisanutto  _et al._ (2008) Paolo E. Trevisanutto, Christine Giorgetti, Lucia Reining, Massimo Ladisa, and Valerio Olevano, “Ab Initio $GW$ Many-Body Effects in Graphene,” Phys. Rev. Lett. 101, 226405 (2008).
  * Feng  _et al._ (2019) Zhenbao Feng, Xiaoyan Zhang, Yoshiharu Sakurai, Zongliang Wang, Hefu Li, and Haiquan Hu, “Compton profile of few-layer graphene investigated by electron energy-loss spectroscopy,” Sci. Rep. 9, 1–5 (2019).
  * Giannozzi  _et al._ (2009) Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, and et al., “Quantum espresso: a modular and open-source software project for quantum simulations of materials,” Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, D Ceresoli, M Cococcioni, N Colonna, I Carnimeo, A Dal Corso, S de Gironcoli, P Delugas, R A DiStasio, A Ferretti, A Floris, G Fratesi, G Fugallo, R Gebauer, U Gerstmann, F Giustino, T Gorni, J Jia, M Kawamura, H-Y Ko, A Kokalj, E Küçükbenli, M Lazzeri, M Marsili, N Marzari, F Mauri, N L Nguyen, H-V Nguyen, A Otero de-la Roza, L Paulatto, S Poncé, D Rocca, R Sabatini, B Santra, M Schlipf, A P Seitsonen, A Smogunov, I Timrov, T Thonhauser, P Umari, N Vast, X Wu, and S Baroni, “Advanced capabilities for materials modelling with quantum ESPRESSO,” Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ (2020) Paolo Giannozzi, Oscar Baseggio, Pietro Bonfà, Davide Brunato, Roberto Car, Ivan Carnimeo, Carlo Cavazzoni, Stefano de Gironcoli, Pietro Delugas, Fabrizio Ferrari Ruffino, Andrea Ferretti, Nicola Marzari, Iurii Timrov, Andrea Urru, and Stefano Baroni, “Quantum espresso toward the exascale,” The Journal of Chemical Physics 152, 154105 (2020), https://doi.org/10.1063/5.0005082 .
  * Perdew  _et al._ (1996) John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
  * Griffin  _et al._ (2021) Sinéad M. Griffin, Katherine Inzani, Tanner Trickle, Zhengkang Zhang, and Kathryn M. Zurek, “Extended calculation of dark matter-electron scattering in crystal targets,” Phys. Rev. D 104, 095015 (2021), arXiv:2105.05253 [hep-ph] .
  * Johansson (2017) F. Johansson, “Arb: efficient arbitrary-precision midpoint-radius interval arithmetic,” IEEE Transactions on Computers 66, 1281–1292 (2017).
  * boo (2022) “The Boost C++ library v1.81.0,”  (2022).
  * Emken (2020) Timon Emken, “libphysica - A C++ library for application in scientific codes [Code],” The code can be found under https://github.com/temken/libphysica. (2020).
  * Emken (2021a) Timon Emken, “obscura: A modular C++ tool and library for the direct detection of (sub-GeV) dark matter via nuclear and electron recoils,” J. Open Source Softw. 6, 3725 (2021a), arXiv:2112.01489 [hep-ph] .
  * Emken (2021b) Timon Emken, “obscura - A C++ library for dark matter detection computations [Code],” The code can be found under https://github.com/temken/obscura. (2021b).
  * Rohatgi  _et al._ (2012) Ankit Rohatgi, Steffen Rehberg, and Zlatan Stanojevic, “WebPlotDigitizer v4.2,”  (2012).
  * Wolfram Research, Inc. (2021) Wolfram Research, Inc., “Mathematica v13.0,”  (2021), champaign, IL, 2021.
  * Saito  _et al._ (1998) R. Saito, G. Dresselhaus, and M.S. Dresselhaus,  _Physical Properties of Carbon Nanotubes_ (Imperial College Press, 1998).