# Stochastic behavior of an $n$-node blockchain under cyber attacks from
multiple hackers with random re-setting times

Xiufeng Xu∗ Liang Hong111Department of Mathematical Sciences, The University
of Texas at Dallas, 800 West Campbell Road, Richardson, TX 75080, USA.

###### Abstract

This paper investigates the stochastic behavior of an $n$-node blockchain
which is continuously monitored and faces non-stop cyber attacks from multiple
hackers. The blockchain will start being re-set once hacking is detected,
forfeiting previous efforts of all hackers. It is assumed the re-setting
process takes a random amount of time. Multiple independent hackers will keep
attempting to hack into the blockchain until one of them succeeds. For
arbitrary distributions of the hacking times, detecting times, and re-setting
times, we derive the instantaneous functional probability, the limiting
functional probability, and the mean functional time of the blockchain.
Moreover, we establish that these quantities are increasing functions of the
number of nodes, formalizing the intuition that the more nodes a blockchain
has the more secure it is.

_Keywords and phrases:_ blockchain; instantaneous functional probability;
limiting functional probability; mean functional times; multiple hackers;
random re-setting time.

##  1 Introduction

The Internet has become an inseparable part of our life. Nowadays most
organizations have an online platform for conducting business with their
customers. The convenience and efficiency provided by this online business
mode come with a serious risk: an organization has become more vulnerable to
cyber attacks that could cost them a fortune. The devastating effect of a
cyber attack is as fresh as in recent memory, as evidenced by the Colonial
Pipeline ransomware attack in 2021, the service attack on many major US
airports in 2022, the cyber attack against DP World Australia Port in 2023,
and the Lurie Children’s Hospital ransomware attack in 2024. Each organization
must protect its information and maintain its online business. One solution is
to back up all its data continuously. Unfortunately, this idea is not
practical for most organizations due to the astronomical amount of information
they have to store. A viable option is to store sensitive data on a
blockchain. A blockchain consists of several identical computers in a network;
each is called a _node_ and stores identical data as other nodes do. All nodes
continuously verify the data according to the majority rule—any piece of data
is considered valid only if it is consistent with its counterpart on at least
half of the nodes. While a hacker can steal the data stored on a node of a
blockchain, the stolen data is often skillfully encrypted and hence is almost
useless to the hacker due to the public key encryption technology (e.g.,
Diffie and Hellman 1976, Rivest et al. 1978, Goldwasser and Micali 1982,
Goldwasser et al. 1988. Therefore, most hackers would either wreak havoc on
the organization by altering the data stored on the blockchain or making a
ransom demand by locking the data on the blockchain. The former type of cyber
attack, called _cyber destructive attack_ , requires a hacker to hack into at
least half of the nodes; the latter type, called _cyber ransom attack_ ,
requires the hacker to hack into all nodes. In view of this, blockchains do
not provide cyber-attack-proof protection, but they reduce an organization’s
likelihood of incurring cyber losses. For a comprehensive review of elements
of blockchains and their applications, see Tama et al. (2017), Zheng et al.
(2018), Casino et al. (2019), Kumar et al. (2020), and references therein.

The security issues of a blockchain have received significant attention from
the scientific community; see, for example, Corbet et al. (2020), Hussain et
al. (2022), Khanum and Mustafa (2022), Khan and Salah (2017), Li et al.
(2071), Meng et al. (2018), and Tsuchiya et al. (2021). However, there are
only a few papers that study the operations-research-theoretic aspect of a
blockchain; see, for instance, Choi et al. (2020), Melo et al. (2021), Xu and
Hong (2014), and Zhang et al. (2020). In particular, Xu and Hong (2024)
investigate the stochastic behavior of an $n$-node blockchain under cyber
attacks from a single hacker. This article is a further study along this line.
For arbitrary distributions of the detecting times, hacking times, and re-
setting times, it derives several key quantities of an $n$-node blockchain
under cyber attacks from multiple hackers when it takes a random amount of
time to re-set the blockchain.

The remainder of the article is organized as follows. Section 2 describes the
problem under consideration and establishes some notational conventions.
Section 3 investigates the stochastic behavior of an $n$-node blockchain in
our setting and gives formulas for the instantaneous functional probability,
limiting functional probability, and mean functional time. It also establishes
that these quantities are increasing functions of the number of nodes. This
substantiates the intuition that the more nodes a blockchain has the harder it
is to hack into it. Section 4 provides several numerical examples. Many
technical details of these examples are delegated to the Appendix. Section 5
concludes the article with some remarks. Python code for all examples is
available at
https://github.com/xuxiufeng/Blockchain_Simulation_Multiple_Hackers_with_Resetting_Times.

##  2 Model setup and notation

We consider a blockchain of $n$ identical nodes where $n\geq 2$. The
blockchain is under continuous monitoring (i.e., 24-hour monitoring) of the IT
department of an organization for potential cyber attacks. To perform a cyber
destructive attack, a hacker must hack into at least $m=\lfloor n/2\rfloor+1$
nodes. We will focus on cyber destructive attacks only since the case of a
cyber ransom attack is similar except $m$ is replaced by $n$. We assume there
are $k$ independent hackers where $k\geq 1$, and each of them continuously
attacks the blockchain until successfully hacking into the blockchain. Once a
hacker hacks into the blockchain, he will immediately start to alter the data
on it. When that happens, we say the blockchain is _dysfunctional_. Hence, the
blockchain will be functional until at least one of the $k$ hackers
successfully hacks into $m$ nodes. The IT department will start re-setting the
blockchain once they detect that the blockchain is under attack. Resetting
takes a random amount of time to finish. After that, the previous effort of
each hacker is forfeited.

Without lost of generality, we label the $n$ identical nodes as node $1$, node
$2$, $\dots$, and node $n$. For $i=1,2,\dots,n$ and $j=1,2,\dots,k$, let
$X_{i}^{j}$ be the time the $j$-th hacker needs to hack into node $i$. For
$j=1,\ldots,k$, we assume that $X_{1}^{j},\ldots,X_{n}^{j}$ are independent
and identically distributed (iid) according to a (cumulative) distribution
function of $F_{X^{j}_{1}}$ with a probability density function
$f_{X^{j}_{1}}(x)$. Let $Y$ be the time the IT department takes to detect any
cyber attack. We assume that $Y$ has a cumulative function $F_{Y}$ with a
probability density function $f_{Y}(y)$. Let $W$ be the time it takes for the
IT department to re-set the blockchain. We assume that $W$ has a distribution
function $F_{W}$ with a probability density function $f_{W}(w)$. Note that
$X_{i}^{j},Y,$ and $W$ are mutually independent for $j=1,\dots,k$. We call
each $X_{i}$ a $hacking$ $time$, $Y$ a $detecting$ $time$, and $W$ a re-set
time. For the $jth$ hacker, if $\sum_{i=1}^{m}X_{i}^{j}<Y$, then the $j$-th
hacker can change the data on the blockchain. If $\sum_{i=1}^{m}X_{i}^{j}\geq
Y$, then the $jth$ hacker is detected by the IT department, and the blockchain
will start being re-set immediately. It takes $W$ amount of time to complete
the re-setting process. Once the re-setting is finished, all the $k$ hackers
will conduct a new round of cyber attacks from the node $1$.

The stochastic behavior of the blockchain can be described as follows. At time
$t=0$, the blockchain starts to operate and all hackers start to work on
hacking into the blockchain. As all nodes are identical, we may assume that
each hacker attacks the $n$ nodes in ascending order. That is, each of them
will first attack node 1. If $Y\leq\min_{1\leq j\leq k}X^{j}_{1}$, the IT
department will detect the cyber attack before any hacker can hack into node
1. In that case, the IT department will immediately start re-setting the
blockchain, rendering all previous efforts of each hacker in vain. The re-
setting process takes $W$ amount of time after which the blockchain operates
as if afresh and all hackers start to attack node 1 again. On the other hand,
if $Y>\min_{1\leq j\leq k}X^{j}_{1}$, that is, at least one of the $k$ hackers
jas successfully bypassed the firewall of node 1, he will immediately move
onto attacking node 2. This process continues until
$Z_{m}=\min_{j=1,\dots,k}\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\bigg{)}<Y$, i.e., one
of the $k$ hackers succeeds in hacking into $m$ nodes.

For given $m$ and $k$, we define the instantaneous functional probability of
the blockchain to be

$P_{mk}(t)=P\\{\text{the blockchain has not been hacked at time $t$}\\}.$ (1)

As $t\rightarrow\infty$, $P_{mk}(t)$ converges to the _limiting functional
probability_ of the blockchain which is defined as

$P_{mk}(\infty)=P\\{\text{the blockchain will never be hacked}\\}.$ (2)

##  3 Stochastic behavior of the proposed blockchain model

###  3.1 Limiting functional probability

To derive the limiting functional probability, we define a $cycle$ to be the
period from the moment the blockchain starts afresh or has just been re-set to
the next moment it starts being re-set. Let $N_{1}$ be the total number of
cycles before one of the $k$ hackers successfully hacks into the blockchain,
i.e., $Z_{m}\leq Y$. Then $N_{1}$ is a geometric random variable222Here we
interpret a geometric random variable with parameter $p$ to be the number of
failures until the first success among a sequence of independent Bernoulli
trails with success probability $p$. with parameter $1-p_{mk}$, where

$\displaystyle p_{mk}$ $\displaystyle=$ $\displaystyle
1-P\bigg{\\{}\sum_{i=1}^{m}X_{i}^{1}>Y,\sum_{i=1}^{m}X_{i}^{2}>Y,\dots,\sum_{i=1}^{m}X_{i}^{k}>Y\bigg{\\}}$
(3) $\displaystyle=$ $\displaystyle
1-P\bigg{\\{}\min_{j=1,\dots,k}\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\bigg{)}>Y\bigg{\\}}$
$\displaystyle=$ $\displaystyle 1-P(Z_{m}>Y)=P(Z_{m}\leq Y)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}F_{Z_{m}}(s)dF_{Y}(s),$

where $F_{Z_{m}}$ is the distribution function of $Z_{m}$. The continuous-time
process seems to be mathematically intractable. However, if we focus only on
those moments the blockchain starts afresh, has just been re-set, or has just
been hacked, then we can identify an embedded three-state discrete-time Markov
chain; see the following diagram.

$1$$2$$3$$1-p_{mk}$$p_{mk}$$1$$1$

Figure 1: Three-state semi-Markov process

In State $1$, the blockchain is functional, i.e., none of the $k$ hackers has
hacked into $m$ nodes. In State $3$, the blockchain is dysfunctional (i.e., it
has been hacked). In State $2$, the blockchain is being re-set. State $2$ is
an absorbing state. The transition matrix of this Markov chain is given by
(4).

$\begin{bmatrix}0&p_{mk}&1-p_{mk}\\\ 0&1&0\\\ 1&0&0\end{bmatrix}$ (4)

Condition on the outcome of the first cycle to obtain

$P_{m}(\infty)=(1-p_{mk})\times 1\times P_{mk}(\infty)+p_{mk}\times 0.$

It follows that $P_{m}(\infty)=0$. Intuitively, this says that if $p_{mk}>0$
then the blockchain will eventually be hacked.

###  3.2 Mean functional time

To derive the mean functional time of blockchain, we identify a renewal
process $\\{N_{2}(t)\\}_{t\geq 0}$ as follows:

$N_{2}(t)=\max\left\\{k:\sum_{i=1}^{k}(Y_{i}+W_{i})\leq t\right\\},$

where $Y_{i}\overset{d}{=}Y\mid Y<Z_{m}$, i.e., $Y_{i}$ and the conditional
random variable $Y\mid Y<Z_{m}$ have the same distribution. We put

$T_{m}=\sum_{i=1}^{N_{1}}(Y_{i}+W_{i})+\left(Z_{m}\middle|\ Z_{m}<Y\right).$
(5)

Then $T_{m}$ is the functional time of the blockchain. Since $N_{1}$ is
independent of $Y_{n+1},Y_{n+2},\ldots$, $N_{1}$ is a stopping time for the
sequence $\\{Y_{i}\\}_{i\geq 1}$. It follows from Wald’s Identity that

$\displaystyle E[T_{m}]$ $\displaystyle=$ $\displaystyle
E[N_{1}](E[Y_{1}]+E[W_{1}])+E\left[Z_{m}\middle|\ Z_{m}<Y\right]$ (6)
$\displaystyle=$ $\displaystyle E[N_{1}]\bigg{\\{}E\left[Y\middle|\
Y<Z_{m}\right]+E[W_{1}]\bigg{\\}}+E\left[Z_{m}\middle|\ Z_{m}<Y\right].$

###  3.3 Instantaneous functional probability

To obtain the instantaneous functional probability, we will establish a
renewal-type equation. To this end, we define $S_{t}$ to be the last moment
the blockchain is re-set before time $t$. That is,

$S_{t}=\sum_{i=1}^{N_{2}(t)}Y_{i}+W_{i}.$

