A Bayesian semiparametric Archimedean copula

TL;DR

This paper introduces a Bayesian semiparametric approach to Archimedean copulas using a spline-based generator, enabling flexible modeling of dependence structures and hypothesis testing in a survival analysis framework.

Contribution

It proposes a novel spline-based semiparametric generator for Archimedean copulas, with Bayesian inference and independence testing capabilities.

Findings

The model covers the full range of Kendall's tau.

Simulation studies validate the model's flexibility and accuracy.

Application to real data demonstrates practical utility.

Abstract

An Archimedean copula is characterised by its generator. This is a real function whose inverse behaves as a survival function. We propose a semiparametric generator based on a quadratic spline. This is achieved by modelling the first derivative of a hazard rate function, in a survival analysis context, as a piecewise constant function. Convexity of our semiparametric generator is obtained by imposing some simple constraints. The induced semiparametric Archimedean copula produces Kendall's tau association measure that covers the whole range $(-1,1)$. Inference on the model is done under a Bayesian approach and for some prior specifications we are able to perform an independence test. Properties of the model are illustrated with a simulation study as well as with a real dataset.

Figures

Tables

Table 1. Table 1: Summary of some parametric Archimedean copulas parametrised such that h ​ ( 0 ) = 1 ℎ 0 1 h(0)=1 , for the first three copulas, and h ​ ( ϵ ) = 1 ℎ italic-ϵ 1 h(\epsilon)=1 , for the Gumbel copula.

Copula	$\Theta$	$\phi \left(t\right)$	${\phi }^{-1}\left(t\right)$	Strict?	h(t)
Product	-	$-\log \left(t\right)$	$e^{-t}$	Yes	$1$
Clayton	$[-1,\infty )$	$\frac{1}{\theta }\left(t^{-\theta }-1\right)$	${\left(1+\theta t\right)}^{-1/\theta }$	If $\theta \ge 0$	$\frac{1}{1+\theta t}$
AMH	$[-1,1)$	$\frac{1}{1-\theta }\log \left(\frac{1-\theta +\theta t}{t}\right)$	$\frac{1-\theta }{e^{\left(1-\theta \right)t}-\theta }$	Yes	$\frac{\left(1-\theta \right)e^{\left(1-\theta \right)t}}{e^{\left(1-\theta \right)t}-\theta }$
Gumbel	$[1,\infty )$	$ϵ{\left(\frac{-\log \left(t\right)}{ϵ\theta }\right)}^{\theta }$	$\exp \left\{-ϵ\theta {\left(\frac{t}{ϵ}\right)}^{1/\theta }\right\}$	Yes	${\left(\frac{ϵ}{t}\right)}^{1-1/\theta }$

Table 2. Table 2: GOF measures, obtained with different partition definitions, for a simulated dataset of size n = 200 𝑛 200 n=200 from the product copula.

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	${\kappa }_{\tau }$	Sup.Norm	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$	$\text{Norm}$

Table 3. Table 3: GOF measures, obtained with different partition definitions, for a simulated dataset of size n = 200 𝑛 200 n=200 from the Clayton copula with θ = − 0.8 𝜃 0.8 \theta=-0.8 .

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	${\kappa }_{\tau }$	Sup.Norm	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$	$\text{Norm}$

Table 4. Table 4: GOF measures, obtained with different partition definitions, for a simulated dataset of size n = 200 𝑛 200 n=200 from the Clayton copula with θ = 1 𝜃 1 \theta=1 .

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	${\kappa }_{\tau }$	Sup.Norm	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$	$\text{Norm}$

Table 5. Table 5: GOF measures, obtained with different partition definitions, for a simulated dataset of size n = 200 𝑛 200 n=200 from the AMH copula with θ = − 0.7 𝜃 0.7 \theta=-0.7 .

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	${\kappa }_{\tau }$	Sup.Norm	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$	$\text{Norm}$

Table 6. Table 6: GOF measures, obtained with different partition definitions, for a simulated dataset of size n = 200 𝑛 200 n=200 from the AMH copula with θ = 0.7 𝜃 0.7 \theta=0.7 .

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	${\kappa }_{\tau }$	Sup.Norm	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$	$\text{Norm}$

Table 7. Table 7: GOF measures, obtained with different partition definitions, for a simulated dataset of size n = 200 𝑛 200 n=200 from the Gumbel copula with θ = 1.4 𝜃 1.4 \theta=1.4 .

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	${\kappa }_{\tau }$	Sup.Norm	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$	$\text{Norm}$

Table 8. Table 8: GOF measures, obtained with different partition definitions, for the real data.

Part.Type	$K$	$Q_{{\kappa }_{\tau }}^{\left(0.025\right)}$	${\hat{\kappa }}_{\tau }$	$Q_{{\kappa }_{\tau }}^{\left(0.975\right)}$	Sample ${\tilde{\kappa }}_{\tau }$	LPML
\csvreader[late after line=
$\text{LPML}$	\Size	$\text{2}5Tau$	$\text{Mtau}$	$\text{9}75Tau$	$\text{RealTau}$

Table 9. Table 9: Independence test for all simulated and real datasets. semiparametric Archimedean copula test (SPAC), contingency table test (CT), mixture model test (MM), and empirical copula test (EC). For the Bayesian tests we report P ​ ( H 1 ∣ data ) P conditional subscript 𝐻 1 data \mbox{P}(H_{1}\mid\mbox{data}) and for the frequentist test we report p-values.

		Bayesian Tests			Freq.Test
Dataset	$\theta$	SPAC	CT	MM	EC
Product	–	0.19	0.04	0.12	0.15
Clayton	$-0.8$	1.00	1.00	0.99	0.00
Clayton	$-0.4$	0.96	0.65	0.83	0.00
Clayton	$0.6$	0.96	0.43	0.79	0.01
Clayton	$1.0$	1.00	0.87	0.75	0.00
AMH	$-0.7$	0.63	0.11	0.21	0.31
AMH	$-0.3$	0.23	0.05	0.08	0.79
AMH	$0.3$	0.22	0.06	0.09	0.81
AMH	$0.7$	0.86	0.38	0.48	0.33
Gumbel	$1.4$	1.00	0.50	0.75	0.00
Gumbel	$2.0$	1.00	0.97	0.93	0.00
Real data	–	0.47	0.07	0.40	0.35

Equations

$$C(u,v)={\phi }^{-1}\left(\phi (u)+\phi (v)\right).$$

$${\kappa }_{\tau }=1+4{\int }_0^1\frac{\phi (t)}{{\phi }^{\prime }(t^{+})}\text{mbox}dt,$$

$$h^{\prime }(t)=\sum _{k=1}^K{\theta }_{k}I({\tau }_{k-1}<t\le {\tau }_k),$$

$$h(t)=\sum _{k=1}^K\left(A_k+{\theta }_{k}t\right)I({\tau }_{k-1}<t\le {\tau }_k),$$

$$H(t)=\sum _{k=1}^K\left(B_k+A_{k}t+\frac{{\theta }_k}{2}{t}^2\right)I({\tau }_{k-1}<t\le {\tau }_k),$$

$${\phi }^{-1}(t)=\exp \{-H(t)\},$$

$$\phi (t)=\sum _{k=1}^K$$

$$-\frac{B_k+\log (t)}{A_k}I({\theta }_k=0))I\left({\phi }^{-1}({\tau }_k)\le t<{\phi }^{-1}({\tau }_{k-1})\right).$$

