Svetlana Gribkova
email: svetlana.gribkova@etu.upmc.fr

Olivier Lopez
email: olivier.lopez0@upmc.fr

Philippe Saint
email: philippe.saint_pierre@upmc.fr

A simplied model for studying bivariate mortality under right-censoring

Keywords: Survival analysis, nonparametric estimation, Kaplan-Meier estimator, bivariate censoring, mortality analysis, copula models, bootstrap Bivariate mortality under right-censoring

In this paper, we provide a nonparametric estimator of the distribution of bivariate censored lifetimes, in a model where the two censoring variables dier only through an additional observed variable. This situation is motivated by a particular application to insurance, where the supplementary variable corresponds to the age dierence between two individuals. Asymptotic results for our estimator are provided. The new tools that we develop are used to perform goodness-of-t tests for survival copula models. The practical performance is illustrated through simulations and a real data analysis.

Introduction

In last survivor insurance, an important issue is to infer on the joint distribution of the lifetimes of two individuals linked through an insurance contract, say (T, U ). One of the diculties in studying such variables comes from the presence of bivariate censoring, with a proportion of censored observations which may be quite high. Therefore, most of the approaches used in this eld are parametric (typically parametric survival copula models, see e.g. [START_REF] Shih | Inferences on the association parameter in copula models for bivariate survival data[END_REF]), while nonparametric tools are rarely used, although they would be required at least to assess the validity of the proposed models. The aim of this paper is to provide a new nonparametric estimator of the joint distribution of two lifetimes under bivariate random censoring, in a framework which is adapted to the study of problems coming from the insurance eld. A specicity of such problems is that an additional variable is generally present, which carries information on the model, this variable being the age dierence between the two individuals under study. Using this information that is often neglected, one can dene a quite simple nonparametric estimator of the distribution of the two lifetimes, which is close to Kaplan-Meier estimator [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF]) and to the estimator of [START_REF] Lin | A simple nonparametric estimator of the bivariate survival function under univariate censoring[END_REF].

Various approaches have been used to perform nonparametric estimation of multivariate lifetimes. Most of them focus on estimating the survival function, without focusing on the joint distribution itself. Therefore, many of them provide consistent estimators of this function, but fail to dene a true distribution. For example, the estimator of the survival function proposed by [START_REF] Campbell | Large-sample properties of nonparametric bivariate estimators with censored data[END_REF] is not monotonic. The nonparametric maximum likelihood (NPMLE) procedure of [START_REF] Hanley | Nonparametric estimation of a multivariate distribution in the presence of censoring[END_REF] leads to an estimator which is sometimes inconsistent for continuous data [START_REF] Tsai | Nonparametric estimation of a bivariate survival function in the presence of censoring[END_REF], while the rate of convergence of a modication of this estimator suggested by [START_REF] Tsai | Nonparametric estimation of a bivariate survival function in the presence of censoring[END_REF] achieves a slow convergence rate (slower that n 1/2 where n denotes the sample size). Another NPMLE approach is proposed by [START_REF] Van Der Laan | Ecient estimation in the bivariate censoring model and repairing NPMLE[END_REF], introducing some modication of the data and using an interval censoring methodology. Although this estimator is shown to be asymptotically ecient for these modied data, the convergence rate is also slower than n 1/2 . On the other hand, the product-limit type estimator proposed by [START_REF] Dabrowska | Kaplan-Meier estimate on the plane[END_REF], which is often used in practice (see e.g. [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF], [START_REF] Fan | Dependence estimation over a nite bivariate failure time region[END_REF], [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF], [START_REF] Gill | Inecient estimators of the bivariate survival function for three models[END_REF], [START_REF] Prentice | Covariance and survivor function estimation using censored multivariate failure time data[END_REF]), assigns negative mass to some points on the plane [START_REF] Pruitt | Strong consistency of self-consistent estimators: general theory and an application to bivariate survival analysis[END_REF]. Nonparametric smoothing techniques have been used by Pruitt (1991a) (but the implicit denition of the estimator leads to diculties and a weak performance according to [START_REF] Van Der Laan | Ecient estimation in the bivariate censoring model and repairing NPMLE[END_REF]), and by [START_REF] Akritas | Estimation of bivariate and marginal distributions with censored data[END_REF] (this last estimator presenting the advantage to dene a true distribution, but, again, with a slower convergence rate, and the necessity of an absolutely continuous censoring variables).

Among the approaches that we mention, each of them suers either from a too slow convergence rate, or from the fact that the corresponding estimators do not provide true probability distributions. The estimator proposed by [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF] does not present these drawbacks, but relies on an assumption on the joint distribution of the censoring variables which may not be reasonable for the particular application we have in mind. Indeed, a specicity of data-sets coming from last-survivor insurance, is that individuals usually quit the study (for some other cause than death) at the same time.

This induces a specic dependence between the two censoring variables involved in the problem. The main idea of the new estimator that we propose consists of using this additional information.

The rest of the paper is organized as follows. In section 2, we present the general censoring framework that we consider. We dene a non parametric estimator of the distribution of (T, U ). In section 3, we provide asymptotic results for estimating quantities of the type E[ϕ(T, U )] for a large class of functions ϕ (the survival function being only a particular case). A bootstrap procedure is proposed to compute the variance of the error in such estimation problems. Application of our nonparametric estimator to goodnessof-t for copula models is considered. Section 4 illustrates our result through simulation studies and a real data analysis.

A simplied model for bivariate right-censoring

We rst present in section 2.1 the general bivariate censoring model that will be considered in the rest of the paper. Estimation of the joint distribution of (T, U ) is introduced in section 2.2.

Bivariate right-censoring

In the following, we consider two lifetimes (T, U ), and i.i.d. replications (T i , U i ) 1≤i≤n of these random variables. In a bivariate right-censoring model, (T i , U i ) 1≤i≤n are not directly observed. Instead, one observes