Let $F_{S_{t}}$ be the distribution function of $S_{t}$. Also, we use
$G(s)=E[{N_{2}(s)}]$ to denote the renewal function of $\\{N_{2}(t)\\}_{t\geq
0}$. Then

$G(s)=\sum_{k=1}^{\infty}P\left(\sum_{i=1}^{k}Y_{i}+W_{i}\leq s\right).$

Moreover, we have

$\displaystyle P(S_{t}=0)$ $\displaystyle=$ $\displaystyle P(Y_{1}+W_{1}>t)$
$\displaystyle dF_{S_{t}}(s)$ $\displaystyle=$ $\displaystyle
P(Y_{1}+W_{1}>t-s)dG(s),\quad\text{$0<s<\infty$}.$

To obtain $P_{m}(t)$, we need to derive the probability that the blockchain is
functional at $t$ and the probability that the blockchain is re-setting at
$t$. For the former, we have

$\displaystyle P\\{\text{the blockchain is functional at t}\\}$
$\displaystyle=$ $\displaystyle P\\{\text{the blockchain is functional at
$t$}\mid S_{t}=0\\}P\\{S_{t}=0\\}$ $\displaystyle+\int_{0}^{t}P\\{\text{the
blockchain is functional at $t$}\mid S_{t}=s\\}dF_{S_{t}}(s)$ $\displaystyle=$
$\displaystyle P\left\\{t<Z_{m}\wedge
Y\right\\}+\int_{0}^{t}P\left\\{t-s<Z_{m}\wedge Y\right\\}dG(s),$

where $dF_{S_{t}}(s)=P\\{t-s<Y_{1}+W_{1}\\}dG(s)$ and
$G(s)=\sum_{k=1}^{\infty}P\left(\sum_{i=1}^{k}Y_{i}+W_{i}\leq s\right)$. For
the latter, we have

$\displaystyle P\\{\text{the blockchain is re-setting at t}\\}$
$\displaystyle=$ $\displaystyle P\\{\text{the blockchain is re-setting at
$t$}\mid S_{t}=0\\}P\\{S_{t}=0\\}$ $\displaystyle+\int_{0}^{t}P\\{\text{the
blockchain is re-setting at $t$}\mid S_{t}=s\\}dF_{S_{t}}(s)$ $\displaystyle=$
$\displaystyle P\left\\{Y\leq t<Y+W_{1},Y\leq Z_{m}\right\\}$
$\displaystyle+\int_{0}^{t}P\left\\{Y\leq t-s<Y+W_{1},Y\leq
Z_{m}\right\\}dG(s).$

Therefore,

$\displaystyle P_{mk}(t)$ $\displaystyle=$ $\displaystyle P\\{\text{the
blockchain has not been hacked at $t$}\mid S_{t}=0\\}P\\{S_{t}=0\\}$
$\displaystyle+\int_{0}^{t}P\\{\text{the blockchain has not been hacked at
$t$}\mid S_{t}=s\\}dF_{S_{t}}(s)$ $\displaystyle=$ $\displaystyle
P\\{\text{the blockchain is functional at $t$}\mid S_{t}=0\\}P\\{S_{t}=0\\}$
(7) $\displaystyle+P\\{\text{the blockchain is re-setting at $t$}\mid
S_{t}=0\\}P\\{S_{t}=0\\}$ $\displaystyle+\int_{0}^{t}P\\{\text{the blockchain
is functional at $t$}\mid S_{t}=s\\}dF_{S_{t}}(s)$
$\displaystyle+\int_{0}^{t}P\\{\text{the blockchain is re-setting at $t$}\mid
S_{t}=s\\}dF_{S_{t}}(s)$ $\displaystyle=$ $\displaystyle
P\left\\{t<Z_{m}\wedge Y\right\\}+P\left\\{Y\leq t<Y+W_{1},Y\leq
Z_{m}\right\\}$ $\displaystyle+\int_{0}^{t}P\left\\{t-s<Z_{m}\wedge
Y\right\\}dG(s)$ $\displaystyle+\int_{0}^{t}P\left\\{Y\leq t-s<Y+W_{1},Y\leq
Z_{m}\right\\}dG(s).$

###  3.4 Monotonicity of $P_{m}(t)$ and $E[T_{m}]$

According to our intuition, the more nodes a blockchain has, the harder it is
to hack into it. Formally, this amounts to saying that $P_{m}(t)$,
$P_{m}(\infty)$, and $E[T_{m}]$ all are increasing functions of $m$. This is
indeed the case, as the next theorem shows.

######  Theorem 3.1.

$P_{mk}(t)$, $P_{mk}(\infty)$, and $E[T_{m}]$ are all increasing functions of
$m$, where $m=\lfloor n/2\rfloor+1$.

###### Proof.

First, the definition of $Z_{m}$ implies

$P\\{t<Y\wedge Z_{m+1}\\}>P\\{t<Y\wedge Z_{m}\\},$

and

$P\\{t-s<Y\wedge Z_{m+1}\\}>P\\{t-s<Y\wedge Z_{m}\\}.$

Similarly, we have

$P\left\\{Y\leq t<Y+W_{1},Y\leq\ Z_{m+1}\right\\}>P\left\\{Y\leq
t<Y+W_{1},Y\leq Z_{m}\right\\}$

and

$P\left\\{Y\leq t-s<Y+W_{1},Y\leq Z_{m+1}\right\\}>P\left\\{Y\leq
t-s<Y+W_{1},Y\leq Z_{m}\right\\}.$

Then it follows from (3.3) that

$P_{m+1,k}(t)>P_{mk}(t),\quad\text{for all positive integer $m$ and $t$},$ (8)

which also implies

$P_{m+1,k}(\infty)=\lim_{t\rightarrow\infty}P_{m+1,k}(t)>\lim_{t\rightarrow\infty}P_{mk}(t)=P_{mk}(\infty).$
(9)

To show $E[T_{m}]$ is increasing in $m$, we write it as

$\displaystyle E[T_{m}]$ $\displaystyle=$
$\displaystyle\frac{1-p_{mk}}{p_{mk}}\times\left\\{\frac{E[Y1_{\\{Y\leq
Z_{m}\\}}]}{1-p_{mk}}+E[W_{1}]\right\\}+\frac{E[Z_{m}1_{\\{Y>Z_{m}\\}}]}{p_{mk}}$
$\displaystyle=$ $\displaystyle\frac{E[Y1_{\\{Y\leq
Z_{m}\\}}]}{p_{mk}}+\frac{E[Z_{m}1_{\\{Y>Z_{m}\\}}]}{p_{mk}}+\frac{1-p_{mk}}{p_{mk}}\times
E[W_{1}].$

Thus,

$\displaystyle E[T_{m+1}]-E[T_{m}]$ $\displaystyle=$
$\displaystyle(p_{mk}p_{m+1,k})^{-1}\bigg{\\{}p_{mk}E[Y1_{\\{Y\leq
Z_{m+1}\\}}]-p_{m+1,k}E[Y1_{\\{Y\leq Z_{m}\\}}]$
$\displaystyle+p_{mk}E\bigg{[}Z_{m+1}1_{\\{Y>Z_{m+1}\\}}\bigg{]}-p_{m+1,k}E\bigg{[}Z_{m}1_{\\{Y>Z_{m}\\}}\bigg{]}$
$\displaystyle+(p_{mk}-p_{m+1,k})E[W_{1}]\bigg{\\}}.$

The first four terms in the brackets can be written as

$\displaystyle
p_{mk}E[Y1_{\\{Y<Z_{m}\\}}]-p_{m+1,k}E[Y1_{\\{Y<Z_{m}\\}}]+p_{mk}E[Y1_{\\{Z_{m}<Y<Z_{m+1}\\}}]$
$\displaystyle+p_{mk}E[(Z_{m+1}-Z_{m})1_{\\{Y>Z_{m+1}\\}}]+p_{mk}E[Z_{m}1_{\\{Y>Z_{m+1}\\}}]$
$\displaystyle-
p_{m+1,k}E[Z_{m}1_{\\{Y>Z_{m+1}\\}}]-p_{m+1,k}E[Z_{m}1_{\\{Z_{m+1}>Y>Z_{m}\\}}]$
$\displaystyle=$
$\displaystyle(p_{mk}-p_{m+1,k})E[Y1_{\\{Y<Z_{m}\\}}]+(p_{mk}-p_{m+1,k})E[Z_{m}1_{\\{Y>Z_{m+1}\\}}]$
$\displaystyle+p_{mk}E[Y1_{\\{Z_{m+1}>Y>Z_{m}\\}}]\
-p_{m+1,k}E[Z_{m}1_{\\{Z_{m+1}>Y>Z_{m}\\}}]$
$\displaystyle+p_{mk}E[(Z_{m+1}-Z_{m})1_{\\{Y>Z_{m+1}\\}}]$
$\displaystyle\geq$
$\displaystyle(p_{mk}-p_{m+1,k})E[Y1_{\\{Y<Z_{m}\\}}]+(p_{mk}-p_{m+1,k})E[Z_{m}1_{\\{Y>Z_{m+1}\\}}]$
$\displaystyle+p_{mk}E[Z_{m}1_{\\{Z_{m+1}>Y>Z_{m}\\}}]\
-p_{m+1,k}E[Z_{m}1_{\\{Z_{m+1}>Y>Z_{m}\\}}]$
$\displaystyle+p_{mk}E[(Z_{m+1}-Z_{m})1_{\\{Y>Z_{m+1}\\}}]$ $\displaystyle=$
$\displaystyle(p_{mk}-p_{m+1,k})E[Y1_{\\{Y<Z_{m}\\}}]+(p_{mk}-p_{m+1,k})E[Z_{m}1_{\\{Y>Z_{m+1}\\}}]$
$\displaystyle+(p_{mk}-p_{m+1,k})E[Z_{m}1_{\\{Z_{m+1}>Y>Z_{m}\\}}]+p_{mk}E[(Z_{m+1}-Z_{m})1_{\\{Y>Z_{m+1}\\}}].$

(3) implies that $p_{mk}$ is decreasing in $m$. It is also clear that
$Z_{m+1}-Z_{m}>0$ almost surely. Hence, $E[T_{m+1}]-E[T_{m}]>0$. ∎

######  Theorem 3.2.

$\lim_{m\rightarrow\infty}P_{m}(t)=1$ and
$\lim_{m\rightarrow\infty}E[T_{m}]=\infty$, where $m=\lfloor n/2\rfloor+1$.

###### Proof.

It is clear from (3) that $\lim_{m\rightarrow\infty}p_{mk}=0$. This is
equivalent to saying that the embedded discrete-time Markov chain will never
be State $2$ almost surely in the limit. Therefore,
$\lim_{m\rightarrow\infty}P_{m}(t)=1$.

Since $\lim_{m\rightarrow\infty}p_{mk}=0$,
$E[N_{1}]=1/p_{mk}-1\rightarrow\infty$ as $m\rightarrow\infty$. It follows
from (6) that

$\displaystyle\lim_{m\rightarrow\infty}E[T_{m}]$ $\displaystyle\geq$
$\displaystyle\lim_{m\rightarrow\infty}E[N_{1}]E\left[Y\middle|\
Y<Z_{m}\right]$ $\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}E[N_{1}]\lim_{m\rightarrow\infty}\left[\frac{E[Y1_{\\{Y<Z_{m}\\}}]}{P\left\\{Y<Z_{m}\right\\}}\right]$
$\displaystyle=$ $\displaystyle\lim_{m\rightarrow\infty}E[N_{1}]E[Y]$
$\displaystyle=$ $\displaystyle\infty.$

∎

Intuitively, the more hackers out there, the more likely the blockchain will
be hacker. The next theorem confirms this intuition.

######  Theorem 3.3.

$P_{mk}(t)$, $P_{mk}(\infty)$, and $E[T_{mk}]$ are all increasing functions of
$k$, where $k$ is a positive integer.

###### Proof.

For $k_{1}>k_{2}$, we have

$\displaystyle
Z_{mk_{2}}=\min_{j=1,\dots,k_{2}}\sum_{i=1}^{m}X_{i}^{j}\geq\min_{j=1,\dots,k_{1}}\sum_{i=1}^{m}X_{i}^{j}=Z_{mk_{1}},$

and hence

$\displaystyle Z_{mk_{2}}\wedge Y\geq Z_{mk_{1}}\wedge Y,$

Then (3.3) implies

$\displaystyle P_{mk_{1}}(t)-P_{mk_{2}}(t)$ $\displaystyle=$ $\displaystyle
P\\{t<Z_{mk_{1}}\wedge Y\\}-P\\{t<Z_{mk_{2}}\wedge Y\\}$
$\displaystyle+\int_{0}^{t}P\\{Y\wedge Z_{mk_{1}}>t-s\\}-\\{Y\wedge
Z_{mk_{2}}>t-s\\}dG(s)$ $\displaystyle+\int_{0}^{t}P\left\\{Y\leq
t-s<Y+W_{1},Y\leq Z_{mk_{1}}\right\\}$ $\displaystyle-P\left\\{Y\leq
t-s<Y+W_{1},Y\leq Z_{mk_{2}}\right\\}dG(s)$ $\displaystyle\leq$ $\displaystyle
0.$