$${\kappa }_{\tau }=-1+2\sum _{k=1}^K{A}_k{\int }_{{\tau }_{k-1}}^{{\tau }_k}\exp \left(-2B_k-2A_{k}t-{\theta }_k{t}^2\right)\text{mbox}dt.$$

$${\phi }^{-1}(t)=\exp \left\{-\sum _{k=1}^K{\theta }_k{w}_k(t)\right\},$$

$$w_{k}(t)=\left\{\begin{array}[]{ll}\Delta_{k},&t>\tau_{k}\\ t-\tau_{k-1}&t\in(\tau_{k-1},\tau_{k}]\\ 0&\mbox{otherwise}\end{array}\right.$$

$$\phi (t)=\sum _{k=1}^K\left\{{\tau }_{k-1}-\frac{1}{{\theta }_k}(\log t+{\vartheta }_{k-1})\right\}I({\vartheta }_{k-1}<-\log t\le {\vartheta }_k),$$

$${\kappa }_{\tau }=1+\sum _{k=1}^K\left\{e^{-2{\vartheta }_k}(1+2{\theta }_k{\tau }_k)-e^{-2{\vartheta }_{k-1}}(1+2{\theta }_k{\tau }_{k-1})\right\}.$$

$$f_C(u,v)={\phi }^{-1(\prime \prime )}\left(\phi (u)+\phi (v)\right){\phi }^{(\prime )}(u){\phi }^{(\prime )}(v),$$

$${\phi }^{-1(\prime \prime )}(t)=\sum _{k=1}^K\left\{(A_k+{\theta }_{k}t{)}^2-{\theta }_k\right\}\exp \left\{-\left(B_k+A_{k}t+\frac{{\theta }_k}{2}{t}^2\right)\right\}I({\tau }_{k-1}<t\le {\tau }_K)$$

$${\phi }^{(\prime )}(t)=-\sum _{k=1}^K\frac{1}{t}{\left(-2{\theta }_k{B}_k+A_k^2-2{\theta }_k\log (t)\right)}^{-1/2}I\left({\phi }^{-1}({\tau }_k)\le t<{\phi }^{-1}({\tau }_{k-1})\right).$$

$$f(\mathbf{\theta })\propto \prod _{k=1}^{K-1}\left\{{\pi }_{0}I({\theta }_k=0)+(1-{\pi }_0)\text{mbox}N({\theta }_k\mid {\mu }_0,{\sigma }_0^2)\right\}I(\mathbf{\theta }\in \Theta ).$$

$$f({\theta }_k\mid {\mathbf{\theta }}_{-k},\text{mbox}data)\propto \text{mbox}lik(\mathbf{\theta }\mid \mathbf{u},\mathbf{v})f(\mathbf{\theta }),$$

$$q({\theta }_k\mid {\mathbf{\theta }}_{-k},{\theta }_k^{(r)})={\pi }_{1}I({\theta }_k=0)+(1-{\pi }_1)\text{mbox}Un({\theta }_k\mid \max \{a_k,{\theta }_k^{(r)}-\delta c_k\},\min \{b_k,{\theta }_k^{(r)}+\delta c_k\})$$

$$p\left({\theta }_k^{\ast },{\theta }_k^{(r)}\right)=\min \left\{1,\frac{f({\theta }_k^{\ast }\mid {\mathbf{\theta }}_{-k},\text{mbox}data)q({\theta }_k^{(r)}\mid {\mathbf{\theta }}_{-k},{\theta }_k^{\ast })}{f({\theta }_k^{(r)}\mid {\mathbf{\theta }}_{-k},\text{mbox}data)q({\theta }_k^{\ast }\mid {\mathbf{\theta }}_{-k},{\theta }_k^{(r)})}\right\}.$$

$${\theta }_k\ge \underset{k\le j\le K-1}{\max }\left\{-\left({\theta }_0+\sum _{i=1,i\ne k}^j({\tau }_i-{\tau }_{i-1}){\theta }_i\right)/({\tau }_k-{\tau }_{k-1})\right\},$$

$${\theta }_k<{\left({\theta }_0+\sum _{i=1}^{k-1}({\tau }_i-{\tau }_{i-1}){\theta }_i\right)}^2.$$

$${\theta }_k>\left(\sqrt{{\theta }_j}-{\theta }_0-\sum _{i=1,i\ne k}^{j-1}({\tau }_i-{\tau }_{i-1}){\theta }_i\right)/({\tau }_k-{\tau }_{k-1})$$

Full text

Keywords: Archimedean copula, Bayes nonparametrics, piecewise constant, survival analysis, quadratic spline.

AMS Classification: 60E05 $\cdot$ 62G05 $\cdot$ 62N86.

1 Introduction

Let $\varphi(\cdot)$ be a continuous, strictly decreasing function from $[0,1]$ to $[0,\infty)$ such that $\varphi(1)=0$. Let $\varphi^{-1}(\cdot)$ be the inverse or the pseudo-inverse of $\varphi$, where the latter is defined as zero for $t>\varphi(0)$. If $\varphi(t)\to\infty$ as $t\to 0^{+}$ the generator is called strict. For instance, $\varphi(t)=-\log(t)$, is an example of a strict generator. An Archimedean copula $C(u,v)$ with generator $\varphi$ is a function from $[0,1]^{2}$ to $[0,1]$ defined as

[TABLE]

A further requirement for (1) to be well defined is that $\varphi$ must be convex (e.g. Nelsen,, 2006).

There are many properties that characterize Archimedean copulas, for instance, they are symmetric, associative and their diagonal section $C(u,u)$ is always less than $u$ for all $u\in(0,1)$. Generators $\varphi(\cdot)$ are usually parametric families defined by a single parameter. Most of them are summarised in (Nelsen,, 2006, Table 4.1) and few of them are also included in Table 1.

Association measures induced by Archimedean copulas are a function of the generator. For instance, Kendall’s tau becomes

[TABLE]

where $\varphi^{\prime}(t^{+})$ denotes the right derivative of $\varphi$ at $t$.

In this work we propose a Bayesian semiparametric generator defined through a quadratic spline. Within a survival analysis context, we model the first derivative of a hazard rate function with a piecewise constant function. The hazard rate and the cumulative hazard functions become linear and quadratic continuous functions, respectively. The induced survival function is used as an inverse generator for an Archimedean copula. Convexity constraints are properly addressed and inference on the model is done under a Bayesian approach.

Other studies on semiparametric generators for Archimedean copulas can be found in Genest and Rivest, (1993) where their model is based on an empirical Kendall’s process. A new approach and extensions of this latter methodology can be found in Genest et al., (2011). In Guillote and Perron, (2015) the model arises from the one-to-one correspondence between an Archimedean generator and a distribution function of a nonnegative random variable. In particular they use a mixture of Pólya trees as a prior for the corresponding distribution function under a Bayesian nonparametric approach. In a work more related to ours, Vandenhende and Lambert, (2005) use the relationship between quantile functions and Archimedean generators to define a semiparametric generator by supplementing a parametric generator with $n+1$ dependence parameters. Differing to their work, our model is not based on any parametric generator and the Kendall’s tau can take values on the whole interval $(-1,1)$.

The contents of the rest of the paper is as follows. In Section 2 we present our proposal and characterise its properties. In Section 3 we provide details of how to make posterior inference under a Bayesian approach. In Section 4 we illustrate the performance of our model with a simulation study as well as with a real data set. We conclude with some remarks in Section 5.

Before proceeding we introduce notation: $\mbox{Un}(a,b)$ denotes a continuous uniform density on the interval $(a,b)$; and $\mbox{N}(\mu,\sigma^{2})$ denotes a normal density with mean $\mu$ and variance $\sigma^{2}$.

2 Model

To elicit our proposal we noticed that $\varphi^{-1}$ is a decreasing function from $[0,\infty)$ to $[0,1]$, so it behaves as a survival function, in a failure time data analysis context (e.g. Klein and Moeschberger,, 2003). The idea is to propose a semi/non parametric form for the inverse generator $\varphi^{-1}$ by using survival analysis ideas. For that we recall some basic definitions.

Let $h(t)$ be a nonnegative function with domain in $[0,\infty)$ such that $H(t)=\int_{0}^{t}h(s)\mbox{d}s$ satisfies $H(t)\to\infty$ as $t\to\infty$. Then $S(t)=\exp\{-H(t)\}$ is a decreasing function from $[0,\infty)$ to $[0,1]$, so it behaves like an inverse generator $\varphi^{-1}(t)$. In a survival analysis context, functions $h(\cdot)$, $H(\cdot)$ and $S(\cdot)$ are the hazard rate, cumulative hazard and survival functions, respectively.

In particular, if $h(t)=\theta$, i.e. constant for all $t$, then $S(t)=e^{-\theta t}$. If we take $\varphi(t)^{-1}=e^{-\theta t}$, then $\varphi(t)=-(\log t)/\theta$. Using (1) we obtain that the resulting copula $C(u,v)=uv$ is the independence copula, and what is interesting, is that it does not depend on $\theta$.

2.1 Main proposal

Using the previous ideas we construct a semiparametric generator in the following way. We first consider a partition of size $K$ of the positive real line, with interval limits given by $0=\tau_{0}<\tau_{1}<\cdots<\tau_{K}=\infty$. Then, we define the first derivative of the hazard rate, as a piecewise constant function of the form

[TABLE]

where $\theta_{K}\equiv 0$. We recover the hazard rate function as $h(t)=\int_{0}^{t}h^{\prime}(s)\mbox{d}s+\theta_{0}$, where $h(0)=\theta_{0}>0$ is an initial condition. Using (3), the hazard rate becomes a piecewise linear function of the form

[TABLE]

where $A_{1}=\theta_{0}$ and $A_{k}=\theta_{0}+\sum_{j=1}^{k-1}(\theta_{j}-\theta_{j+1})\tau_{j}$, for $k=2,\ldots,K$.

Integrating now the hazard function (4), the cumulative hazard is a piecewise quadratic function given by

[TABLE]

where $B_{1}=0$ and $B_{k}=\sum_{j=2}^{k}(\theta_{j}-\theta_{j-1})\tau_{j-1}^{2}/2$, for $k=2,\ldots,K$.

We therefore define a semiparametric inverse generator as the induced survival function, which can be written as

[TABLE]

where $H(t)$ is given in (5). We now study some properties of this inverse generator.

Proposition 1

Consider the semiparametric inverse generator $\varphi^{-1}(t)$, given in (6), and assume that $\{\theta_{k},k=0,1,\ldots,K\}$ are such that $\theta_{0}>0$, $\theta_{K}=0$ and satisfy conditions (C1) and (C2) given by

1. (C1)

$A_{k}+\theta_{k}t>0$, for $t\in(\tau_{k-1},\tau_{k}]$ and for all $k=1,\ldots,K$. 2. (C2)

$(A_{k}+\theta_{k}t)^{2}>\theta_{k}$, for $t\in(\tau_{k-1},\tau_{k}]$ and for all $k=1,\ldots,K$.

Then,

1. (i)

$\varphi^{-1}(t)$* is a continuous and injective function of $t$,* 2. (ii)

$\varphi^{-1}(t)$* is a convex function,* 3. (iii)

$\varphi^{-1}(t)$* induces a strict generator.*

Proof For (i) we know that $h^{\prime}(t)$, as in (3), is a piecewise constant discontinuous function, however, function $h(t)$, as in (4), is continuous. To see this, for each $k=1,\ldots,K$, the limit from the left is $\lim_{t\to\tau_{k}^{-}}h(t)=\lim_{t\to\tau_{k}^{-}}A_{k}+\theta_{k}t=A_{k}+\theta_{k}\tau_{k}$, and the limit from the right becomes $\lim_{t\to\tau_{k}^{+}}h(t)=\lim_{t\to\tau_{k}^{+}}A_{k+1}+\theta_{k+1}t=A_{k+1}+\theta_{k+1}\tau_{k}$. Since $A_{k+1}=A_{k}+(\theta_{k}-\theta_{k+1})\tau_{k}$, then both limits coincide. For the second part of (i), we have that $H(t)$ is a monotonous function whose derivative is strictly positive, due to condition (C1), therefore $H(t)$ is injective and invertible on its image (Rudin,, 1987). For (ii) we take the second derivative of $\varphi^{-1}(t)$ which becomes $\varphi^{-1(\prime\prime)}(t)=\{h(t)\}^{2}\exp\{-H(t)\}-h^{\prime}(t)\exp\{-H(t)\}$, this is positive if and only if $\{h(t)\}^{2}-h^{\prime}(t)>0$. For this to happen we require condition $(C2)$. For (iii), $\varphi^{-1}(t)$ must be a proper survival function, that is, $h(t)$ must be nonnegative, which is achieved by imposing condition $(C1)$. Furthermore, we need $\lim_{t\to\infty}\varphi^{-1}(t)=0$, which is equivalent to prove that $\lim_{t\to\infty}H(t)=\lim_{t\to\infty}\left(B_{K}+A_{K}t+\theta_{K}t^{2}/2\right)=\infty$. This is true since $B_{K}$ is a finite constant, by definition $\theta_{K}=0$, and this together with $(C2)$ imply $A_{K}>0$, so the linear part goes to infinity when $t\to\infty$. $\diamond$

By property (i) in Proposition 1, we can invert equation (6) to obtain an expression for the generator. This is given by

[TABLE]

The value $K$ controls the flexibility of the generator, and thus of the copula. If $K=1$, the induced Archimedean copula is the independence copula, whereas for larger $K$, the generator, and the induced copula, become semiparametric. Potentially $K$ could be infinite implying a nonparametric model. We now discuss some association properties of our semiparametric generator.

To see the kind of association induced by our proposal, we computed the Kendall’s tau using expression (2) with generator (7). This is given in the following result.

Proposition 2

The Kendall’s tau obtained by the Archimedean copula with semiparametric generator (7) is given by

[TABLE]

Moreover, this $\kappa_{\tau}\in(-1,1)$.

Proof Rewriting expression (2) in terms of the inversed generator we obtain $\kappa_{\tau}=1-4\int_{0}^{\infty}t\{\varphi^{-1(\prime)}(t)\}^{2}\mbox{d}t$. Computing the derivative we get $\varphi^{-1(\prime)}(t)=-\sum_{k=1}^{K}(A_{k}+\theta_{k}t)\times$ $\exp\{-(B_{k}+A_{k}t+\theta_{k}t/2)\}I(\tau_{k-1}<t\leq\tau_{k})$. Doing the integral we obtain the expression. To obtain the range of possible values of $\kappa_{\tau}$ it is easier to re-write $\kappa_{\tau}$ in terms of $h(t)$ and $H(t)$. This becomes $\kappa_{\tau}=-2\int_{0}^{\infty}th^{\prime}(t)\exp\{-H(t)\}\,dt$. Here it is straightforward to see that the sign of $\kappa_{\tau}$ is determined by the sign of $h^{\prime}(t)$, therefore $h^{\prime}(t)>0$ for all $t$ implies $-1<\kappa_{\tau}<0$ and $h^{\prime}(t)\leq 0$ implies $0\leq\kappa_{\tau}<1$. $\diamond$

The expression for $\kappa_{\tau}$ tells us that the concordance induced by our semiparametric copula is a function of both, the parameters $\{\theta_{k}\}$, as well as of the partition limits $\{\tau_{k}\}$. It depends on a definite integral and can be evaluated numerically. What is more important is that $\kappa_{\tau}$ covers the whole range from $-1$ to $1$, showing that our proposal is very flexible.

To illustrate the flexibility of our model we define a partition of the positive real line of size $K=10$, such that $\tau_{k}=-\log(1-k/10)$ for $k=0,1,\ldots,10$. We consider two scenarios for the values of the parameters $\{\theta_{k}\}$. The first scenario is defined by $\theta_{k}<0$ for all $k<K$, whereas the second scenario contains $\theta_{k}>0$ for all $k<K$. Conditions $(C1)$ and $(C2)$ were satisfied in both cases. Figure 1 contains functions $h^{\prime}(t)$, $H(t)$ and $\varphi^{-1}(t)$ for two different scenarios, the solid (blue) line corresponds to the first scenario and the dotted (red) line to the second scenario. In the first case the corresponding hazard function (middle panel) is decreasing, whereas for the second case the hazard function is increasing. The induced concordance values are $\kappa_{\tau}=0.368$ and $\kappa_{\tau}=-0.202$, respectively.

As a second example, we consider a partition of size $K=50$, such that $\tau_{k}=-\log(1-k/50)$ for $k=0,1,\ldots,50$. We consider three different scenarios for the parameters $\{\theta^{(i)}_{k}\}$ with $i=1,2,3$, respectively. In the first scenario we assume $\theta_{1}^{(1)}\sim\mbox{Un}(-1,1)$, in the second $\theta_{1}^{(2)}\sim\mbox{Un}(-50,0)$ and in the third $\theta_{1}^{(3)}\sim\mbox{Un}(0,1)$. Posteriorly, we define sequentially $\theta_{k}^{(i)}\sim\mbox{Un}(a_{k}^{(i)},b_{k}^{(i)})$ with $a_{k}^{(i)}$ and $b_{k}^{(i)}$ constants such that constraints $(C1)$ and $(C2)$ are satisfied, for $k=2,\ldots,K-1$ and $i=1,2,3$. We repeated sampling from these distributions a total of 5,000 times, and for each repetition we computed $\kappa_{\tau}$. The induced histogram densities for the three scenarios are presented in Figure 2. For the first scenario, the values of $\kappa_{\tau}$ range from $-0.3$ to $0.4$, showing that our model can capture both negative and positive concordance measures. For the second scenario, the values of $\kappa_{\tau}$ are all positive and the distribution is right skewed, and for the third scenario the values of $\kappa_{\tau}$ are all negative showing a left skewed distribution.

According to Nelsen, (2006), new generators can be defined if we apply a scale transformation of the form $\phi^{-1}(t)=\varphi^{-1}(\beta t)$ if and only if $\phi(t)=\varphi(t)/\beta$, for $\beta>0$, where $\phi(t)$ becomes a new Archimedean copula generator. More recently, Di Bernardino and Rullière, (2013) realised that the new generator $\phi(t)$ induces exactly the same copula (1) as that obtained with $\varphi(t)$. To see this we have that $C_{\phi}(u,v)=\phi^{-1}(\phi(u)+\phi(v))=\varphi^{-1}\left(\beta\left\{\frac{1}{\beta}\varphi(u)+\frac{1}{\beta}\varphi(v)\right\}\right)=C_{\varphi}(u,v)$. In other words, an Archimedean copula generator is not unique.

Moreover, in terms of the hazard rate functions, $h_{\phi}(t)$ and $h_{\varphi}(t)$, induced by generators $\phi$ and $\varphi$, respectively, the relationship becomes $h_{\phi}(t)=\beta h_{\varphi}(\beta t)$. In order to make our semiparametric generator identifiable, without loss of generality, we impose the new constraint

1. (C3)

$\theta_{0}=1$.

This constraint is equivalent to imposing $h(0)=1$ in definition (4).

2.2 Alternative construction

Instead of starting with a piecewise function for the derivative of a hazard rate, we could start by defining a piecewise constant function for the hazard rate itself. That is $h(t)=\sum_{k=1}^{K}\theta_{k}I(\tau_{k-1}<t\leq\tau_{k})$ with $\theta_{k}>0$, and $\{\tau_{k}\}$ a partition of the positive real line. In this case the cumulative hazard function becomes $H(t)=\sum_{j=1}^{k}\theta_{j}\Delta_{j}+\theta_{k}(t-\tau_{k-1})$, for $t\in(\tau_{k-1},\tau_{k}]$, with $\Delta_{j}=\tau_{j}-\tau_{j-1}$. The inverse generator is then a linear spline of the form

[TABLE]

with

[TABLE]

and the corresponding Archimedean generator has the form

[TABLE]

with $\vartheta_{k}=\sum_{j=1}^{k}\theta_{j}\Delta_{j}$. To ensure convexity of the generator we further require $\theta_{1}\geq\theta_{2}\geq\cdots\geq\theta_{K}$. Furthermore, the Kendall’s tau has a simpler expression

[TABLE]

However, it can be shown that this expression for the Kendall’s tau only allows positive values, constraining the possible associations captured by the model. Therefore, in the remainder of the paper we will concentrate on our main proposal defined in Section 2.1.

3 Posterior inference

The copula density $f_{C}(u,v)$, of an Archimedean copula, can be obtained by taking the second crossed derivatives with respect to $u$ and $v$ in expression (1). In terms of the generator and its inverse this density becomes

[TABLE]

where the single and double primes denote first and second derivatives, respectively, and are given by

[TABLE]

and

[TABLE]

Let $(U_{i},V_{i})$, $i=1,\ldots,n$ be a bivariate sample of size $n$ from $f_{C}(u,v)$ defined in (8). With this we can construct the likelihood for $\boldsymbol{\theta}=(\theta_{0},\theta_{1},\ldots,\theta_{K})$ as $\mbox{lik}(\boldsymbol{\theta}\mid{\bf u},{\bf v})=\prod_{i=1}^{n}f_{C}(u_{i},v_{i}\mid\boldsymbol{\theta})$, where we have made explicit the dependence on $\boldsymbol{\theta}$ in the notation of the copula density. Recall that the parameter space $\Theta$ contains the values of $\boldsymbol{\theta}$ that satisfy several conditions, $(C1)$ and $(C2)$ given in Proposition 1, $(C3)$ to make our generator unique, and $\theta_{K}=0$.

We assume a joint prior distribution for $\boldsymbol{\theta}$ of the form

[TABLE]

Note that we explicitly allow the $\theta_{k}$’s, for $k=1,\ldots,K-1$ to be zero with positive probability $\pi_{0}$. This prior choice is useful to define an independence test. Specifically, we consider the hypotheses $H_{0}:U\mbox{ and }V$ independent, which is equivalent to $H_{0}:\theta_{1}=\cdots=\theta_{K-1}=0$, versus the alternative $H_{1}:U\mbox{ and }V$ dependent, which is equivalent to $H_{1}:\theta_{k}\neq 0$ for at least one $k=1,\ldots,K-1$. To perform the test we can compute the posterior probabilities of $H_{0}$ and $H_{1}$ and make the decision using decision theory (DeGroot,, 2004), or use the corresponding Bayes factor (Kass and Raftery,, 1995), which in case that $\mbox{P}(H_{0})=\mbox{P}(H_{1})$, this becomes the odds in favour of $H_{1}$, that is, $B_{10}=\mbox{P}(H_{1}\mid\mbox{data})/\mbox{P}(H_{0}\mid\mbox{data})$. Here we follow the approach of Filippi et al., (2016) and report $\mbox{P}(H_{1}\mid\mbox{data})$ as an evidence in favour of dependence.

The posterior distribution of $\boldsymbol{\theta}$ is simply given by the product of expressions (8) and (9), up to a proportionality constant. It is easier to characterize the posterior distribution by implementing a Gibbs sampler (Smith and Roberts,, 1993) and sampling from the conditional posterior distributions

[TABLE]

for $k=1,\ldots,K-1$. However, sampling from conditional distributions (10) is not trivial since the parameter $\theta_{k}$ appears everywhere in the likelihood, the parameter space is complex and no closed expression can be obtained for the normalising constant, we therefore propose a Metropolis-Hastings step (Tierney,, 1994) by sampling $\theta_{k}^{*}$ at iteration $(r+1)$ from a random walk proposal distribution

[TABLE]

where the interval $(a_{k},b_{k})$ represents the conditional support of $\theta_{k}$, $c_{k}=b_{k}-a_{k}$ is its length, with $a_{k}=\max_{k\leq j\leq K-1}\left\{\left(\sqrt{\theta_{j+1}}I(\theta_{j+1}\geq 0)-\theta_{0}-\sum_{i=1,i\neq k}^{j}(\tau_{i}-\tau_{i-1})\theta_{i}\right)/(\tau_{k}-\tau_{k-1})\right\}$, for $k=1,\ldots,K-1$, $b_{k}=\left(\theta_{0}+\sum_{j=1}^{k-1}(\tau_{j}-\tau_{j-1})\theta_{j}\right)^{2}$, for $k=2,\ldots,K-1$, and $b_{1}=1$. The justification of these bounds obeys the inclusion of constraints $(C1)$ and $(C2)$ and their derivations are given in Appendix Appendix. The parameters $\pi_{1}$ and $\delta$ are tuning parameters that control the acceptance rate.

Therefore, at iteration $r+1$ we accept $\theta_{k}^{*}$ with probability

[TABLE]

This Metropolis-Hastings within Gibbs procedure to obtain posterior inference of our model was implemented in Python and the code is available upon request from the first author.

To perform the independent test, posterior probability of $H_{0}$ can be approximated via Monte Carlo by using the MCMC posterior draws of the vector $\boldsymbol{\theta}$ and computing the relative frequency of the event $\theta_{k}=0$ for all $k=1,\ldots,K-1$, we therefore obtain posterior probability of $H_{1}$ by computing the complement.

4 Numerical studies

We illustrate the performance of our model in two ways, through a simulation study, and with a real data set.

To define the partition $\{\tau_{k}\}$ of the positive real line, inspired by the generator of the product copula, we consider a Log-$\alpha$ partition defined by $\tau_{k}=-\alpha\log(1-k/K)$ for $k=0,\ldots,K-1$, with $\alpha>0$. This partition is the result of transforming a uniform partition in the interval $[0,1]$ via a convex function. Larger values of $\alpha$ increase the spread of the partition along the positive real line.

4.1 Simulation study

We generated simulated data from four parametric Archimedean copulas, namely the product, Clayton, Ali-Mikhail-Haq (AMH) and Gumbel copulas. Their features are summarised in Table 1, where we include the parameter space, the generator, the inverse generator, an indicator whether the copula is strict or not and the induced $h(t)$ function obtained through inversion of relationship (6).

Note that, due to the nonunicity of an Archimedean generator, an equivalent constraint to $(C3)$ has to be imposed to the parametric generators that we are going to compare to. That is we set $h(0)=1$ for the product, Clayton and AMH copulas, and $h(\epsilon)=1$ for the Gumbel copula, for say $\epsilon=0.01$. The difference in the latter case is because, for a Gumbel copula, $h(t)\to\infty$ when $t\to 0$. These conditions are already included in the parametrisation used in Table 1.

For each parametric copula we took a sample of size $n=200$. To specify the copulas we took particular values in the parametric space that induce negative and positive dependence. In particular we set $\theta\in\{-0.4,-0.8,0.6,1\}$ for the Clayton copula, $\theta\in\{-0.3,-0.7,0.3,0.7\}$ for the AMH copula, and $\theta\in\{1.4,2.0\}$ for the Gumbel copula. For the partition size we compared $K\in\{10,20\}$ and tried values $\alpha\in\{0.3,0.5,0.9,1,2,\ldots,10\}$.

For the prior distributions (9) we took $\pi_{0}=0,\mu_{0}=-1$ and $\sigma_{0}^{2}=10$. We implemented a MH step within the Gibbs sampler where the proposal distributions were specified by $\pi_{1}=0$ and $\delta=0.25$. The acceptance rate attained with these specifications are around 30%, which according to Robert and Casella, (2010) are optimal for random walks. Finally, the chains were ran for 20,000 iterations with a burn-in of 2,000 and keeping one of every 5th iteration to produce posterior estimates. Convergence of the chains was assessed informally by looking at the trace and ergodic means plots. Computational times using an intel core i7 microprocessor averaged 6 and 15 minutes for the partition sizes $K=10$ and $K=20$, respectively.

To assess goodness of fit (GOF) we computed several statistics. The logarithm of the pseudo marginal likelihood (LPML), originally suggested by Geisser and Eddy, (1979), to assess the fitting of the model to the data. The supremum norm, defined by $\sup_{t}|\varphi^{-1}(t)-\widehat{\varphi}^{-1}(t)|$ to assess the discrepancy between our posterior estimate (posterior mean) $\widehat{\varphi}^{-1}(t)$ from the true inverse generator $\varphi^{-1}(t)$. We also computed the Kendall’s tau coefficient and compare the point (posterior mean) and 95% interval estimates with the true value. Additionally, as a graphical aid to see the performance of our model, we compare the posterior estimates, point (posterior mean) and 95% pointwise credible intervals, of functions $h(t)$ and $\varphi^{-1}(t)$ with the true ones. In general, the idea of our model is to properly estimate the joint density of a particular dataset, say $f(u,v)$, but in Archimedean copulas such a density is characterised by the generator, like in (8). This is why we concentrate on comparing the inverse generator and its associated hazard function.

To avoid overwhelming the reader with many tables and graphs, we only show results for some of the simulated datasets to illustrate, the performance of our model in the other datasets not shown is analogous. The GOF statistics are shown in Tables A Bayesian semiparametric Archimedean copula to A Bayesian semiparametric Archimedean copula. Although we fitted our model with all values of $\alpha$ mentioned above, we only show results for those around the best fitting model in the tables. Posterior estimates of the functions are depicted in Figures 3 to 8. Here we only show estimates with the best fitting model.

For the product copula the GOF measures are presented in Table A Bayesian semiparametric Archimedean copula. With exception of the partition Log-$3$ for $K=20$, for all settings considered, the true $\kappa_{\tau}$ lies inside the 95% credible intervals. The LPML chooses the model with Log-$1$ partitions of size $K=10$, and corresponds to the second smallest value of the supremum norm. Posterior estimates of functions $h(t)$ and $\varphi^{-1}(t)$ are shown in Figure 3. In both cases the true function lies inside the 95% credible intervals.

For the Clayton copula we have two choices of $\theta$, $-0.8$ and $1$. The first choice, $\theta=-0.8$, corresponds to a generator that is not strict, that is, $\varphi^{-1}(t)>0$ for $t\in[0,5/4]$, and $\varphi^{-1}(t)=0$ for $t>5/4$. This is an interesting challenge because our model defined only strict generators. The settings with smallest supremum norm, Log-$0.5$ with $K=10$, produces the 95% credible interval for $\kappa_{\tau}$ closest to the true value, however it does not achieve the largest LPML. The inconsistency of the GOF measures might be due to the non strictness feature of the true generator. Moreover, if we look at the graphs of the posterior estimates of $h(t)$ and $\varphi^{-1}(t)$ (Figure 4), for larger values of $t$ the true functions lie outside of our posterior estimates. For $\theta=1$, the best model is obtained with a Log-$6$ partition of size $10$. In this case, posterior estimates of functions $h(t)$ and $\varphi^{-1}(t)$ with the best fitting (Figure 5), contain the true functions.

For the AMH copula we have two values of $\theta$, $-0.7$ and $0.7$. The best fitting chosen by at least two of the three GOF criteria is obtained with a Log-$6$ and Log-$1$ partitions of size $K=10$, respectively for the two values of $\theta$. Posterior estimates of functions $h(t)$ and $\varphi^{-1}(t)$ with the best fitting are shown in Figures 7 and 6, respectively. In all cases the true functions lie within the 95% credible intervals.

For the Gumbel copula with $\theta=1.4$ we have an interesting behaviour. The true $h(t)$ function has the feature that $h(0)=\infty$. This represents a challenge for our model since we have imposed the constraint $(C3)$ which is equivalent to $h(0)=1$. The highest LPML value is obtained with a Log-$7$ partition of size $K=10$, however the posterior 95% credible interval for $\kappa_{\tau}$ does not contain the true value. On the other hand, the second best value of LPML is obtained with an Log-$3$ partition of size $10$, and in this case the 95% credible interval for $\kappa_{\tau}$ does contain the true value. We select this latter as the best fitting. Posterior estimates of functions $h(t)$ and $\varphi^{-1}(t)$ are shown in Figure 8. Recalling that the true hazard function goes asymptotically to infinity when $t\to 0$, therefore, for values close to zero the true $h(t)$ lies outside our posterior credible intervals, something similar happens in the estimates of the inverse generator. Apart from this, our posterior estimates are very good for $t>\epsilon$.

An important learning from the previous examples is that increasing the partition size does not necessarily imply better fitting.

4.2 Real data analysis

In public health it is important to study the factors that determine the birth weight of a child. Low birth weight is associated with high perinatal mortality and morbility (e.g. Stevens-Simon and Orleans,, 2001). We study the dependence structure between the age of a mother ($X$) and the weight of her child ($Y$), and concentrated on mothers of 35 years old and above. The dataset was obtained from the General Hospital of Mexico through the open data platform that can be accessed at https://datos.gob.mx/busca/dataset/perfiles-metabolicos-neonatales/resource/4ab603eb-b73a-498f-8c56-0dc6d21930e8. It contains $n=208$ records from the first sample of the neonatal metabolic profile of male babies registered in the year 2017 in Mexico City.

The marginal distributions for variables $U$ and $V$, induced by copula (1), are uniform. In practice, copulas are used to model the dependence for any pair of random variables regardless of their marginal distributions. Let $X$ and $Y$ be two random variables with marginal cumulative distributions $F(x)$ and $G(y)$ respectively. Then the joint cumulative distribution function for $(X,Y)$ is obtained as (Sklar,, 1959), $H(x,y)=C(F^{-1}(x),G^{-1}(y))$, where $C$ is given in (1).

Since we are just interested in modelling the dependence between $X$ and $Y$, it is common in practice to transform the original data, $(X_{i},Y_{i})$, $i=1,\ldots,n$, to the unit interval via a modified rank transformation (Deheuvels,, 1979) in the following way. Let ${\bf X}^{\prime}=(X_{1},\ldots,X_{n})$ and ${\bf Y}^{\prime}=(Y_{1},\ldots,Y_{n})$ then $U_{i}=\mbox{rank}(i,{\bf X})/n$ and $V_{i}=\mbox{rank}(i,{\bf Y})/n$ are the transformed data, where $\mbox{rank}(i,{\bf X})=k$ if and only if $X_{i}=X_{(k)}$ for $i,k=1,\ldots,n$. This is based on the probability integral transform using the empirical cumulative distribution function of each coordinate.

In Figure 9 we show a dispersion diagram of the original data (left panel) and the rank transformed data (right panel). To avoid problems due to ties in the original data, we first include a perturbation to the data by adding a uniform random variable $\mbox{Un}(-0.0,0.01)$ to each coordinate. The sample Kendall’s tau value for the transformed data is $\tilde{\kappa}_{\tau}=-0.1162$.

We fitted our model to the transformed data with the following specifications. To define the partitions we took values $\alpha\in\{0.3,0.5,0.9,1,2,\ldots,10\}$ with sizes $K\in\{10,20\}$. For the prior we took $\pi_{0}=0$, $\mu_{0}=-1$ and $\sigma_{0}^{2}=10$. The MCMC specifications were the same as those used for the simulated data.

The GOF measures computed were the LPML and the posterior estimates (point and 95% credible interval) of $\kappa_{\tau}$. The results are reported in Table A Bayesian semiparametric Archimedean copula. The best fitting model according to LPML is that obtained with a partition of size $K=10$ and Log-$10$. The sample concordance $\tilde{\kappa}_{\tau}$ is included in our posterior 95% credible interval estimate $\kappa_{\tau}\in(-0.213,-0.098)$.

The estimated hazard rate function $h(t)$ and the inverse generator $\varphi^{-1}(t)$, with the best fitting model, are included in the top row in Figure 10. The solid thick line corresponds to the point estimates and the solid thin lines to the 95% credible intervals. For a visual comparison, the blue dotted line corresponds to the functions of the independence (product) copula. Additionally, we include an estimate of the joint density as well as the corresponding contour plots (bottom row in Figure 9. These estimates suggest that there is a negative (weak) dependence between the age of the mother and the birth weight of the child. The older the mother, the less weight of the child. This finding could potentially help the policy makers to focus campaigns to help the awareness of future mothers.

4.3 Independence test

As mentioned in Section 2, we can use our model to undertake an independence test. For that we choose the prior distribution for the $\theta_{k}$’s, as in (9), such that the prior probability of $H_{0}:\theta_{1}=\cdots=\theta_{K-1}=0$ is $1/2$, in other words, we want $\mbox{P}(H_{0})=\pi_{0}^{K-1}=1/2$. Particularly, for a partition of size $K=10$ we need to specify $\pi_{0}=0.9258$. In order to get a point of mass proposal in the MH step we consider $\pi_{1}=0.3$. We re-ran our model using these values with the other specifications left unchanged and performed the test for all simulated and real datasets.

To place our test in context, we compare our results with the recently proposed independence test of Filippi et al., (2016), based on Dirichlet process mixture models. These authors actually proposed two tests, one based on a contingency table approach (CT) and another based on a mixture model approach (MM). Additionally, we implemented a frequentist test based on the empirical copula (EC) given in Deheuvels, (1979). For the three Bayesian tests, ours (SPAC) and the other two competitors, we report $\mbox{P}(H_{1}\mid\mbox{data})$, whereas for the frequentist test we report the p-values. All these values are included in Table 9.

We first mention that the values $\mbox{P}(H_{1}\mid\mbox{data})$ from the Bayesian tests have to be calibrated with respect to that obtained for the product (independent) dataset. The three tests assign small evidence of dependence to the product dataset, as it should be, whereas the frequentist test assigns a p-value of 0.15 to the same product dataset, which is large enough to not to reject the null hypothesis of independence.

For the Clayton and Gumbel datasets, all four test are consistent giving enough evidence to dependence. For the AMH datasets we have mixed decisions. None of the four tests are able to detect dependence for the cases of $\theta=-0.3$ and $\theta=0.3$. This is understandable since the AMH copula produces data that look similar to the product copula for values of $\theta$ close to zero. For the other two values, $\theta=-0.7$ and $\theta=0.7$, the frequentist test EC does not detect dependence, however the Bayesian tests give more evidence of dependence, being our SPAC test the one that best supports dependence for these two datasets.

Finally, for the real dataset, we also have mixed decisions. Tests CT and EC do not detect any dependence, however, our new test SPAC and MM give enough support to dependence, which is also consistent to the estimated generator obtained with our model and presented in the right panel of Figure 10.

5 Concluding remarks

We have proposed a semiparametric Archimedean copula that is flexible enough to capture the behaviour of several families of parametric Archimedean copulas. Our model is capable of modelling positive and negative dependence. The number of parameters in the model to produce a good estimation of the dependence in the data should not be extremely high. For most of the examples considered here ten parameters are enough.

Defining an appropriate partition to analyse real data sets is not trivial. We suggest to try different values of $\alpha$ in a wide range and compare using a GOF criteria like the LPML we used here.

Our proposal is also suitable to perform an independent test, which compares favourable with alternative independence tests. For the datasets considered here, our proposal assigned the largest evidence of dependence for the dependent datasets.

In the exposition and in examples considered here, we concentrated on bivariate copulas, however extensions to more than two dimensions is also possible, say $C(u_{1},\ldots,u_{m})=\varphi^{-1}\left(\varphi(u_{1})+\cdots+\varphi(u_{m})\right)$. Performance of our semiparametric copula in this multivariate setting is worth studying.

Our model is motivated by semiparametric proposals for survival analysis functions (Nieto-Barajas and Walker,, 2002) and appropriately modified to satisfy the properties of an Archimedean generator. The semiparametric generator presented here turned out to be based on quadratic splines, however, alternative proposals are possible as the one described in Section 2.2.

Although, the motivation of our proposal lies within a survival analysis context, the inclusion of right censored observations into the analysis is not straightforward. The likelihood contribution would involve the cumulative distribution function of the bivariate density induced by the copula, and this is not available in closed form. A data augmentation technique, like those in Tanner, (1991), would be the way to proceed.

Acknowledgements

This research was done while the first author was doing a post doctorate at the Department of Statistics, ITAM, and was supported by Asociación Mexicana de Cultura, A.C. The authors are grateful to three anonymous referees for their insightful comments.

Appendix

Derivation of posterior conditional support of $\theta_{k}$.

In order to satisfy constraint $(C1)$, we consider first the case $\theta_{k}\leq 0$. Therefore

$\min_{t\in(\tau_{k-1},\tau_{k}]}A_{k}+\theta_{k}t=A_{k}+\theta_{k}\tau_{k}$. This implies the following constraint for $\theta_{k}$,

[TABLE]

for $k=1,\ldots,K-1$, where we define the empty sum as zero.

On the order hand, if $\theta_{k}>0$ we have $\min_{t\in(\tau_{k-1},\tau_{k}]}A_{k}+\theta_{k}t=A_{k}+\theta_{k}\tau_{k-1},$ and we get, from condition $(C2)$, the following restriction

[TABLE]

This defines the upper bound $b_{k}$, for $k=2,\ldots,K-1$, and $b_{1}=1$.

Because the term $\theta_{k}$ appears on the right side of the previous inequality for $j=k+1,\ldots,K-1$, we need to consider the following restriction

[TABLE]

if $\theta_{j}\geq 0$. Combining this with the constraint when $\theta_{k}\leq 0$ above, we get the lower bound $a_{k}$ for $k=1,\ldots,K-1$.