{ Y i = inf(T i , C i ), and δ i = 1 T i ≤C i Z i = inf(U i , D i ), and γ i = 1 U i ≤D i ,
where (C i , D i ) 1≤i≤n consists of i.i.d. replications of a random bivariate censoring vector (C, D), and (δ i , γ i ) 1≤i≤n are indicator functions allowing the distinction between censored and uncensored observations.

In many applications, such as last survivor insurance, there exists some clear relationship between the two censoring variables (C, D). In this particular case, if T (resp. U ) denotes the total lifetime of the husband (resp. his wife), C (resp. D) denotes the age at which the husband (resp. the wife) stops being under observation for any other cause than death.

Usually, censoring causes are twofold: the end of the observation period (if the person is not dead at this time, his/her lifetime is not observed), or the surrender of the contract.

In both situations, one can observe that, for many cases, this event which stops observation occurs at the same time for both members of the couple. Taking the example of a pension contract with a reversion clause, one can see that surrendering the contract will automatically remove the two members of the couple from the database (unless one of them possesses additional contracts that could allow the company to keep some track on him/her, which is usually not the case due to the complexity of such a tracking process).

If ε denotes the age dierence between the two members of the couple, then D = C + ε.

Moreover, the variables (ε i ) 1≤i≤n are observed for all couples.

To summarize, in such a framework, observations are made of

(Y i , Z i , ε i , δ i , γ i ) 1≤i≤n ,
where the random variables ε i represents some age dierence between the two persons observed. We now state some identiability assumptions that will hold throughout this paper.

Assumption 1 Assume that 1. D = C + ε.

(T, U

) is independent from ε, and from C, and

P(T = C) = P(U = C + ε) = 0.
3. C is independent from ε.

In Assumption 1, points 2 and 3 are a direct multivariate extension of the classical identiability assumption required to ensure consistency of Kaplan-Meier estimator in the univariate case (see [START_REF] Stute | The strong law under random censorship[END_REF]). In section 2.2 below, we show how one can estimate nonparametrically the distribution of (T, U ) under Assumption 1.

Nonparametric estimation of the distribution of (T, U )

In this section, we show how to estimate quantities of the type E[ϕ(T, U )] for some function ϕ. A particular case is the joint survival function S F (t, u) = P(T > t, U > u), which corresponds to ϕ(T, U ) = 1 T >t,U >u . In absence of censoring, the answer to this estimation problem consists of using empirical means. Dening F (t, u) = P(T ≤ t, U ≤ u), one can rewrite E[ϕ(T, U )] = ∫ ϕ(t, u)dF (t, u), which can be consistently estimated by ∫ ϕ(t, u)dF emp (t, u), where F emp (t, u) = n -1 ∑ n i=1 1 T i ≤t,U i ≤u denotes the empirical distribution function. In our framework, the empirical distribution in unfortunately unavailable, since (T i , U i ) 1≤i≤n are not directly observed.

We propose to rely on an estimator of the type

F (t, u) = n ∑ i=1 δ i γ i W n (Y i , Z i , ε i )1 Y i ≤t,Z i ≤u , (2.1)
to generalize the empirical distribution function to our framework. Using such type of estimators, one can straightforwardly dene an estimator of

E[ϕ(T, U )] = ∫ ϕ(t, u)dF (t, u) by ∫ ϕ(t, u)d F (t, u) = n ∑ i=1 δ i γ i W n (Y i , Z i , ε i )ϕ(Y i , Z i ). (2.2)
The idea is similar to [START_REF] Lopez | A generalization of kaplan-meier estimator for analyzing bivariate mortality under right-censoring and left-truncation with applications to model-checking for survival copula models[END_REF] and [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF]: instead of assigning the same n -1 -mass to each observation (as it is the case when considering the empirical distribution function), one assigns mass at doubly uncensored observations (since only these observations are completely relevant to understand the dependence structure between T and U ), while the mass

W n (Y i , Z i , ε i ) is designed to compensate censoring.
Dening S G (t) = P(C > t), (with G(t) = P(C ≤ t)) and

F * (t, u) = ∑ n i=1 δ i γ i W * (Y i , Z i , ε i )1 Y i ≤t,Z i ≤u , with W * (y, z, e) = n -1 S G (max(y, z -e)-)
E[δγψ(Y, Z)] = E [ E [ 1 max(T,U -ε)≤C |T, U, ε ] ψ(T, U ) ] = E [S G (max(T, U -ε)-)ψ(T, U )] .
Unfortunately, this ideal estimator F * can not be computed in practice, since it relies on the unknown survival function S G . Nevertheless, it is possible to estimate this function

S G . Dene η i = 1 -δ i γ i , and A i = max(T i , U i -ε i ). The variable C i is observed as long as C i < A i (that is η i = 1
). Hence, C can also be considered as a right-censored variable (censored by variables A), provided that the event {C i = A i } has probability 0. This is actually the case, from point 2 in Assumption 1, where we assumed that P(T

= C) = P(U -ε = C) = 0.
Therefore, the distribution of C can be estimated by a Kaplan-Meier estimator based on the censored sample (B i , η i ) 1≤i≤n , where

B i = inf(C i , A i ). Moreover, it
is important to notice that Assumption 1 ensures consistency of Kaplan-Meier estimator, that is, following [START_REF] Stute | The strong law under random censorship[END_REF], P(A = C) = 0 and A independent from C.

Therefore, dening the Kaplan-Meier estimator ŜG of S G , that is,

ŜG (t) = ∏ k:B k ≤t ( 1 - d Ĥ0 (B k ) Ĥ(B k )
)

,

where

Ĥ0 (t) = n -1 ∑ n i=1 η i 1 B i ≤t , and Ĥ(t) = n -1 ∑ n i=1 1 B i ≥t , a natural choice of a function W n in (2.1) is W n (y, z, e) = 1 n ŜG (max(y, z -e)-) . (2.3)
This estimator is close to the estimator proposed by [START_REF] Lin | A simple nonparametric estimator of the bivariate survival function under univariate censoring[END_REF]. The dierence, in our approach, is the presence of the additional random variable ε i corresponding to the age dierence.

Asymptotic theory

The present section is devoted to the asymptotic results on the nonparametric estimator dened in section 2.2. A Central Limit Theorem for (2.2) is provided in section 3.1. As a corollary of this result, we deduce asymptotic convergence results for estimating Kendall's τ coecient, which is a classical dependence measure. Section 3.2 provides a bootstrap procedure in order to investigate the estimation error. In section 3.3, we derive theoretical results that may be used to perform goodness-of-t tests for survival copula models.

An asymptotic representation Theorem for estimator (2.1)

We aim to obtain an asymptotic representation for quantities of the type (2.1). Instead of considering a single function ϕ, we focus on obtaining results that hold uniformly for functions ϕ ∈ F, F denoting a class of functions. This uniformity result is required if we wish to obtain, for example, uniform consistency results for the estimation of the distribution function. Considering this problem, the natural class of functions to be considered is F 1 = {(t, u) → 1 t≤x,u≤y : x ∈ T , y ∈ U}, where T and U denote the support of the distribution of each marginal.

In the following, we consider a class of functions F, with envelope Φ, satisfying Assumptions 2 and 3 below. Since our proof will rely on empirical processes theory, Assumption 2 consists of assuming that a class of functions related to F is Donsker, that is a class with an uniform central limit theorem property (see [START_REF] Van Der Vaart | Weak convergence and empirical processes with applications to statistics[END_REF] for a precise denition of Donsker classes).

Assumption 2 Let G denote the class of positive, monotonic functions bounded by 1, and

χ(T, U, C, ε) = δγS G (max(T, U -ε)-) -2 . For any (t 0 , u 0 ) in R 2 such that S F (t 0 , u 0 ) > 0, dene H t 0 ,u 0 = {(T, U, C, ε) → χ(T, U, C, ε)f (T, U )g(max(T, U -ε)-)1 T ≤t 0 ,U ≤u 0 , f ∈ F , g ∈ G},
and assume that H t 0 ,u 0 is a Donsker class of functions.

Assumption 3 is required only if we wish to obtain consistency on the whole support of (T, U ). It automatically holds if one considers bounded functions with compact support strictly included in the support of the distribution. It can be understood as an assumption on the tail of the distributions of T and U. Similar assumptions have been used in [START_REF] Gill | Large sample behaviour of the product-limit estimator on the whole line[END_REF], [START_REF] Stute | Distributional convergence under random censorship when covariables are present[END_REF], or [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF].

Assumption 3 Assume that E [Φ(T, U ) 2 S G (max(T, U -ε)-) -1 ] < ∞. Moreover, dene C(y) = ∫ y -∞ dG(t) [1 -F (t)][1 -G(t-)] 2 , and assume that E[Φ(T, U )C 1/2+υ (max(T, U -ε)-)S G (max(T, U -ε)-) -1 ] < ∞, for some υ > 0 (arbitrary small).
If we consider the particular case of F 1 , Assumption 2 automatically holds, provided that the moment conditions of Assumption 3 hold. Generally, this will also be the case for parametric class of functions, or suciently smooth classes of functions using permanence properties of Donsker classes.

We now state the main theoretical result of this section.

Theorem 3.1 Recall that

B i = inf(A i , C i ), where A i = max(T i , U i -ε i ), and that η i = 1 C i ≤A i , and let F A (t) = P(A ≤ t) and H(t) = P(B > t). Under Assumptions 1 to 3, ∫ ϕ(t, u)d( F -F * )(t, u) = 1 n n ∑ i=1 ψ ϕ (T i , U i , ε i ) + R n (ϕ), (3.1)
where sup ϕ∈F |R n (ϕ)| = o P (n -1/2 ), and

ψ ϕ (Y i , Z i , ε i ) = ∫ { η i S G (B i ∨ a) H(B i ) - ∫ 1 B i ≥u S G (u ∨ y)dF A (u) H(u)F A (u) + (1 -η i )1 B i >y F A (B i ) -S G (y) } ϕ(t, u)dP (T,U,C,ε) (t, u, c, e) S G (a-) ,
where we used a = max(t, u -e) to shorten the notation, and where P (T,U,C,ε) denotes the true law of (T, U, C, ε). As a consequence, since E[ψ ϕ (Y, Z, ε)] = 0, we have, for all

ϕ ∈ F , n 1/2 (∫ ϕ(t, u)d F (t, u) -E[ϕ(T, U )] ) =⇒ N (0, σ 2 ϕ ), (3.2)
with

σ 2 ϕ = E [ { δ i γ i ϕ(Y i , Z i ) S G (B i -) -E[ϕ(T, U )] + ψ ϕ (Y i , Z i , ε i ) } 2 ] ,
and =⇒ denotes the weak convergence.

The proof of this result is postponed to the Appendix section. Equation (3.2) can be used to compute asymptotic condence intervals, provided that one is able to consistently estimate the asymptotic variance σ 2 ϕ . This can be done by replacing all the unknown distributions functions involved in the expression of σ 2 ϕ by their empirical counterparts. However, this approximation may be too rough in practice, and bootstrap procedures seem to be more appropriate if one wishes to compute condence interval. This bootstrap method is shown in section 3.2.

In addition to the application of Theorem 3.1 to the class F 1 (corresponding to the estimation of the joint survival function), we show how this result may be used to provide asymptotic results for the estimation of Kendall's τ coecient. Kendall's τ coecient is a classical dependence measure which can be dened in the following way. For two random variables (T, U ),

τ = P((T 1 -U 1 )(T 2 -U 2 ) > 0)-P((T 1 -U 1 )(T 2 -U 2 ) < 0), where (T 1 , U 1 )
and (T 2 , U 2 ) are independent replications of (T, U ). There exists a relationship between τ and the distribution function [START_REF] Nelsen | An introduction to copulas[END_REF]. Therefore, a natural estimator of τ is

F, that is τ = 4 ∫ F (x, y)dF (x, y) -1, see e.g.
τ = 4 ∫ F (x, y)d F (x, y) -1. (3.3)
As it is shown in Wang and Wells (2000a), censoring may cause this estimator not to be consistent in some particular situations. Indeed, dening S H (y, z) = P(Y > y, Z > z) the survival function of the observed times, S 1 = {(t, u) : S H (t, u) > 0}, and S 2 = {(t, u) :

S F (t, u) > 0}, we can see that some part of the distribution, namely S 2 -S 1 is never observed, since the corresponding observations are always censored. If this dierence of sets is empty, this does not introduce bias in the estimation of τ. In other situations, some bias will arise, which can be evaluated according to the method of Wang and Wells (2000a).

Corollary 3.2 below shows that this estimator admits an asymptotic representation.

Corollary 3.2 Let ψ F denote function ψ as dened in Theorem 3.1 for function ϕ = F.

To shorten the notation, we will denote ψ t,u the function corresponding to ϕ(Y, Z) =

1 Y ≤t,Z≤u . Assume that S 2 -S 1 = ∅. Then, τ -τ = 4 {∫ F (t, u)d[F * -F ](t, u) + ∫ [F * (t, u) -F (t, u)]dF (t, u) + 1 n n ∑ i=1 { ψ F (Y i , Z i , ε i ) + ∫ ψ t,u (Y i , Z i , ε i )dF (t, u) } } + o P (n -1/2 ). (3.4)
In the representation of Corollary (3.2), each term is a sum of i.i.d. quantities with zero expectation and nite variance. Therefore, Corollary (3.2) shows that τ is asymptotically Gaussian. Its asymptotic variance (which as a complex form) can be deduced from this representation. Nevertheless, we do not emphasize this variance, since we recommend using bootstrap procedures to investigate the law of τ (see section 4.1).

Let us also mention that, if the assumption S 2 -S 1 = ∅ does not hold, the result is still true, but with τ replaced by 4 ∫ S 1 F (x, y)dF (x, y) -1.

Proof of Corollary 3.2. Write

τ -τ = 4 {∫ F (t, u)d[ F -F ](t, u) + ∫ [ F (t, u) -F (t, u)]dF (t, u) + ∫ [ F (t, u) -F (t, u)]d[ F -F ](t, u) } . (3.5)
Applying Theorem 3.1 to function F , the rst term of (3.5) can be expanded as

∫ F (t, u)d[F * -F ](t, u) + 1 n n ∑ i=1 ψ F (Y i , Z i , ε i ) + o P (n -1/2 ).
Moreover, again from Theorem 3.1,

∫ [ F (t, u) -F (t, u)]dF (t, u) = ∫ [F * (t, u) -F (t, u)]dF (t, u) + 1 n n ∑ i=1 ∫ ψ t,u (Y i , Z i , ε i )dF (t, u) + o P (n -1/2 ).
The third term of (3.5) can be rewritten as

∫ [F * (t, u) -F (t, u)]d[F * -F ](t, u) + ∫ [ F (t, u) -F * (t, u)]d[F * -F ](t, u) + ∫ [ F (t, u) -F * (t, u)]d[ F -F * ](t, u) := T 1 + T 2 + T 3 .
The term T 1 is a second order degenerate U -statistics and is therefore of order O P (n -1 ).

To study T 2 , apply Theorem 3.1 to the class of indicator functions 1 T ≤t,U ≤u to obtain

T 2 = ∫ ψ(T i , U i , ε i )d[F * -F ](t, u) + o P (n -1/2 ).
The integral in this decomposition is zero, since E[ψ ϕ (T, U, ε)] = 0, and since

∫ ϕd[F * - F ](t, u) = 0.
Finally, observe that T 3 can be rewritten as

T 3 = 1 n n ∑ i=1 δ i γ i [ F (Y i , Z i ) -F * (Y i , Z i )][ ŜG (B i -) -S G (B i -)] S G (B i -) ŜG (B i -) ,
which is bounded by

|T 3 | ≤ sup t,u | F (t, u) -F * (t, u)| sup b | ŜG (b) -S G (b)| sup b S G (b) ŜG (b) × 1 n n ∑ i=1 η i S G (B i -) 2 .
The rst supremum is O P (n -1/2 ) from Theorem 3.1, the second one is o P (1) from the uniform consistency of Kaplan-Meier estimator (see [START_REF] Stute | The strong law under random censorship[END_REF]), while the third one is O P (1) (see [START_REF] Gill | Large sample behaviour of the product-limit estimator on the whole line[END_REF]). Moreover, the empirical mean on the right-hand side is O P (1) provided that each term of the sum has nite expectation, which is the case from Assumption 3. Combining these facts leads to T 3 = o P (n -1/2 ) and concludes the proof.

Bootstrap procedure

As we already mentioned, the asymptotic results of Theorem 3.1 may be dicult to use when it comes to approximate the law of ∫ ϕ(t, u)d F (t, u). The problem comes from the complex from of the asymptotic variance, and therefore from the diculty to estimate it in an accurate way. Therefore, the aim of the present section is to propose a bootstrap procedure that allows to circumvent this problem.

Under univariate censoring, two main methodologies have been proposed in the literature to perform bootstrap, see [START_REF] Efron | Censored data and the bootstrap[END_REF] and [START_REF] Reid | Estimating the median survival time[END_REF]. The methodology of Efron (1981) consists of using the nonparametric estimators of the distribution of the lifetime and of the censoring to resimulate samples. [START_REF] Akritas | Bootstrapping the Kaplan-Meier estimator[END_REF] showed that only Efron's methodology was consistent. We therefore propose to adapt this strategy.

The basic idea consists of simulating variables (T, U ) according to the estimated distribution dened by F (and renormalized in order to ensure that the total mass is equal to one). The censoring can be simulated similarly using Ĝ, while ε is simulated according to its empirical distribution Fε (t) = n -1 ∑ n i=1 1 ε i ≤t . The procedure is summarized below. To compute B bootstrap n-samples, repeat for b = 1, ..., B the following simulation scheme, 1. Simulate independent variables (T b i , U b i ) 1≤i≤n under the probability distribution F / F (R 2 ).

2. Simulate independent variables (ε b i ) 1≤i≤n under the probability distribution Fε .

3. Simulate independent variables (C b i ) 1≤i≤n under the probability distribution Ĝ/ Ĝ(R 2 ).

The

b-th bootstrap sample is composed of (Y b i , Z b i , δ b i , γ b i , ε i ) 1≤i≤n , where Y b i = inf(T b i , C b i ), Z b i = inf(U b i , C b i + ε b i ), δ b i = 1 T b i ≤C b i , γ b i = 1 U b i ≤C b i +ε b i .

Application to survival copula inference

Survival copula models are a common tool to model dependence between two lifetimes (T, U ). Indeed, the bivariate survival function S F (t, u) = P(T > t, U > u) of the random vector (T, U ) admits, by Sklar's Theorem [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF]), a copula representation, that is

S F (t, u) = C(S T (t), S U (u)),
where S T (t) = P(T > t) and S U (u) = P(U > u), and where C is a survival copula function (see e.g. [START_REF] Nelsen | An introduction to copulas[END_REF]). To understand the dependence between T and U, which is represented by the copula function C, it is natural to search for an estimator of C, usually based on a parametric or semiparametric model, see e.g. [START_REF] Shih | Inferences on the association parameter in copula models for bivariate survival data[END_REF]. Nonparametric inference is then required to assess the validity of the model. [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] proposed to extend the methodology of [START_REF] Genest | Statistical inference procedures for bivariate Archimedean copulas[END_REF] in presence of censoring. This approach relies on the estimation of the function

v → K(v) = P(S F (T, U ) ≤ v).
If we consider the particular case of Archimedean copula families (that is copulas dened as

C(u, v) = ϕ -1 (ϕ(u) + ϕ(v))
where the generator ϕ is a convex function satisfying the conditions of Theorem 4.3.4 in [START_REF] Nelsen | An introduction to copulas[END_REF]), there exists a one-to-one correspondance between the generator ϕ and the function K, through the relationship

K(v) = v - ϕ(v) ϕ ′ (v) . (3.6)
The basic idea of goodness-of-t procedures based on K consists of comparing a parametric estimator, (based on an estimator ϕ θ depending on the parametric model and on θ, the association parameter estimated from the data) to a nonparametric one. [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] used an estimator based on the nonparametric estimator of [START_REF] Dabrowska | Kaplan-Meier estimate on the plane[END_REF].

The nonparametric estimator that we propose to use is dened as

K(v) = ∫ 1 ŜF (t,u)≤v d F (t, u). (3.7)
In Proposition 3.3, we show that the process n 1/2 ( K(v) -K(v)) converges towards a

Gaussian process. This kind of result is essential to legitimate goodness-of-t techniques that will be fully discussed in section 4.2. Nevertheless, we do not focus on the estimation of the asymptotic covariance process. In practice, since its computation seems rather delicate, it is preferable to rely on bootstrap procedure (see section 4.2).

Proposition 3.3 Assume that:

1. The distribution function K(v) admits a continuous bounded derivative k(v).

2. Given S F (t, u) = v, there exists a version of the conditional distribution of (Y, Z) and a countable family P of partitions E on I (where I denotes the support of (T, U )) into a nite number of Borel sets satisfying inf E∈P max E∈E diam(E) = 0, such that, for all E ∈ E, the mapping

v → µ v (E) = k(v)P ((T, U ) ∈ E|S(T, U ) = v) is continuous.
Then, there exists a zero-mean Gaussian process such that,

n 1/2 ( K(•) -K(•) ) =⇒ - ∫ ∫ 1 S F (t,u)>• dW (t, u) - ∫ ∫ W (t, u)dµ • (t, u).
Proof. From Theorem 3.1, one can deduce

n 1/2 ( ŜF (t, u) -S F (t, u)) =⇒ W (t, u),
where W is a Gaussian process with mean zero. Consequently, Theorem 1 in [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] applies.

Simulations and real data example

In this section, we investigate the nite sample size behaviour of our procedure. This investigation is done through simulation studies and illustrated on a real data example.

The data that we consider has been initially studied by [START_REF] Frees | Annuity valuation with dependent mortality[END_REF], and was studied by [START_REF] Carriere | Bivariate survival models for coupled lives[END_REF], [START_REF] Youn | Statistical aspects of joint life insurance pricing[END_REF], [START_REF] Youn | Pricing practices for joint last survivor insurance[END_REF] and [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF]. We refer to [START_REF] Frees | Annuity valuation with dependent mortality[END_REF] for a more detailed description of this dataset, containing lifetimes of two members of a couple who subscribed an insurance contract. The dataset concerns 14947 contracts from a large Canadian insurer, observed between Decembre 29th, 1988 and Decembre 31th, 1993 . After elimination of samesex contracts and of couples with more than one policies (for which we only keep one policy), 11454 contracts remain. In addition to bivariate censoring, observations are subject to left truncation. Nevertheless, we do not consider left-truncation in the approach that we develop in the present paper. Neglecting left-truncation will lead to a slight over-estimation of the lifetimes, which, from the prospective of an insurer who wishes to evaluate his liabilities in the case of a pension contract, represents a cautious approach.

In section 4.1, we discuss the problem of estimating Kendall's τ coecient, illustrating the theoretical results of Corollary 3.2. In section 4.2, we study the practical implementation of a goodness-of-t procedure for copula models, based on the process K dened in (3.7).

Nonparametric estimation of Kendall's τ coecient

Real data example. Using τ dened in equation ( 3.3), we nd an estimated value of Kendall's τ coecient which is τ = 0.6696, which is roughly of the same order as the values obtained by other authors on the same data-set (for example, for a specic generation, [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF] obtained an estimation which is 0.6039). We used the nonparametric bootstrap procedure described in section 3.2 to approximate the law of τ .

Through B = 1000 bootstrap replications, we obtain an estimation of the distribution of τ which is represented in Figure 4.1 below.

The author wishes to thank the Society of Actuaries, through the courtesy of Edward J. We can observe on Figure 4.1 that the distribution of τ obtained using the bootstrap procedure does not seem to be Gaussian. Therefore, it legitimates to rely rather on this bootstrap procedure than on normal approximation to investigate uncertainty in estimating τ.

Simulation study. To illustrate the convergence of τ , we present some results of a simulation study. The random lifetimes (T, U ) are simulated from a Clayton copula model (see Table 2 in section 4.2 for a precise denition) with association parameter θ = 2 (which corresponds to a value τ = 0.5), with marginals following a Weibull distribution. Weibull distribution is parametrized through a shape parameter α and a scale parameter β, and admits a density

f (t) = α β ( t β ) α-1 exp ( - t α β α
) ,

for t ≥ 0. We consider the case α = 10 and β = 1.7. The censoring variables C are simulated according to an exponential distribution with parameter λ (with mean λ -1 ).

Dierent values of λ are considered in order to change the average proportion of doubly uncensored observations. Random variables ε i are simulated according to an exponential distribution with parameter µ = 50.

For each considered value of the parameter λ, we generate n-samples for dierent values of n. We repeat N = 1000 times the simulation scheme in order to estimate the bias E[τ -τ ], the variance V ar(τ ), and the mean-squared error E[(τ -τ ) 2 ]. Results are presented in Table 1 below. Table 1: Estimation of the mean-squared error and related quantities for the estimation of Kendall's τ coecient.

Model

Goodness-of-t for semiparametric copula models

A goodness-of-t procedure based on K. Consider a parametric family of Archimedean survival copulas F C = {C θ : θ ∈ Θ}. We will denote ϕ θ the Archimedean generator of copula C θ . We describe how to extend the procedure proposed by [START_REF] Genest | Statistical inference procedures for bivariate Archimedean copulas[END_REF] to test

H 0 : C ∈ F C , against H 1 : C / ∈ F C .
The principle of the test consists of computing an estimator θ (assuming that H 0 holds) from the data, and then use ϕ θ and (3.6) to compute a parametric estimator K θ of function K. Next, considering some distance d between curves, the test statistic is

T n = d( K, Kθ ),
where K is dened in (3.7). H 0 is rejected when T n > s α , where s α is a critical value that ensures that the procedure achieves level α. In the following, we will consider the

particular case d( K, K θ ) = [ ∫ 1 0 ( K(v) -K θ (v) ) 2 dv] 1/2 .
To estimate θ, one can either rely on a semiparametric maximum likelihood procedure, as it is done in [START_REF] Shih | Inferences on the association parameter in copula models for bivariate survival data[END_REF], or take θ = arg min θ∈Θ d( K, K θ ), which has been done in [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF] and seems more natural in our framework. Therefore, we will use this second approach, and our test statistic may be rewritten as T n = min θ∈Θ d( K, K θ ).

To compute the critical values, a bootstrap procedure is required. In our framework, [START_REF] Wang | Model selection and semiparametric inference for bivariate failure-time data[END_REF] proposed a bootstrap methodology, which has been shown to fail to be consistent by [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF]. Therefore, we prefer to adapt the consistent resampling plan dened in [START_REF] Genest | Goodness-of-t procedures for copula models based on the probability integral transformation[END_REF] to the presence of censoring. This results on using the bootstrap procedure dened in section 3.2, but replacing Step 1 by 1'. Simulate independent variables (T b i , U b i ) 1≤i≤n under the distribution dened by C θ and with marginal distributions dened by the Kaplan-Meier estimators (univariate) of T and U, which corresponds to an approximation of the law of (T, U ) under H 0 . Alternatively, in a full parametric modelling of the distribution of (T, U ), parametric distributions may be used instead of the nonparametric Kaplan-Meier estimators. Based on B bootstrap replications of T n , the critical value s α can be determined in order to ensure a level α of the procedure.

Real data example. We consider three copula models that have been used by [START_REF] Luciano | Modelling stochastic mortality for dependent lives[END_REF] to study the mortality of a particular generation in the data-set (in the present paper, we do not distinguish between generations). These models are Clayton, Frank copula models, and a copula called Nelsen 4.2.20 (corresponding to the copula dened in formula 4.2.20 in [START_REF] Nelsen | An introduction to copulas[END_REF]). Denition of these three Archimedean families in recalled in Table 2 below. Estimators θ are computed by minimization of the distance d( K, K θ ).

Model

ϕ θ (t) C θ (u, v) θ Clayton θ -1 (t -θ -1) (u -θ + v -θ + 1) -1/θ 4.8991
Frank -log

( exp(-θt)-1 exp(-θ)-1 ) -θ -1 log ( 1 + (exp(-θu)-1)(exp(-θv)-1) (exp(-θ)-1) ) 11.4115 4.2.20 exp(t -θ ) -e log ( exp(u -θ ) + exp(v -θ ) -e ) -1/θ 1.338
Table 2: Expression of the dierent copula families considered. The column θ presents the estimated association parameter on the data-set, minimizing distance d.

The graphical comparison between K and K θ is presented in Figure 2 below. Table 3 presents the results of the test procedure described above, comparing the value of the test statistic T n to the quantiles of the distribution of T n under the null hypothesis (computed using our bootstrap procedure). 4.2.20 case. However, in this last model, the corresponding p-value is small, while it is not the case for the two other models. Graphically, it seems that all the models that we consider have diculties to capture the behaviour of function K for values of v between 0.2 and 0.4.

Conclusion

The estimator that we considered in this paper is designed for applications in which the censoring times for both individuals only diers through an observed random variable. In the example that we consider, this observed variable represents the age dierence. We only considered the case of two random lifetimes (T, U ), but the procedure can easily be generalized to more lifetimes. The main diculty of this extension comes from the fact that the procedure requires to put mass only at fully observed observations. To obtain a sucient number of such observations when the number of lifetimes is high, one would need to have a large value of the sample size n. In this paper, we focused on applying our results to the study of copula models. Other applications of this technique could be considered, as regression models (see [START_REF] Lopez | Bivariate censored regression relying on a new estimator of the joint distribution function[END_REF] for related problems), or to the evaluation of various dependence measures, see [START_REF] Fan | Dependence estimation over a nite bivariate failure time region[END_REF] or Hougaard (2000).

6 Appendix 6.1 Proof of Theorem 3.1

Let (t 0 , u 0 ) denote some point in R 2 such that S F (t 0 , u 0 ) > 0.

First case:

ϕ(t, u) = 0 for t ≥ t 0 or u ≥ u 0 . One can write ∫ ϕ(t, u)d( F -F * )(t, u) = 1 n n ∑ i=1 δ i γ i [ ŜG (B i -) -S G (B i -)]ϕ(Y i , Z i ) S G (B i -) 2 + R 1n (ϕ), with |R 1n (ϕ)| ≤ sup a≤a 0 | ŜG (a) -S G (a)| 2 S G (a) ŜG (a) 2 × ( 1 n n ∑ i=1 δ i γ i Φ(Y i , Z i ) ) , (6.1)
where a 0 is some point in R such that S G (a 0 ) > 0. It follows from the uniform n 1/2 -consistency of ŜG (see Gill (1983)) that the right-hand side in (6.1), which does not depend on ϕ, is O P (n -1/2 ). Moreover, let us observe that the functions

f n (Y i , Z i , ε i , δ i , γ i ) = δ i γ i ŜG (B i -)ϕ(Y i , Z i )[S G (B i -)] -2 , f (Y i , Z i , ε i , δ i , γ i ) = δ i γ i S G (B i -)ϕ(Y i , Z i )[S G (B i -)] -2 ,
are two elements of the Donsker class H dened in Assumption 3. Moreover, using the uniform convergence rate of ŜG , one obtains that ∥f n -f ∥ ∞ → 0. Therefore, the asymptotic equicontinuity of Donsker classes (see Lemma 19.24 in van der Vaart (1998))

ensures that ∫ ϕ(t, u)d( F -F * )(t, u) = ∫ [ ŜG (a-) -S G (a-)]ϕ(t, u)dP (T,U,C,ε) (t, u, c, e) S G (a-) + R 2n (ϕ), (6.2)
where sup ϕ∈F |R 2n (ϕ)| = o P (n -1/2 ). Next, the representation (3.1) follows from [START_REF] Stute | Distributional convergence under random censorship when covariables are present[END_REF] or [START_REF] Gijbels | Almost sure asymptotic representation for a class of functionals of the Kaplan-Meier estimator[END_REF], since

ŜG (a) -S G (a) = 1 n n ∑ i=1 { η i S G (B i ∨ a) H(B i ) - ∫ 1 B i ≥u S G (u ∨ y)dF A (u) H(u)F A (u) } + { (1 -η i )1 B i >y F A (B i ) -S G (y) } + R(y),
where sup y≤y 0 |R(y)| = o P (n -1/2 ), where y 0 is some point such that P(B > y 0 ) > 0 (by assuming that ϕ is zero for large values of t or u, we ensure that, in (6.2), only terms with b smaller than such a y 0 appear). To deduce equation (3.2), it suces to observe that

∫ ϕ(t, u)d[ F -F ](t, u) = 1 n n ∑ i=1 δ i γ i ϕ(Y i , Z i ) S G (B i -) -E [ϕ(T, U )]+ 1 n n ∑ i=1 ψ ϕ (Y i , Z i , ε i )+o P (n -1/2 ).
Each of these two i.i.d. sums have zero mean, and the Central Limit Theorem applies.

General case: the general case follows from a combination of the rst case and of Lemma 6.1, in order to make (t 0 , u 0 ) tend to innity. Point 1 in Lemma 6.1 is easily checked by observing that the terms in the i. in [START_REF] Gill | Large sample behaviour of the product-limit estimator on the whole line[END_REF]). Moreover, dening