This shows that $P_{mk}(t)$ is an increasing function of $k$.

Next, (6) implies

$\displaystyle E[T_{mk}]$ $\displaystyle=$ $\displaystyle\frac{E[Y1_{\\{Y\leq
Z_{mk}\\}}]}{p_{mk}}+\frac{E[Z_{mk}1_{\\{Y>Z_{mk}\\}}]}{p_{mk}}+\frac{1-p_{mk}}{p_{mk}}\times
E[W_{1}].$

Thus,

$\displaystyle E[T_{mk_{1}}]-E[T_{mk_{2}}]$ $\displaystyle=$
$\displaystyle(p_{mk_{2}}p_{mk_{1}})^{-1}\bigg{\\{}p_{mk_{2}}E[Y1_{\\{Y<Z_{mk_{1}}\\}}]-p_{mk_{1}}E[Y1_{\\{Y<Z_{mk_{2}}\\}}]$
$\displaystyle+p_{mk_{2}}E[Z_{mk_{1}}1_{\\{Y>Z_{mk_{1}}\\}}]-p_{mk_{1}}E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{2}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[W_{1}]\bigg{\\}}.$

The first four terms in the brackets can be written as

$\displaystyle
p_{mk_{2}}E[Y1_{\\{Y<Z_{mk_{1}}\\}}]-p_{mk_{1}}E[Y1_{\\{Y<Z_{mk_{1}}\\}}]-p_{mk_{1}}E[Y1_{\\{Z_{mk_{1}}<Y<Z_{mk_{2}}\\}}]$
$\displaystyle+p_{mk_{2}}E[(Z_{mk_{1}}-Z_{mk_{2}})1_{\\{Y>Z_{mk_{1}}\\}}]+p_{mk_{2}}E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle-
p_{mk_{1}}E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]+p_{mk_{1}}E[Z_{mk_{2}}1_{\\{Z_{mk_{2}}>Y>Z_{mk_{1}}\\}}]$
$\displaystyle=$
$\displaystyle(p_{mk_{2}}-p_{mk_{1}})E[Y1_{\\{Y<Z_{mk_{1}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle+p_{mk_{1}}E[Z_{mk_{2}}1_{\\{Z_{mk_{2}}>Y>Z_{mk_{1}}\\}}]\
-p_{mk_{1}}E[Y1_{\\{Z_{mk_{1}}<Y<Z_{mk_{2}}\\}}]$
$\displaystyle+p_{mk_{2}}E[(Z_{mk_{1}}-Z_{mk_{2}})1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle\leq$
$\displaystyle(p_{mk_{2}}-p_{mk_{1}})E[Y1_{\\{Y<Z_{mk_{1}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle+p_{mk_{1}}E[Z_{mk_{2}}1_{\\{Z_{mk_{2}}>Y>Z_{mk_{1}}\\}}]\
-p_{mk_{1}}E[Z_{mk_{1}}1_{\\{Z_{mk_{1}}<Y<Z_{mk_{2}}\\}}]$
$\displaystyle+p_{mk_{2}}E[(Z_{mk_{1}}-Z_{mk_{2}})1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle=$
$\displaystyle(p_{mk_{2}}-p_{mk_{1}})E[Y1_{\\{Y<Z_{mk_{1}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle+p_{mk_{1}}E[(Z_{mk_{2}}-Z_{mk_{1}})1_{\\{Z_{mk_{2}}>Y>Z_{mk_{1}}\\}}]$
$\displaystyle+p_{mk_{2}}E[(Z_{mk_{1}}-Z_{mk_{2}})1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle\leq$
$\displaystyle(p_{mk_{2}}-p_{mk_{1}})E[Y1_{\\{Y<Z_{mk_{1}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle+p_{mk_{1}}E[(Z_{mk_{2}}-Z_{mk_{1}})1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle+p_{mk_{2}}E[(Z_{mk_{1}}-Z_{mk_{2}})1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle=$
$\displaystyle(p_{mk_{2}}-p_{mk_{1}})E[Y1_{\\{Y<Z_{mk_{1}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[Z_{mk_{2}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle+(p_{mk_{2}}-p_{mk_{1}})E[(Z_{mk_{1}}-Z_{mk_{2}})1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle=$
$\displaystyle(p_{mk_{2}}-p_{mk_{1}})E[Y1_{\\{Y<Z_{mk_{1}}\\}}]+(p_{mk_{2}}-p_{mk_{1}})E[Z_{mk_{1}}1_{\\{Y>Z_{mk_{1}}\\}}]$
$\displaystyle\leq$ $\displaystyle 0,$

Besides, $(p_{mk_{2}}-p_{mk_{1}})E[W_{1}]\leq 0$. Therefore, $E[T_{mk}]$ is an
increasing function of $k$. ∎

###  3.5 Cost-benefit analysis

Suppose $C_{1}(m)$ is the cost per unit of time for re-setting the blockchain,
$C_{2}(m)$ is the cost per unit of time for running the blockchain, and $R(m)$
is the revenue earned per unit of time by running the blockchain. Then the
total net revenue (i.e., the total profit), denoted as $TNR$, can be written
as

$TNR=\left\\{\sum_{i=1}^{N_{1}}Y_{i}+\left(Z_{m}\middle|\
Z_{m}<Y\right)\right\\}(R(m)-C_{2}(m))-C_{1}(m)\sum_{i=1}^{N_{1}}W_{i}.$ (10)

It follows from Wald’s Identity that the total expected net revenue is given
by the following equation:

$E_{m}[TNR]=\left\\{E[N_{1}]E\left[Y\middle|\
Y<Z_{m}\right]+E\left[Z_{m}\middle|\
Z_{m}<Y\right]\right\\}(R(m)-C_{2}(m))-E[N_{1}]E[W_{1}]C_{1}(m).$ (11)

Therefore, the expected net revenue per unit of time, denoted as $E[NR]$, can
be written as

$E_{m}[NR]=\frac{E[TNR]}{E[T_{m}]}.$ (12)

##  4 Examples

In this section, we provide three numerical examples. Derivations of all
formulas are provided in the Appendix.

######  Example 1 (Exponential hacking, exponential re-set, and exponential
detecting times).

Suppose $X_{i}^{j}\overset{iid}{\sim}{\sf Exp}(\lambda_{j})$ for $i=1,\dots,n$
and $j=1,\dots,k$, $Y\sim{\sf Exp}(\delta)$, and $W\sim{\sf Exp}(\eta)$ such
that $X_{1}^{j}\dots X_{n}^{j}$, $Y$, and $W$ are mutually independent, where
${\sf Exp}(\lambda_{j})$, ${\sf Exp}(\delta)$, and ${\sf Exp}(\eta)$ denote
the exponential distributions with rate $\lambda_{j}>0$, $\delta>0$, and
$\eta>0$, respectively. We have $X\sim{\sf Exp}(\lambda_{j})$ if and only if
its probability density function is given by

$f(x)=\lambda_{j}e^{-\lambda_{j}x},\quad x>0.$

Clearly, $\sum_{i=1}^{m}X_{i}^{j}\sim{\sf Gamma}(m,\lambda_{j})$, where ${\sf
Gamma}(m,\lambda_{j})$ denotes the gamma distribution with shape parameter $m$
and rate parameter $\lambda_{j}$. That is, $X\sim{\sf Gamma}(m,\lambda_{j})$
if and only if its probability density function and cumulative density
function are given by

$f(x)=\frac{\lambda_{j}^{m}}{\Gamma(m)}x^{m-1}e^{-\lambda_{j}x},\quad x>0,$

and

$F(x)=\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}x),$

where $\Gamma(m)=\int_{0}^{\infty}t^{m-1}e^{-t}dt$ is the Gamma function. The
cumulative distribution function of $Z_{m}$ is given by

$\displaystyle F_{Z_{m}}(z)$ $\displaystyle=$ $\displaystyle
1-\prod_{j=1}^{k}\bigg{[}1-P\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\leq
z\bigg{)}\bigg{]}=1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}z)\bigg{]}.$

We have

$\displaystyle E[T_{m}]=$ $\displaystyle\frac{\int_{0}^{\infty}y\delta
e^{-\delta
y}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]}dy}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$
$\displaystyle+\frac{(1-\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta y}ds)\int_{0}^{\infty}w\eta e^{-\eta
w}dw}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}s(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]})^{\prime}e^{-\delta
s}ds}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds},$

and

$\displaystyle P_{mk}(t)=$
$\displaystyle\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]}-(1-e^{-\delta
t})\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]}$
$\displaystyle+(\frac{P(Y\leq t)}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq
Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{\infty}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}(t-s))\bigg{]}$
$\displaystyle-(1-e^{-\delta(t-s)})\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}(t-s))\bigg{]}dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{P(Y\leq t-s)}{P(Y\leq Z_{m})}1_{\\{0\leq
t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}-\frac{P\\{Y+W\leq t-s\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t-s\leq W+Z_{m}\\}}$
$\displaystyle-1_{\\{t-s>W+Z_{m}\\}})\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}dydG(s).$

We perform a simulation study to estimate $E[T_{m}]$ and $P_{m5}(t)$ as
function of $m$ and $t$ for $k=5$. For $m=1,2,\dots,40$, we generate $T_{m}$
for the $N=30,000$ iterations. Then we estimate $E[T_{m}]$ as
$\frac{1}{N}\sum_{i=1}^{N}T_{mi}$. Similarly, for $P_{m5}(t)$, we also
performed $30,000$ iterations to estimate each $P_{m5}(t)$ for $m=1,\dots,40$.
Figure 1 displays the values of $E[T_{m}]$ and $P_{m5}(3)$ under exponential
hacking, exponential re-setting, and exponential detecting times for
$m=1,2,\dots,40$. It shows that $E[T_{m}]$ and $P_{m5}(3)$ are both increasing
function of $m$.

(a) $E(T_{m})$ as a function of $m$

(b) $P_{m5}(3)$ as a function of $m$

Figure 1: Figure 1: Plots of $E[T_{m}]$ and $P_{m5}(3)$ under exponential
hacking, exponential re-set, and exponential detecting times.

Figure 2 provides the $3D$ version of $P_{m5}(t)$ and expected net revenue
$E_{m}[NR]$, based on $N=30,000$ iterations. Panel $(a)$ of Figure 2 plots
$P_{m5}(t)$ as a function of $m$ and $t$. The blue dotted line in $P_{m5}(t)$
is a $3D$ version of $P_{m5}(3)$. Panel $(b)$ of Figure 2 demonstrates the
expected net revenues per unit of time for various values of $m$ when
$R(m)=0.2m,C_{1}(m)=2m^{0.2},$ and $C_{2}(m)=2m^{0.3}$. It is also based on
$30,000$ iterations. The red vertical line shows that the expected net revenue
achieves its maximum value of $1.87$ million dollars at $m=40$ for $1\leq
m\leq 40$.

(a) $P_{m5}(t)$ as a bivariate function of $m$ and $t$

(b) $E_{m}(NR)$ as a function of $m$

Figure 2: Figure 2: Plots of $P_{m5}(t)$ and $E_{m}(NR)$ under exponential
hacking, exponential re-set, and exponential detecting times.

######  Example 2 (Gamma hacking, Gamma re-set, and Gamma detecting times).

Suppose $X_{1}^{j},\ldots,$
$X_{n}^{j}\overset{iid}{\sim}{\sf Gamma}(\eta_{j},\delta_{j})$, $Y\sim{\sf
Gamma}(\alpha,\beta)$, and $W\sim{\sf Gamma}(\theta,\tau)$ such that
$X_{1}\dots X_{n}$, $Y$, and $W$ are mutually independent, where ${\sf
Gamma}(\eta,\delta)$, ${\sf Gamma}(\alpha,\beta)$, and ${\sf
Gamma}(\theta,\tau)$ denote the gamma distributions with shape parameters
$\eta>0$, $\alpha>0$, and $\theta>0$ and rate parameters $\delta>0$,
$\beta>0$, and $\tau>0$, respectively. Then $\sum_{i=1}^{m}X_{i}\sim{\sf
Gamma}(m\eta_{j},\delta_{j})$. It follows that

$\displaystyle E[T_{m}]$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}e^{-\beta
y}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]}dy}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}\int_{0}^{\infty}\bigg{(}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}w\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau
w}dsdw}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})^{\prime}(s-\frac{s}{\Gamma(\alpha)}\gamma(\alpha,\beta
s))ds}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds},$

and

$\displaystyle P_{mk}(t)$ $\displaystyle=$
$\displaystyle\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)$
$\displaystyle+\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}\bigg{)}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)$ $\displaystyle+(\frac{\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq
Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
y)(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{t}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]}-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}-\frac{P\\{Y+W\leq
t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq t-s\leq W+Z_{m}\\}}$
$\displaystyle-1_{\\{t-s>W+Z_{m}\\}})\int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
y)(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dydG(s).$

As in Example 1, we performed $N=30,000$ iterations to estimate $E[T_{m}]$ for
$m=1,2,\dots,40$. Figure 3 displays plots of $E[T_{m}]$ and $P_{m4}(2)$ for
various values of $m$. Clearly, $P_{m4}(2)$ and $E[T_{m}]$ are increasing
functions of $m$. As $m$ increases, $P_{m4}(2)$ approaches $1$.

(a) $E(T_{m})$ as a function of $m$

(b) $P_{m4}(2)$ as a function of $m$

Figure 3: Figure 3: Plots of $E[T_{m}]$ and $P_{m4}(2)$ under gamma hacking,
gamma re-set, and gamma detecting times.

(a) $P_{m4}(t)$ as a bivariate function of $m$ and $t$

(b) $E_{m}(NR)$ as a function of $m$

Figure 4: Figure 4: Plots of $P_{m4}(t)$ and $E_{m}(NR)$ under gamma hacking,
gamma re-set, and gamma detecting times.

Panel (a) of Figure 4 provides a 3D plot of $P_{m4}(t)$ as a function of $m$
and $t$ for $k=4$ based on $30,000$ iteration. The blue dotted line is the 3D
version of $P_{m4}(2)$. Panel (b) of Figure 4 demonstrates $E_{m}[NR]$ for
$R(m)=0.3m+4,C_{1}(m)=2.5m^{2},$ and $C_{2}(m)=2m^{0.3}$, based on $50,000$
iterations. It shows that $E_{m}[NR]$ achieves the maximum value of $2.08$
million dollars when $m=1$.

######  Example 3 (Gamma hacking, Gamma re-set, and Weibull detecting times).

The setup here is the same as the one in Example 2 except that here we have
$Y\sim\mathsf{Weibull}(\alpha,\beta)$. That is, the probability density
function of $Y$ is given by

$\displaystyle
f(y)=\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(y/\alpha\right)^{\beta}},\quad
y>0.$

In this case, we have

$\displaystyle E[T_{m}]$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(y/\alpha\right)^{\beta}}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}y)\bigg{]}dy}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle+\frac{(1-\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})f_{Y}(s)ds)\int_{0}^{\infty}w\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau
w}dw}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})^{\prime}e^{-(s/\alpha)^{\beta}}ds}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$

and

$\displaystyle P_{mk}(t)$ $\displaystyle=$
$\displaystyle\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}-1+e^{-(t/\alpha)^{\beta}}$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]})(1-e^{-(t/\alpha)^{\beta}})$
$\displaystyle+(\frac{1-e^{-(t/\alpha)^{\beta}}}{P(Y\leq Z_{m})}1_{\\{0\leq
t\leq Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-(y/\alpha)^{\beta}})(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{t}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]}-1+e^{-((t-s)/\alpha)^{\beta}}$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})(1-e^{-((t-s)/\alpha)^{\beta}})dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{1-e^{-((t-s)/\alpha)^{\beta}}}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}$
$\displaystyle-\frac{P\\{Y+W\leq t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq
t-s\leq W+Z_{m}\\}}-1_{\\{t-s>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-(y/\alpha)^{\beta}})(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dydG(s).$

Here we performed $N=20,000$ iterations to estimate $E[T_{m}]$ for
$m=1,2,\dots,40$. Figure 5 displays plots of $E[T_{m}]$ and $P_{m7}(7)$ as $m$
increases. It shows that $P_{m7}(7)$ and $E[T_{m}]$ are increasing functions
of $m$.

(a) $E(T_{m})$ as a function of $m$

(b) $P_{m7}(7)$ as a function of $m$

Figure 5: Figure 5: Plots of $E[T_{m}]$ and $P_{m7}(7)$ under gamma hacking,
gamma re-set, and Weibull detecting times.

(a) $P_{m7}(t)$ as a bivariate function of $m$ and $t$

(b) $E_{m}(NR)$ as a function of $m$

Figure 6: Figure 6: Plots of $P_{m7}(t)$ and $E_{m}(NR)$ under gamma hacking,
gamma re-set, and Weibull detecting times.

Figure 6 demonstrates the $3D$ version of $P_{m7}(t)$ and $E_{m}[NR]$. Panel
$(a)$ of Figure 6 provides a plot of $P_{m7}(t)$ as a function of $m$ and $t$,
based on $N=20,000$ iterations. The blue dotted line is the $3D$ version of
$P_{m7}(4)$. Panel $(b)$ of Figure 6 demonstrates $E_{m}[NR]$ for different
values of $m$ when $R(m)=1.3m^{1.2},C_{1}(m)=0.5m^{2},$ and
$C_{2}(m)=0.4m^{1.5}$, based on $N=20,000$ iterations. The red vertical line
shows that the expected net revenue achieved its maximum value $7.47$ millions
dollars at $m=24$.

##  5 Concluding remarks

This article studied the stochastic behavior of an $n$-node blockchain under
cyber attacks from multiple hackers when re-setting the blockchain takes a
random amount of time. Under the assumption that the detecting times, hacking
times, and re-setting times follow arbitrary distributions, we have derived
the instantaneous functional probability, limiting functional probability, and
mean functional time. Moreover, we have shown that these three quantities are
all increasing functions of the number of nodes, providing theoretical backup
for the intuition that the more nodes a blockchain has, the safer it is.

## Reference

Choi, T., Guo, S., Liu, N., and Shi, X. (2020). Optimal pricing in on-demand-
service-platform-operations with hired agents and risk-sensitive customers in
the blockchain era. _European Journal of Operational Research_ 284(3),
1031–1042.

Casino, F., Dasaklis, T. K. and Patsakis, C. (2019). A systematic literature
review of blockchain-based applications: current status, classification and
open issues. _Telematics and Informatics_ 36, 55–81.

Corbet, S., Cumming, D. , Lucey, B. , Peat, M. , Vigne, S. (2020). The
destabilising effects of cryptocurrency cybercriminality. _Economics Letters_
volume 191.

Diffie, W. and Hellman, M.E. (1976). New directions in cryptography. _IEEE
Transactions on Information Theory_ 22(5), 644–654.

Goldwasser, S. and Micali, S. (1982). Probabilistic encryption. _Journal of
Computer System and Sciences_ 28(2), 270–299.

Goldwasser, S., Micali, S. and Rivest, R.L. (1988). A digital signature scheme
secure against adaptive chosen-message attacks. _SIAM Journal on Computing_
18(1), 283–308.

Hussain, S., Sivakumar, T.B., and Khang, A. (2022). Cryptocurrency
Methodologies and Techniques. _The Data-Driven Blockchain Ecosystem_ page 9.

Khan, M.A. and Salah, K. (2017). IoT security: review, blockchain solutions,
and open challenges. _Future Generation Computer Systems._ 82, 395–411.

Khanum, S. and Mustafa, K. (2022). A systematic literature review on sensitive
data protection in blockchain applications. _Concurrency and Computation
Practice and Experience._ 35(1).

Kumar, A., Abhishek, K., Ghalib, M.R., Nerurkar, P., Bhirud, S., Alnumay, W.,
Kumar, S.A., Chatterjee, P. and Ghosh, U. (2020). Securing logistics system
and supply chain using Blockchain. _Applied stochastic models in business and
industry_ 37(3), 413–428.

Li, C. and Zhang, L.J. (2017). A blockchain based new secure multi-layer
network model for internet of things. _Proceedings–2017 IEEE 2nd Intertional
Congress on Internet of Things_ , 33–41.

Melo, C., Dantas, J. , Pereira, P. and Maciel, P. (2021). Distributed
application provisioning over Ethereum-based private and permissioned
blockchain: availability modeling, capacity, and costs planning. _The Journal
of Supercomputing_ 77, 9615–9641.

Meng, W., Tischhauser, E. , Wang, Q. , Wang, Y. and Han, J. (2017). When
instruction detection meets blockchain technology: a review. _IEEE Access_ 6,
10179–10188.

Xu, X., and Hong, L. (2024). Instantaneous and limiting behavior of an n-node
blockchain under cyber attacks from a single hacker. Business Analytics in
Practice: Enhancing Decision Making (Edited by Emrouznejad et al.) _Lecture
Notes in Operations Research_ , Springer, New York, to appear.

Rivest, R.L. , Shamir, A. and Adelman, L. (1978). A method for obtaining
digital signatures and public-key cryptosystems. _Communications of the ACM_
21(2), 120–126.

Tama, B.A., Kweka, B.J., Park, Y. and Rhee, K.H. (2017). A critical review of
blockchain and its current applications. _2017 International Conference on
Electrical Engineering and Computer Science IEEE_ , 109–113.

Tsuchiya, Y. and Hiramoto, N. (2021). How cryptocurrency is laundered: Case
study of Coincheck hacking incident. _Forensic Science International: Reports_
volume 4.

Zhang, Z., Zargham., and M., V. (2020). On modeling blockchain-enabled
economic networks as stochastic dynamical systems. _Applied Network Science_
5, 19.

Zheng, Z., Xie, S. , Dai, H.N., and Wang, H. (2018). Blockchain challenges and
opportunities: a survey. _International Journal of Web and Grid Services_
14(4), 352–375.

## Appendix

### Some general formulas

Since $Z_{m}=\min_{j=1,\dots,k}\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\bigg{)}$, we
have

$\displaystyle F_{Z_{m}}(z)$ $\displaystyle=$ $\displaystyle
P\bigg{(}\min_{j=1,\dots,k}\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\bigg{)}\leq
z\bigg{)}$ $\displaystyle=$ $\displaystyle
1-P\bigg{(}\min_{j=1,\dots,k}\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\bigg{)}>z\bigg{)}$
$\displaystyle=$ $\displaystyle
1-P\bigg{(}\sum_{i=1}^{m}X_{i}^{1}>z,\sum_{i=1}^{m}X_{i}^{2}>z,\dots,\sum_{i=1}^{m}X_{i}^{k}>z\bigg{)}$
$\displaystyle\stackrel{{\scriptstyle ind}}{{=}}$ $\displaystyle
1-\prod_{j=1}^{k}P\bigg{(}\sum_{i=1}^{m}X_{i}^{j}>z\bigg{)}$ $\displaystyle=$
$\displaystyle 1-\prod_{j=1}^{k}\bigg{[}1-P\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\leq
z\bigg{)}\bigg{]}.$

Equation (5) in Hong and Sarkar (2013) implies

$\displaystyle E[Y_{1}]$ $\displaystyle=$ $\displaystyle
E\bigg{[}Y\bigg{|}Y<Z_{m}\bigg{]}=\int_{0}^{\infty}yf_{Y|Y<Z_{m}}(y)dy$
$\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}\int_{y}^{\infty}yg(y,s)dsdy}{P(Y<Z_{m})}=\frac{\int_{0}^{\infty}\int_{y}^{\infty}yf_{Y}(y)f_{Z_{m}}(s)dsdy}{1-\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds},$

where the joint probabiliyt density function of $Y$ and $Z_{m}$, by
independence of $Y$ and $Z_{m}$, is given by

$\displaystyle g(y,s)=f_{Y}(y)f_{Z_{m}}(s),\quad y>0,s>0.$

It is clear that (3) implies that

$p_{mk}=P\bigg{\\{}Y>Z_{m}\bigg{\\}}=\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds.$

Hence,

$\displaystyle E[N_{1}][Y_{1}]$ $\displaystyle=$
$\displaystyle\frac{1-p_{mk}}{p_{mk}}E\bigg{[}Y\bigg{|}Y<Z_{m}\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}\int_{y}^{\infty}yf_{Y}(y)f_{Z_{m}}(s)dsdy}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}$
$\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}yf_{Y}(y)\int_{y}^{\infty}f_{Z_{m}}(s)dsdy}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}$
$\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}yf_{Y}(y)(1-F_{Z_{m}}(y))dy}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}.$

Similarly,

$E[N_{1}][W_{1}]=\frac{1-p_{mk}}{p_{mk}}E[W]=\frac{(1-\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds)\int_{0}^{\infty}wf_{W}(w)dw}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds},$

and

$\displaystyle E\bigg{[}Z_{m}\bigg{|}Z_{m}<Y\bigg{]}$
$\displaystyle=\frac{\int_{0}^{\infty}\int_{s}^{\infty}sf_{Y}(y)f_{Z_{m}}(s)dyds}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}$
$\displaystyle=\frac{\int_{0}^{\infty}sf_{Z_{m}}(s)\int_{s}^{\infty}f_{Y}(y)dyds}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}$
$\displaystyle=\frac{\int_{0}^{\infty}sf_{Z_{m}}(s)(1-F_{Y}(s))ds}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}.$

Therefore,

$\begin{split}E[T_{m}]=&\frac{\int_{0}^{\infty}yf_{Y}(y)(1-F_{Z_{m}}(y))dy+(1-\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds)\int_{0}^{\infty}wf_{W}(w)dw}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}\\\
&+\frac{\int_{0}^{\infty}sf_{Z_{m}}(s)(1-F_{Y}(s))ds}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}.\end{split}$
(13)

Next, we shift our attention to the terms in (3.3). We have