Γ n (t 0 , u 0 ) = n -1 n ∑ i=1 δ i γ i C 1/2+υ (B i -)1 Y i ≥t 0 ,Z i ≥u 0 Φ(Y i , Z i )S G (B i -) -2 ,
we see that conditions 4 and 5 in Lemma 6.1 hold, thanks to Assumption 3.

A technical Lemma

Lemma 6.1 Let I (t 0 ,u 0 ) = {t ≤ t 0 , u ≤ u 0 }, and let I denote the support of (Y, Z), and

(t 1 , u 1 ) the upper bound of the support. Let F be a class of functions. Let P n (t, u, ϕ) be a process on I (t 0 ,u 0 ) × F. Dene, for any (t, u) ∈ I, R n (t, u, ϕ) = P n (t 1 , u 1 , ϕ) -P n (t, u, ϕ).

Assume that for all (t 0 , u 0 ) such that I (t 0 ,u 0 ) is strictly included in the interior part of I, (P n (t, u, ϕ)) (t,u)∈I (t 0 ,u 0 ) ,ϕ∈F =⇒ (W (V ϕ (t))) t∈I (t 0 ,u 0 ) ,ϕ∈F ,

where W (V ϕ (t)) is a Gaussian process with covariance function V ϕ , and =⇒ denotes the weak convergence. Assume that the following conditions hold, 1. lim (t 0 ,u 0 )→(t 1 ,u 1 ) V ϕ (t 0 , u 0 ) = V ϕ (t 1 , u 1 ), with sup ϕ∈F |V ϕ (t 1 , u 1 )| < ∞, 4. Γ n (t 0 , u 0 ) → Γ(t 0 , u 0 ) in probability, 5. lim (t 0 0 )→(t 1 ,u 1 ) Γ(t 0 , u 0 ) = 0.

Then P (t 1 , u 1 , ϕ) =⇒ N (0, V ϕ (t 1 , u 1 )).

Proof. From Theorem 13.5 in [START_REF] Billingsley | of probability measures[END_REF]

Figure 1 :

 1 Figure 1: Histogram of the distribution of τ using the bootstrap procedure.

Figure 2 :

 2 Figure2: Graphical comparison between K and K θ for the three copula models considered.

  i.d. sum in the right-hand side of decomposition (3.1) have nite variance. To show that points 2-4 hold, deneI (t 0 ,u 0 ) = {(t, u) : t ≤ t 0 , u ≤ u 0 }, and observe that n 1/2 ∫ ϕ(t, u)1 (t,u) / ∈I (t 0 ,u 0 ) d[F -F * ](t, u) ≤ sup a n 1/2 ŜG (a) -S G (a) i C 1/2+υ (B i -)1 Y i ≥t 0 ,Z i ≥u 0 Φ(Y i , Z i ) S G (B i -) 2 . Dening M n = sup a n 1/2 |[ ŜG (a) -S G (a)][C(a) -1/2-υ ]| × sup a |S G (a) ŜG (a) -1 |, we get M n = O P (1) from Theorem 2.1 in[START_REF] Gill | Large sample behaviour of the product-limit estimator on the whole line[END_REF], since ∫ C(a) -1-2υ dC(a) < ∞ (see condition (2.1)

2.

  |R n (t ′ , u ′ , ϕ)| ≤ M n × Γ n (t 0 , u 0 ), for all τ ′ ∈ I -I (t 0 ,u 0 ) , 3. M n = O P (1),

  ′ ,u ′ )>(t,u),ϕ∈F |R n (t ′ , u ′ , ϕ)| > ε > M ) = 0,using the fact that M n = O P (1) (condition 3).

	As a consequence,
	lim (t,u)→(t 1 ,u 1 )	lim sup n→∞	P	( (t ) sup	≤ lim M →∞	lim sup
							)
							= 0.	(6.3)
	Using condition 2 in the Lemma, the probability in equation (6.3) is bounded, for all
	M > 0, by
		P(|Γ (6.4)
	Moreover, from condition 4,
		lim sup
	Since Γ(t, u) → 0 (condition 5), we can deduce that
		lim (t,u)→τ H	n→∞ lim sup	P(|Γ n (t, u) -Γ(t, u)| > ε/M -Γ(t, u)) = 0.
	Hence,
	lim

and from condition 1, it suces to show that, for all ε > 0, lim

(t 0 ,u 0 )→(t 1 ,u 1 ) lim sup n→∞ P ( sup (t,u)∈I-I (t 0 ,u 0 ) ,ϕ∈F |R n (t, u, ϕ)| > ε n (t, u) -Γ(t, u)| > ε/M -Γ(t, u)) + P(M n > M ). n→∞ P(|Γ n (t, u) -Γ(t, u)| > ε/M -Γ(t, u)) = 1 ε/M -Γ(t,u)≥0 . (t,u)→(t 1 ,u 1 ) lim sup n→∞ P ( sup (t ′ ,u ′ )>(t,u),ϕ∈F |R n (t ′ , u ′ , ϕ)| > ε ) ≤ lim sup n→∞ P(M n > M ). n→∞ P (M n

Acknowledgements: This work is supported by French Agence Nationale de la Recherche (ANR) ANR Grant ANR-09-JCJC-0101-01.


G Sandell

H Wiesemeyer

M A Requena-Torres

S Heyminck

R Güsten

J Stutzki

R Simon

U U Graf

GREAT [C II] and CO observations of the BD+40 • 4124 region

Keywords: ISM, molecules -Stars, circumstellar matter -Stars, pre-main sequence -Stars, variables, T Tauri, Herbig Ae/Be -(ISM, ) photon-dominated region (PDR)

The BD+40 • 4124 region was observed with high angular and spectral resolution with the German heterodyne instrument GREAT in CO J = 13 → 12 and [C II] on SOFIA. These observations show that the [C II] emission is very strong in the reflection nebula surrounding the young Herbig Ae/Be star BD+40 • 4124. A strip map over the nebula shows that the [C II] emission approximately coincides with the optical nebulosity. The strongest [C II] emission is centered on the B2 star and a deep spectrum shows that it has faint wings, which suggests that the ionized gas is expanding. We also see faint CO J = 13 → 12 at the position of BD+40 • 4124, which suggests that the star may still be surrounded by an accretion disk. We also detected [C II] emission and strong CO J = 13 → 12 toward V 1318 Cyg. Here the [C II] emission is fainter than in BD+40 • 4124 and appears to come from the outflow, since it shows red and blue wings with very little emission at the systemic velocity, where the CO emission is quite strong. It therefore appears that in the broad ISO beam the [C II] emission was dominated by the reflection nebula surrounding BD+40 • 4124, while the high J CO lines originated from the adjacent younger and more deeply embedded binary system V 1318 Cyg.

Introduction

The BD+40 • 4124, at a distance of 980 pc (Shevchenko et al. 1991), forms a a small pre-main-sequence cluster (Herbig 1960;Hillenbrand et al. 1995). Herbig (1960) included BD+40 • 4124 in his original list of Herbig Ae/Be (HAEBE) stars and noted that it had three companions: LkHα 224 (V 1686 Cyg), LkHα 225 (V 1318 Cyg), and LkHα 226 within 35 ′′ of BD+40 • 4124. Both BD+40 • 4124 and V 1686 Cyg are HAEBE stars with spectral types of B2 Ve and A7 Ve, respectively (van den Ancker et al. 1998;van den Ancker, Wesselius & Tielens 2000). Hillenbrand et al. (1995), however, classified V 1686 Cyg as B5 Ve. V 1318 Cyg is a heavily obscured binary system with a separation of ∼ 5 ′′ (Aspin, Sandell & Weintraub 1994). The southern member of this binary, V 1318 Cyg S, is more deeply embedded with a visual extinction of 10 mag or more. Several studies indicate that V 1318 Cyg may have a hidden companion about 1 ′′ to the NE that powers the outflow and the H 2 O maser (Davies et al. 2001;Marvel 2005). The spectral type of V 1318 Cyg is very uncertain, but based on analysis of the spectral energy distribution, both Aspin, Sandell & Weintraub (1994) and van den Ancker, Wesselius & Tielens (2000) argued that at least one of the stars must be an intermediate mass star. Hillenbrand et al. (1995) gave spectral classifications of mid-A to Fe for both components. A K-band image of the BD+40 • 4124 field is shown in Fig. 1. BD+40 • 4124 illuminates a bright reflection nebula with a radius of ∼ 30 ′′ (Herbig 1960; Loren 1977; Hillenbrand et al. 1995), while V 1318 Cyg excites an H 2 O maser and drives a bipolar molecular outflow (Palla et al. 1995). The BD+40 • 4124 region was studied in the infrared by several groups using both the ISO SWS and LWS spectrometers (van den Ancker, Wesselius & Tielens 2000;Creech-Eakman et al. 2002;Lorenzetti et al. 2002). ISO LWS spectra were obtained toward BD+40 • 4124, V 1686 Cyg, and V 1318 Cyg, even though they are all covered with one pointing, since the beam width for LWS is ∼ 80 ′′ . Creech-Eakman et al. (2002) and Lorenzetti et al. (2002) attributed all the [C II] and high transition CO (J = 14 → 13 to J = 16 → 15) emission to BD+40 • 4124, while van den Ancker, Wesselius & Tielens (2000) assigned about the same amount of [C II] emission to BD+40 • 4124 and V 1318 Cyg, and no emission at all to V 1686 Cyg.

In this letter we revisit the BD+40 • 4124 region with the GREAT1 heterodyne instrument on SOFIA to see where the [C II] and high J CO emission really comes from. Is it dominated by the bright reflection nebulosity illuminated by BD+40 • 4124 or does the emission come from the younger, embedded star V 1318 Cyg? Given the low spatial resolution of the ISO data, it is clear that higher resolution data are needed to determine where the [C II] and CO emission comes from and which physical mechanisms generate this emission.

Observations

BD+40 • 4124 and V 1318 Cyg were observed with GREAT on SOFIA on April 6, 2011 on a 53 minute leg at an altitude of 43,000 ft. GREAT is a modular heterodyne instrument, with two (2001). BD+40 • 4124 is the bright star in the northern part of the image, while the second bright star 26 ′′ to the south east is V 1686 Cyg. V 1318 Cyg is the double star system at the south eastern edge of the image. The southern star is the active one, showing strong H 2 emission and powers a massive bipolar molecular outflow. channels that are both used simultaneously. For a more complete description of the instrument, see Heyminck et al. (2012). For this flight, which was the first science flight with GREAT, the configuration was set to the low-frequency channel 1b (L1#b), which covers the frequency range 1.42 -1.52 THz, and the lowfrequency channel 2 (L2), covering 1.82 -1.92 THz. The L1#b channel was tuned to CO J = 13 → 12 (1496.923 GHz) in the lower sideband and the L2 channel was centered on the [C II] 2 P 3/2 → 2 P 1/2 (1900.5369 GHz) in the upper sideband. The half power beam width (HPBW) is ∼ 21 ′′ at 1.5 THz and 16 ′′ at 1.9 THz. As backend we used the fast Fourier transform spectrometers (AFFTS), which have a bandwidth of 1.5 GHz and 8192 channels and provide a frequency resolution of 0.212 MHz. The measured system temperatures were ∼ 4,500 K for [C II] and 3,000 K for CO = 13 → 12. All observations were made in dual beam switch mode using a 150 ′′ chop amplitude in equatorial coordinates at a position angle of 135 • . We took "long" integration spectra toward both BD+40 • 4124 and V1318 Cyg with total integration times of 2.5 and 3.7 minutes, respectively. We also obtained a strip map over BD+40 in RA with 15 ′′ spacing going from -15 ′′ ,0 ′′ to +45 ′′ ,0 ′′ to explore how extended the [C II] emission was. The integration time for these spectra was 0.6 minutes/position. The spectra were analyzed in CLASS, where we coadded individual spectra and removed a linear baseline. All spectra are calibrated in T *

A and corrected for a forward scattering efficiency of 0.95.

Results

BD+40 • 4124

The [C II] emission is quite strong toward BD+40 • 4124 with an antenna temperature of T A ∼ 13.5 K (Fig. 2, Table 1). In Fig. 3 we show the positions of the strip map that we obtained of BD+40 • 4124 marked as crosses on the 850 µm SCUBA image from Sandell, Weintraub & Hamidouche (2011). The strip map (Fig. 4) shows that the [C II] emission peaks within errors on Table 1. Gaussian fits to long integration spectra Source Line the star and slowly decreases to the east. There does not appear to be any interaction between the [C II] emission and the dense molecular ridge (see Fig. star. The long integration spectrum (Fig. 2) shows a broad, ∼ 6 km s -1 wide pedestal (Table 1), indicating that the ionized gas is slowly expanding outward.

T * A dV T * A ∆V V LSR [K km s -1 ] [K] [km s -1 ] [km s -1 ] BD+40 • 41247 [C
The CO J = 13 → 12 emission, which was observed simultaneously with [C II], was also detected, but only at position (0 ′′ , 0 ′′ ) toward the star itself (Figure 2, Table 1). If the line indeed had a similar strength in the surrounding reflection nebula, it should have been detected when we took an average of the three spectra east of the star in the strip-map. However, there is no hint of CO emission in these positions in the strip-map, confirming that the CO emission comes from BD+40 • 4124 proper, which therefore must be surrounded by hot, molecular gas, presumably in an accretion disk.

V 1318 Cyg

We only obtained one long integration spectrum toward V 1318 Cyg (Fig. 2). CO J = 13 → 12 was easily detected with T * A ∼ 1.7 K. The emission is dominated by faint, broad emission extending from ∼ -5 km s -1 on the blueshifted side to ∼ +22 km s -1 on the redshifted side. A two-component Gaussian fit gives a FWHM of 3.3 km s -1 for the narrow component and 11.7 km s -1 for the wide pedestal (Table 1). The broad component is clearly associated with the outflow powered by V 1318 Cyg (Palla et al. 1995) and dominates the CO emission, i.e., the integrated line emission is ∼1.5 times higher than the emission from the narrow component. The radial velocity of the narrow component, 8.1 km s -1 , agrees well with the systemic velocity of the star, 7.9 km s -1 determined from high spatial resolution CO J = 2 → 1 observations [START_REF] Looney | [END_REF]) and is likely to originate in the accretion disk or the surrounding envelope.

In contrast, the [C II] emission looks quite different. At first glance it looks as if the spectrum is self-absorbed. However, there is no reason why the [C II] would be affected by selfabsorption. Instead it appears that the emission is completely dominated by the outflow. Even though the CO J = 13 → 12 line profile suggests that the blueshifted wing is somewhat stronger than the redshifted one, the difference is quite marginal, while the blueshifted emission completely dominates in [C II]. Most of the emission is at low velocities and we see no "high-velocity" Fig. 2. Deep [C II] and CO J = 13 → 12 spectra toward BD+40 • 4124 and V 1318 Cyg. [C II] is plotted in black and the CO line profiles in red. For BD+40 • 4124, where the CO J = 13 → 12 line is quite faint, we multiplied the spectrum with a factor of four. The systemic velocity is marked by a vertical gray line.

gas in [C II] . Since we have not mapped the emission, we cannot say how extended it is. There is probably some [C II] emission throughout the outflow lobes, although the [C II] is likely to be stronger close the exciting star, see the discussion section below.

Discussion

[C II] emission

Our observations of [C II] in the BD+40

• 4124 region show that most of the [C II] emission comes from the reflection nebulosity surrounding the B2 star. Based on our strip map, we estimate the [C II] emission to have a FWHM of ∼ 60 ′′ similar to that of the reflection nebula. The integrated line inten- (2000) found for BD+40 • 4124, although they assigned a similar line intensity to V 1318 Cyg. Even though we did not map the region around V 1318 Cyg, it is quite clear from our [C II] spectrum of V 1318 Cyg (Figure 2, Table 1) that the emission is dominated by the outflow. The line intensity toward V 1318 Cyg is only about a third of what we see toward BD+40 • 4124.

There are very few observations of outflows in [C II] that have both sufficient spatial and/or spectral resolution to investigate the source of the [C II] emission. The jet-like low-mass outflow HH 46 was observed with the integral field spectrometer PACS centered on the low-mass protostar (van Kempen et al. 2010). These observations show that the [C II] emission is about twice as strong in the blueshifted outflow lobe compared to the redshifted outflow lobe or the central protostar. Since there is no velocity information in these data, it is not clear whether the emission comes from low-velocity UV-heated gas in the cavity walls or whether [C II] is mixed with the molecular high-velocity gas. Unfortunately, there were no observations far away from the central source, which makes it difficult to check whether [C II] requires direct UV-excitation or whether Cshocks in the outflow would be equally efficient. The similar line strength throughout the red-and blueshifted outflow lobes suggests that in a low-mass outflow we mostly see emission from the cavity walls, which would agree with what we see toward V 1318 Cyg with GREAT, where we only see [C II] at low velocities. HIFI observations of the massive DR 21 outflow, however, show [C II] emission that is as broad as the molecular line emission (Ossenkopf et al. 2010). DR 21, however, is a high-mass star formation region with intense UV-emission.

The size of the outflow powered by V 1318 Cyg is ∼ 55 ′′ (Palla et al. 1995). If we assume that the [C II] emission is similar in strength as what is seen toward V 1318 Cyg, we would expect something like 15 × 10 -15 W m -2 . This confirms that the reflection nebula surrounding BD+40 • 4124 dominates the [C II] emission in the large ISO beam with a minor contribution from the V 1318 Cyg outflow.

CO J = 13 → 12 emission

The detection of CO J = 13 → 12 emission towards BD+40 • 4124 is intriguing. Sandell, Weintraub & Hamidouche (2011) did not detect any submillimeter continuum emission towards the star with an upper limit of 35 mJy beam -1 at 850 µm, suggesting that the star may already have dispersed its disk. Skinner, Brown, & Stewart (1993), however, reported a marginal detection of free-free emission at 3.6 and 6 cm, which, if real, could mean that there is still some remnant disk around the star, since free-free emission from early B-stars is generally believed to originate from photoionization of circumstellar disks (Hollenbach, Yorke, & Johnstone 2000). Because we only see hot CO toward the star, it is unlikely that it would come from the surrounding envelope and it is therefore more plausible that it originates in a circumstellar disk that surrounds the star. Such a disk would have moderate inclination, because the extinction toward the star is low, A V ∼ 3 mag (van den Ancker et al. 1998), which agrees with the relatively narrow linewidth of the CO 13 → 12 line, 2.5 km s -1 (Table 1). Looking for high J CO emission may therefore be a good way to probe the gas in disks around early B-stars, since the gas is expected to be much hotter than in low-mass stars, owing to the much higher UV-field that illuminates and heats the disk. It would therefore be very valuable to search for hot CO in the few early B-stars that show strong evidence for a circumstellar disk, such as HD200775, which illuminates the reflection nebula NGC 7023 (Okamoto et al. 2009), or MWC 349, which has a largely ionized disk, see e.g. Sandell, Weintraub & Hamidouche (2011).

In contrast, the CO J = 13 → 12 emission from the deeply embedded V 1318 Cyg is much stronger than for BD+40 • 4124, and dominated by the hot outflow, not by the disk. The CO line has a broad pedestal with a FWHM of ≥ 10 km s -1 , which comes from the outflow and contributes to more than half of the line intensity. There is also a narrow component, which could come from the accretion disk, although it could also originate in the dense envelope that surrounds the star.

Conclusions

Observations with GREAT on SOFIA toward the BD+40 • 4124 group show that the [C II] emission seen by ISO is dominated by emission from the reflection nebula illuminated by BD+40 • 4124. However, we do also see [C II] emission from the outflow powered by V 1318 Cyg, which is a deeply embedded young HAEBE star ∼ 35 ′′ southeast of BD+40 • 4124. On the other hand, the high J CO emission observed by ISO is completely dominated by hot gas in the outflow from V 1318 Cyg. We also detected faint CO J = 13 → 12 emission toward BD+40 • 4124, which suggests that the stars is still surrounded by an accretion disk.

Fig. 1 .

 1 Fig. 1. A K-band image of the BD+40 • 4124 field from Davies et al.(2001). BD+40 • 4124 is the bright star in the northern part of the image, while the second bright star 26 ′′ to the south east is V 1686 Cyg. V 1318 Cyg is the double star system at the south eastern edge of the image. The southern star is the active one, showing strong H 2 emission and powers a massive bipolar molecular outflow.

  3), which extends from V 1318 Cyg to the north past BD+40 • 4124(Sandell, Weintraub & Hamidouche 2011;[START_REF] Looney | [END_REF]. There are no changes in the [C II] brightness as it crosses over the ridge, nor is there any change in radial velocity or line width, suggesting that BD+40 • 4124 is in the foreground, and not directly connected to the molecular ridge. The extent of the [C II] emission with a full width half maximum (FWHM) of ∼ 60 ′′ , or about the same size as the reflection nebulosity surrounding the star, confirms that the line emission originates from the reflection nebulosity, and is almost certainly photo-ionized by the UV-emission from the young B2

Fig. 3 .

 3 Fig. 3. 850 µm continuum image of the BD+40 • 4124 region from (Sandell, Weintraub & Hamidouche 2011). The positions observed in [C II] strip map are marked with black crosses. BD+40 • 4124 is indicated in the figure. The two star symbols at the peak of the 850 µm contours are V 1318 Cyg N and S, while the one to the west of them is V 1686 Cyg.

Fig. 4 .

 4 Fig. 4. Strip map in RA over BD+40 • 4124 in [C II] with 15 ′′ spacing. The map suggests that [C II] peaks slightly west of the star, although the offset could be caused by the pointing.

GREAT is a development by the MPI für Radioastronomie and the KOSMA/ Universität zu Köln, in cooperation with the MPI für Sonnensystemforschung and the DLR Institut für Planetenforschung.

Acknowledgements. Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy. SOFIA Science Mission Operations are conducted jointly by the Universities Space Research Association, Inc., under NASA contract NAS2-97001, and the Deutsches SOFIA Institut under DLR contract 50 OK 0901. We also thank Hans Zinnecker for a critical reading of the paper.