$\displaystyle P\bigg{\\{}t<Y\wedge Z_{m}\bigg{\\}}$
$\displaystyle=1-P\bigg{\\{}Z_{m}\leq t\bigg{\\}}-P\\{Y\leq
t\\}+P\bigg{\\{}Z_{m}\leq t,Y\leq t\bigg{\\}}.$

Similarly,

$\displaystyle P\bigg{\\{}Y\wedge Z_{m}>t-s\bigg{\\}}=$ $\displaystyle
1-P\bigg{\\{}Z_{m}\leq t-s\bigg{\\}}-P\\{Y\leq t-s\\}$
$\displaystyle+P\bigg{\\{}Z_{m}\leq t-s,Y\leq t-s\bigg{\\}}.$

Also,

$\displaystyle P\left\\{Y\leq t<Y+W_{1},Y\leq Z_{m}\right\\}$ $\displaystyle=$
$\displaystyle P\left\\{Y\leq t<Y+W_{1}\mid Y\leq Z_{m}\right\\}P\left\\{Y\leq
Z_{m}\right\\}$ $\displaystyle=$ $\displaystyle P\left\\{Y_{1}\leq
t<Y_{1}+W_{1}\right\\}P\left\\{Y\leq Z_{m}\right\\}$ $\displaystyle=$
$\displaystyle\left[1-P\\{Y_{1}>t\text{ or }Y_{1}+W_{1}\leq
t\\}\right]P\left\\{Y\leq Z_{m}\right\\}$ $\displaystyle=$
$\displaystyle[1-P\\{Y_{1}>t\\}-P\\{Y_{1}+W_{1}\leq
t\\}+P\\{Y_{1}>t,Y_{1}+W_{1}\leq t\\}]P\left\\{Y\leq Z_{m}\right\\}$
$\displaystyle=$ $\displaystyle[1-P\\{Y_{1}>t\\}-P\\{Y_{1}+W_{1}\leq
t\\}]P\left\\{Y\leq Z_{m}\right\\}$ $\displaystyle=$
$\displaystyle[P\\{Y_{1}\leq t\\}-P\\{Y_{1}+W_{1}\leq t\\}]P\left\\{Y\leq
Z_{m}\right\\},$

where $Y_{1}\overset{d}{=}Y\mid Y\leq Z_{m}$. Similarly,

$\displaystyle P\left\\{Y\leq t-s<Y+W_{1},Y\leq Z_{m}\right\\}$
$\displaystyle=$ $\displaystyle[P\\{Y_{1}\leq t-s\\}-P\\{Y_{1}+W_{1}\leq
t-s\\}]P\left\\{Y\leq Z_{m}\right\\}.$

Therefore,

$\begin{split}P_{mk}(t)=&1-F_{Z_{m}}(t)-F_{Y}(t)+F_{Z_{m}}(t)F_{Y}(t)\\\
&+(F_{Y_{1}}(t)-F_{Y_{1}+W_{1}}(t))\int_{0}^{\infty}F_{Y}(y)f_{Z_{m}}(y)dy\\\
&+\int_{0}^{t}1-F_{Z_{m}}(t-s)-F_{Y}(t-s)+F_{Z_{m}}(t-s)F_{Y}(t-s)dG(s)\\\
&+\int_{0}^{t}(F_{Y_{1}}(t-s)-F_{Y_{1}+W_{1}}(t-s))\int_{0}^{\infty}F_{Y}(y)f_{Z_{m}}(y)dydG(s),\end{split}$
(14)

where $G(s)=\sum_{n=1}^{\infty}F_{\sum_{i=1}^{n}(Y_{i}+W_{i})}(s)$. Also, we
have

$F_{Y_{1}}(t)=\frac{P\bigg{\\{}Y\leq t,Y\leq Z_{m}\bigg{\\}}}{P\bigg{\\{}Y\leq
Z_{m}\bigg{\\}}}=\begin{cases}1&\text{if $t>Z_{m}$},\\\ \frac{P\\{Y\leq
t\\}}{P\\{Y\leq Z_{m}\\}}&\text{if $0\leq t\leq Z_{m}$},\\\ 0&\text{if
$t<0$},\end{cases}$ (15)
$\begin{split}F_{Y_{1}+W_{1}}(t)&=\frac{P\bigg{\\{}Y+W\leq t,Y\leq
Z_{m}\bigg{\\}}}{P\bigg{\\{}Y\leq Z_{m}\bigg{\\}}}\\\
&=\begin{cases}1&\text{if $t>W+Z_{m}$},\\\ \frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}&\text{if $0\leq t\leq W+Z_{m}$},\\\ 0&\text{if
$t<0$},\end{cases}\end{split}$ (16)

and

$\begin{split}F_{\sum_{i=1}^{n}(Y_{i}+W_{i})}(s)&=\frac{P\bigg{\\{}\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq
s,Y^{\star}_{i}\leq Z_{mi}\bigg{\\}}}{P\bigg{\\{}Y^{\star}_{i}\leq
Z_{mi}\bigg{\\}}}\\\ &=\begin{cases}1&\text{if
$t>\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}$},\\\
\frac{P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq t\\}}{P\\{Y\leq
Z_{m}\\}}&\text{if $0\leq t\leq\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}$},\\\
0&\text{if $t<0$},\end{cases}\end{split}$ (17)

where $Y^{\star}_{i}\overset{iid}{\sim}F_{Y}(y)$ and $Z_{mi}$ denotes the
minimum time spend on hacking into the $m$ node during the $ith$ cycle. Thus,
(15), (16), and (17) can also be written as

$\displaystyle F_{Y_{1}}(t)$ $\displaystyle=$ $\displaystyle\frac{P(Y\leq
t)}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq Z_{m}\\}}+1_{\\{t>Z_{m}\\}},$
$\displaystyle F_{Y_{1}+W_{1}}(t)$ $\displaystyle=$
$\displaystyle\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq t\leq
W+Z_{m}\\}}+1_{\\{t>W+Z_{m}\\}},$ $\displaystyle
F_{\sum_{i=1}^{n}(Y_{i}+W_{i})}(s)$ $\displaystyle=$
$\displaystyle\frac{P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq
t\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq
t\leq\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}\\}}$
$\displaystyle+1_{\\{t>\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}\\}}.$

### From Example 1

Since $Y\sim{\sf Exp}(\delta)$ and $W\sim{\sf Exp}(\eta)$, (Some general
formulas) implies

$\displaystyle p_{mk}=P\bigg{\\{}Y>Z_{m}\bigg{\\}}$
$\displaystyle=\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\bigg{)}\delta
e^{-\delta s}ds,$

where $\gamma(m,\lambda_{j}s)=\int_{0}^{\lambda_{j}s}t^{m-1}e^{-t}dt$ is the
lower incomplete gamma function. By (13), we have

$\displaystyle E[T_{m}]=$
$\displaystyle\frac{\int_{0}^{\infty}yf_{Y}(y)(1-F_{Z_{m}}(y))dy+(1-\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds)\int_{0}^{\infty}wf_{W}(w)dw}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}$
$\displaystyle+\frac{\int_{0}^{\infty}sf_{Z_{m}}(s)(1-F_{Y}(s))ds}{\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds}$
$\displaystyle=$ $\displaystyle\frac{\int_{0}^{\infty}y\delta e^{-\delta
y}(1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]}))dy}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$

$\displaystyle+\frac{(1-\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta y}ds)\int_{0}^{\infty}w\eta e^{-\eta
w}dw}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}s(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]})^{\prime}(1-(1-e^{-\delta
s}))ds}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\delta e^{-\delta
y}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]}dy}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$
$\displaystyle+\frac{(1-\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta y}ds)\int_{0}^{\infty}w\eta e^{-\eta
w}dw}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}s(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]})^{\prime}e^{-\delta
s}ds}{\int_{0}^{\infty}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\delta
e^{-\delta s}ds}.$

By (14), we have

$\displaystyle P_{mk}(t)$ $\displaystyle=$ $\displaystyle
1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]})-(1-e^{-\delta
t})$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]})(1-e^{-\delta
t})$ $\displaystyle+(\frac{1-e^{-\delta t}}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq
Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{\infty}1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}(t-s))\bigg{]})-(1-e^{-\delta(t-s)})$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}(t-s))\bigg{]})(1-e^{-\delta(t-s)})dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{1-e^{-\delta(t-s)}}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}-\frac{P\\{Y+W\leq
t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq t-s\leq W+Z_{m}\\}}$
$\displaystyle-1_{\\{t-s>W+Z_{m}\\}})\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}dydG(s)$
$\displaystyle=$
$\displaystyle\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]}-(1-e^{-\delta
t})\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]}$
$\displaystyle+(\frac{1-e^{-\delta t}}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq
Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{\infty}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}(t-s))\bigg{]}-(1-e^{-\delta(t-s)})\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}(t-s))\bigg{]}dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{1-e^{-\delta(t-s)}}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}-\frac{P\\{Y+W\leq
t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq t-s\leq W+Z_{m}\\}}$
$\displaystyle-1_{\\{t-s>W+Z_{m}\\}})\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}dydG(s),$

where

$\displaystyle\gamma(m,\lambda_{j}s)=\int_{0}^{\lambda_{j}s}z^{m-1}e^{-z}dz,$
$\displaystyle G(s)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq
t\\}}{\int_{0}^{\infty}(1-e^{-\delta
y})(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}t)\bigg{]})^{\prime}dy}1_{0\leq
t\leq\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}}$
$\displaystyle+1_{t>\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}},$

and

$\displaystyle
f_{Z_{m}}(y)=\frac{F_{Z_{m}}(y)}{dy}=(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}y)\bigg{]})^{\prime}.$

Note that $P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq t\\}$ is the
distribution function of the sum of two independent gamma random variables,
$\sum_{i=1}^{n}Y_{i}^{\star}\sim{\sf Gamma}(n,\delta)$ and
$\sum_{i=1}^{n}W_{i}\sim{\sf Gamma}(n,\eta)$. Thus, Moschopoulos (1985)
implies

$\displaystyle P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq
t\\}=C\sum_{k=0}^{\infty}\delta_{l}\int_{0}^{t}\frac{y^{2n+k-1}e^{-y(\delta\vee\eta)}(\delta\vee\eta)^{2n+l}}{\Gamma(2n+l)}dy,$

where

$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{(\delta\vee\eta)^{2n}}{(\delta\eta)^{n}},$
$\displaystyle\delta_{l+1}$ $\displaystyle=$
$\displaystyle\frac{1}{l+1}\sum_{i=1}^{l+1}i\gamma_{i}\delta_{l+1-i},$

for $l=0,1,2,\dots$, and

$\displaystyle\gamma_{j}=\frac{n}{j}((1-\frac{\delta}{\delta\vee\eta})^{j}+(1-\frac{\eta}{\delta\vee\eta})^{j}),$

for $j=1,2,\dots$. Moreover,

$\displaystyle P(Y+W\leq t)=1-\frac{\eta}{\eta-\delta}e^{-\delta
t}+\frac{\delta}{\eta-\delta}e^{-\eta t},$

and

$\displaystyle P(Y\leq Z_{m})$ $\displaystyle=$ $\displaystyle
1-\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m)}\gamma(m,\lambda_{j}s)\bigg{]}\bigg{)}\delta
e^{-\delta s}ds.$

### From Example 2

$\sum_{i=1}^{m}X_{i}^{j}\sim{\sf Gamma}(m\eta_{j},\delta_{j})$, $Y\sim{\sf
Gamma}(\alpha,\beta)$, and $W\sim{\sf Gamma}(\theta,\tau)$. It is clear that

$\displaystyle F_{Z_{m}}(z)$ $\displaystyle=$ $\displaystyle
1-\prod_{j=1}^{k}\bigg{[}1-P\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\leq
z\bigg{)}\bigg{]}$ $\displaystyle=$ $\displaystyle
1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]},$
$\displaystyle
f_{Y}(y)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}e^{-\beta y},$
$\displaystyle F_{Y}(y)=\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta y),$
$\displaystyle
f_{W}(w)=\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau w},$

and

$\displaystyle F_{W}(w)=\frac{1}{\Gamma(\theta)}\gamma(\theta,\tau w).$

By (13), we have

$\displaystyle E[T_{m}]$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}e^{-\beta
y}(1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]}))dy}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$
$\displaystyle+\frac{(1-\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds)\int_{0}^{\infty}w\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau
w}dw}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}s(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})^{\prime}(1-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
s))ds}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}e^{-\beta
y}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]}dy}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}\int_{0}^{\infty}\bigg{(}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}w\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau
w}dsdw}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})^{\prime}(s-\frac{s}{\Gamma(\alpha)}\gamma(\alpha,\beta
s))ds}{\int_{0}^{\infty}\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]}\bigg{)}\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds}.$ $\displaystyle P_{mk}(t)$ $\displaystyle=$ $\displaystyle
1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]})-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)$
$\displaystyle+\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}\bigg{)}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)$ $\displaystyle+(\frac{\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq
Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
y)(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{t}1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}-\frac{P\\{Y+W\leq
t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq t-s\leq W+Z_{m}\\}}$
$\displaystyle-1_{\\{t-s>W+Z_{m}\\}})\int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
y)(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dydG(s)$
$\displaystyle=$
$\displaystyle\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)$
$\displaystyle+\bigg{(}1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}\bigg{)}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)$ $\displaystyle+(\frac{\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
t)}{P(Y\leq Z_{m})}1_{\\{0\leq t\leq
Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
y)(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{t}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]}-\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta(t-s))}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}-\frac{P\\{Y+W\leq
t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq t-s\leq W+Z_{m}\\}}$
$\displaystyle-1_{\\{t-s>W+Z_{m}\\}})\int_{0}^{\infty}\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta
y)(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dydG(s),$

where

$\displaystyle P\\{Y+W\leq
t\\}=C\sum_{k=0}^{\infty}\delta_{k}\int_{0}^{t}\frac{y^{\alpha+\theta+k-1}e^{-y(\beta\vee\tau)}(\beta\vee\tau)^{\alpha+\theta+k}}{\Gamma(\alpha+\theta+k)}dy,$

and

$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{(\beta\vee\tau)^{2n}}{(\beta\tau)^{n}},$
$\displaystyle\delta_{k+1}$ $\displaystyle=$
$\displaystyle\frac{1}{k+1}\sum_{i=1}^{k+1}i\gamma_{i}\delta_{k+1-i},$

for $k=0,1,2,\dots$,

$\displaystyle\gamma_{j}=\frac{\alpha}{j}(1-\frac{\beta}{\beta\vee\tau})^{j}+\frac{\theta}{j}(1-\frac{\tau}{\beta\vee\tau})^{j}$

for $j=1,2,\dots$, and

$\displaystyle G(S)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq
t\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq
t\leq\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}\\}}$
$\displaystyle+1_{\\{t>\sum_{i=1}^{n}W_{i}+\sum_{i=1}^{n}Z_{mi}\\}}.$

Since $\sum_{i=1}^{n}Y_{i}^{\star}\sim{\sf Gamma}(n\alpha,\beta)$ and
$\sum_{i=1}^{n}W_{i}\sim{\sf Gamma}(n\theta,\tau)$,

$\displaystyle P\\{\sum_{i=1}^{n}(Y^{\star}_{i}+W_{i})\leq
t\\}=C\sum_{k=0}^{\infty}\delta_{k}\int_{0}^{t}\frac{y^{n\alpha+n\theta+k-1}e^{-y(\beta\vee\tau)}(\beta\vee\tau)^{n\alpha+n\theta+k}}{\Gamma(n\alpha+n\theta+k)}dy,$

where

$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{(\beta\vee\tau)^{2n}}{(\beta\tau)^{n}},$
$\displaystyle\delta_{k+1}$ $\displaystyle=$
$\displaystyle\frac{1}{k+1}\sum_{i=1}^{k+1}i\gamma_{i}\delta_{k+1-i},$

for $k=0,1,2,\dots$, and

$\displaystyle\gamma_{j}=\frac{n\alpha}{j}(1-\frac{\beta}{\beta\vee\tau})^{j}+\frac{n\theta}{j}(1-\frac{\tau}{\beta\vee\tau})^{j},$

for $j=1,2,\dots$.Moreover,

$\displaystyle P(Y\leq Z_{m})$ $\displaystyle=$ $\displaystyle
1-\int_{0}^{\infty}F_{Z_{m}}(s)f_{Y}(s)ds$ $\displaystyle=$ $\displaystyle
1-\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta^{\alpha}}{\Gamma(\alpha)}s^{\alpha-1}e^{-\beta
s}ds.$

### From Example 3

Since $\sum_{i=1}^{m}X_{i}^{j}\sim{\sf Gamma}(m\eta_{j},\delta_{j})$,
$Y\sim\mathsf{Weibull}(\alpha,\beta)$, and $W\sim{\sf Gamma}(\theta,\tau)$. We
have

$\displaystyle F_{Z_{m}}(z)$ $\displaystyle=$ $\displaystyle
1-\prod_{j=1}^{k}\bigg{[}1-P\bigg{(}\sum_{i=1}^{m}X_{i}^{j}\leq
z\bigg{)}\bigg{]}$ $\displaystyle=$ $\displaystyle
1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]},$
$\displaystyle
f_{Y}(y)=\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(y/\alpha\right)^{\beta}},$
$\displaystyle F_{Y}(y)=1-e^{-(y/\alpha)^{\beta}},$ $\displaystyle
f_{W}(w)=\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau w},$

and

$\displaystyle F_{W}(w)=\frac{1}{\Gamma(\theta)}\gamma(\theta,\tau w).$

We have

$\displaystyle E[T_{m}]$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(y/\alpha\right)^{\beta}}(1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}y)\bigg{]}))dy}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle+\frac{(1-\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})f_{Y}(s)ds)\int_{0}^{\infty}w\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau
w}dw}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})^{\prime}(1-(1-e^{-(s/\alpha)^{\beta}}))ds}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}y\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(y/\alpha\right)^{\beta}}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}y)\bigg{]}dy}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle+\frac{(1-\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})f_{Y}(s)ds)\int_{0}^{\infty}w\frac{\tau^{\theta}}{\Gamma(\theta)}w^{\theta-1}e^{-\tau
w}dw}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}$
$\displaystyle+\frac{\int_{0}^{\infty}(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})^{\prime}e^{-(s/\alpha)^{\beta}}ds}{\int_{0}^{\infty}(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}s)\bigg{]})\frac{\beta}{\alpha}\left(\frac{y}{\alpha}\right)^{\beta-1}e^{-\left(s/\alpha\right)^{\beta}}ds}.$

Also,

$\displaystyle P_{mk}(t)$ $\displaystyle=$ $\displaystyle
1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]})-(1-e^{-(t/\alpha)^{\beta}})$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]})(1-e^{-(t/\alpha)^{\beta}})$
$\displaystyle+(\frac{1-e^{-(t/\alpha)^{\beta}}}{P(Y\leq Z_{m})}1_{\\{0\leq
t\leq Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-(y/\alpha)^{\beta}})(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{t}1-(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})-(1-e^{-((t-s)/\alpha)^{\beta}})$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})(1-e^{-((t-s)/\alpha)^{\beta}})dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{1-e^{-((t-s)/\alpha)^{\beta}}}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}$
$\displaystyle-\frac{P\\{Y+W\leq t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq
t-s\leq W+Z_{m}\\}}-1_{\\{t-s>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}1-e^{-(y/\alpha)^{\beta}}(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dydG(s)$
$\displaystyle=$
$\displaystyle\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]}-1+e^{-(t/\alpha)^{\beta}}$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}t)\bigg{]})(1-e^{-(t/\alpha)^{\beta}})$
$\displaystyle+(\frac{1-e^{-(t/\alpha)^{\beta}}}{P(Y\leq Z_{m})}1_{\\{0\leq
t\leq Z_{m}\\}}+1_{\\{t>Z_{m}\\}}-\frac{P\\{Y+W\leq t\\}}{P\\{Y\leq
Z_{m}\\}}1_{\\{0\leq t\leq W+Z_{m}\\}}-1_{\\{t>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-(y/\alpha)^{\beta}})(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dy$
$\displaystyle+\int_{0}^{t}\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]}-1+e^{-((t-s)/\alpha)^{\beta}}$
$\displaystyle+(1-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}(t-s))\bigg{]})(1-e^{-((t-s)/\alpha)^{\beta}})dG(s)$
$\displaystyle+\int_{0}^{t}(\frac{1-e^{-((t-s)/\alpha)^{\beta}}}{P(Y\leq
Z_{m})}1_{\\{0\leq t-s\leq Z_{m}\\}}+1_{\\{t-s>Z_{m}\\}}$
$\displaystyle-\frac{P\\{Y+W\leq t-s\\}}{P\\{Y\leq Z_{m}\\}}1_{\\{0\leq
t-s\leq W+Z_{m}\\}}-1_{\\{t-s>W+Z_{m}\\}})$
$\displaystyle\times\int_{0}^{\infty}(1-e^{-(y/\alpha)^{\beta}})(-\prod_{j=1}^{k}\bigg{[}1-\frac{1}{\Gamma(m\eta_{j})}\gamma(m\eta_{j},\delta_{j}z)\bigg{]})^{\prime}dydG(s).$


Denoting the above without the last ${\rm o}_{\rm P}(1)$ as
$v(\bm{u}^{(1)},\bm{u}^{(2)})$, since the random function in (31) is convex
with respect to $(\bm{u}^{(1)},\bm{u}^{(2)})$, the argmin lemma in Hjort and
Pollard (2011) verifies that the minimum point of the random function
converges in probability to $\mathop{\rm
argmin}_{(\bm{u}^{(1)},\bm{u}^{(2)})}v(\bm{u}^{(1)},\bm{u}^{(2)})$. If we are
in the setting of doubly robust estimation where either the modeling of
$e^{[h]}(\bm{\alpha})$ or the modeling of ${\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}$ is correct, the first
expectation, which is the coefficient of $\bm{u}^{(1)}$, becomes ${\rm
E}\\{\sum_{h=1}^{H}\partial\ell^{[h]}(\bm{\beta}_{\lambda}^{*})/\partial\bm{\beta}^{(1)}\\}$,
and we see that each of its components has an absolute value less than
$\lambda$ from (10). Therefore, the minimum point of
$v(\bm{u}^{(1)},\bm{u}^{(2)})$ with respect to $\bm{u}^{(1)}$ is $\bm{0}$, and
then (18) holds. On the other hand, the minimum point for $\bm{u}^{(2)}$ is

$\displaystyle{\rm
E}\Bigg{(}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\Bigg{[}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\Bigg{)}^{-1}$
$\displaystyle\Bigg{\\{}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}-\sqrt{N}\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})$
$\displaystyle\phantom{\Bigg{\\{}}+{\rm
E}\Bigg{(}\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\alpha}^{{\rm
T}}}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}[\ell^{[h]}(\bm{\beta}_{\lambda}^{*})-{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}]\Bigg{)}\sqrt{N}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})$
$\displaystyle\phantom{\Bigg{\\{}}+{\rm
E}\Bigg{[}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial^{2}}{\partial\bm{\gamma}^{{\rm
T}}\partial\bm{\beta}^{(2)}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\sqrt{N}(\hat{\bm{\gamma}}-\bm{\gamma}^{*})\Bigg{\\}},$

and the first expectation for which the inverse matrix is taken is rewritten
as ${\rm
E}\\{\sum_{h=1}^{H}\partial^{2}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})/\allowbreak\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}\\}$ in the setting of doubly robust estimation. Accordingly,
$N\hat{\bm{\beta}}_{\lambda}^{(1)}$ converges in probability to $\bm{0}$ and
$\sqrt{N}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})$ is
asymptotically equivalent to (22). Also, if the modeling of
$e^{[h]}(\bm{\alpha})$ and the modeling of ${\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}$ are both correct, the
minimum point for $\bm{u}^{(2)}$ reduces to

$\displaystyle{\rm
E}\Bigg{\\{}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}^{-1}$
$\displaystyle\Bigg{(}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}-\sqrt{N}\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{)}$

and, thus, (19) holds.

###  A.2 Performance comparison of criteria in case of continuous outcome

Similarly to Section 5.2, the explanatory variable
$\bm{x}=(x_{1},x_{2},\ldots,x_{p})^{{\rm T}}$ follows ${\rm
N}(\bm{0}_{p},\bm{I}_{p})$ and the confounding variable $z$ follows ${\rm
N}(0,1)$ independently of $\bm{x}$. On the other hand, the conditional
distribution of $y^{[h]}$ given the explanatory and confounding variables is
Gaussian, as in Section 5.1. For the distributions of $t^{[h]}$ and $y^{[h]}$
given $z$, which are necessary for the estimation, we suppose a multinomial
logit model and a linear regression model defined by

$\displaystyle{\rm E}(t^{[h]}\mid
z)=\frac{\exp(z\alpha^{[h]})}{\sum_{h^{\prime}=1}^{H}\exp(z\alpha^{[h^{\prime}]})},\qquad
y^{[h]}\mid z\sim{\rm N}(z\gamma,\sigma^{[h]2})$

with $\alpha^{[h]}$ and $\gamma$ as parameters ($h\in\\{1,2,\ldots,H\\}$). As
a natural flow from these models, letting
$\bm{z}\equiv(z,z^{2}/\sqrt{2}-1/\sqrt{2})^{{\rm T}}$,
$\bm{\alpha}\equiv(\alpha_{1}^{[h]},\alpha_{2}^{[h]})^{{\rm T}}$ and
$\bm{\gamma}\equiv(\gamma_{1},\gamma_{2})^{{\rm T}}$, we generate data from

$\displaystyle{\rm E}(t^{[h]}\mid\bm{z})=\frac{\exp(\bm{z}^{{\rm
T}}\bm{\alpha}^{[h]})}{\sum_{h^{\prime}=1}^{H}\exp(\bm{z}^{{\rm
T}}\bm{\alpha}^{[h^{\prime}]})},\qquad y^{[h]}\mid(\bm{x},\bm{z})\sim{\rm
N}(\bm{x}^{{\rm T}}\bm{\beta}^{[h]}+\bm{z}^{{\rm T}}\bm{\gamma},1).$

The true values used here, are $\bm{\beta}^{[1]}=(2\beta^{*},0)^{{\rm T}}$,
$\bm{\beta}^{[2]}=(0,\beta^{*})^{{\rm T}}$,
$\bm{\beta}^{[3]}=(-\beta^{*},0)^{{\rm T}}$ and
$\bm{\beta}^{[4]}=(0,-2\beta^{*})^{{\rm T}}$ when $(p,H)=(2,4)$; when
$(p,H)=(2,6)$, $\bm{\beta}^{[5]}=(2\beta^{*},0)^{{\rm T}}$ and
$\bm{\beta}^{[6]}=(0,\beta^{*})^{{\rm T}}$ are added to these, and when
$(p,H)=(2,8)$, $\bm{\beta}^{[7]}=(-\beta^{*},0)^{{\rm T}}$ and
$\bm{\beta}^{[8]}=(0,-2\beta^{*})^{{\rm T}}$ are further added. When
$(p,H)=(4,4)$, we use the above acted on by $\otimes\bm{1}_{2}$ as
$\bm{\beta}^{[j]}$ ($j\in\\{1,2,3,4\\}$). The true value of $\bm{\alpha}$ is
$(\alpha^{*},0)$ or $(\alpha^{*},\alpha^{*})$, and the true value of
$\bm{\gamma}$ is $(\gamma^{*},0)$ or $(\gamma^{*},\gamma^{*})$. When the true
value of $\bm{\alpha}$ is $(\alpha^{*},\alpha^{*})$, the model of the
assignment variable is misspecified, and when the true value of $\bm{\gamma}$
is $(\gamma^{*},\gamma^{*})$, the model of the outcome variable is
misspecified. In addition, $\alpha^{*}$ is 0.1 or 0.2, $\beta^{*}$ is 0.2 or
0.4, $\gamma^{*}$ is 0.5 or 1.0, and the sample size $N$ is 200 or 400. As for
the distribution of $y^{[h]}$ in the causal inference model, considering that
$\bm{z}\sim{\rm N}(\bm{0}_{2},\bm{I}_{2})$, we use the one obtained by
marginalizing ${\rm N}(\bm{x}^{{\rm T}}\bm{\beta}^{[h]}+\bm{z}^{{\rm
T}}\bm{\gamma},1)$) with respect to $\bm{z}$. Letting
$\bm{x}^{[h]}\equiv\bm{x}\otimes\bm{e}^{[h]}$ and
$\bm{\beta}\equiv(\bm{\beta}^{[1]{\rm T}},\bm{\beta}^{[2]{\rm
T}},\ldots,\bm{\beta}^{[H]{\rm T}})^{{\rm T}}$, it holds that $\bm{x}^{{\rm
T}}\bm{\beta}^{[h]}=\bm{x}^{[h]{\rm T}}\bm{\beta}$, so this model can also be
expressed as (3).

Table 6 shows whether the asymptotic bias in (16) given by Theorem 2 and the
asymptotic bias in (27) given by Theorem 3 are accurate approximations of the
true bias. The number of Monte Carlo iterations is set to 200. It confirms
that we can evaluate the bias somewhat accurately, except in the lower two
settings, where neither the doubly robust estimation nor the bias evaluation
in this paper is theoretically guaranteed. Since QICw evaluates all these
values at $pH/4=2$, the approximation accuracy of the asymptotic bias is at
least better than that of QICw.

Table 6: Bias evaluation in causal inference models with Gaussian noise. The upper entry of the parameters is $(p,H,N)$, the middle entry is $(\beta^{*},\alpha^{*},\gamma^{*})$, and the lower entry represents the misspecified model. The $j\ (\in\\{1,2,\ldots,8\\})$ column lists the bias when the number of elements in the active set is $\hat{p}=j$ minus the bias when $\hat{p}=j-1$. On the other hand, the True1 and True2 rows show the Monte Carlo evaluated true values for the inverse-probability-weighted and doubly robust estimation, while the IPIC and DRIC rows show the average of the unbiased estimators for True1 and True2 based on the respective criteria. parameters |  | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
(2, 4, 200) | True1 | 5.28 | 4.12 | 5.93 | 6.85 | 6.20 | 6.08 | 3.78 | 8.98
(0.2, 0.2, 0.5) | IPIC | 9.82 | 9.01 | 7.55 | 7.09 | 7.21 | 7.22 | 4.63 | 10.73
none | True2 | 5.80 | 3.97 | 4.89 | 4.99 | 4.30 | 5.13 | 2.64 | 5.94
DRIC | 8.93 | 8.12 | 6.88 | 6.50 | 6.36 | 6.68 | 3.96 | 8.20
(2, 4, 200) | True1 | 9.02 | 7.46 | 9.31 | 11.33 | 11.04 | 8.58 | 6.87 | 13.77
(0.2, 0.2, 1.0) | IPIC | 10.30 | 9.03 | 7.38 | 7.98 | 7.45 | 7.19 | 5.46 | 10.24
none | True2 | 6.80 | 4.41 | 5.57 | 5.97 | 5.12 | 5.91 | 2.12 | 8.01
DRIC | 5.96 | 5.51 | 5.47 | 4.81 | 4.84 | 4.70 | 2.58 | 6.61
(2, 4, 200) | True1 | 8.16 | 6.68 | 8.18 | 7.46 | 9.27 | 9.24 | 5.80 | 13.23
(0.2, 0.1, 1.0) | IPIC | 8.46 | 7.93 | 7.50 | 6.89 | 6.87 | 7.43 | 4.68 | 10.36
none | True2 | 6.28 | 4.49 | 4.37 | 4.92 | 4.58 | 4.82 | 2.48 | 7.76
DRIC | 5.67 | 5.27 | 4.68 | 4.76 | 4.59 | 3.93 | 2.26 | 6.75
(2, 4, 200) | True1 | 9.38 | 8.44 | 9.67 | 10.17 | 10.38 | 8.69 | 6.54 | 13.94
(0.4, 0.2, 1.0) | IPIC | 12.97 | 10.30 | 7.59 | 5.71 | 6.45 | 7.16 | 3.94 | 10.76
none | True2 | 11.01 | 6.61 | 5.30 | 3.95 | 3.40 | 4.79 | 2.01 | 6.47
DRIC | 9.62 | 7.01 | 4.95 | 3.98 | 3.40 | 3.57 | 1.61 | 5.88
(2, 4, 400) | True1 | 7.85 | 9.24 | 8.16 | 11.08 | 11.39 | 8.77 | 8.49 | 10.76
(0.2, 0.2, 1.0) | IPIC | 9.74 | 9.31 | 8.58 | 8.23 | 8.51 | 8.68 | 5.95 | 11.32
none | True2 | 6.01 | 6.02 | 4.54 | 4.82 | 7.39 | 4.79 | 4.08 | 5.41
DRIC | 5.91 | 5.71 | 4.95 | 5.01 | 5.19 | 5.07 | 2.93 | 6.83
(2, 4, 200) | True1 | 20.31 | 18.32 | 9.84 | 11.00 | 11.45 | 9.03 | 0.10 | 22.65
(0.2, 0.2, 1.0) | IPIC | 30.37 | 2.91 | 6.17 | 5.49 | 6.39 | 7.41 | 2.62 | 10.75
treatment | True2 | 10.81 | 12.99 | 4.63 | 3.67 | 6.72 | 5.63 | 3.31 | 7.45
DRIC | 5.81 | 5.02 | 5.06 | 4.41 | 4.85 | 4.47 | 2.48 | 6.46
(2, 4, 200) | True1 | 33.42 | 26.09 | 24.23 | 24.23 | 24.11 | 23.16 | 14.60 | 33.78
(0.2, 0.2, 1.0) | IPIC | 41.28 | 27.97 | 18.35 | 16.38 | 17.94 | 16.48 | 11.14 | 23.05
outcome | True2 | 31.36 | 22.93 | 18.70 | 21.77 | 21.32 | 18.41 | 6.82 | 27.42
DRIC | 41.23 | 23.29 | 17.58 | 16.88 | 17.24 | 15.03 | 6.51 | 22.37
(2, 4, 200) | True1 | 14.37 | 12.23 | 8.63 | 10.99 | 10.28 | 8.88 | 4.00 | 16.99
(0.2, 0.1, 0.5) | IPIC | 31.97 | 13.01 | 10.55 | 10.79 | 11.84 | 9.62 | 5.24 | 16.75
both | True2 | 11.30 | 11.61 | 8.29 | 7.92 | 6.39 | 5.52 | -2.10 | 16.97
DRIC | 14.22 | 12.42 | 10.43 | 10.79 | 9.69 | 10.29 | 4.44 | 14.92
(2, 4, 200) | True1 | 137.77 | 80.40 | 35.99 | 26.20 | 24.13 | 20.48 | -38.72 | 94.26
(0.2, 0.2, 1.0) | IPIC | 167.58 | -31.90 | 15.08 | 8.56 | 9.21 | 12.73 | -6.21 | 31.79
both | True2 | 133.76 | 173.72 | 16.90 | -61.39 | 28.37 | -4.32 | -139.65 | 236.77
DRIC | 39.05 | 18.78 | 12.71 | 14.87 | 12.29 | 11.78 | 1.82 | 19.01

Table 7 compares the proposed criteria IPIC and DRIC with the existing
criterion QICw. The mean squared error of $\hat{\bm{\beta}}$ is the main
indicator of the comparison, and the number of elements in the active set,
$\hat{p}=|\hat{\mathcal{J}}_{\lambda}^{(2)}|$, is also considered to be an
important indicator. In this table, the number of Monte Carlo iterations is
set to 200. Looking at the main indicators, IPIC outperforms QICw and DRIC
outperforms IPIC in all settings. As expected, as $H$ increases, such as $H=6$
and $H=8$, IPIC improves QICw significantly. The difference between IPIC and
DRIC is clearer compared with what is shown in Section 5.2. When there is a
model misspecification, the difference between the proposed and existing
criteria appears to widen, while the difference between IPIC and DRIC is not
as large as intended. Looking at $\hat{p}$, the proposed criteria narrow down
the explanatory variables considerably as expected, suggesting that the
estimation results vary substantially among these criteria.

Table 7: Performance comparison of QICw, IPIC and DRIC in causal inference models with Gaussian noise. The upper entry of the parameters is $(p,H,N)$, the middle entry is $(\beta^{*},\alpha^{*},\gamma^{*})$, and the lower entry represents the misspecified model. The $\hat{p}_{1}$ and $\sqrt{{\rm MSE}_{1}}$ columns list the dimension of selection and the square root of the mean squared error ($\times 10$) for nonzero parameters, whereas the $\hat{p}_{2}$ and $\sqrt{{\rm MSE}_{2}}$ columns show those for zero parameters and the $\hat{p}$ and $\sqrt{{\rm MSE}}$ columns show those for all parameters. parameters |  | $|\hat{p}_{1}|$ | $\sqrt{{\rm MSE}_{1}}$ | $\hat{p}_{2}$ | $\sqrt{{\rm MSE}_{2}}$ | $\hat{p}$ | $\sqrt{{\rm MSE}}$
---|---|---|---|---|---|---|---
(2, 4, 200) | QICw | 3.5 [0.7] | 4.02 [1.62] | 2.9 [1.1] | 3.66 [1.81] | 6.4 [1.5] | 5.71 [1.72]
(0.2, 0.2, 1.0) | IPIC | 2.4 [1.3] | 4.35 [1.47] | 1.3 [1.3] | 2.11 [2.23] | 3.8 [2.4] | 5.30 [1.55]
none | DRIC | 2.8 [1.2] | 3.78 [1.38] | 1.4 [1.4] | 1.58 [1.76] | 4.3 [2.4] | 4.46 [1.38]
(2, 4, 200) | QICw | 3.9 [0.3] | 4.12 [1.70] | 3.1 [1.0] | 3.82 [1.78] | 7.0 [1.1] | 5.87 [1.76]
(0.4, 0.2, 1.0) | IPIC | 3.6 [0.7] | 4.51 [1.91] | 1.9 [1.3] | 2.85 [2.16] | 5.6 [1.6] | 5.76 [1.87]
none | DRIC | 3.7 [0.6] | 4.57 [1.74] | 1.9 [1.4] | 2.16 [1.83] | 5.7 [1.8] | 5.42 [1.61]
(2, 4, 200) | QICw | 3.6 [0.6] | 3.93 [1.45] | 2.8 [1.0] | 3.34 [1.55] | 6.5 [1.3] | 5.37 [1.54]
(0.2, 0.1, 1.0) | IPIC | 2.5 [1.4] | 4.37 [1.43] | 1.3 [1.3] | 1.98 [1.97] | 3.9 [2.4] | 5.18 [1.44]
none | DRIC | 2.8 [1.2] | 3.80 [1.32] | 1.4 [1.4] | 1.63 [1.64] | 4.3 [2.3] | 4.46 [1.27]
(2, 4, 200) | QICw | 3.6 [0.6] | 3.14 [1.23] | 2.9 [1.1] | 2.82 [1.32] | 6.5 [1.4] | 4.42 [1.27]
(0.2, 0.2, 0.5) | IPIC | 3.0 [1.1] | 3.42 [1.33] | 1.6 [1.3] | 1.92 [1.60] | 4.7 [2.0] | 4.23 [1.31]
none | DRIC | 3.0 [1.1] | 3.45 [1.31] | 1.5 [1.3] | 1.60 [1.63] | 4.5 [2.1] | 4.12 [1.37]
(2, 4, 400) | QICw | 3.7 [0.5] | 2.82 [1.08] | 2.9 [1.1] | 2.58 [1.27] | 6.7 [1.3] | 4.01 [1.16]
(0.2, 0.2, 1.0) | IPIC | 3.1 [1.0] | 3.15 [1.29] | 1.6 [1.3] | 1.74 [1.54] | 4.8 [2.1] | 3.93 [1.22]
none | DRIC | 3.4 [0.8] | 2.70 [0.95] | 1.8 [1.4] | 1.43 [1.12] | 5.3 [1.9] | 3.28 [0.83]
(2, 6, 200) | QICw | 5.3 [0.8] | 6.32 [2.27] | 4.9 [1.2] | 5.93 [2.46] | 10.2 [1.6] | 8.91 [2.62]
(0.2, 0.2, 1.0) | IPIC | 2.8 [2.1] | 6.18 [1.58] | 1.9 [2.0] | 3.06 [3.25] | 4.8 [3.8] | 7.45 [2.25]
none | DRIC | 2.9 [2.0] | 5.34 [1.51] | 1.7 [1.9] | 2.01 [2.27] | 4.8 [3.7] | 6.12 [1.58]
(2, 8, 200) | QICw | 7.3 [0.9] | 8.54 [2.41] | 6.9 [1.1] | 8.14 [2.53] | 14.2 [1.6] | 12.02 [2.60]
(0.2, 0.2, 1.0) | AIC | 3.7 [2.6] | 7.95 [1.58] | 2.7 [2.4] | 3.97 [3.58] | 6.5 [4.8] | 9.47 [2.14]
none | DRIC | 3.8 [2.6] | 7.38 [1.57] | 2.4 [2.5] | 2.52 [2.54] | 6.4 [4.9] | 8.20 [1.59]
(4, 4, 200) | QICw | 7.1 [0.9] | 6.02 [1.58] | 6.2 [1.5] | 5.42 [1.64] | 13.4 [2.0] | 8.25 [1.63]
(0.2, 0.2, 1.0) | IPIC | 5.5 [2.2] | 6.20 [1.54] | 3.6 [2.3] | 3.64 [2.31] | 9.2 [4.1] | 7.55 [1.50]
none | DRIC | 5.7 [2.1] | 5.46 [1.45] | 3.2 [2.3] | 2.59 [1.84] | 9.1 [4.0] | 6.36 [1.27]
(2, 4, 200) | QICw | 3.5 [0.7] | 4.12 [1.52] | 2.9 [1.1] | 3.53 [1.84] | 6.5 [1.5] | 5.64 [1.82]
(0.2, 0.1, 1.0) | IPIC | 2.4 [1.3] | 4.48 [1.51] | 1.4 [1.3] | 2.10 [2.08] | 3.9 [2.3] | 5.33 [1.65]
treatment | DRIC | 2.8 [1.2] | 3.86 [1.26] | 1.5 [1.4] | 1.62 [1.50] | 4.3 [2.3] | 4.46 [1.19]
(2, 4, 200) | QICw | 3.5 [0.7] | 4.35 [1.82] | 2.9 [1.0] | 4.06 [2.15] | 6.4 [1.5] | 6.22 [2.16]
(0.2, 0.2, 1.0) | IPIC | 2.7 [1.3] | 4.38 [1.73] | 1.7 [1.4] | 2.85 [2.54] | 4.5 [2.4] | 5.70 [2.06]
treatment | DRIC | 2.6 [1.3] | 4.00 [1.39] | 1.4 [1.4] | 1.69 [1.71] | 4.1 [2.4] | 4.69 [1.30]
(2, 4, 200) | QICw | 3.6 [0.6] | 4.04 [1.53] | 3.0 [1.0] | 3.91 [2.06] | 6.7 [1.3] | 5.83 [2.03]
(0.2, 0.2, 0.5) | IPIC | 2.4 [1.4] | 4.45 [1.51] | 1.4 [1.3] | 2.29 [2.42] | 3.8 [2.5] | 5.48 [1.78]
outcome | DRIC | 2.2 [1.4] | 4.58 [1.43] | 1.1 [1.3] | 1.86 [2.47] | 3.5 [2.5] | 5.39 [1.88]
(2, 4, 200) | QICw | 3.6 [0.7] | 6.35 [2.66] | 3.2 [0.9] | 6.41 [3.72] | 6.9 [1.2] | 9.37 [3.82]
(0.2, 0.2, 1.0) | IPIC | 1.5 [1.5] | 6.05 [1.86] | 1.0 [1.3] | 2.75 [4.09] | 2.4 [2.6] | 7.34 [3.21]
outcome | DRIC | 1.4 [1.3] | 5.88 [1.62] | 0.8 [1.2] | 2.27 [4.12] | 2.3 [2.4] | 6.93 [3.37]

###  A.3 SURE theory for elastic net

Here, we extend SURE theory so that it can be applied to the elastic net
proposed in Zou and Hastie (2005). Specifically, we deal with the estimator
$\hat{\bm{\theta}}_{\lambda^{(1)},\lambda^{(2)}}$ with two regularization
parameters, defined in (30), multiplied by some constant. The following is
essentialy the same procedure as in Zou et al. (2007), and where the content
is actually almost identical, only the facts will be mentioned. For a while,
we will fix the outcome variable $\bm{y}=(y_{1},y_{2},\ldots,y_{N})$. Let
$\mathcal{J}_{\lambda^{(1)},\lambda^{(2)}}$ be the active set of $j$ such that
$\hat{\theta}_{\lambda^{(1)},\lambda^{(2)},j}$ is not 0, let
$s_{\lambda^{(1)},\lambda^{(2)},j}$ be the sign of
$\hat{\theta}_{\lambda^{(1)},\lambda^{(2)},j}$ if
$j\in\mathcal{J}_{\lambda^{(1)},\lambda^{(2)}}$ and $0$ if
$j\notin\mathcal{J}_{\lambda^{(1)},\lambda^{(2)}}$, and let
$\bm{s}_{\lambda^{(1)},\lambda^{(2)}}$ be the $p$-dimensional vector with
making $s_{\lambda^{(1)},\lambda^{(2)},j}$ the $j$-th component. Then, the
same type of equation as in (7) is established using
$\hat{\bm{\theta}}_{\lambda^{(1)},\lambda^{(2)}}$ and
$\bm{s}_{\lambda^{(1)},\lambda^{(2)}}$. Now, let us consider a partition of
the domain as follows:

For each mutually exclusive subdomain $\Lambda_{k}$ of
$(\lambda^{(1)},\lambda^{(2)})$ where
$\bigcup_{k}\Lambda_{k}=\mathbb{R}^{2}_{+}$, the active set
$\mathcal{J}_{k}=\mathcal{J}_{\lambda^{(1)},\lambda^{(2)}}$ and the sign
function $\bm{s}_{\mathcal{J}_{k}}=\bm{s}_{\lambda^{(1)},\lambda^{(2)}}$ are
constant.

The points on the boundary $\bigcup_{k}\partial\Lambda_{k}$ of these
subdomains are called transition points. In view of the continuous variation
of $\hat{\bm{\theta}}_{\lambda^{(1)},\lambda^{(2)}}$ with respect to
$(\lambda^{(1)},\lambda^{(2)})$, we can suppose that in the neighboring
$\Lambda_{k}$ and $\Lambda_{h}$, one of the active sets has one element larger
than the other. Without loss of generality, let $\Lambda_{k}$ be the larger,
and let the explanatory variable for the larger portion have the subscript
$j^{*}$ in $\mathcal{J}_{k}$. In this case,
$(\hat{\bm{\theta}}_{\mathcal{J}_{k}})_{j^{*}}=0$ on the boundary
$\partial\Lambda_{k}\cap\partial\Lambda_{h}$. In an expression similar to (7)
above, taking the limit
$(\lambda^{(1)},\lambda^{(2)})\to\partial\Lambda_{k}\cap\partial\Lambda_{h}$
in $\Lambda_{k}$. Then, from
$(\hat{\bm{\theta}}_{\mathcal{J}_{k}})_{j^{*}}=0$, we can see that a linear
equation with respect to $(\lambda^{(1)},\lambda^{(2)})$ holds on the
boundary. This is summarized as follows.

######  Lamma 2.

If $\Lambda_{k}$ and $\Lambda_{h}$ are neighboring regions, then the points in
$\partial\Lambda_{k}\cap\partial\Lambda_{h}$ satisfy a linear equation with
respect to $(\lambda^{(1)},\lambda^{(2)})$.

From now on, we will move $\bm{y}$. In addition, let us suppose that
$(\lambda^{(1)},\lambda^{(2)})$, which was originally considered fixed, is a
transition point. That is, we suppose that there are $k$ and $h$ such that
$(\lambda^{(1)},\lambda^{(2)})\in\partial\Lambda_{k}\cup\partial\Lambda_{h}$.
In this case, $\bm{y}$ satisfies the equation in Lemma 2. In other words, only
when $\bm{y}$ satisfies the equation, $(\lambda^{(1)},\lambda^{(2)})$ can be a
transition point on $\partial\Lambda_{k}\cup\partial\Lambda_{h}$. We can
easily check that the set of $\bm{y}$ satisfying the equation is of measure 0.
Here, considering all possible combinations of $\mathcal{J}_{k}$,
$\mathcal{J}_{h}$, $j^{*}$, $\bm{s}_{\mathcal{J}_{k}}$ and
$\bm{s}_{\mathcal{J}_{h}}$, we define
$\mathcal{N}_{\lambda^{(1)},\lambda^{(2)}}$ as the set of all $\bm{y}$ such
that $(\lambda^{(1)},\lambda^{(2)})$ is a transition point somewhere. Note
that since the number of combinations is finite,
$\mathcal{N}_{\lambda^{(1)},\lambda^{(2)}}$ is of measure 0. This fact is
summarized as follows.

######  Lamma 3.

For any $(\lambda^{(1)},\lambda^{(2)})\in\mathbb{R}^{2}_{+}$, there exists a
set ${\cal N}_{\lambda^{(1)},\lambda^{(2)}}$ of measure $0$ in
$\mathbb{R}^{N}$ such that $(\lambda^{(1)},\lambda^{(2)})$ is not a transition
point if
$\bm{y}\in\mathcal{G}_{\lambda^{(1)},\lambda^{(2)}}\equiv\mathbb{R}^{N}\setminus{\cal
N}_{\lambda^{(1)},\lambda^{(2)}}$.

Using the Bolzano-Weierstrass theorem, we obtain the following.

######  Lamma 4.

For any $(\lambda^{(1)},\lambda^{(2)})\in\mathbb{R}^{2}_{+}$,
$\hat{\bm{\theta}}_{\lambda^{(1)},\lambda^{(2)}}$ is continuous with respect
to $\bm{y}$.

In the following, we will move $\bm{y}$ just locally. Supposing a sequence
$\\{\bm{y}_{t}\\}$ which converges to a fixed point $\bm{y}^{*}$, we advance
the discussion using Lemma 4 and the Karush-Kuhn-Tucker condition and obtain
the following.

######  Lamma 5.

For any $(\lambda^{(1)},\lambda^{(2)})\in\mathbb{R}^{2}_{+}$ and
$\bm{y}^{*}\in\mathcal{G}_{\lambda^{(1)},\lambda^{(2)}}$, the active set
$\mathcal{J}_{\lambda^{(1)},\lambda^{(2)}}$ and the sign function
$\bm{s}_{\lambda^{(1)},\lambda^{(2)}}$ are locally constant with respect to
$\bm{y}$ at $\bm{y}^{*}$.

Let $\bm{y}\in\mathcal{G}_{\lambda^{(1)},\lambda^{(2)}}$. From Lemma 5, there
exists a sufficiently small positive constant $\varepsilon$ such that the
active set and the sign function for $\bm{y}+\Delta\bm{y}$ do not change if
$\|\Delta\bm{y}\|<\varepsilon$. That is,
$\hat{\bm{\theta}}_{\lambda^{(1)},\lambda^{(2)}}$ is uniform Lipschitz
continuous for $\bm{y}\in\mathcal{G}_{\lambda^{(1)},\lambda^{(2)}}$, and
furthermore, from Lemma 3, $\hat{\bm{\theta}}_{\lambda^{(1)},\lambda^{(2)}}$
is almost differentiable for $\bm{y}\in\bm{R}^{N}$. Thus, Stein’s lemma can be
used, and we can say that SURE theory is valid.